Wavelet-based parallel dynamic mesh adaptation for magnetohydrodynamics in the AMROC framework

TL;DR

This paper introduces a wavelet-based parallel mesh adaptation technique for magnetohydrodynamics simulations within the AMROC framework, enhancing computational efficiency and precision for space physics applications.

Contribution

It presents a novel wavelet-based mesh adaptation method integrated into a parallel MHD solver in the AMROC framework, improving efficiency and accuracy.

Findings

Enhanced mesh adaptation criteria using wavelet techniques

Improved computational efficiency in MHD simulations

Effective control of simulation precision with dynamic meshes

Abstract

Computational magnetohydrodynamics (MHD) for space physics has become an essential area in understanding the multiscale dynamics of geophysical and astrophysical plasma processes, partially motivated by the lack of space data. Full MHD simulations are typically very demanding and may require substantial computational efforts. In particular, computational space-weather forecasting is an essential long-term goal in this area, motivated for instance by the needs of modern satellite communication technology. We present a new feature of a recently developed compressible two- and three-dimensional MHD solver, which has been successfully implemented into the parallel AMROC (Adaptive Mesh Refinement in Object-oriented C++) framework with improvements concerning the mesh adaptation criteria based on wavelet techniques. The developments are related to computational efficiency while controlling the precision using dynamically adapted meshes in space-time in a fully parallel context.

Tables

Table 1. Table 1 : MR and SG: accuracy and number of cells.

Method	$ϵ$	$L_{1,\mathrm{AMR}}$	# of cells
SG	$0.20$	$2.30$	$411,392$
MR	$0.05$	$2.15$	$333,064$

Table 2. Table 2 : Comparison of the CPU time (s) for the two-dimensional magnetic cloudy test case corresponding to a refined mesh 1 , 024 2 1 superscript 024 2 1,024^{2} , ϵ = 0.01 italic-ϵ 0.01 \epsilon=0.01 , η = 0.80 𝜂 0.80 \eta=0.80 , and base mesh 16 2 superscript 16 2 16^{2} .

	Number of Processors
$2$D	$1$	$2$	$4$	$8$	$16$
Carmen-MHD	$17,520$	-	-	-	-
AMROC	$1,889$	$1021$	$635$	$478$	$447$

Table 3. Table 3 : AMROC 3 3 3 D adaptive computations of the magnetic cloudy test case corresponding to ϵ = 0.025 italic-ϵ 0.025 \epsilon=0.025 , η = 0.80 𝜂 0.80 \eta=0.80 , and base mesh 32 2 superscript 32 2 32^{2} .

	Corresponding refined mesh
	${256}^3$	${512}^3$	$1,{024}^3$
# cells	$5,917,288$	$27,444,592$	$122,076,544$
% cells	$35.3$	$20.4$	$11.4$
# blocks	$148$	$894$	$2,502$
CPU time (s)	$1,468$	$11,350$	$117,731$

Equations

$$\frac{\partial \rho }{\partial t}+\nabla \cdot \left(\rho \mathbf{u}\right)$$

$$\frac{\partial \left(\rho \mathbf{u}\right)}{\partial t}+\nabla \cdot \left[\rho \mathbf{u}{\mathbf{u}}^{\mathbf{T}}+\mathbf{I}\left(p+\frac{B.B}{2}\right)-\mathbf{B}{\mathbf{B}}^{\mathbf{T}}\right]$$

$$\frac{\partial \mathbf{B}}{\partial t}+\nabla \cdot \left[\mathbf{u}{\mathbf{B}}^{\mathbf{T}}-\mathbf{B}{\mathbf{u}}^{\mathbf{T}}+\psi \mathbf{I}\right]$$

$$\frac{\partial E}{\partial t}+\nabla \cdot \left[\left(E+p+\frac{B.B}{2}\right)\mathbf{u}-\mathbf{B}\left(u.B\right)\right]$$

$$\frac{\partial \psi }{\partial t}+c_h^2\nabla \cdot \mathbf{B}$$

$$E=\frac{p}{\gamma -1}+\rho \frac{\mathbf{u}\cdot \mathbf{u}}{2}+\frac{\mathbf{B}\cdot \mathbf{B}}{2},$$

$${\bar{\mathbf{Q}}}^{\ell +1}\stackrel{projection}{\underset{prediction}{\rightleftharpoons }}{\bar{\mathbf{Q}}}_{\text{MR}}^{\ell +1}=\{{\bar{\mathbf{Q}}}^{\ell }\}\cup \{d^{\ell }\},$$

$${\bar{\mathbf{Q}}}_i^{\ell }=P_{\ell +1\to \ell }\left({\bar{\mathbf{Q}}}_{2i}^{\ell +1},{\bar{\mathbf{Q}}}_{2i+1}^{\ell +1}\right)=\frac{1}{2}\left({\bar{\mathbf{Q}}}_{2i}^{\ell +1}+{\bar{\mathbf{Q}}}_{2i+1}^{\ell +1}\right),$$

$$Projection:P^{i}_{\ell+1\rightarrow\ell}:{\bar{\mathbf{Q}}}^{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\ell+1}}\rightarrow{\bar{\mathbf{Q}}}^{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\ell}}$$

$$Prediction:P^{i}_{\ell\rightarrow\ell+1}:{\bf Q}^{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\ell}}\rightarrow\tilde{{\bf Q}}^{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\ell+1}}$$

$${\tilde{\mathbf{Q}}}_{2i}^{\ell +1}=P_{\ell \to \ell +1}^0\left({\bar{\mathbf{Q}}}_{i-1}^{\ell },{\bar{\mathbf{Q}}}_i^{\ell },{\bar{\mathbf{Q}}}_{i+1}^{\ell }\right)={\bar{\mathbf{Q}}}_i^{\ell }-\frac{1}{8}({\bar{\mathbf{Q}}}_{i+1}^{\ell }-{\bar{\mathbf{Q}}}_{i-1}^{\ell })$$

$${\tilde{\mathbf{Q}}}_{2i+1}^{\ell +1}=P_{\ell \to \ell +1}^1\left({\bar{\mathbf{Q}}}_{i-1}^{\ell },{\bar{\mathbf{Q}}}_i^{\ell },{\bar{\mathbf{Q}}}_{i+1}^{\ell }\right)={\bar{\mathbf{Q}}}_i^{\ell }+\frac{1}{8}({\bar{\mathbf{Q}}}_{i+1}^{\ell }-{\bar{\mathbf{Q}}}_{i-1}^{\ell }),$$

$$d_i^{\ell }={\bar{\mathbf{Q}}}_i^{\ell }-{\tilde{\mathbf{Q}}}_i^{\ell }.$$

$${ϵ}^{\ell }=\frac{ϵ}{|\Omega |}{2}^{d\left(\ell -\mathbb{L}\right)},1\le \ell \le \mathbb{L},$$

$$|{\rho }_{j+1,k}-{\rho }_{j,k}|>{ϵ}^{\rho },|{\rho }_{j,k+1}-{\rho }_{j,k}|>{ϵ}^{\rho },|{\rho }_{j+1,k+1}-{\rho }_{j,k}|>{ϵ}^{\rho },$$

Full text

††journal: Computers & Fluids

1 Introduction

Space weather forecasting is concerned with predicting disturbances in the Earth upper atmosphere and the magnetic field that can have dramatic consequences for modern technology, especially satellites and communication electronics. A combination of theoretical studies, observations, numerical simulations including data assimilation, and forecasting is the key for achieving success in this strategic area [1]. For numerical predictions of space plasma physics, the considerable range of dynamically active space-time scales is a major obstacle. Accordingly, fast, robust, and efficient numerical models that merge physics-based, accurate simulation with timely observations are of fundamental importance. A particularly successful computational approach is the magnetohydrodynamic model [2]. Roughly speaking, the magnetohydrodynamic model consists of a system of eight nonlinear partial differential equations describing the dynamics of a compressible, inviscid, and perfectly electrically conducting fluid interacting with a magnetic field, combining thus the Euler equations of hydrodynamics with the Maxwell equations of electrodynamics. The latter yields an evolution equation for the magnetic field, the so-called induction equation, and an incompressibility constraint of the magnetic field using Gauss’s law [3]. Numerical approaches to solve these systems are computationally costly and a full mesh approach to compute them is prohibitive in most cases of interest in space physics [2]. Therefore, adaptive techniques have been combined with on-the-fly mesh refinement [2]. More recently, multiresolution techniques have also been used [4, 5]. In such approaches, the mesh is refined locally just in the regions where structures or discontinuities are or can be present in short integration time [6].

Adaptive techniques reduce the computing time significantly while preserving the high accuracy of the numerical solutions. In a previous article [7] we presented a comparison of these mesh adaptation techniques, which we have recently extended in [8]. In these works, we used the generic open-source framework for patch-structured adaptive mesh refinement AMROC [9, 10]. A core finding of these publications is that the multiresolution approach is mathematically more rigorous and leads to a more faithful mesh adaptation; yet, the patch-based adaptation approach – thanks to advantageous data structures – reduces the overall computation time drastically. Hence, in this work we present the multiresolution approach in our MHD solver [11] within the AMROC framework.

The organisation of the remainder of the manuscript is as follows. In Section 2, we contextualise the space environment of the Earth and introduce the primary phenomena of interest for simulation test cases. In Section 3 we describe the governing equations, their discretisation using finite volumes, and outline the main ideas of the patch-based adaptive mesh refinement (AMR) approach briefly as implemented in AMROC. In Section 4 we develop numerical experiments and discuss the accuracy and efficiency of our implementation. Finally, in Section 5 we draw some conclusions.

2 Earth space environment

The Sun is the source of several phenomena that affect the sidereal bodies, even the human artefacts, such as interplanetary probes, existing in the heliosphere [12]. The three main solar agents are electromagnetic radiation, high energetic corpuscular radiation, and magnetised plasma structures evolving in the solar wind.

As the Earth has an atmosphere and an intrinsic magnetic field, a particular situation occurs; the solar radiation ionises the upper atmosphere of the Earth, located above an altitude of $70~{}\mathrm{km}$. Also, a magnetic field of inner origin imprecates the whole atmosphere. This conjunction of ionisation and magnetic field creates a kind of shield to the solar plasma displacement, i.e., the fully ionised and magnetised solar wind plasma cannot mix with the terrestrial plasma, which establishes a geomagnetic field domain surrounding the planet, named magnetosphere. As a consequence, the expanding solar plasma deviates from its original direction involving the Earth’s domain [13].

Incident upon this obstacle, the solar plasma wind, which moves with supersonic speed, creates a bow shock involving the terrestrial domain. Immediately after the shock the flux presents a subsonic speed, thermalised particles, and an intensified interplanetary magnetic field, characterising a region designated as magnetosheath. The plasma of this region compresses the region dominated by the geomagnetic field in a process that defines an interface between the two physical media, the magnetopause. As a manifestation of a tangential discontinuity, this region is a surface of total pressure equilibrium between the solar wind-magnetosheath plasma and the geomagnetic field confined in the magnetosphere. Inside the magnetosphere, several processes establish distinct regions of plasma, energetic particle distributions, and a sophisticated building of electrical current systems. All the features of this real scenario are consequences of the electrodynamical interaction between the incident solar wind and the Earth atmosphere [14].

An enormous amount of studies provided by experimental, observational research, such as the one obtained by satellite data-set analyses, and theoretical approaches, such as the propositions of phenomena from the magnetohydrodynamic formalism, has contributed to an in-depth interpretation of the space environment [15, 16]. The efforts of numerical simulations can significantly help the scientific concepts under development and provide more realistic process descriptions. The comprehensive view and understanding of important space plasma processes and their geoeffective events depend on much more realistic performances of the currently available magnetohydrodynamic models. However, to be utilised in emerging space weather programs, competitive-in-time codes are increasingly demanded nowadays.

In this work, we aim at developing an up-to-date numerical tool for modelling magnetohydrodynamics. We address two fundamental cases, prominent in the investigation of applied space sciences. The first case is the Orszag–Tang vortex, originally introduced in [17] for incompressible MHD, a canonical model problem for testing the transition to supersonic two-dimensional MHD turbulence. The second is the magnetic shock cloud [15, 18], a common occurrence evolving superimposed with the solar wind through the interplanetary medium.

2.1 GLM-MHD model

In numerical simulations of the ideal magnetohydrodynamics equations, including the divergence-free constraint of the magnetic field, typical methodologies consider an additional correction term to facilitate the enforcement of this physical property. There are different kinds of well-known methods in order to minimise this effect, as described in [19] and references therein. In the context of this study, having in mind the application of the multiresolution method using explicit time integration, we adopt the approach proposed in [20], with the non-dimensional adjustment added by Mignone and Tzeferacos in [21]. Namely, the divergence-free constraint is treated by the introduction of a new variable $\psi$ and a corresponding balance equation is added to the ideal MHD equations. This process leads to the well-known Generalised Lagrangian Multiplier (GLM) hyperbolic conservation system

[TABLE]

where $\rho$ represents density, $p$ is the pressure, ${\bf u}$ is the fluid velocity vector, ${\bf B}$ is the magnetic field vector, ${\bf I}$ is the identity tensor of order $2$, and the superscript symbol $T$ indicates the transposed matrix. The parameter $c_{h}$ is defined as $c_{h}=\max\left[\nu\frac{\Delta h}{\Delta t},\,\max\left(|u_{i}|\pm c_{f}\right)\right],$ where $\Delta h$ is the minimal value of the mesh sizes in each direction, $\nu$ the Courant number, $u_{i}$ is the velocity of the $i$-th component, and $c_{f}$ is the fast magneto-acoustic wave of the MHD model. The $c_{p}$ value is defined in terms of the parameter $\alpha_{p}=\Delta h\,\dfrac{c_{h}}{c_{p}^{2}}$, where $\alpha_{p}\in[0,1]$, as described in [21]. The total energy density $E$ is given by the constitutive law

[TABLE]

in which $\gamma$ is the adiabatic constant with $\gamma>1$. Moreover, the above MHD system is completed by suitable initial and boundary conditions as presented in the numerical experiments section. These equations have been rewritten in non-dimensional form such that the magnetic permeability yields the identity, i.e., $\mu=1$. We also consider the divergence control parameter $\alpha_{p}=0.4$, and $\psi\equiv 0$ in the initial conditions.

In near Earth space, the governing equations of MHD modelling can in particular develop discontinuities, i.e., shocks and contact waves. Therefore, we use finite volume shock-capturing methods that are constructed to properly handle this behaviour in a robust and oscillation-free way, as discussed in detail in [22]. Moreover, we are interested in studying the development of MHD instabilities that are very local and can present complex local multiscale behaviour. For this reason, it can be undoubtedly beneficial to have an economical, accurate, and efficient mesh representation of these features.

3 Numerical discretisation

As reference discretisation for these equations in the conservation form, the numerical solution is represented by the quantity vector $\mathbf{Q}$ of the approximated cell averages on a uniform mesh of the computation domain. For space discretisation, a finite volume method is chosen, which results in a system of ordinary differential equations with a vector of numerical flux function differences with respect to each cell. In all numerical schemes throughout this paper enhanced numerical flux functions with comparable second-order-accurate reconstruction and flux limiting are used. For time integration, we adopt an explicit second order Runge–Kutta scheme.

The adaptive mesh refinement method (AMR) [23, 24, 25] follows a patch-oriented refinement approach, where non-overlapping rectangular submeshes $G_{\ell,m}$ define the domain $G_{\ell}:=\bigcup_{m=1}^{M_{\ell}}G_{\ell,m}$ of an entire level $\ell=0,\dots,L$. As the construction of refinement proceeds recursively, a hierarchy of submeshes successively contained within the next coarser level domain is created. The recursive nature of the algorithm only allows the addition of one new level in each refinement operation. The patch-based approach does not require special coarsening operations; submeshes are simply removed from the hierarchy. The coarsest possible resolution is thereby restricted to the level [math] mesh. Typically, it is assumed that all mesh widths on level $\ell$ are $r_{\ell}$-times finer than on the level $\ell-1$, i.e., $\Delta t_{\ell}=\Delta t_{\ell-1}/r_{\ell}$ and $\Delta x_{n,\ell}=\Delta x_{n,\ell-1}/r_{\ell}$, with $r_{\ell}\in\mathbb{N},r_{\ell}\geq 2$ for $\ell>0$ and $r_{0}=1$. This ensures that a time-explicit finite volume scheme remains stable under a CFL-type condition on all levels of the hierarchy. In our MR implementations we always use $r_{\ell}=2$ here.

The numerical update is applied on the level $\ell$ by calling a single-mesh routine implementing the finite volume scheme in a loop over all the submeshes $G_{\ell,m}$. The regularity of the input data allows a straightforward implementation of the scheme and furthermore permits optimisation to take advantage of high-level caches, pipelining, etc. New refinement meshes are initialised by interpolating the vector of conservative quantities ${\bf Q}$ from the next coarser level. However, data in an already refined cell are copied directly from the previous refinement patches. Ghost cells around each patch are used to decouple the submeshes computationally. Ghost cells outside of the root domain $G_{0}$ are used to implement physical boundary conditions. Ghost cells in $G_{\ell}$ have a unique interior cell analogue and are set by copying the data value from the patch where the interior cell is contained (synchronisation). For $\ell>0$, internal boundaries can also be used. If recursive time step refinement is employed, ghost cells at the internal refinement boundaries on the level $\ell$ are set by time-space interpolation from the two previously calculated time steps of level $\ell-1$. Otherwise, spatial interpolation from the level $\ell-1$ is sufficient.

One feature of the AMR algorithm is that refinement patches overlay coarser mesh data structures instead of being embedded, again preventing data fragmentation. Values of cells covered by finer submeshes are subsequently overwritten by averaged fine mesh values, which, in general, would lead to a loss of conservation on the coarser mesh. A remedy to this problem is to replace the coarse mesh numerical fluxes at refinement boundaries with the sum of fine mesh fluxes along the corresponding coarse cell boundary.

3.1 MR refinement indicator

The principle of MR methods is the transformation of the cell averages given by the finite volume discretisation into a multiscale representation. A detailed review can be found in [26, 27] and references therein.

We consider a discrete solution of the discretisation as initial cell average data $\bar{\mathbf{Q}}^{\ell+1}$ at level $\ell+1$. Then for instance, in one decomposition level a two-level MR transformation can be written as follows,

[TABLE]

where $d^{\ell}$ contains the information between the two consecutive levels $\ell$ and $\ell+1$; and $\bar{\mathbf{Q}}^{\ell}$ stores a smooth version of the original numerical solution $\bar{\mathbf{Q}}^{\ell+1}$.

These ideas are a natural extension of the work [28, 29]. The numerical solution at the finest resolution level is transformed into a set of coarser scale approximations plus a series of prediction errors corresponding to wavelet coefficients. These coefficients describe the difference between subsequent resolutions. The main principle is then to use the decay of the wavelet coefficients to estimate the local regularity of the solution [30, 31]. In regions where the solution is smooth these coefficients are small, while they have large magnitude in regions of steep gradients or discontinuities.

In order to perform the MR method, some operations for projection and prediction are required. For the MR scheme with finite volumes, where the cell values are local averages, a coarser cell $\Omega^{\ell}_{i}$ has its value estimated using the smaller scale values and a unique projection operator $P_{\ell+1\rightarrow\ell}\,:\,\bar{\mathbf{Q}}^{\ell+1}\,\mapsto\,\bar{\mathbf{Q}}^{\ell}$. In this scheme, the projection operator to obtain the solution on a coarser cell is given by the average value of its children. For the one-dimensional case (cf. Fig. 1, top), the projection is performed by

[TABLE]

where $\bar{\mathbf{Q}}^{\ell}_{i}$ is the average value of the cell $\Omega^{\ell}_{i}$.

The prediction operators are used to perform the opposite path of the projection operators and allow to obtain the values of the finer cells using the values of the coarser ones (cf. Fig. 1, bottom). For each child cell $i$ to be predicted, there is a different prediction operator represented by $P^{i}_{\ell\rightarrow\ell+1}\,:\,\bar{\mathbf{Q}}^{\ell}\,\mapsto\,\bar{\mathbf{Q}}^{\ell+1}$ for the one-dimensional case. These operators yield a non-unique approximation of $\bar{\mathbf{Q}}^{\ell+1}_{i}$ by interpolation.

We use polynomial interpolation of second degree on the cell-averages as proposed by Harten [28], which yields third-order accuracy. For the one-dimensional case, it follows that

[TABLE]

[TABLE]

where $\tilde{\mathbf{Q}}^{\ell}_{i}$ is an approximation of the value $\bar{\mathbf{Q}}^{\ell}_{i}$. With this choice, the operator satisfies the locality and the consistency properties, namely: the interpolation into a child cell is computed from the cell-averages of its parent and its nearest uncle cells in each direction; and prediction and projection operator are consistent, i.e. $P_{\ell+1\rightarrow\ell}\circ P_{\ell\rightarrow\ell+1}=\mbox{Identity}$.

The prediction operator is used to obtain the wavelet coefficients $d^{\ell}_{i}$ of the finer cells. The wavelet coefficients are then given by the difference between the values on the finer level and the predicted values as

[TABLE]

The values $d^{\ell}_{i}$ are also used for reconstructing the finest levels without errors due to their property of being the interpolation error. Their norm yields the local approximation error. Moreover, the information of the cell-average value of the two children is equivalent to the knowledge of the cell-average value of the parent and one independent detail.

The same idea can be extended to higher-dimensional cases. For instance, for two dimensions the information of the cell-averages of four children is equivalent to the knowledge of three wavelet coefficients in the different directions and the nodal value on the coarser mesh.

In the process to obtain the adaptive meshes, we flag all cells in which the associate wavelet coefficients $d^{\ell}$ are larger than a threshold. We can select $d^{\ell}$ for only a scalar value from the MHD quantities or adopt other combinations between the variables. For the examples presented below we consider only a scalar value for density or pressure.

Choice of the threshold

There are different possible choices for the threshold, which enable the identification of the retained wavelet coefficients having magnitude above the threshold. In practice, its value should be chosen such that the perturbation related to the thresholding and the discretisation errors is of the same order. Moreover, it is possible, for instance, to use a constant threshold value $\epsilon$ for all levels. However, in the context of finite volumes, and based on our own experiences, [32, 5], we usually follow Harten’s thresholding strategy, i.e,

[TABLE]

in order to control the $L_{1}$-norm. In this case, $\mathbb{L}$ is the finest scale level added to the base mesh level in AMROC, the dimension parameter is $d=2$ or $3$ according to the dimension used, and $|\Omega|$ is the cell area.

3.2 Clustering algorithm

After evaluating the refinement indicators and flagging cells for refinement, a special clustering algorithm [25] is used to create new refinement patches until the ratio between all cells and flagged ones in every new submesh is above a prescribed value $0<\eta\leq 1$. Central to the patch-based mesh refinement approach is the utilisation of a dedicated algorithm to create blocks (or patches) from individual cells tagged for refinement by any of the above criteria. We use a recursive algorithm proposed by Bell et al. [33]. This method, inspired by techniques used in image detection, counts the number of flagged cells in each row and column on the entire domain. The sums $\Upsilon$ are called signatures. First, cuts into new boxes are placed on all edges where $\Upsilon$ is equal to zero. In the second step, cuts are placed at zero crossings of the discrete second derivative $\Delta=\Upsilon_{\nu+1}-2\,\Upsilon_{\nu}+\Upsilon_{\nu-1}$. The algorithm starts with the steepest zero crossing and uses recursively weaker ones, until the ratio between all cells and flagged ones in every new mesh are above the prescribed value $\eta$. An illustration of the general clustering procedure is given in Figure 2.

4 Computational experiments and discussions

In order to compare and assess different refinement criteria quantitatively, we compute the $L_{1}$ error between the adaptive solution and its corresponding uniform mesh solution related to the same maximal resolution used in the adaptive computation, as discussed in [8]. We denote this error by $L_{1,\mathrm{AMR}}$.

Computations are run in parallel on nodes of a recent GNU/LINUX compute cluster that provides $20$ cores with shared memory per node. The AMROC system is parallelised through the MPI library with dynamic re-partitioning. Load balancing is carried out after each level-0 time step in the adaptive cases. AMROC is pursuing a rigorous domain decomposition strategy, in which the increased computational expense on higher refinement levels in the patch-based AMR algorithm is considered in evaluating parallel workload; however, only units of smallest resolution corresponding to a cell on level zero are utilised [34]. This approach simplifies the implementation and reduces the expense of the partitioning algorithm, but it can lead to slight load imbalances on deep hierarchies. The algorithm used for partitioning is always a multi-dimensional space filling curve [9, 10].

In all experiments we use a Cartesian mesh with HLLD numerical flux introduced in [35] and a MinMod limiter as discussed in [36].

4.1 Orszag–Tang vortex

Nowadays, this classical experiment is also used to test how robust a compressible MHD code is at handling the formation of MHD shocks [37], shock-shock interactions, discontinuities, and other structure formations, as presented in [38]. Therefore, this test is also frequently used for code verification and comparison, and it can also demonstrate how significant magnetic monopoles, i.e the $\nabla\cdot B=0$ condition, affect the numerical solution, as discussed in [39] and references therein. Here we use this test to verify how the wavelet-based criteria handle shock formations. We also present a comparison with the Scaled Gradient (SG) refinement criteria already implemented in the AMROC framework [10, 8, 40]. For SG criteria in two-dimensions, a cell at a position $(j,k)$ is flagged for refinement if at least one of the following relations is satisfied for mass density, for instance,

[TABLE]

where the constant $\epsilon^{\,\rho}$ denotes the prescribed refinement limit; in subsequent parts we simply denote it by $\epsilon=\epsilon^{\,\rho}$.

Computational set-up

We consider as initial conditions the density $\rho=\gamma^{2}$, the pressure $p=\gamma$, and the periodic velocities with $u_{x}=-\sin(y)$, and $u_{y}=\sin(x)$ together with the magnetic field $B_{x}=-\sin(y)$, $B_{y}=\sin(2x)$, using the parameters $\gamma=5/3$ and CFL $\nu=0.3$. The computational domain is $0~{}\leq~{}x~{}\leq~{}1$, $0~{}\leq~{}y~{}\leq~{}1$ with periodic boundary conditions. This magnetic field is constructed using a periodic vector potential to guarantee vanishing divergence of the magnetic field.

Numerical results

In this experiment we have observed that the wavelet-based MR criteria detect the shock formations reproducing the expected structures. Moreover, we have found that it produces less refinement of additional features at the maximal level compared to the gradient criterion and additionally stronger coarsening. Therefore, MR prevents unnecessary over-refinement while preventing improper coarsening. These representative effects are the reason why the MR criterion with hierarchical thresholding achieves a smaller error than the SG criterion. An example of these effects is presented in Figure 3 (b, and c, in the right panel). The adaptive mass density solution is presented in Panel (a) for the uniform mesh and it is represented by isolines. Pseudo-colours are used to identify the levels. These results are computed for $L=5$ maximum refinement level, coarsest mesh $32^{2}$, i.e., corresponding to a uniform mesh of $512^{2}$ cells, at final time $t_{e}=\pi$.

To compare the number of cells we have used in Table 1 the value $\eta=0.80$ and threshold parameters which lead to similar errors for SG and MR. Therefore, even with slightly smaller errors, MR needs less cells to obtain the expected representation, resulting in a gain of $1.3$ considering the ratio between the number of cells and the error.

This result is highlighted in panels (b) and (c) on the right of Figure 3. The SG method presents refinement from levels $2\!-\!5$, whereas the MR method from $3\!-\!5$. We have also observed that SG uses much more refinement at the highest level. On the left panel, we present the processor distribution for the three cases. All show a good distribution on the processors. Especially for the adaptive cases, the distribution pattern follows the data organisation on the level refinements. As the MR case has less blocks in the highest level and more regions at the same level, its distribution is slightly less fragmented than for the SG one.

The accuracy of the adaptive solution depends on the choice of the MR threshold parameter with a computational cost determined by the number of cells used in the representation. In our case the number of cells is directly related to the CPU time. Figure 4 illustrates the accuracy behaviour measured by the $L_{\mathrm{1,AMR}}$ error, which is related to the choice of the threshold (vertical axis at right). Therefore, we can estimate roughly the error when choosing the threshold based on this behaviour. On the left vertical axis, we present the number of cells. The number of cells in the adaptive mesh increases as the error decreases. Again, we can deduce a rough estimator how to proceed for choosing a certain threshold to obtain a certain number of cells based on the desired accuracy.

4.2 Magnetic shock cloud

This experiment models a disruption of a high-density magnetic cloud by a strong shock wave, as described in [41]. The initial condition of this Riemann problem defines a region of the advancing plasma – which causes the shock – and a stationary state where the shock advances. We also define a spherical magnetic cloud as a high density region, like a plasmoid, in hydrostatic equilibrium with the surrounding plasma.

We compare the performance obtained with the serial multiresolution Carmen-MHD code [42, 4, 5, 43], which is the prototype solver we had implemented initially during this development project. The detailed software design aspects of both software systems are discussed in [8]. This experiment also provides the performance of the wavelet patch-based method when dealing with high speed flows.

Computational set-up

The cloud region has its centre at $(0.25,0.5,0.5)$ and radius $r_{0}=0.15$ with density $\rho=10$ inside the cloud and $\rho=1$ otherwise. We have adopted outlet boundaries and these regions are limited by the domain boundaries and a plane parallel to the $yz$ plane at $x=0.05$. The advancing plasma initial condition is given by $\rho=3.86859$, $p=167.345$, $u_{x}=11.2536$, $u_{y}=u_{z}=B_{x}=0$, $B_{y}=2.1826182$, and $B_{z}=-B_{y}$, with a stationary state given by $p=1$, $\mathbf{u}=0$, $B_{y}=B_{z}=0.56418958$. The computational domain is $[0,1]^{d}$, where $d=2$ or $3$ respectively stand for the two- and the three-dimensional case. We use the following parameters: CFL $\nu=0.3$, $\gamma=5/3$, the MR threshold is applied to the pressure variable and we run the simulations until the final time $t_{e}=0.06$.

Numerical results

For two-dimensional simulations the wavelet-based adaptive solutions present similar behaviour for both Carmen-MHD and AMROC solvers. Both capture well all the expected relevant physical structures, especially the bow shock, and moreover the symmetry of the solution is almost perfectly preserved.

Figure 5 presents pressure solutions and their respective adaptive meshes at $t_{e}$ for both environments. In the Carmen-MHD code, the cell-based structure requires many less cells in the representation than in the patch-based one used in the AMROC framework. At the final time $t_{e}$, AMROC needs $311,612$ cells clustered in $528$ blocks with $\eta=0.99$ and Carmen-MHD $40,329$ cells. On the other hand, we obtain an improvement in CPU time of about $9$ times in AMROC compared to Carmen-MHD, considering mono-processor runs. Similar results were obtained for the hydrodynamic Euler solvers in both environments [8].

With the additional parallelisation in the AMROC framework, the improvement becomes even larger, with a speed-up of $39$ using $16$ processors, as reflected in the CPU times in Table 2. AMROC simulations are performed with seven refinement levels over a base mesh $16^{2}$, using $2^{n}$ processors with $n$ from one to four. We observe that in these experiments the maximum AMROC scalability is near $1.8$ for one to two processors and it reduces to a factor of around $1.3$ for the four and eight processor cases. Considering eight and sixteen processors, this factor is smaller due to a reduction of the problem size, as expected.

Considering the $3$D case with computations using one processor and a MR threshold $\epsilon=0.01$ until time $t_{e}$, the Carmen-MHD code with a corresponding $128^{3}$ refined mesh needs a CPU time of $16,816$ seconds, while AMROC needs only $2,620$ seconds considering a base mesh $32^{2}$, and $\eta=0.80$. Therefore, we achieve a speed-up of a factor $6$. Furthermore, if $60$ processors are considered, the speed-up increases to around $84$, as this AMROC computation uses only $201$ seconds, in this case with $\epsilon=0.025$ and the same $\eta=0.80$.

Two-dimensional cuts of the pressure solution on the adaptive mesh at time $t_{e}=0.06$ are presented in Fig. 6 (left panel), for $\epsilon=0.025$, $\eta=0.80$, base mesh $32^{2}$, and corresponding fine mesh $1024^{3}$. We can again observe that the symmetry is almost perfectly preserved and the adaptive mesh refines all relevant structures, particularly the bow shock, the front cap, and the internal structures. Additionally, new features could be developed in the three-dimensional scenario, especially, in the tail, and in front of the cloudy structure. Conjointly, the mesh adaptation at time $t_{e}=0.06$ fits very well the main features of the solution (Fig. 6, right panel). Furthermore, the adaptive mesh distribution is well balanced considering the difference in computational costs required by the adaptive meshes (Fig. 6, right panel), as it can be observed in the $2$D cuts in Fig. 7. The processor distribution (not shown) of the $x-z$ plane presents similar results as for the $x-y$ plane. Using Table 3, we show the $3$D performance of AMROC related to three corresponding locally refined meshes, and $60$ processors (distributed on $3$ nodes) at time $t_{e}=0.06$. We observe that the adaptation rises proportionally with the mesh size as we expand the corresponding refined mesh by a factor of two. Basically, a factor of $10$ increases the CPU time as we expand by a factor $2^{3}$ the fully refined mesh. However, the effective number of cells on the adaptive mesh growths only by a factor of $5$.

5 Conclusions

We presented and benchmarked a parallel solver with dynamic mesh adaptation for magnetohydrodynamics, implemented into the MPI-parallel distributed memory framework AMROC. The GLM-MHD model is discretised using finite volumes in two- and three-dimensional Cartesian geometries with explicit time integration. MR criteria are employed for triggering the mesh refinement. We considered two classical MHD benchmarks, the Orszag–Tong vortex and the magnetic shock cloud configuration. The accuracy and CPU time of the developed code were assessed and parallelisation issues including load balancing were analysed. We showed that in comparison with the scaled gradient criterion the new MR implementation yields much better results in terms of accuracy and memory compression. Moreover, our implementation presented a significant improvement compared to its MR-MHD base serial code.

This improvement is in agreement with the conclusions drawn in [8] for hydrodynamic problems where the total CPU time was the primary concern and it was found that patch-based hierarchical data structures yield a better choice. These data structures preserve some memory coherence on the computing data and using auxiliary data avoids repeated generation of topological and numerical procedures. However, this advantage comes at the cost of a more complex implementation. Yet, by using AMROC, including the verified patch-based AMR algorithm implementation plus parallelisation [10], and incorporating the MHD-GLM method as a patch integrator and implementing MR mesh adaptation as a refinement criterion, a performance enhanced two- and three-dimensional adaptive parallel MHD solver based on multiresolution principles has been realised in minimal time.

Acknowledgements

The authors thank the FAPESP SPRINT – University of Southampton (Grant:16/ 50016-9), FAPESP (Grant: 2015/ 25624-2), CNPq (Grants: 307083/2017-9, 306038/2015-3, 140626/ 2014-0, 141741/2013-9), and FINEP (Grant: 0112052700) for financial support of this research. Also ML is grateful to CAPES and University of Southampton (SOU) for financial support for his PhD internships and SOU and I2M, Aix-Marseille Université for kind hospitality. KS acknowledges financial support from the ANR-DFG, grant AIFIT (Grant 15-CE40-0019), and support by the French Research Federation for Fusion Studies (FRCM) within the framework of the European Fusion Development Agreement (EFDA). We are indebted to Eng. V. E. Menconi for his valuable computational assistance.