ABSTRACT. In weakly coupled plasmas, the integro-differential Fokker-Planck collision operator describes the dynamical collisional relaxation of the plasma distribution function in the velocity space. The operator is quadratically nonlinear, features exact conservation of mass, momentum, and energy, satisfies the Boltzmann H-theorem, preserves positivity of the distribution function, and features a unique equilibrium solution (i.e., a drifting Maxwellian). Additionally, the collisional velocity space transport coefficients are obtained from integrals of the distribution function (the so-called Rosenbluth potentials). The integral formulation, together with the above-mentioned properties make the numerical solution aspects of this system exceedingly challenging to deal with.

Often, collision times, τcol, are much shorter than dynamical timescales of interest (e.g., in ICF problems, often >104), and the use of implicit methods (i.e., Newton, accelerated Picard [1], etc.) is warranted. Fully implicit methods are advocated here, owing for the need to exactly preserve conservation laws, and to ensure asymptotic convergence to the Maxwellian when appropriate [2]. However, fully implicit methods demand nonlinear iterative solvers, which developing an effective one is challenging in this context owing to the non-locality of the formulation (via the Rosenbluth potentials), leading to dense linear systems.

To effectively deal with the numerical challenges arising from non-locality, we explore a multiscale iterative strategy based on high-order/low-order (HOLO) convergence accelerator scheme [3,4] for the full nonlinear Rosenbluth-Fokker-Planck collision operator. HOLO employ an LO (fluid) moment system to accelerate the convergence of the HO (kinetic) system. In turn, the LO system is closed self-consistently with the HO system. The LO system is comprised of equations for plasma number density, drift velocity, and temperature. These quantities, in turn, inform an LTE approximation for the Rosenbluth potentials plus a perturbation term [of O(Δt/τcol)<< 1] computed from the HO solution. This reformulation shifts the non-local contributions through the Rosenbluth potentials from the HO system to the LO one, which lives in a low-dimensional space and where they can be dealt with efficiently. Numerical experiments in challenging applications in ICF (i.e., shocks and implosions) will demonstrate the enabling capabilities of the HOLO scheme when Δt >> τcol.

[1] Anderson J. Assoc. Comput. Mach., 12, 547-560 (1965). [2] Taitano et al., J. Comp. Phys., 297, 357-380 (2015). [3] Taitano et al., J. Comp. Phys., 284, 737-757 (2015). [4] Chacón et al., J. Comp. Phys., 330, 21-45 (2017).

ABSTRACT. Much progress has been made in recent years demonstrating the viability of fully implicit particle-in-cell (PIC) algorithms [1-4]. When applied to plasma models supporting stiff time scales, fully implicit schemes are known to produce stable simulations at substantially lower resolutions than would be required for explicit discretizations, resulting in large computational gains. In the context of magnetically confined fusion plasmas, a fully implicit discretization scheme can resolve a well-known issue, the so called “cancellation problem”, in which an inexact cancellation of large terms due to different numerical representations introduces large inaccuracies in the simulation results. In this talk, we describe our efforts to apply the work of G. Chen and L. Chacon [4] to a gyrokinetic ion, drift kinetic electron electromagnetic model in the fusion plasma code XGC1 [5]. Since particles are coupled only through the electric and magnetic fields, the nonlinear system of equations resulting from the implicit discretization can be formulated in terms of a residual from the field equations, requiring far less solver memory than a residual in terms of the full particle system. A preconditioned Picard iteration scheme is then applied to this formulation. An effective preconditioner is derived from an electron fluid model, which accurately captures the fast time scale physics in the kinetic model. The convergence rate can be further improved using Anderson acceleration.

[1] G. Chen, L. Chacon, and D.C. Barnes, J. Comput. Phys. 230 (2011) 7018. [2] S. Markidis, G. Lapenta, J. Comput. Phys. 230 (2011) 7037. [3] G. Chen, L. Chacon, C. Liebs, D. Knoll, and W. Taitano, J. Comput. Phys. 258 (2014) 555. [4] G. Chen, L. Chacon, Comput. Phy. Comm. 197, (2015) 73-87. [5] S. Ku, C.S. Chang, and P.H. Diamond, Nucl. Fusion 49 (2009) 115021.

ABSTRACT. Bridging the inherent multi-scale nature of many important problems in plasma physics requires a combination of model reduction and algorithmic innovation. The hybrid model with full-orbit kinetic ions and fluid electrons is a promising approach to describe a wide range of space and laboratory plasmas. In particular, it has been recently demonstrated that the hybrid model is the minimum sufficient model to correctly capture the rate and global evolution of a reconnecting system for arbitrary guide magnetic field [1,2,3], with important consequences for space weather modeling. Here we focus on the quasi-neutral hybrid model using macro-particles to model the kinetic ion species, which avoids the need to solve for the distribution function in 6D (3D-3V) phase space, but is subject to discrete particle noise that scales as sqrt(N) for the number of macro-particles N. The majority of existing algorithms employ explicit time-stepping, and much of the development has focused on the accuracy and stability of the time integration schemes. Key unresolved issues with these approaches involve the numerical instability of cold ion beams moving through the spatial grid due to aliasing errors [4], and the quadratic CFL condition (dt ~ dx^2) associated with fast dispersive Whistler waves. Recent work [5,6] has explored the use of fully implicit methods to step over such timescales in a stable manner, but have not addressed the topic of momentum or energy conservation. Here we present a novel fully implicit hybrid algorithm that features discrete mass, momentum and energy conservation, as well as discrete preservation of the solenoidal condition [7]. The algorithm features sub-cycling and orbit averaging of the ions, with cell-centered finite differences and implicit midpoint time advance. To reduce numerical noise, the algorithm allows arbitrary-order shape functions and conservative smoothing on the gather and scatter operations. We discuss the implementation of the algorithm into a Jacobian-Free Newton Krylov solver framework, preconditioned with a multilevel moment-based (i.e., extended MHD) solver, and then verify it for a number of test problems. These demonstrate the correctness of our implementation, the unique conservation properties, its favorable stability properties, and the effectiveness of our preconditioning approach. References: [1] A. Stanier, W. Daughton, Andrei N. Simakov, L. Chacón, A. Le, H. Karimabadi, J. Ng, and A. Bhattacharjee, Phys. Plasmas 24, 022124 (2017). [2] J. Ng, Y.-M. Huang, A. Hakim, A. Bhattacharjee, A. Stanier, W. Daughton, L. Wang, K. Germaschewski, Phys. Plasmas 22, 112104 (2015). [3] A. Stanier, W. Daughton, L. Chacón, H. Karimabadi, J. Ng, Y.-M. Huang, A. Hakim, and A. Bhattacharjee, Phys. Rev. Lett. 115, 175004 (2015). [4] P. W. Rambo, Journal of Computational Physics, 118, 152-158 (1995). [5] B. Sturdevant, S. E. Parker, Y. Chen and B. B. Hause, Journal of Computational Physics, 316, 519 (2016). [6] J. Cheng, S. Parker, Y. Chen, D. Uzdensky, Journal of Computational Physics, 245, 364 (2013). [7] A. Stanier, L. Chacón, G. Chen, Journal of Computational Physics, 376, 597-616 (2019)

ABSTRACT. Accurate modeling of radiation-hydrodynamics problems is important for many high-energy density physics (HEDP) applications such as astrophysics and Inertial confinement fusion (ICF) research. In the HEDP regime, the material and radiation energies become similar in the magnitude, and thermal radiative transfer (TRT) becomes a dominant energy transfer mechanism.

The stiff coupling between the radiation and material energies via absorption-emission and a large-dimensionality makes the numerical solution of the TRT system a computationally demanding task. To remedy the computational demand of the TRT solver, we have recently developed a moment-based acceleration algorithm (a.k.a., High-Order Low-Order, or HOLO). Our HOLO algorithm employs a discretely-consistent low-dimensional system to accelerate the solution of the stiff collisional physics. Furthermore, in order to properly model a more realistic ICF problems, we have extended algorithm to a Lagrangian radiation-hydrodynamics system.

Within our HOLO algorithm, a majority of numerical stiffness exists in the low-order (LO) system, which resembles advection-diffusion equations. The linear/nonlinear solver robustness for this LO system is a key. Unlike Eulerian radiation-hydrodynamics problems, Lagrangian radiation-hydrodynamics schemes involve mesh-motions, which often cause highly-stretched elements in radiation-driven problems (e.g., ICF implosion simulations). This highly-stretched mesh adds a difficulty in a linear/nonlinear solver convergence. We use a multigrid preconditioned Newton-Krylov method with a nonlinear elimination to solve this LO system. In this talk, we present our recent findings of effective and robust preconditioning strategies and multigrid options for solving stiff Lagrangian radiation-hydrodynamics problems.

ABSTRACT. We develop and study a new framework of multiscale finite element method (MsFEM) for solving the convective-diffusive problems, Helmholtz equations, and heterogeneous elliptic problems. In the proposed framework, which we call the multiscale finite element method with adaptive bubble function enrichment (or MsFEM_bub), we decompose the solution function space into two parts. The first part is the one that can be resolved by the multiscale basis functions and the second part is an unsolved part that is taken care by a set of bubble-shaped functions. These bubble functions are defined similarly to the way to construct to multiscale basis functions. We perform the local-global information exchanges through updated local boundary condition for these bubble function functions. Once the multiscale solution is recovered from the solution of global numerical formulation on coarse grids, which couples these multiscale basis functions, it provides feedback for updating the local boundary conditions on each coarse element. As the approach iterates, the quality of MsFEM_bub solutions gets improved, since these multiscale basis functions with bubble function enrichment are expected to be able to capture the multiscale feature of the approximate solution more accurately. We demonstrate the advantage of the proposed method through some numerical examples for both steady and unsteady convective-diffusive, Helmholtz, and porous media flow benchmark problems.

ABSTRACT. Many problems in porous media science and geophysics comprise interactions of processes, and are typically formulated as a system of coupled PDEs. In most cases these systems are transient and often also non-linear. Developing efficient solvers is a delicate task, since one must to combine suitable schemes for (i) time integration, (ii) linearization, and (iii) (geometric and/or algebraic) multilevel solvers, finally being employed in a (iv) parallel computing environment. In this presentation, we take an application oriented approach, and discuss problems from poroelasticity and density-driven-flow. For these two classes, we describe a unified framework, using linearly-implicit time integration and parallel adaptive multilevel solvers, and provide numerical results.

ABSTRACT. Problems arising in the earth's mantle convection involve finding the solution to Finite Element Stokes systems with large viscosity contrasts. These systems contain localized features which, even with adaptive mesh refinement, result in linear systems that can be on the order of 100+ million unknowns. One common approach while preconditioning the velocity space of these systems is to apply an Algebraic Multigrid (AMG) v-cycle (as is done in the ASPECT software, for example), however, with AMG, robustness can often be difficult with respect to problem size and number of parallel processes. Additionally, we have seen an increase in iteration counts during steps adaptive refinement when using AMG. In contrast, the Geometric Multigrid (GMG) method, by using information about the geometry of the problem, should offer a more robust option, i.e., should have convergence properties independent of the mesh size and be equivalent across processor counts.

Here we present both matrix-based and matrix-free variants of the GMG v-cycle which will work on adaptively refined, distributed meshes, and we will compare it against the current AMG preconditioner used in the ASPECT software. We will demonstrate the robustness of GMG with respect to problem size and show scaling up to 1024 cores, as well as show other advantages of matrix-free computations.

All the computations are performed using the open source, C++, Finite Element library deal.II, which contains functionality for adaptive mesh refinement, multigrid methods, and matrix-free computations.

ABSTRACT. The mathematical characterization of the sound of a musical instrument still follows Schumann’s laws. According to this theory, the resonances of the instrument body, “the formants”, filter the oscillations of the sound generator (e.g., strings) and produce the characteristic “timbre” of an instrument. This is a strong simplification of the actual situation. It applies to a point source and can be easily performed by a loudspeaker, disregarding the three dimensional structure of music instruments. To describe the effect of geometry and material of the instruments, we set up a 3d model and simulate it using the simulation system UG4.

In this work, we present FEM based numerical simulations of a harpsichord. We aim to capture the oscillation behavior of eigenfrequencies of its soundboard. We resolve the complicated geometry by several unstructured 3d grids and take into account the anisotropy of wood. The eigenvalue problem is solved using the PINVIT method with an efficient GMG preconditioner. The latter allows us to resolve the harpsichord with a high resolution grid, which is required to capture fine modes of the simulated eigenfrequencies. To verify our results, we compare them with measurement data obtained from an experimental modal analysis of a real reference harpsichord.

ABSTRACT. Motivated by uncertainty quantification for the neutron diffusion criticality problem, we will study an elliptic eigenvalue problem with coefficients that depend on infinitely many stochastic parameters. The stochasticity in the coefficients causes the eigenvalues and eigenfunctions to also be stochastic, and so our goal is to compute the expectation of the minimal eigenvalue. In practice, to approximate this expectation one must: 1) truncate the stochastic dimension; 2) discretise the eigenvalue problem in space (e.g., by finite elements); and 3) apply a quadrature rule to estimate the expected value.

In this talk, we will present a multilevel Monte Carlo method for approximating the expectation of the minimal eigenvalue, which is based on a hierarchy of finite element meshes and truncation dimensions. To improve the efficiency over Monte Carlo we will also use a quasi-Monte Carlo rule to generate the sampling points. Quasi-Monte Carlo rules are deterministic (or quasi-random) quadrature rules that are well-suited to high-dimensional integration and can converge at a rate of $1/N$, which is faster than the rate of $1/\sqrt{N}$ for Monte Carlo. Also, to make each eigenproblem solve on a given level more efficient, we utilise the two-grid scheme from [Xu \& Zhou 1999] to obtain the eigenvalue on the fine mesh from the coarse eigenvalue (and eigenfunction) with a single linear solve.