In this thesis, we consider the problem of solving large and sparse linear systems of saddle point type stemming from optimization problems. The focus of the thesis is on iterative methods, and new preconditioning srategies are proposed, along with novel spectral estimtates for the matrices involved.

This Ph.D thesis focuses on iterative regularization methods for regularizing linear and nonlinear ill-posed problems. Regarding linear problems, three new stopping rules for the Conjugate Gradient method applied to the normal equations are proposed and tested in many numerical simulations, including some tomographic images reconstruction problems. Regarding nonlinear problems, convergence and convergence rate results are provided for a Newton-type method with a modified version of Landweber iteration as an inner iteration in a Banach space setting.

In this thesis we present a new and accurate series of computation methods for compressible multi-phase flows with capillary effects based upon the full seven-equation Baer-Nunziato model. For that reason, there are some numerical methods to obtain high accuracy solutions, which will be shown here. First, a high resolution shock capturing Total Variation Diminishing (TVD) finite volume scheme is used on both Cartesian and unstructured triangular grids. Regarding the TVD finite volume scheme on the unstructured grid, time-accurate local time stepping (LTS) is applied to compute the solutions of the governing PDE system, in which the results are also compared with time-accurate global time stepping. Second, we propose a novel high order accurate numerical method for the solution of the seven equation Baer-Nunziato model based on ADER discontinuous Galerkin (DG) finite element schemes combined with a posteriori subcell finite volume limiting and adaptive mesh refinement (AMR). In multi-phase flows, the difficulty is to design accurate numerical methods for resolving the phase interface, as well as the simulation of the phenomena occurring at the interface, such as surface tension effects, heat transfer and friction. This is because of the interactions of the fluids at the phase interface and its complex geometry. So the accurate simulation of compressible multi-phase flows with surface tension effects is currently still one of the most challenging problems in computational fluid dynamics (CFD). In this work, we present a novel path-conservative finite volume discretization of the continuum surface force method (CSF) of Brackbill et al. to account for the surface tension effect due to curvature of the phase interface. This is achieved in the context of a diffuse interface approach, based on the seven equation Baer-Nunziato model of compressible multi-phase flows. Such diffuse interface methods for compressible multi-phase flows including capillary effects have first been proposed by Perigaud and Saurel. Regarding the high order accuracy of a diffuse interface approach, the interface is captured and allowed to travel across one single possibly refined cell, and is computed in the context of multi-dimensional high accurate space/time DG schemes with AMR and a posteriori sub-cell stabilization. The surface tension terms of the CSF approach are considered as a part of the non-conservative hyperbolic system. We propose to integrate the CSF source term as a non-conservative product and not simply as a source term, following the ideas on path conservative finite volume schemes put forward by Castro and Parés. For the validation of the current numerical methods, we compare our numerical results with those published previously in the literature.

In this work the numerical discretization of the partial differential governing equations for compressible and incompressible flows is dealt within the discontinuous Galerkin (DG) framework along space-time adaptive meshes. Two main research fields can be distinguished: (1) fully explicit DG methods on collocated grids and (2) semi-implicit DG methods on edge-based staggered grids. DG methods became increasingly popular in the last twenty years mainly because of three intriguing properties: i) non-linear L2 stability has been proven; ii) arbitrary high order of accuracy can be achieved by simply increasing the polynomial order of the chosen basis functions, used for approximating the state-variables; iii) high scalability properties make DG methods suitable for large-scale simulations on general unstructured meshes. It is a well known fact that a major weakness of high order DG methods lies in the difficulty of limiting discontinuous solutions, which generate spurious oscillations, namely the so-called ’Gibbs phenomenon’. Over the years, several attempts have been made to cope with this problem and different kinds of limiters have been proposed. Among them, a rather intriguing paradigm has been defined in the work of [71], in which the nonlinear stabilization of the scheme is sequentially and locally introduced only for troubled cells on the basis of a multidimensional optimal order detection (MOOD) criterion. In the present work the main benefits of the MOOD paradigm, i.e. the computational robustness even in the presence of strong shocks, are preserved and the numerical diffusion is considerably reduced also for the limited cells by resorting to a proper sub-grid. In practice the method first produces a so-called candidate solution by using a high order accurate unlimited DG scheme. Then, a set of numerical and physical detection criteria is applied to the candidate solution, namely: positivity of pressure and density, absence of floating point errors and satisfaction of a discrete maximum principle in the sense of polynomials. Then, in those cells where at least one of these criteria is violated the computed candidate solution is detected as troubled and is locally rejected. Next, the numerical solution of the previous time step is scattered onto cell averages on a suitable sub-grid in order to preserve the natural sub-cell resolution of the DG scheme. Then, a more reliable numerical solution is recomputed a posteriori by employing a more robust but still very accurate ADER-WENO finite volume scheme on the sub-grid averages within that troubled cell. Finally, a high order DG polynomial is reconstructed back from the evolved sub-cell averages. Moreover, handling typical multiscale problems, dynamic adaptive mesh refinement (AMR) and adaptive polynomial order methods are probably the two main ways of preserving accuracy and efficiency, and saving computational effort. The here adopted AMRapproach is the so called ’cell by cell’ refinement because of its formally very simple tree-type data structure. In the here-presented ’cell-by-cell’ AMR every single element is recursively refined, from a coarsest refinement level l0 = 0 to a prescribed finest (maximum) refinement level lmax, accordingly to a refinement-estimator function X that drives step by step the choice for recoarsening or refinement. The combination of the sub-cell resolution with the advantages of AMR allows for an unprecedented ability in resolving even the finest details in the dynamics of the fluid. First, the Euler equations of compressible gas dynamics and the magnetohydrodynamics (MHD) equations have been treated [281]. Then, the presented method has been readily extended to the special relativistic ideal MHD equations [280], but also the the case of diffusive fluids, i.e. fluid flows in the presence of viscosity, thermal conductivity and magnetic resistivity [116]. In particular, the adopted formalism is quite general, leading to a novel family of adaptive ADER-DG schemes suitable for hyperbolic systems of partial differential equations in which the numerical fluxes also depend on the gradient of the state vector because of the parabolic nature of diffusive terms. The presented results show clearly that the shock-capturing capability of the news schemes are significantly enhanced within the cell-by-cell Adaptive Mesh Refinement (AMR) implementation together with time accurate local time stepping (LTS). The resolution properties of the new scheme have been shown through a wide number of test cases performed in two and in three space dimensions, from low to high Mach numbers, from low to high Reynolds regimes. In particular, concerning MHD equations, the divergence-free character of the magnetic field is taken into account through the so-called hyperbolic ’divergence-cleaning’ approach which allows to artificially transport and spread the numerical spurious ’magnetic monopoles’ out of the computational domain. A special treatment has been followed for the incompressible Navier-Stokes equations. In fact, the elliptic character of the incompressible Navier-Stokes equations introduces an important difficulty in their numerical solution: whenever the smallest physical or numerical perturbation arises in the fluid flow then it will instantaneously affect the entire computational domain. Thus, a semi-implicit approach has been used. The main advantage of making use of a semi-implicit discretization is that the numerical stability can be obtained for large time-steps without leading to an excessive computational demand [117]. In this context, we derived two new families of spectral semi-implicit and spectral space-time DG methods for the solution of the two and three dimensional Navier-Stokes equations on edge-based staggered Cartesian grids [115], following the ideas outlined in [97] for the shallow water equations. The discrete solutions of pressure and velocity are expressed in the form of piecewise polynomials along different meshes. While the pressure is defined on the control volumes of the main grid, the velocity components are defined on edge-based dual control volumes, leading to a spatially staggered mesh. In the first family, high order of accuracy is achieved only in space, while a simple semi-implicit time discretization is derived by introducing an implicitness factor theta in [0.5, 1] for the pressure gradient in the momentum equation. The real advantages of the staggering arise after substituting the discrete momentum equation into the weak form of the continuity equation. In fact, the resulting linear system for the pressure is symmetric and positive definite and either block penta-diagonal (in 2D) or block hepta-diagonal (in 3D). As a consequence, the pressure system can be solved very efficiently by means of a classical matrix-free conjugate gradient method. Moreover, a rigorous theoretical analysis of the condition number of the resulting linear systems and the design of specific preconditioners, using the theory of matrix-valued symbols and Generalized Locally Toeplitz (GLT) algebra has been successfully carried out with promising results in terms of numerical efficiency [102]. The resulting algorithm is stable, computationally very efficient, and at the same time arbitrary high order accurate in both space and time. The new numerical method has been thoroughly validated for approximation polynomials of degree up to N = 11, using a large set of non-trivial test problems in two and three space dimensions, for which either analytical, numerical or experimental reference solutions exist. Moreover, the here mentioned semi-implicit DG method has been successfully extended to a novel edge-based staggered ’cell-by-cell’ adaptive meshes [114].

The main topic of the thesis is the synthesis of control laws for interacting agent-based dynamics and their mean-field limit. In particular, after a general introduction, in the second chapter a linearization-based approach is used for the computation of sub-optimal feedback laws obtained from the solution of differential matrix Riccati equations. Quantification of dynamic performance of such control laws leads to theoretical estimates on suitable linearization points of the nonlinear dynamics. Subsequently, the feedback laws are embedded into a nonlinear model predictive control framework where the control is updated adaptively in time according to dynamic information on moments of linear mean-field dynamics. The performance and robustness of the proposed methodology is assessed through different numerical experiments in collective dynamics. In the other chapters of the thesis I present some related projects, robustness of systems with uncertainties, a proximal gradient approach for sparse control and an application in crowd evacuation dynamics.

IIn this work we present a new class of well-balanced, arbitrary high order accurate semi-implicit discontinuous Galerkin methods for the solution of the shallow water and incompressible Navier-Stokes equations on staggered unstructured curved meshes. Isoparametric finite elements are used to take into account curved domain boundaries. Regarding two-dimensional shallow water equations, the discrete free surface elevation is defined on a primal triangular grid, while the discrete total height and the discrete velocity field are defined on an edge-based staggered dual grid. Similarly, for the two-dimensional incompressible Navier-Stokes case, the discrete pressure is defined on the main triangular grid and the velocity field is defined on the edge-based staggered grid. While staggered meshes are state of the art in classical finite difference approximations of the incompressible Navier-Stokes equations, their use in the context of high order DG schemes is novel and still quite rare. High order (better than second order) in time can be achieved by using a space-time finite element framework, where the basis and test functions are piecewise polynomials in both space and time. Formal substitution of the discrete momentum equation on the dual grid into the discrete continuity equation on the primary grid yields a very sparse system for the scalar pressure involving only the direct neighbor elements, so that it becomes a block four-point system in 2D and a block five-point system for 3D tetrahedral meshes. The resulting linear system is conveniently solved with a matrix-free GMRES algorithm. Note that the same space-time DG scheme on a collocated grid would lead to ten non-zero blocks per element in 2D and seventeen non-zero blocks in 3D, since substituting the discrete velocity into the discrete continuity equation on a collocated mesh would involve also neighbors of neighbors. From numerical experiments we find that our linear system is well-behaved and that the GMRES method converges quickly even without the use of any preconditioner, which is a unique feature in the context of high order implicit DG schemes. A very simple and efficient Picard iteration is then used in order to derive a space-time pressure correction algorithm that achieves also high order of accuracy in time, which is in general a non-trivial task in the context of high order discretizations for the incompressible Navier-Stokes equations. The special case of high order in space low order in time allows us to recover further regularity about the main linear system for the pressure, such as the symmetry and the positive semi-definiteness in the general case. This allows us to use a very fast linear solver such as the conjugate gradient (CG) method. The flexibility and accuracy of high order space-time DG methods on curved unstructured meshes allows to discretize even complex physical domains with very coarse grids in both space and time. We will further extend the previous method to three-dimensional incompressible Navier-Stokes system using a tetrahedral main grid and a corresponding face-based hexaxedral dual grid. The resulting dual mesh consists in non-standard 5-vertex hexahedral elements that cannot be represented using tensor products of one dimensional basis functions. Indeed a modal polynomial basis will be used for the dual mesh. This new family of numerical schemes is verified by solving a series of typical numerical test problems and by comparing the obtained numerical results with available exact analytical solutions or other numerical reference data. Furthermore, the comparison with available experimental results will be presented for incompressible Navier-Stokes equations.

In this work we develop a new class of high order accurate Arbitrary-Lagrangian-Eulerian (ALE) one-step finite volume schemes for the solution of nonlinear systems of conservative and non-conservative hyperbolic partial differential equations. The numerical algorithm is designed for two and three space dimensions, considering moving unstructured triangular and tetrahedral meshes, respectively. As usual for finite volume schemes, data are represented within each control volume by piecewise constant values that evolve in time, hence implying the use of some strategies to improve the order of accuracy of the algorithm. In our approach high order of accuracy in space is obtained by adopting a WENO reconstruction technique, which produces piecewise polynomials of higher degree starting from the known cell averages. Such spatial high order accurate reconstruction is then employed to achieve high order of accuracy also in time using an element-local space-time finite element predictor, which performs a one-step time discretization. Specifically, we adopt either the continuous Galerkin (CG) predictor, which does not allow discontinuities in time and is suitable for smooth time evolutions, or the discontinuous Galerkin (DG) predictor which can handle stiff source terms that might produce jumps in the local space-time solution. Since we are dealing with moving meshes the elements deform while the solution is evolving in time, hence making the use of a reference system very convenient. Therefore, within the space-time predictor, the physical element is mapped onto a reference element using a high order isoparametric approach, where the space-time basis and test functions are given by the Lagrange interpolation polynomials passing through a predefined set of space-time nodes. The computational mesh continuously changes its configuration in time, following as closely as possible the flow motion. The entire mesh motion procedure is composed by three main steps, namely the Lagrangian step, the rezoning step and the relaxation step. In order to obtain a continuous mesh configuration at any time level, the mesh motion is evaluated by assigning each node of the computational mesh with a unique velocity vector at each timestep. The node solver algorithm preforms the Lagrangian stage, while we rely on a rezoning algorithm to improve the mesh quality when the flow motion becomes very complex, hence producing highly deformed computational elements. A so-called relaxation algorithm is finally employed to partially recover the optimal Lagrangian accuracy where the computational elements are not distorted too much. We underline that our scheme is supposed to be an ALE algorithm, where the local mesh velocity can be chosen independently from the local fluid velocity. Once the vertex velocity and thus the new node location has been determined, the old element configuration is connected with the new one at the future time level with straight edges to represent the local mesh motion, in order to maintain algorithmic simplicity. The final ALE finite volume scheme is based directly on a space-time conservation formulation of the governing system of hyperbolic balance laws. The nonlinear system is reformulated more compactly using a space-time divergence operator and is then integrated on a moving space-time control volume. We adopt a linear parametrization of the space-time element boundaries and Gaussian quadrature rules of suitable order of accuracy to compute the integrals. In our algorithm either a simple and robust Rusanov-type numerical flux or a more sophisticated and less dissipative Osher-type numerical flux is employed. We apply the new high order direct ALE finite volume schemes to several hyperbolic systems, namely the multidimensional Euler equations of compressible gas dynamics, the ideal classical and relativistic magneto-hydrodynamics (MHD) equations and the non-conservative seven-equation Baer-Nunziato model of compressible multi-phase flows with stiff relaxation source terms. Numerical convergence studies as well as several classical test problems will be shown to assess the accuracy and the robustness of our schemes. Furthermore we focus on the following issues to improve the algorithm efficiency: the time evolution, the numerical flux computation across element boundaries and the high order WENO reconstruction procedure. First, a time-accurate local time stepping (LTS) algorithm for unstructured triangular meshes is derived and presented, where each element can run at its own optimal time step, given by a local CFL stability condition. Then, we propose a new and efficient quadrature-free formulation for the flux computation, in which the space-time boundaries of each element are split into simplex sub-elements. This leads to space-time normal vectors as well as Jacobian matrices that are constant within each sub-element, hence allowing the flux integrals to be evaluated on the space-time reference control volume once and for all analytically during a preprocessing step. Finally, we consider the very new a posteriori MOOD paradigm, recently proposed for the Eulerian framework, to overcome the expensive WENO approach on moving meshes. The MOOD technique requires the use of only one central reconstruction stencil because the limiting procedure is carried out a posteriori instead of a priori, as done in the WENO formulation.

Model order reduction is a technique to reduce computational times of parameterized PDEs while maintaining good accuracy of the approximated solution. Reduced basis methods (RB) are the most common algorithms for reducing the complexity of parameterized PDEs and nowadays they are widely applied and very actively researched in numerous fields. We propose two ideas to further enhance model reduction: the Fundamental Order Reduction Method (FOR) and offline error estimators for RB methods. The FOR method uses nonlinear combinations of the solutions to build the reduced model and use simple affine evaluations to execute the online stage. On the other hand, offline estimators are a class of estimators that move a-posteriori operations to the offline stage, reducing in this way the load of computations in the online stage. We apply these two ideas to an EEG equation which is useful for detecting the position where an epilepsy seizure begins inside the brain. We present two known ways to solve this equation: direct approach and subtraction approach, and show theoretical and numerical results of the application of the RB and FOR methods. We prove that is not feasible to apply model reduction in the direct approach but show that it is possible in the subtraction approach. Afterwards we solve the inverse problem associated with the EEG equation using a combination of the FOR method and neural networks.

Recent advances in medical science regarding the role of the venous system in the development of neurological conditions has renewed the attention of researchers in this district of the cardiovascular system. The main goal of this thesis is to perform a theoretical study of Chronic CerebroSpinal Venous Insufficiency (CCSVI), a venous pathology that has been associated to Multiple Sclerosis. CCSVI is a condition in which main cerebral venous drainage pathways are obstructed. Its impact in cerebral hemodynamics and its connection to Multiple Sclerosis is subject of current debate in the medical community. In order to perform a credible study of the haemodynamical aspects of CCSVI, a sufficiently accurate mathematical model of the problem under investigation must be used. The venous system has not received the same attention as the arterial counterpart by the medical community. As a consequence, the mathematical modeling and numerical simulation of the venous system lies far behind that of the arterial system. The venous system is a low-pressure system, formed by very thin-walled vessels, if compared to arteries, that are likely to collapse under the action of gravitational or external forces. These properties set special requirements on the mathematical models and numerical schemes to be used. In this thesis we present a closed-loop multi-scale mathematical model of the cardiovascular system, where medium to large arteries and veins are represented as one-dimensional (1D) vessels, whereas the heart, the pulmonary circulation, capillary beds and intracranial pressure are modeled as lumped parameter models. A characteristic feature of our closed-loop model is the detailed description of head and neck veins. Due to the large inter-subject variability of the venous system, we perform a patient-specific characterization of major veins of the head and neck using MRI data collected in collaboration with the Magnetic Resonance Research Facility of the Wayne State University, Detroit (USA). Computational results are carefully validated using published data for the arterial system and most regions of the venous system. For head and neck veins validation is carried out through a detailed comparison of simulation results against patient-specific Phase-Contrast MRI flow quantification data. Regarding the development of novel numerical schemes, we construct high-order accurate, robust and efficient numerical schemes for 1D blood flow in elastic and viscoelastic vessels, as well as a solver for vessel networks. The solver is validated in the context of an in vitro network of vessels for which experimental and numerical results are available. After validation of both, the mathematical model and the numerical methodology, we use our theoretical tool to study the influence of different CCSVI patterns on cerebral hemodynamics. CCSVI patterns are defined by the medical literature as combinations of venous obstructions at different locations. Here we used two strategies. First, we take a venous configuration corresponding to a healthy control and explore the effect of different CCSVI patterns by modifying this network. Then, we characterize our venous network with the geometry of a real CCSVI patient and compare results with the ones obtained for the healthy control. The presented model provides a powerful tool to study still unresolved aspects of cerebral blood flow physiology, as well as several venous pathologies. Furthermore, it constitutes an ideal platform for improving currently used algorithms and for integrating fundamental physiological processes, such as detailed hemodynamics, regulatory mechanisms and transport of substances.

In the present work we solve systems of partial differential equations (PDE) for hyperbolic conservation laws using semi-implicit numerical methods on staggered meshes applied both to the class of finite volume (FV) and both to the family of high order Discontinuous Galerkin (DG) finite elements schemes. In particular, we want to show that these new semi-implicit schemes can be applied in several fields of applied sciences, such as in geophysical flows and compressible fluids in compliant tubes. Inside this thesis we distinguish two big parts. First, we consider staggered semi-implicit schemes for compressible viscous fluids flowing in elastic pipes. This topic is very important in several practical applications of civil, environmental, industrial and biomedical engineering. Here, we analyse the accuracy and the computational efficiency of fully explicit and semi-implicit 1D and 2D finite volume schemes for the simulation of highly unsteady viscous compressible flows in laminar regime in axially symmetric rigid and elastic pipes. We consider two families of differential models that can be used to predict the pressure and velocity distribution along the tube. One is the so called 2Dxr PDE model which is derived from the full compressible Navier-Stokes equations under the assumptions of a hydrostatic pressure and an axially symmetric geometry. The second family is a simple 1D non-conservative PDE system based on the cross-sectionally averaged version of the Navier-Stokes equations in cylindrical coordinates. In this last case, the use of a simple steady friction model is not enough to simulate the wall friction phenomena in highly transient regime. As a consequence the wall friction model has to be frequency dependent and, following previous studies present in the literature, we consider the classes of convolution integral (CI) models and instantaneous acceleration (IA) models. We carry out a rather complete analysis of the previously-mentioned methods for the simulation of flows characterized by fast transient regime in rigid and compliant tubes. The numerical results show that the convolution integral models are clearly better than instantaneous acceleration models concerning accuracy. Moreover, for CI models, instead of computing the convolution integrals, which is very time- and memoryconsuming, we express these methods via a set of additional ODEs for appropriate auxiliary variables. This trick improves the computational efficiency of these methods substantially, since it avoids the direct computation of the convolution integral. In addition, semi-implicit finite volume methods are significantly superior to classical explicit finite volume schemes in terms of computational efficiency, however, providing the same level of accuracy. We then proceed by extending the finite volume discretization of the 1D and 2Dxr PDE models to arbitrary high-order of accuracy in space introducing a new SIDG scheme on staggered meshes. Both models include the effects of the viscosity and of the wall motion. The nonlinear convective terms are discretized explicitly by using a classical RKDG scheme of arbitrary high-order of accuracy in space and third order of accuracy in time. The continuity equation is integrated over the elements that belong to the main grid, while the momentum equation is integrated over the control volumes of the edgebased staggered dual grid. Inserting the discrete momentum equation into the discrete continuity leads to a mildly nonlinear algebraic system for the degrees of freedom of the pressure, which is solved by using the (nested) Newton method of Brugnano, Casulli and Zanolli. We use the -method in order to get second order of accuracy in time for the implicit part of the scheme. In addition, the schemes have to obey only a mild CFL condition based on the fluid velocity and not based on the sound speed; consequently these schemes work also in the low Mach number regime and even in the incompressible limit of the Navier-Stokes equations. This is a very important property, which is the so-called asymptotic preserving (AP) property of the scheme. We carry out several numerical tests in order to validate this novel family of numerical methods against available exact solutions and experimental data. We also report numerical convergence tables in order to show that the new schemes indeed achieve high order of accuracy in space. In the second part of the thesis, we present a new class of a posteriori sub-cell finite volume limiters for spatially high order accurate semi-implicit discontinuous Galerkin schemes on staggered Cartesian grids for the solution of the 1D and 2D shallow water equations (SWE) and of the Euler equations both expressed in conservative form. Here, the starting point is the unlimited arbitrary high order accurate staggered SIDG scheme proposed by Dumbser and Casulli (2013). For this metho d, the mass conservation equation and the momentum equations are integrated using a discontinuous finite element strategy on staggered control volumes, where the discrete free surface elevation is defined on the main grid and the discrete momentum is defined on edge-based staggered dual control volumes. According to the semi-implicit approach, pressure terms are discretized implicitly, while the nonlinear convective terms are discretized explicitly. Inserting the momentum equations into the discrete continuity equation leads to a well conditioned block diagonal linear system for the free surface elevation which can be efficiently solved with modern iterative methods. Furthermore, the staggered SIDG is also extended to the Euler equations of compressible gasdynamics. Here, the governing PDE are rewritten using a flux vector splitting technique. The convective terms are updated using an explicit Runge-Kutta DG integrator. Then, the discrete momentum equation, which is integrated again on the dual grid, is coupled with the discrete energy equation that is discretized on the control volumes of the main grid. The pressure is efficiently obtained solving a linear system combined with an iterative Picard iteration procedure.