131 research outputs found
Stabilized mixed approximation of axisymmetric Brinkman flows
This paper is devoted to the numerical analysis of an augmented finite element approximation of the axisymmetric Brinkman equations. Stabilization of the variational formulation is achieved by adding suitable Galerkin least-squares terms, allowing us to transform the original problem into a formulation better suited for performing its stability analysis. The sought quantities (here velocity, vorticity, and pressure) are approximated by Raviart−Thomas elements of arbitrary order k ≥ 0, piecewise continuous polynomials of degree k + 1, and piecewise polynomials of degree k, respectively. The well-posedness of the resulting continuous and discrete variational problems is rigorously derived by virtue of the classical Babuška–Brezzi theory. We further establish a priori error estimates in the natural norms, and we provide a few numerical tests illustrating the behavior of the proposed augmented scheme and confirming our theoretical findings regarding optimal convergence of the approximate solutions
A Computational Framework for Axisymmetric Linear Elasticity and Parallel Iterative Solvers for Two-Phase Navier–Stokes
This dissertation explores ways to improve the computational efficiency of linear elasticity and the variable density/viscosity Navier--Stokes equations. While the approaches explored for these two problems are much different in nature, the end goal is the same - to reduce the computational effort required to form reliable numerical approximations.\\
The first topic considered is the axisymmetric linear elasticity problem. While the linear elasticity problem has been studied extensively in the finite-element literature, to the author\u27s knowledge, this is the first study of the elasticity problem in an axisymmetric setting. Indeed, the axisymmetric nature of the problem means that a change of variables to cylindrical coordinates reduces a three-dimensional problem into a decoupled one-dimensional and two-dimensional problem. The change of variables to cylindrical coordinates, however, affects the functional form of the divergence operator and the definition of the inner products. To develop a computational framework for the linear elasticity problem in this context, a new projection operator is defined that is tailored to the cylindrical form of the divergence and inner products. Using this framework, a stable finite-element quadruple is derived for . These computational rates are then validated with a few computational examples.\\
The second topic addressed in this work is the development of a new Schur complement approach for preconditioning the two-phase Navier--Stokes equations. Considerable research effort has been invested in the development of Schur complement preconditioning techniques for the Navier--Stokes equations, with the pressure-convection diffusion (PCD) operator and the least-squares commutator being among the most popular. Furthermore, more recently researchers have begun examining preconditioning strategies for variable density / viscosity Stokes and Navier--Stokes equations. This work contributes to recent work that has extended the PCD Schur complement approach for single phase flow to the variable phase case. Specifically, this work studies the effectiveness of a new two-phase PCD operator when applied to dynamic two-phase simulations that use the two-phase Navier--Stokes equations. To demonstrate the new two-phase PCD operators effectiveness, results are presented for standard benchmark problems, as well as parallel scaling results are presented for large-scale dynamic simulations for three-dimensional problems
A Mixed Method for Axisymmetric Div-Curl Systems
We present a mixed method for a three-dimensional axisymmetric div-curl system reduced to a two-dimensional computational domain via cylindrical coordinates. We show that when the meridian axisymmetric Maxwell problem is approximated by a mixed method using the lowest order Nédélec elements (for the vector variable) and linear elements (for the Lagrange multiplier), one obtains optimal error estimates in certain weighted Sobolev norms. The main ingredient of the analysis is a sequence of projectors in the weighted norms satisfying some commutativity properties
Multigrid methods for Maxwell\u27s equations
In this work we study finite element methods for two-dimensional Maxwell\u27s equations and their solutions by multigrid algorithms. We begin with a brief survey of finite element methods for Maxwell\u27s equations. Then we review the related fundamentals, such as Sobolev spaces, elliptic regularity results, graded meshes, finite element methods for second order problems, and multigrid algorithms. In Chapter 3, we study two types of nonconforming finite element methods on graded meshes for a two-dimensional curl-curl and grad-div problem that appears in electromagnetics. The first method is based on a discretization using weakly continuous P1 vector fields. The second method uses discontinuous P1 vector fields. Optimal convergence rates (up to an arbitrary positive epsilon) in the energy norm and the L2 norm are established for both methods on graded meshes. In Chapter 4, we consider a class of symmetric discontinuous Galerkin methods for a model Poisson problem on graded meshes that share many techniques with the nonconforming methods in Chapter 3. Optimal order error estimates are derived in both the energy norm and the L2 norm. Then we establish the uniform convergence of W-cycle, V-cycle and F-cycle multigrid algorithms for the resulting discrete problems. In Chapter 5, we propose a new numerical approach for two-dimensional Maxwell\u27s equations that is based on the Hodge decomposition for divergence-free vector fields. In this approach, an approximate solution for Maxwell\u27s equations can be obtained by solving standard second order scalar elliptic boundary value problems. We illustrate this new approach by a P1 finite element method. In Chapter 6, we first report numerical results for multigrid algorithms applied to the discretized curl-curl and grad-div problem using nonconforming finite element methods. Then we present multigrid results for Maxwell\u27s equations based on the approach introduced in Chapter 5. All the theoretical results obtained in this dissertation are confirmed by numerical experiments
Electrohydrodynamic Simulations of Capsule Deformation Using a Dual Time-Stepping Lattice Boltzmann Scheme
Capsules are fluid-filled, elastic membranes that serve as a useful model for synthetic and biological membranes. One prominent application of capsules is their use in modeling the response of red blood cells to external forces. These models can be used to study the cell’s material properties and can also assist in the development of diagnostic equipment. In this work we develop a three dimensional model for numerical simulations of red blood cells under the combined influence of hydrodynamic and electrical forces. The red blood cell is modeled as a biconcave-shaped capsule suspended in an ambient fluid domain. Cell deformation occurs due to fluid motion and electrical forces that arise due to differences in the electrical properties between the internal fluid, external fluid, and cell membrane. The electrostatic equations are solved using the immersed interface method. A finite element method is used to compute the membrane’s elastic forces and the membrane’s bending resistance is described by the Helfrich bending energy functional. The membrane forces are coupled to the fluid equations through the immersed boundary method, where the elastic, bending, and electric forces appear as force densities in the Navier-Stokes equations. The fluid equations are solved using a novel dual time-stepping (DTS) lattice Boltzmann method (LBM), which decouples the fluid and capsule discretizations. The computational efficiency of the DTS scheme is studied for capsules in shear flow where it is found that the newly proposed scheme decreases computational time by a factor of 10 when compared to the standard LBM capsule model. The method is then used to study the dynamics of spherical and biconcave capsules in a combined shear flow and DC electric field. For spherical capsules the effect of field strength, shear rate, membrane capacitance, and membrane conductance are studied. For biconcave capsules the effect of the electric field on the tumbling and tank-treading modes of biconcave capsules is discussed
Improved Numerical Methods for Elliptic Problems in Astrophysics
Les ePDEs (elliptic partial differential equations, en anglès) apareixen en una àmplia varietat d'àrees de les matemàtiques, la física i l'enginyeria. Són de particular interès en Astrofísica on apareixen, per exemple, quan es calcula el potencial gravitacional, en la solució de l'equació de Grad-Shafranov per magnetosferes lliures de forces, o d'imposar lligadures de divergència zero en la integració numèrica de les equacions MHD (magnetohydrodynamics, en anglès). En general, les ePDEs s'han de resoldre numèricament, establint una demanda cada vegada més gran d'algoritmes eficients i altament paral·lels per abordar la seua resolució computacional.
El SRJ (scheduled relaxation Jacobi, en anglès) pertany a una prometedora classe de mètodes, atípic per la combinació de senzillesa i eficàcia, que s'ha introduït recentment per resoldre ePDEs lineals de tipus Poisson. És una extensió del mètode iteratiu clàssic de Jacobi utilitzat per resoldre sistemes d'equacions lineals del tipus Au = b. Hereta, d'entre altres, la seua robustesa. La seua metodologia es basa en el càlcul d'uns paràmetres apropiats per a una aproximació multinivell amb l'objectiu de minimitzar el nombre d'iteracions necessàries per a reduir el residual per davall d'una tolerància especificada.
L'eficiència en la reducció del residual augmenta amb el nombre de nivells emprats en l'algoritme. Tanmateix, l'aplicació de la metodologia original per calcular els paràmetres d'estos esquemes SRJ òptims més enllà de 5 nivells és enormement dificultosa. Això és degut fonamentalment a la presència d'un sistema mixt algebraic-diferencial (no lineal) d'equacions el qual es torna cada vegada més rígid a mesura que augmenta el nombre de nivells.
D'una banda, hem trobat una nova metodologia per a l'obtenció dels paràmetres dels esquemes òptims de l'algoritme SRJ que supera les limitacions de la metodologia original i proporciona els paràmetres per a estos esquemes amb un nombre elevat de nivells, fóra bo fins a 15, i per a resolucions de fins a 215 punts per dimensió. Això dóna lloc a factors d'acceleració de diversos centenars respecte del mètode de Jacobi en el cas de resolucions típiques i de milers en alguns casos amb altes resolucions. La major part de l'èxit en la recerca d'estos esquemes òptims amb més de 10 nivells es basa en una reducció analítica de la complexitat del sistema d'equacions abans esmentat. A més, s'estén l'algoritme original per aplicar-lo a certs sistemes d'equacions el·líptiques no lineals.
D'altra banda, en un esquema típic SRJ, s'empra l'anterior conjunt de paràmetres en cicles de M iteracions consecutives fins que s'arriba a la tolerància prescrita. Presentem la forma analítica del conjunt òptim de factors de relaxació per al cas en què tots ells són estrictament diferents, i veiem que l'algoritme resultant és equivalent al mètode no estacionari de Richardson generalitzat, en el que es precondiciona la matriu del sistema d'equacions multiplicant per D = diag(A). El nostre mètode per estimar els pesos té l'avantatge que el càlcul explícit dels valors propis mínim i màxim de la matriu A (o la matriu d'iteració corresponent de l'esquema de Jacobi amb pes subjacent) es substitueix pel càlcul (molt més fàcil) de les freqüències mínima i màxima derivades de l'anàlisi d'estabilitat de von Neumann de l'operador el·líptic continu. Este conjunt de pesos també és l'òptim per al problema general, la qual cosa ens dóna la convergència més ràpida de tots els possibles esquemes SRJ per una estructura de malla donada. Ens referirem a ell com el mètode de Chebyshev-Jacobi. El factor d'amplificació del mètode es pot trobar analíticament i permet l'estimació exacta del nombre d'iteracions necessàries per a assolir la tolerància desitjada. També mostrem que a partir del conjunt de pesos calculats per l'esquema SRJ òptim per a una mida de cicle fix és possible calcular numèricament el valor òptim del factor de relaxació del mètode SOR (successive overrelaxation, en anglès) en alguns casos.
Demostrem amb exemples pràctics, d'aplicació en Astrofísica, que el nostre mètode també funciona molt bé per als problemes de tipus Poisson en els que es fa servir una discretització d'alt ordre de l'operador Laplacià (per exemple, discretitzacions de 9- o 17- punts). Això té molt d'interès, ja que estes discretitzacions no produeixen matrius CO (consistently ordered, en anglès) i, per tant, la teoria de Young no es pot utilitzar per calcular el valor òptim del paràmetre de relaxació òptim de SOR. D'altra banda, els esquemes SRJ òptims deduïts ací són avantatjoses respecte a les implementacions existents per SOR pel que fa a discretitzacions d'alt ordre de l'operador Laplacià en la mesura que no cal recórrer als esquemes multicolors per a la seua execució en paral·lel.
Presentem el mètode de Chebyshev-Jacobi fent servir una implementació purament MPI i una implementació híbrida OpenMP/MPI, ambdues sobre màquines de memòria compartida i de memòria distribuïda. Mostrem el seu rendiment i com escalen. També mostrem com arribar a velocitats de convergència notables amb execucions en paral·lel sobre GPUs quan la resolució d'equacions en derivades parcials el·líptiques amb diferències finites es fa utilitzant de manera conjunta el mètode de Chebyshev-Jacobi i les discretitzacions d'alt ordre.
Finalment, tractar d'aplicar els nostres mètodes més enllà de l'àmbit de l'Astrofísica. En particular, abordem el problema de trobar els modes normals de vibració de l'ull humà. Este problema es pot resoldre amb una variant millorada de la metodologia que ací es presenta. La millora consisteix a estendre el càlcul del conjunt òptim de paràmetres al cas de matrius no definides positives. Les nostres idees sobre com procedir en este camp s'esbossen en el treball futur d'esta tesi.Elliptic partial differential equations appear in a wide variety of areas of mathematics, physics and engineering. They are of particular interest in Astrophysics and appear, e.g., when computing the gravitational potential, in the solution of the Grad-Shafranov equation for force-free magnetospheres, to impose divergence free constraints in the numerical integration of MHD equations or when solving the constraint equations in General Relativity. Typically, elliptic equations must be solved numerically, which sets an ever-growing demand for efficient and highly parallel algorithms to tackle their computational solution.
The Scheduled Relaxation Jacobi is a promising class of methods, atypical for combining simplicity and efficiency, that has been recently introduced for solving linear Poisson-like elliptic equations.
It is an extension of the classical Jacobi iterative method to solve linear systems of equations (Au=b). It also inherits its robustness. Its methodology relies on computing the appropriate parameters of a multilevel approach with the goal of minimizing the number of iterations needed to cut down the residuals below specified tolerances.
The efficiency in the reduction of the residual increases with the number of levels employed in the algorithm. Applying the original methodology to compute the algorithm parameters with more than 5 levels notably hinders obtaining optimal schemes, as the mixed (non-linear) algebraic-differential system of equations from which they result become notably stiff.
On one hand, we present a new methodology for obtaining the parameters of Scheduled Relaxation Jacobi schemes that overcomes the limitations of the original algorithm and provides parameters for these schemes with up to 15 levels and resolutions of up to 215 points per dimension, allowing for acceleration factors larger than several hundreds with respect to the Jacobi method for typical resolutions and, in some high resolution cases, close to 1000. Most of the success in finding these optimal schemes with more than 10 levels is based on an analytic reduction of the complexity of the previously mentioned system of equations. Furthermore, we extend the original algorithm to apply it to certain systems of non-linear elliptic equations.
On the other hand, in a typical Scheduled Relaxation Jacobi scheme, the former set of factors is employed in cycles of M consecutive iterations until a prescribed tolerance is reached. We present the analytic form for the optimal set of relaxation factors for the case in which all of them are strictly different, and find that the resulting algorithm is equivalent to a non-stationary generalized Richardson's method where the matrix of the system of equations is preconditioned multiplying it by D=diag(A). Our method to estimate the weights has the advantage that the explicit computation of the maximum and minimum eigenvalues of the matrix A (or the corresponding iteration matrix of the underlying weighted Jacobi scheme) is replaced by the (much easier) calculation of the maximum and minimum frequencies derived from a von Neumann analysis of the continuous elliptic operator. This set of weights is also the optimal one for the general problem, resulting in the fastest convergence of all possible Scheduled Relaxation Jacobi schemes for a given grid structure. We refer to it as the Chebyshev-Jacobi method. The amplification factor of the method can be found analytically and allows for the exact estimation of the number of iterations needed to achieve a desired tolerance.
We also show that with the set of weights computed for the optimal SRJ scheme for a fixed cycle size it is possible to estimate numerically the optimal value of the relaxation factor in the successive overrelaxation method in some cases.
We demonstrate with practical examples, with application in Astrophysics, that our method also works very well for Poisson-like problems in which a high-order discretization of the Laplacian operator is employed (e.g., a 9- or 17-points discretization). This is of interest since the former discretizations do not yield consistently ordered A matrices and, hence, the theory of Young cannot be used to predict the optimal value of the SOR parameter. Furthermore, the optimal SRJ schemes deduced here are advantageous over existing SOR implementations for high-order discretizations of the Laplacian operator in as much as they do not need to resort to multi-coloring schemes for their parallel implementation.
We present the implementation of the Chebyshev-Jacobi method using a purely MPI implementation, an openMP / MPI hybrid implementation over both shared memory machines and distributed memory machines. They show ideal speed up. We also show how to reach a remarkable speed up when solving elliptic partial differential equations with finite differences thanks to the joint use of the Chebyshev-Jacobi method with high order discretizations and its parallel implementation over GPUs.
Finally, we have tried to apply our methods beyond the realm of Astrophysics with limited success though. Specially, we have addressed the problem of finding normal modes of human eyeballs. This problem is ready for being solved with an improved variant of the methodology here presented. The improvement consists on extending the calculation of the optimal set of parameters to non positive-definite matrices. Our ideas on how to proceed in this field are sketched in the outlook of this thesis
- …