Elliptic Solver Research Articles

Reconstructing the Coronal Magnetic Field: The Role of Cross-field Currents in Solution Uniqueness

Abstract We present a new 3D magnetohydrostatic (MHS) direct elliptic solver for extrapolating the coronal magnetic field from photospheric boundary conditions in a manner consistent with an assumed plasma distribution. We use it to study the uniqueness of the reconstructed magnetic field as a function of how significant the plasma forcing is on the force balance of the magnetic field. To this end, we consider an analytic MHS model as ground truth. The model uses two free parameters to decompose the current into two parts: a magnetic-field-aligned component and a cross-field component. We perform a comprehensive study of the 2D parameter space to understand under what conditions the ground truth can be reproduced uniquely. We find that current oriented perpendicular to the magnetic field has a smaller solution space than the same amount of current oriented parallel to the magnetic field, and so MHS regimes with larger proportions of plasma-related forcing may be a promising avenue toward finding unique magnetic field reconstructions.

A fully implicit, scalable, conservative nonlinear relativistic Fokker–Planck 0D-2P solver for runaway electrons

Upon application of a sufficiently strong electric field, electrons break away from thermal equilibrium and approach relativistic speeds. These highly energetic ‘runaway’ electrons (∼ MeV) play a significant role in tokamak disruption physics, and therefore their accurate understanding is essential to develop reliable mitigation strategies. For this purpose, we have developed a fully implicit solver for the 0D-2P (i.e., including two momenta coordinates) relativistic nonlinear Fokker–Planck equation (rFP). As in earlier implicit rFP studies (NORSE, CQL3D), electron–ion interactions are modeled using the Lorentz operator, and synchrotron damping using the Abraham–Lorentz–Dirac reaction term. However, our implementation improves on these earlier studies by (1) ensuring exact conservation properties for electron collisions, (2) strictly preserving positivity, and (3) being scalable algorithmically and in parallel. Key to our proposed approach is an efficient multigrid preconditioner for the linearized rFP equation, a multigrid elliptic solver for the Braams–Karney potentials (Braams and Karney, 1987), and a novel adaptive technique to determine the associated boundary values. We verify the accuracy and efficiency of the proposed scheme with numerical results ranging from small electric-field electrical conductivity measurements to the accurate reproduction of runaway tail dynamics when strong electric fields are applied.

A fourth-order kernel-free boundary integral method for implicitly defined surfaces in three space dimensions

The kernel-free boundary integral (KFBI) method is a finite difference version of the traditional boundary integral method for elliptic and parabolic partial differential equations on complex domains. It evaluates boundary or volume integrals involved in the solution of boundary integral equations (BIEs) by solving equivalent but simple interface problems on regular grids, so that the integral kernel or Green function is never needed or computed. This is the essential difference of the KFBI method from the traditional ones. It takes advantage of the well-conditioning property of discrete BIEs so that the number of Krylov subspace iterations is essentially independent of discretization parameter or system dimension. This paper presents a fourth-order kernel-free boundary integral method for second-order elliptic partial differential equations on complex domains in three space dimensions, whose boundaries are given by implicitly defined surfaces. It represents the domain boundary and discretizes data on it by intersection points of the surface with Cartesian grid lines. The approach has a variety of advantages. The current work solves simple interface problems with corrected 27-point compact finite difference schemes and calculates the discrete equations with a fast Fourier transform based elliptic solver. Numerical examples show that the proposed method is efficient as well as accurate.

A comparison of the shared-memory parallel programming models OpenMP, OpenACC and Kokkos in the context of implicit solvers for high-order FEM

We consider the application of three performance-portable programming models in the context of a high-order spectral element, implicit time-stepping solver for the Navier–Stokes equations. We aim to evaluate whether the use of these models allows code developers to deliver high-performance solvers for computational fluid dynamics simulations that are capable of effectively utilising both many-core CPU and GPU architectures. Using the core elliptic solver for the Navier–Stokes equations as a benchmarking guide, we evaluate the performance of these models on a range of unstructured meshes and give guidelines for the translation of existing codebases and their data structures to these models.

A geometric multigrid library for quadtree/octree AMR grids coupled to MPI-AMRVAC

We present an efficient MPI-parallel geometric multigrid library for quadtree (2D) or octree (3D) grids with adaptive refinement. Cartesian 2D/3D and cylindrical 2D geometries are supported, with second-order discretizations for the elliptic operators. Periodic, Dirichlet, and Neumann boundary conditions can be handled, as well as free-space boundary conditions for 3D Poisson problems, for which we use an FFT-based solver on the coarse grid. Scaling results up to 1792 cores are presented. The library can be used to extend adaptive mesh refinement frameworks with an elliptic solver, which we demonstrate by coupling it to MPI-AMRVAC. Several test cases are presented in which the multigrid routines are used to control the divergence of the magnetic field in magnetohydrodynamic simulations.

Solving hyperbolic-elliptic problems on singular mapped disk-like domains with the method of characteristics and spline finite elements

A common strategy in the numerical solution of partial differential equations is to define a uniform discretization of a tensor-product multi-dimensional logical domain, which is mapped to a physical domain through a given coordinate transformation. By extending this concept to a multi-patch setting, simple and efficient numerical algorithms can be employed on relatively complex geometries. The main drawback of such an approach is the inherent difficulty in dealing with singularities of the coordinate transformation.This work suggests a comprehensive numerical strategy for the common situation of disk-like domains with a singularity at a unique pole, where one edge of the rectangular logical domain collapses to one point of the physical domain (for example, a circle). We present robust numerical methods for the solution of Vlasov-like hyperbolic equations coupled to Poisson-like elliptic equations in such geometries. We describe a semi-Lagrangian advection solver that employs a novel set of coordinates, named pseudo-Cartesian coordinates, to integrate the characteristic equations in the whole domain, including the pole, and a finite element elliptic solver based on globally C1 smooth splines (Toshniwal et al., 2017 [27]). The two solvers are tested both independently and on a coupled model, namely the 2D guiding-center model for magnetized plasmas, equivalent to a vorticity model for incompressible inviscid Euler fluids. The numerical methods presented show high-order convergence in the space discretization parameters, uniformly across the computational domain, without effects of order reduction due to the singularity. Dedicated tests show that the numerical techniques described can be applied straightforwardly also in the presence of point charges (equivalently, point-like vortices), within the context of particle-in-cell methods.

Scaling to the stars – a linearly scaling elliptic solver for p-multigrid

High-order methods gain increased attention in computational fluid dynamics. However, due to the time step restrictions arising from the semi-implicit time stepping for the incompressible case, the potential advantage of these methods depends critically on efficient elliptic solvers. Due to the operation counts of operators scaling with the polynomial degree p times the number of degrees of freedom nDOF, the runtime of the best available multigrid solvers scales with O(p⋅nDOF). This scaling with p significantly lowers the applicability of high-order methods to high orders. While the operators for residual evaluation can be linearized when using static condensation, Schwarz-type smoothers require their inverses on fixed subdomains. No explicit inverse is known in the condensed case and matrix-matrix multiplications scale with p⋅nDOF. This paper derives a matrix-free explicit inverse for the static condensed operator in a cuboidal, Cartesian subdomain. It scales with p3 per element, i.e. nDOF globally, and allows for a linearly scaling additive Schwarz smoother, yielding a p-multigrid cycle with an operation count of O(nDOF). The resulting solver uses fewer than four iterations for all polynomial degrees to reduce the residual by ten orders and has a runtime scaling linearly with nDOF for polynomial degrees at least up to 48. Furthermore the runtime is less than one microsecond per unknown over wide parameter ranges when using one core of a CPU, leading to time-stepping for the incompressible Navier-Stokes equations using as much time for explicitly treated convection terms as for the elliptic solvers.

Parallel Implementation of a PETSc-Based Framework for the General Curvilinear Coastal Ocean Model

The General Curvilinear Coastal Ocean Model (GCCOM) is a 3D curvilinear, structured-mesh, non-hydrostatic, large-eddy simulation model that is capable of running oceanic simulations. GCCOM is an inherently computationally expensive model: it uses an elliptic solver for the dynamic pressure; meter-scale simulations requiring memory footprints on the order of 10 12 cells and terabytes of output data. As a solution for parallel optimization, the Fortran-interfaced Portable–Extensible Toolkit for Scientific Computation (PETSc) library was chosen as a framework to help reduce the complexity of managing the 3D geometry, to improve parallel algorithm design, and to provide a parallelized linear system solver and preconditioner. GCCOM discretizations are based on an Arakawa-C staggered grid, and PETSc DMDA (Data Management for Distributed Arrays) objects were used to provide communication and domain ownership management of the resultant multi-dimensional arrays, while the fully curvilinear Laplacian system for pressure is solved by the PETSc linear solver routines. In this paper, the framework design and architecture are described in detail, and results are presented that demonstrate the multiscale capabilities of the model and the parallel framework to 240 cores over domains of order 10 7 total cells per variable, and the correctness and performance of the multiphysics aspects of the model for a baseline experiment stratified seamount.

Simulation of glioblastoma growth using a 3D multispecies tumor model with mass effect.

In this article, we present a multispecies reaction-advection-diffusion partial differential equation coupled with linear elasticity for modeling tumor growth. The model aims to capture the phenomenological features of glioblastoma multiforme observed in magnetic resonance imaging (MRI) scans. These include enhancing and necrotic tumor structures, brain edema and the so-called "mass effect", a term-of-art that refers to the deformation of brain tissue due to the presence of the tumor. The multispecies model accounts for proliferating, invasive and necrotic tumor cells as well as a simple model for nutrition consumption and tumor-induced brain edema. The coupling of the model with linear elasticity equations with variable coefficients allows us to capture the mechanical deformations due to the tumor growth on surrounding tissues. We present the overall formulation along with a novel operator-splitting scheme with components that include linearly-implicit preconditioned elliptic solvers, and a semi-Lagrangian method for advection. We also present results showing simulated MRI images which highlight the capability of our method to capture the overall structure of glioblastomas in MRIs.

Fluctuating hydrodynamics of electrolytes at electroneutral scales

At mesoscopic scales electrolyte solutions are modeled by the fluctuating generalized Poisson-Nernst-Planck (PNP) equations [J.-P. P\'eraud et al., Phys. Rev. F, 1(7):074103, 2016]. However, at length and time scales larger than the Debye scales, electrolytes are effectively electroneutral, and the charged-fluid PNP equations become too stiff to solve numerically. Here we formulate the isothermal incompressible equations of fluctuating hydrodynamics for reactive multispecies mixtures involving charged species in the electroneutral limit, and design a numerical algorithm to solve these equations. Our model does not assume a dilute electrolyte solution but rather treats all species on an equal footing, accounting for cross-diffusion and non-ideality using Maxwell-Stefan theory. By enforcing local electroneutrality as a constraint, we obtain an elliptic equation for the electric potential that replaces the Poisson equation in the fluctuating PNP equations. We develop a second-order midpoint predictor-corrector algorithm to solve either the charged-fluid or electroneutral equations with only a change of the elliptic solver. We use the electroneutral algorithm to study a gravitational fingering instability, triggered by thermal fluctuations, at an interface where an acid and base react to neutralize each other. Our results demonstrate that, because the four ions diffuse with very different coefficients, one must treat each ion as an individual species, and cannot treat the acid, base, and salt as neutral species. This emphasizes the differences between electrodiffusion and classical Fickian diffusion, even at electroneutral scales.

Non-oscillatory forward-in-time integrators for viscous incompressible flows past a sphere

A non-oscillatory forward-in-time (NFT) integrator is developed to provide solutions of the Navier–Stokes equations for incompressible flows. Simulations of flows past a sphere are chosen as a benchmark representative of a class of engineering flows past obstacles. The methodology is further extended to moderate Reynolds number, stably stratified flows under gravity, for Froude numbers that typify the characteristic regimes of natural flows past distinct isolated features of topography in weather and climate models. The key elements of the proposed method consist of the Multidimensional Positive Definite Advection Transport Algorithm (MPDATA) and a robust non-symmetric Krylov-subspace elliptic solver. The solutions employ a finite volume spatial discretisation on unstructured and hybrid meshes and benefit from a collocated arrangement of all flow variables while being inherently stable. The development includes the implementation of viscous terms with the detached-eddy simulation (DES) approach employed for turbulent flows. Results demonstrate that the proposed methodology enables direct comparisons of the numerical solutions with corresponding laboratory studies of viscous and stratified flows while illustrating accuracy, robustness and flexibility of the NFT schemes. The presented simulations also offer a better insight into stably stratified flows past a sphere.

FVM 1.0: a nonhydrostatic finite-volume dynamical core for the IFS

Abstract. We present a nonhydrostatic finite-volume global atmospheric model formulation for numerical weather prediction with the Integrated Forecasting System (IFS) at ECMWF and compare it to the established operational spectral-transform formulation. The novel Finite-Volume Module of the IFS (henceforth IFS-FVM) integrates the fully compressible equations using semi-implicit time stepping and non-oscillatory forward-in-time (NFT) Eulerian advection, whereas the spectral-transform IFS solves the hydrostatic primitive equations (optionally the fully compressible equations) using a semi-implicit semi-Lagrangian scheme. The IFS-FVM complements the spectral-transform counterpart by means of the finite-volume discretization with a local low-volume communication footprint, fully conservative and monotone advective transport, all-scale deep-atmosphere fully compressible equations in a generalized height-based vertical coordinate, and flexible horizontal meshes. Nevertheless, both the finite-volume and spectral-transform formulations can share the same quasi-uniform horizontal grid with co-located arrangement of variables, geospherical longitude–latitude coordinates, and physics parameterizations, thereby facilitating their comparison, coexistence, and combination in the IFS. We highlight the advanced semi-implicit NFT finite-volume integration of the fully compressible equations of IFS-FVM considering comprehensive moist-precipitating dynamics with coupling to the IFS cloud parameterization by means of a generic interface. These developments – including a new horizontal–vertical split NFT MPDATA advective transport scheme, variable time stepping, effective preconditioning of the elliptic Helmholtz solver in the semi-implicit scheme, and a computationally efficient implementation of the median-dual finite-volume approach – provide a basis for the efficacy of IFS-FVM and its application in global numerical weather prediction. Here, numerical experiments focus on relevant dry and moist-precipitating baroclinic instability at various resolutions. We show that the presented semi-implicit NFT finite-volume integration scheme on co-located meshes of IFS-FVM can provide highly competitive solution quality and computational performance to the proven semi-implicit semi-Lagrangian integration scheme of the spectral-transform IFS.

Maxwell\u2013Stefan model of multicomponent ion transport inside a monolayer Nafion membrane for intensified chlor-alkali electrolysis

A mathematical model based on a generalized Maxwell–Stefan equation has been developed to describe multicomponent ion and water transport inside a cation-exchange membrane. This model has been validated using experimental data and has been used to predict concentration profiles, membrane potential drop, and transport numbers of ions and water for the chlor-alkali process at increased current densities. Several improvements have been made to previously developed Maxwell–Stefan models. In our model, the generalized Maxwell–Stefan equation is written in terms of concentration instead of mole fraction and the fixed group (membrane) concentration is assumed to be constant. We have adapted the Augmented matrix method using the built-in partial differential equation parabolic elliptic (pdepe) solver in Matlab®, and both the concentration and the electrical potential gradients have been solved simultaneously. The boundary conditions are determined with the Donnan equilibrium at the membrane–solution interface. We have also employed semi-empirical correlations to define the Maxwell–Stefan diffusivities inside the membrane. For the bulk diffusivities, we applied the correlations for the concentrated solution instead of the values at infinite dilution. With the diffusivities presented in this work, the model shows a better fit to the experimental data than with previously reported fitted diffusivities. Prediction of the sodium transport number and water transport number is generally good, whereas the deviations with regard to membrane potential might also be related to issues with the experimental data. The model predicts an increase in both sodium and water transport numbers at increased current density operation of chlor-alkali production.Graphical abstract

Toward computing gravitational initial data without elliptic solvers

Two new methods have been proposed for solving the gravitational constraints without using elliptic solvers by formulating them as either an algebraic–hyperbolic or parabolic–hyperbolic system. Here, we compare these two methods and present a unified computational infrastructure for their implementation as numerical evolution codes. An important potential application of these methods is the prescription of initial data for the simulation of black holes. This paper is meant to support progress and activity in that direction.

Cost reduction of stochastic Galerkin method by adaptive identification of significant polynomial chaos bases for elliptic equations

One widely used and computationally efficient method for uncertainty quantification using spectral stochastic finite element is the stochastic Galerkin method. Here the solution is represented in polynomial chaos expansion, and the residual of the discretized governing equation is projected on the polynomial chaos bases. This results in a system of deterministic algebraic equations with the polynomials chaos coefficients as unknown. However, one impediment for its large scale applications is the curse of dimensionality, that is, the exponential growth of the number of polynomial chaos bases with the stochastic dimensionality and degree of expansion. Here, for a stochastic elliptic problem, an adaptive selection of polynomial chaos bases is proposed. Accordingly, during the first few iterations in the preconditioned conjugate gradient method for solving the system of linear algebraic equations, the chaos bases with maximal contribution –in an appropriately defined metric – to the solution are first identified. Subsequently, only these bases are retained for further iterations until convergence is achieved. Using numerical studies a three times cost saving over the existing method is observed. Furthermore, for enhancing the computational cost gain, the stochastic Galerkin method is reformulated as a generalized Sylvester equation. This step allowed efficient usage of the sparsity of moments of product of polynomial chaos bases. Through numerical studies on problems with large stochastic dimensionality, an additional cost saving of up to one order of magnitude –twenty times –is observed. This amounts to sixty times speedup over the existing method, when adaptive selection and generalized Sylvester equation formulation are used together. The proposed methodology can be easily incorporated in an existing standard stochastic Galerkin method solver for elliptic problems.

Parallelization of elliptic solver for solving 1D Boussinesq model

In this paper, a parallel implementation of an elliptic solver in solving 1D Boussinesq model is presented. Numerical solution of Boussinesq model is obtained by implementing a staggered grid scheme to continuity, momentum, and elliptic equation of Boussinesq model. Tridiagonal system emerging from numerical scheme of elliptic equation is solved by cyclic reduction algorithm. The parallel implementation of cyclic reduction is executed on multicore processors with shared memory architectures using OpenMP. To measure the performance of parallel program, large number of grids is varied from 28 to 214. Two test cases of numerical experiment, i.e. propagation of solitary and standing wave, are proposed to evaluate the parallel program. The numerical results are verified with analytical solution of solitary and standing wave. The best speedup of solitary and standing wave test cases is about 2.07 with 214 of grids and 1.86 with 213 of grids, respectively, which are executed by using 8 threads. Moreover, the best efficiency of parallel program is 76.2% and 73.5% for solitary and standing wave test cases, respectively.

Framework of Multiblock Structured Grid Smoothing Algorithms for Aircraft Configurations

This paper presents the development of a new framework of multiblock structured elliptic grid smoothing for full-aircraft configurations. The framework is fully automated and applies different solvers on blocks based on their properties. Blocks with a solid wall are initialized by algebraic orthogonal grids and solved with a line SOR implicit elliptic solver. The remaining blocks are initialized by transfinite interpolation grids and solved with an SOR explicit elliptic solver. Both solvers apply spacing and curvature source terms. This combination results in orthogonal grids close to walls while imposing elliptic spacing and curvature control using variable relaxation factors in the field, including across the block boundaries. The orthogonality close to the solid walls is improved by using the properties of the initial algebraic orthogonal grids. The code is parallelized using Message Passing Interface for acceleration. The performance of the framework is demonstrated in a first step on 2D problems: a canonical test case and a multi-element high-lift airfoil. In a second step, the performance is analyzed on a full three-dimensional aircraft configurations (wing/body/engine) case at Bombardier Aerospace. The paper examines the effect of the proposed approach through various assessments of grid metrics and flow solver behaviors. In particular, the study highlights the improvements in reducing spurious drag through a far-field drag analysis.

Iterative Solvers for Elliptic Problems with Arbitrary Anisotropy Strengths

In this paper new asymptotic-preserving schemes are introduced for the iterative resolution of anisotropic elliptic equations arising in magnetized plasma simulations. These methods overcome the resolution of the saddle point problems systematically implemented in precedent realizations. This allows us to easily derive iterative solvers for these asymptotic-preserving schemes. It brings a leap forward in the computational efficiency of the method for three-dimensional problems. This is conclusively outlined thanks to three-dimensional serial computations carrying out tens of millions of unknowns. The gains are substantial in terms of memory as well as computational requirements compared to sparse direct solvers, which represent the only alternative successfully operated so far.

Building blocks for a leading edge high‐order flow solver

AbstractAn important trend in Computational Fluid Dynamics is towards high‐order methods, as they offer a substantially lower discretization error for the same number of degrees of freedom (DOF). Examples are the Spectral‐Element Methods (SEM) and Discontinuous Galerkin (DG) methods. Unfortunately, with most implementations the work load of such solvers increases drastically with the number of DOF, for example, when increasing the polynomial degree of the approximation. This issue gets particular pressing for elliptic solvers which are a vital building block in the time‐stepping of the incompressible Navier‐Stokes equations, resulting from pressure projection methods or implicit treatment of viscous terms. So far, this drastic increase of resources has hampered the use of SEM for higher polynomial degrees, such as 16 or more.The present contribution is located at this particular “Frontier of CFD” and proposes an SEM with linear scaling in the number of degrees of freedom. It is achieved independent of the polynomial degree, independent of the aspect ratio of elements, and involves constant iteration count when increasing the number of elements. The method is based on combining static condensation with block‐Jacobi preconditioning and iterative substructuring. The latter two lead to a constant iteration count, while the former warrants linear operator complexity.In the presentation, the scheme is described in detail and applied to the construction of a Helmholtz solver. This solver is extremely fast and able to solve the Helmholtz equation with 1010 unknowns on 240 cores in acceptable time. It thus enables competitive high‐order SEM simulations even on small clusters and in this way expands the frontier of CFD towards highly accurate results requiring comparatively modest resources. (© 2017 Wiley‐VCH Verlag GmbH & Co. KGaA, Weinheim)

Array-based, parallel hierarchical mesh refinement algorithms for unstructured meshes

In this paper, we describe an array-based hierarchical mesh refinement capability through uniform refinement of unstructured meshes for efficient solution of PDE’s using finite element methods and multigrid solvers. A multi-degree, multi-dimensional and multi-level framework is designed to generate the nested hierarchies from an initial coarse mesh that can be used for a variety of purposes such as in multigrid solvers/preconditioners, to do solution convergence and verification studies and to improve overall parallel efficiency by decreasing I/O bandwidth requirements (by loading smaller meshes and in-memory refinement). We also describe a high-order boundary reconstruction capability that can be used to project the new points after refinement using high-order approximations instead of linear projection in order to minimize and provide more control on geometrical errors introduced by curved boundaries.The capability is developed under the parallel unstructured mesh framework “Mesh Oriented dAtaBase” (MOAB Tautges et al. (2004)). We describe the underlying data structures and algorithms to generate such hierarchies in parallel and present numerical results for computational efficiency and effect on mesh quality. We also present results to demonstrate the applicability of the developed capability to study convergence properties of different point projection schemes for various mesh hierarchies and to a multigrid finite-element solver for elliptic problems.