Ideal and Resistive Magnetohydrodynamic Two-Dimensional Simulation of the Kelvin-Helmholtz Instability in the Context of Adaptive Multiresolution Analysis

ABSTRACT

This work is concerned with the numerical simulation of the Kelvin-Helmholtz instability using an ideal and resistive two-dimensional magnetohydrodynamics model in the context of an adaptive multiresolution approach. The Kelvin-Helmholtz instabilities are caused by a velocity shear and normally expected in a layer between two fluids with different velocities. Due to its complexity, this kind of problem is a well-known test for numerical schemes and it is important for the verification of the developed code. The aim of this paper is to verify the implemented numerical model with the well-known astrophysics FLASH code.

Keywords:
magnetohydrodynamics; Kelvin-Helmholtz instability; adaptive multiresolution analysis; numerical simulation; scientific computing

RESUMO

Este trabalho refere-se à simulação numérica da instabilidade de Kelvin-Helmholtz usando um modelo magneto-hidrodinâmico ideal e resistivo bidimensional no contexto de uma abordagem de multirresolução adaptativa. As instabilidades de Kelvin-Helmholtz são causada por uma velocidade de cisalhamento, e normalmente é esperada em uma camadaentre dois fluidos com diferentes velocidades. Devido à sua complexidade, esse tipo de problema é um teste bem conhecido para esquemas numéricos, sendo importante para a verificação do código desenvolvido. O objetivo principal deste trabalho é verificar a implementação numérica do modelo utilizando o conhecido código FLASH, aplicado ao contexto das instabilidades de Kelvin-Helmholtz.

Palavras-chave:
magneto-hidrodinâmica; instabilidade de Kelvin-Helmholtz; análise multirre solução adaptativa; simulação numérica; computação científica

1 INTRODUCTION

The magnetohydrodynamics (MHD) theory describes the dynamics of a conducting fluid inpresence of magnetic fields and constitutes an important tool to study the macroscopic behavior of plasmas ¹¹. In this context, the Kelvin-Helmholtz instability, which is commonly expected in boundary layers separating two fluids, is an important and a complex physical problem that can be studied with the MHD models, and should often occur in both astrophysical and space geophysical environments ⁶. On the discretization of the MHD system, we use a finite volume (FV) method combined with an adaptive multiresolution (MR) approach to create a computational mesh refined where local structures are presented. The FV method is based on the integral form of conservation laws and guarantees the conservation of the model quantities. The MR for cell-averages was firstly introduced by Ami Harten ¹⁰. The idea is to represent a set of data in different levels of resolution by using a wavelet transform. In the C ⁺⁺ code named CARMEN ¹⁵ the MR algorithm is implemented for compressible Navier-Stokes and five more system of equations. The ideal MHD equations were added later to the CARMEN code ^{3), (8), (9}, and it is employed herein.

We use the FLASH code ⁷, developed by the Flash Center in University of Chicago, well-known in astrophysics and space geophysics, to create a reference MHD solution to our results, since it is not possible to obtain a exact solution. The goal of this work is to verify the numerical results of CARMEN code for the Kelvin-Helmholtz instability problem by comparing them with the reference solution, which is obtained in a regular Cartesian mesh.

The content is organized as follows. In Section 2, we briefly present both the MHD model and the MR approach we use to simulate the Kelvin-Helmholtz instabilities; in Section 3, the numerical methodology and implementation; in Section 4, the results and discussion. The final remarks are presented in Section 5.

2 THE MHD MODEL

The ideal model describes the behavior of a perfectly conducting fluid under the influence of a magnetic field. By adding a resistive term to the MHD system, there is no magnetic flux conservation anymore, which can lead to a more diffusive behavior. The resistivity is associated to the parameter η, which comes from the Ohm’s law. For η ≠ 0 it can trigger a different behavior to be studied in a plasma. It is important to note that when η → 0 the resistive MHD model becomes the ideal MHD equations, which describe the conservation of mass, energy, momentum and magnetic flux. In this work we consider η as a constant, but it is also possible to choose a scalar function. In this context, we introduce the resistive MHD equations.

$\frac{\partial \rho }{\partial t}+\nabla ·\left(\rho u\right)=\mathrm{0,}$ (2.1a)

$\frac{\partial e}{\partial t}+\nabla ·\left[\left(e+p+\frac{|B{|}^2}{2}\right)u-B\left(u·B\right)\right]=\nabla ·\left[B\times \eta \left(\nabla \times B\right)\right],$ (2.1b)

$\frac{\partial \left(\rho u\right)}{\partial t}+\nabla ·\left[\rho u^{t}u+I\left(p+\frac{|B{|}^2}{2}\right)-B^{t}B\right]=0,$ (2.1c)

$\frac{\partial B}{\partial t}+\nabla ·\left[u^{t}B-B^{t}u\right]=-\nabla \times \left(\eta \nabla \times B\right),$ (2.1d)

where ρ represents density, p the pressure, u = (_{^{ux, uy, uz}} ) the velocity vector, B = (_{^{Bx, By, Bz}} ) the magnetic field vector, I the identity tensor of order 2, and γ the ratio of specific heats (γ > 1). The pressure is given by the constitutive law

$p=\left(\gamma -1\right)\left(e-\rho \frac{{\left|u\right|}^2}{2}-\frac{{\left|B\right|}^2}{2}\right),$

where e is the energy density. For the magnetic field, we have the Gauss’ law for magnetism ∇ · B = 0, which means there is no magnetic monopole in the solution of the MHD model. This is an initial condition for the model. From the Faraday’s law ∇ × E = $\frac{\partial B}{\partial t}$ (by applying the divergence operator on the equation), we obtain $\frac{\partial \left(\nabla ·B\right)}{\partial t}$ = 0, i.e., this means there is no variation of the divergence of B over time.

In numerical simulation, the divergence of B does not always vanish. Then it becomes necessary to implement a correction (or divergence cleaning) scheme so the solution will not lead to unphysical behavior or unwanted instabilities. In the next section, we present the numerical methodology for the simulation, including the divergence correction scheme.

3 NUMERICAL APPROACH

To introduce the MHD simulation numerical methodology we firstly present the initial value problem for conservation laws of the form

$\frac{\partial U}{\partial t}+\nabla ·F\left(U\right)=\mathcal{S}\left(U\right),$ (3.1)

$U\left(x,y,t=0\right)=U_0\left(x,y\right),\left(x,y\right)\in \Omega ,$ (3.2)

where U = (ρ, _{^{e, ux, uy, uz, Bx, By, Bz}} ) is the vector of conservative variables, F = F(U) the flux tensor, 𝒮 = 𝒮(U) the source term vector, Ω the domain and t the time. Using the definition of Equation (2.1), we have

$\mathbf{F}\left(\mathbf{U}\right)=\left(\begin{array}{c}\rho \mathbf{u}\\ \left(e+p+\frac{|\mathbf{B}{|}^2}{2}\right)\mathbf{u}-\mathbf{B}\left(\mathbf{u}·\mathbf{B}\right)\\ \rho {\mathbf{u}}^t\mathbf{u}+\mathbf{I}\left(p+\frac{|\mathbf{B}{|}^2}{2}\right)-{\mathbf{B}}^t\mathbf{B}\\ \mathbf{}{\mathbf{u}}^t\mathbf{}\mathbf{B}-\mathbf{}{\mathbf{B}}^t\mathbf{u}\end{array}\right),\mathcal{S}\left(\mathbf{U}\right)=\left(\begin{array}{c}0\\ \nabla ·\left[\mathbf{B}\times \eta (\nabla \times \mathbf{B})\right]\\ \mathbf{0}\\ -\nabla \times (\eta \nabla \times \mathbf{B})\end{array}\right).$ (3.3)

In this section we describe the numerical methodology we use for the MHD simulation in the context of adaptive multiresolution approach.

Divergence Cleaning approach. As discussed in the previous section, ∇ · B = 0 is not satisfied numerically and it can lead to unphysical behavior in the numerical solution of the MHD model. We use the parabolic-hyperbolic correction ¹, in which the errors are propagated and dissipated, by introducing a new scalar variable ψ to the model, adding the term ∇ψ to the right-hand size of Equation (2.1d) and defining a new equation to the ideal MHD system

$\frac{\partial \psi }{\partial t}+c_h^2\nabla ·B=-\frac{c_h^2}{c_p^2}\psi ,$ (3.4)

where _^cp and _^ch are the parabolic-hyperbolic parameters, with _^ch > 0, defined as

$c_h=c_h\left(t\right):={\nu }_{CFL}\frac{\min \left\{\Delta x,\Delta y\right\}}{\Delta t},$

where the Courant number ν_CFL ∈ (0, 1), ∆x and ∆y are the space steps and ∆t is the time step. We also consider the parameter α = ∆h$c_h/c_p^2$, where ∆h = min(∆x, ∆y) ¹³. The presented scheme prevents magnetic monopoles from occurring and preserves the topology of the solution. The model used to perform the simulation consists of the resistive MHD model and the respective additional divergence correction terms or equation.

Finite Volume discretization. The discretization of the resistive MHD model is performed by a finite volume method, which is based on the integral form of the conservation laws and ensures the conservation of the system. The domain is partitioned into cells (or volumes) that are associated to indexes i, j, with i, j ∈ {1, …, N}, where N is the number of cells in each x or y directions. A cell-average at _^tn , n ∈ ℕ, is associated to each cell _^Ci,j of the domain, obtained by the integral of a certain quantity U ⁿ = (_{^{x, y, tn}} ) over the cell, i.e.,

${\overline{\mathbf{U}}}_{i,j}^n=\frac{1}{\left|C_{i,j}\right|}\int {\int }_{C_{i,j}}{\mathbf{U}}^n\mathrm{d}x\mathrm{d}y.$ (3.5)

For simplicity, the superscript n is omitted in the following text. By applying the integral operator to Equation (3.1) in terms of dx, dy for _{^{t = tn}} , and the Divergence Theorem, we obtain

$\frac{\partial }{\partial t}{\overline{\%\mathbf{U}}}_{i,j}=-\frac{1}{\left|C_{i,j}\right|}{\int }_{\partial C_{i,j}}\mathbf{F}·{\mathbf{n}}_k\mathrm{d}{x}_k+{\overline{\mathbf{S}}}_{i,j},$ (3.6)

where ∂_^Ci,j is the boundary of cell _^Ci,j , and k is related to the flux direction. From the first right-hand term we conclude that the flux has to be evaluated on the interfaces of the cells. The values of the flux are based on the centered cell-averages, so it is necessary to approximate the values of F through the cell interfaces. Since the interfaces of a cell _^Ci,j in x direction are located at I ± 1/2, j, and i, j ± 1/2 in y direction, we have

$\begin{array}{l}\frac{\partial }{\partial t}{\overline{U}}_{i,j}=-\frac{1}{\Delta x}\left({\mathcal{F}}_{i+1/\mathrm{2,}j}-{\mathcal{F}}_{i-1/\mathrm{2,}j}\right)\text{}\\ \begin{array}{ccc}& & \end{array}-\frac{1}{\Delta y}\left({\mathcal{G}}_{i,j+1/2}-{\mathcal{G}}_{i,j-1/2}\right)+{\overline{S}}_{i,j},\end{array}$ (3.7)

where ℱ, 𝒢 are the estimated fluxes (or numerical fluxes) on the boundaries. The numerical scheme we use is the Harten-Lax-van Leer Discontinuities (HLLD) Riemann solver ¹² due to its efficiency to resolve isolated discontinuities, and the monotonized central (MC) reconstruction ¹⁶ to ensure the 2^nd -order in space.

Time evolution. The time evolution of Equation (3.7) is performed by a compact second-order Runge-Kutta explicit scheme. The system is completed by suitable initial and boundary conditions. The two-dimensional form of this system is considered.

Adaptive multiresolution analysis. The MR method is based on an adaptive cell average approach as discussed in ^{2), (10}. To briefly present the adaptive MR analysis, we should first introduce the MR representation. The main idea of MR is to decompose the data into several scaled versions of it. In the context of our work, our data is an MHD solution represented by a set of cell averages (mesh) ⁹. The refinement levels are associated to a multiresolution mesh structure that creates dyadic embedded cell meshes. Those meshes have different numbers of cells according to the level they belong to.

Let ℓ be a given level of refinement such as 0 ≤ ℓ ≤ L, where L is the most refined level, and ${\overline{U}}^{\mathcal{l}}$ is the cell-averages at level ℓ. The computational mesh for each level ℓ has 2^2ℓ cells. Since we aim to decompose/reconstruct the data into the refinement levels and also navigate between them, we define the projection and prediction operators

$P^{\mathcal{l}+1\to \mathcal{l}}:{\overline{U}}^{\mathcal{l}+1}\to {\overline{U}}^{\mathcal{l}},$ (3.8)

$P^{\mathcal{l}\to \mathcal{l}+1}:{\overline{U}}^{\mathcal{l}}\to {\overline{U}}^{\mathcal{l}+1}.$ (3.9)

To obtain coarser data from a refined set, for instance ℓ + 1 → ℓ, we use the projection operator defined in Equation (3.8), which is exact and unique. Otherwise, to obtain data in a more refined level from coarse data, the prediction operator is used to perform an estimation of these values. The approximation errors obtained with the prediction are called details or wavelet coefficients, denoted by d ^ℓ = [$d_{i,j}^{\mathcal{l}}$, i ∈ {1, …, 2^ℓ}, j ∈ {1, …, 2^ℓ}] and they give the information about the local regularity of the data, i.e., if the solution is locally smooth or not. As a consequence of these processes, we obtain a relation between each two adjacent levels ${\overline{U}}^{\mathcal{l}+1}\leftrightarrow \left\{d^{\mathcal{l}},{\overline{U}}^{\mathcal{l}}\right\}$, which means data in a more refined level can be acquired from the data in the coarser level along with its wavelet coefficients. Analogously, by extending this concept, it is possible to create a relation between the coarsest and most refined levels as following

${\overline{U}}^L\leftrightarrow \left\{d^{L-1},d^{L-2},\dots ,d^0,{\overline{U}}^0\right\}.$ (3.10)

Such relation is the basis of the multiresolution transform. This transform allows us to navigate through the levels and obtain the data at each level of refinement. From this representationcomes the adaptive MR.

The idea of adaptivity starts from the wavelet coefficients, which can measure the local regularity of the data according to a given threshold parameter є^ℓ = є(є⁰, ℓ), where ℓ denotes the cell scale level and є⁰ is the initial threshold parameter. In this work, we take into account thelevel-dependent and constant threshold parameters. The first varies with the level ℓ, in other words ^{2), (10}, at each level of refinement we compute є^ℓ for the two-dimensional case as

${ϵ}^{\mathcal{l}}=\frac{{ϵ}_0}{\left|\Omega \right|}{2}^{2\left(\mathcal{l}-L+1\right)},0\le \mathcal{l}\le L-\mathrm{1,}$ (3.11)

where |Ω| is the area of the domain. The constant parameter remains unchanged for every ℓ, with 0 ≤ ℓ ≤ L - 1. When the details are larger than є^ℓ, the computational mesh needs to be more refined locally; otherwise the mesh can remain coarser. This methodology allows the computational mesh to be more refined only where it is required. For simplicity, we denote the level dependent and constant threshold parameters as є⁰ and є, respectively.

3.1 CARMEN code

The original CARMEN code was developed by Olivier Roussel ¹⁵ in C ⁺⁺ to simulate theAdvection-diffusion, Burgers-diffusion, Flame front, Flame ball, Flame-curl interaction andNavier-Stokes equations with the finite volume method in the context of the adaptive multiresolution analysis for cell-averages. The ideal MHD model implementation in two dimensions along with two Riemann solvers have being done since 2012 ^{3), (8}. After that, some modifications were made to improve the CPU time and boundary conditions, reconstruct the variables, fix MHD waves evaluation and more ⁹. Subsequently, the three-dimensional model, finite volume approach for MHD, and resistive MHD model is implemented, which is the most recent modification. In this section we present the numerical methodology implemented in this code.

3.1.1 Notes on implementation

The additional source terms in Equations (2.1b) and (2.1d) require some modifications in the code. First we present these terms in a less simplified form separately, and consider the current density J = ∇ · B, with J = (_{^{Jx, Jy, Jz}} ). In two-dimensions, we obtain:

$J=\left(\frac{\partial B_z}{\partial y},-\frac{\partial B_z}{\partial x},\frac{\partial B_y}{\partial x}-\frac{\partial B_x}{\partial y}\right),$ (3.12)

as there is no variation in z direction.

Energy flux term. This term is part of the energy density Equations (2.1b), and it is denoted by ∇ · (B × ηJ). Since we have already computed a divergence operator for the flux evaluation, it is possible to use the same computation for the term B × ηJ. Thus, the flux functions ℱ and 𝒢 in Equation (3.7) for the variable ℰ can be rewritten as

$\hat{\mathcal{F}}\left(\mathcal{E}\right)=\mathcal{F}\left(\mathcal{E}\right)-\eta (B_y{J}_z-B_z{J}_y)$ (3.13)

$\hat{\mathcal{G}}\left(\mathcal{E}\right)=\mathcal{G}\left(\mathcal{E}\right)-\eta (B_z{J}_x-B_x{J}_z),$ (3.14)

by adding the components of B × ηJ. We should recall that the fluxes $\widehat{\mathcal{F}}$ and $\widehat{\mathcal{G}}$ are computed on the interfaces of the cells in x and y directions, respectively.

Ohmic resistivity term. In this case, we denote the Ohmic resistivity in terms of J and we obtain

$\nabla \times \left(\eta J\right)=\eta \left(\frac{\partial J_z}{\partial y},\frac{\partial J_z}{\partial x},\frac{\partial J_y}{\partial x},\frac{\partial J_x}{\partial y}\right),$ (3.15a)

$=\frac{\partial }{\partial x}(0,-\eta J_z,\eta J_y)+\frac{\partial }{\partial y}(\eta J_z,0,-\eta J_x).$ (3.15b)

Here we also take advantage of the divergence operator. Since we have computed a derivative in x for ℱ and in y for 𝒢, the vector components in Equation (3.15b) are added to the fluxes according to its space derivative, i.e., we add the x derivatives to ℱ and, similarly, the y derivatives to 𝒢, always taking into account that each component of this vector is related to an equation of the magnetic field. Thereby, we update the fluxes as following

$\widehat{\mathcal{F}}\left(B_x\right)=\mathcal{F}\left(B_x\right),$ (3.16)

$\hat{\mathcal{F}}\left(B_y\right)=\mathcal{F}\left(B_y\right)-\eta J_z,$ (3.17)

$\hat{\mathcal{F}}\left(B_z\right)=\mathcal{F}\left(B_z\right)+\eta J_y,$ (3.18)

and for y direction

$\hat{\mathcal{G}}\left(B_x\right)=\mathcal{G}\left(B_x\right)+\eta J_z,$ (3.19)

$\widehat{\mathcal{G}}\left(B_y\right)=\mathcal{G}\left(B_y\right),$ (3.20)

$\hat{\mathcal{G}}\left(B_z\right)=\mathcal{G}\left(B_z\right)-\eta J_x,$ (3.21)

which are calculated on the cells interfaces. Thus Equation (3.7) becomes

$$\begin{array}{l}\frac{\partial }{\partial t}{\overline{U}}_{i,j}=-\frac{1}{\Delta x}\left({\widehat{\mathcal{F}}}_{i+1/\mathrm{2,}j}-{\widehat{\mathcal{F}}}_{i-1/\mathrm{2,}j}\right)\text{}\\ \begin{array}{ccc}& & \end{array}-\frac{1}{\Delta y}\left({\widehat{\mathcal{G}}}_{i,j+1/2}-{\widehat{\mathcal{G}}}_{i,j-1/2}\right).\end{array}$$ (3.22)

It is important to note that the approach we use here is not unique, so other forms to evaluate these terms may be used.

3.2 Code verification: reference solution

The FLASH code is a modular, parallel metaphysics code developed in FORTRAN90 and C in the Flash Center of the University of Chicago, well-known in the space geophysics and astrophysics fields. The code includes the ideal and resistive MHD models with the finite volume method and it is also possible to use an adaptive mesh refinement or non-adaptive approaches for the simulation. Since FLASH is a verified FV code, we use it in simulations as a reference solution.

The similarities between these codes in the context of this work are the following: a Harten-Lax-Van Leer-Discontinuities Riemann solver and monotonized central variable reconstruction. Moreover the following configurations are used: For the time evolution, FLASH performs a one-step Hancock; FLASH also uses the 8-wave based divergence correction ¹⁴, which stabilizes the numerical method by adding source terms proportional to ∇ · B on the right-hand side of the system, i.e. the source term vector 𝒮 = (0, -B∇ · B, -u∇ · B, -(u · B)∇ · B) is added to the MHD equations.

FLASH code can also be adaptive, but we are interested only in its non-adaptive form as our reference. We do not compare the mesh refinement and MR adaptivity, however that is possible, as discussed in ⁴.

In particular, for FLASH code the source terms in Equations (2.1b) and (2.1d) are added to ℰ and B quantities after the time evolution at each timestep.

4 NUMERICAL EXPERIMENTS

In this section we present the results obtained with CARMEN code for L = 9 and the respective reference computed with FLASH with the same level i.e., 512 × 512 mesh. We reinforce here that the reference results allow us to compare the solution qualitatively and quantitatively, then it is possible to measure how the CARMEN code solution is converging to the expected results.

In fluids or plasmas, the Kelvin-Helmholtz instability is triggered by a velocity shear and creates a vortex as discussed in ⁶. That occurs in several phenomena in nature and space, in the clouds and magnetosphere, for instance ⁵. In our simulations, we consider the initial conditions presented in Table 1, with

$$\begin{array}{l}{u}_x^0=5\left(\tanh \left(20\left(y+0.5\right)\right)-\left(\tanh \left(20\left(y-0.5\right)\right)+1\right)\right),\text{}\\ u_y^0=0.25\sin \left(2\pi x\right)\left(\exp \left[-100{\left(y+0.5\right)}^2\right]-\exp \left[-100{\left(y-0.5\right)}^2\right]\right),\end{array}$$

and periodic boundary conditions everywhere. The computational domain is Ω = [0, 1.0] × [-1.0, 1.0] and the finest scale is L = 9. We also define as parameters the Courant number ν = 0.4, γ = 1.4, and the physical time t = 0.5. The choice of resistivity η = 0.02 is made according to astrophysics uses, as discussed in ⁵ and for the ideal case we consider η = 0. We have also tested other values of η as [0.005, 0.05]. We consider the threshold parameters є⁰ = 0.1 and 0.01 both for the ideal and resistive cases. For the divergence-cleaning, we adopt α = 0.4, as discussed in ⁹.

Table 1:
Kelvin-Helmholtz instability initial condition.

ρ	p	ux	uy	uz	Bx	By	Bz
1.0	50.0	$u_x^0$	$u_y^0$	0.0	1.0	0.0	0.0

The reference (left) and CARMEN (right) code solutions for ideal and resistive MHD are presented for the variables density ρ, velocity y-component _^uy , and magnetic field x-component _^Bx in Figures 1 and 2. In both simulations the obtained results have the same maximum and minimum values. The larger gradients values in the solution are presented where the vortexes are located, near to y = -0.5 and y = 0.5. For the variable density, for instance, it is possible to find gradients of one order of magnitude for the ideal case, and even not using the entire simulation mesh these structures are well captured. The resistive results are clearly more diffusive when compared to the ideal case, as expected, since it is the effect of the additional resistive terms. Figures 3 and 4 present the cuts for the variables ρ, _^uy , and _^Bx , in x = 0.5 and y = 0.5 at t = 0.5, respectively, of reference and CARMEN code solutions for ideal and resistive MHD. These cuts are located where the most critical structures in solution are placed, due to many discontinuities, and they remain very close for each case with both codes.

Figure 1:
Ideal MHD variables ρ, _^uy and _^Bx (from top to bottom) obtained at t = 0.5 and L = 9, with reference and CARMEN code solutions for є⁰ = 0.1.

Figure 2:
Resistive MHD variables ρ, _^uy and _^Bx (from top to bottom) obtained at t = 0.5 and L = 9, with reference and CARMEN code solutions for є⁰ = 0.1.

Figure 3:
Ideal MHD cuts of variables ρ, _^uy and _^Bx (from top to bottom) obtained at t = 0.5 and L = 9, with reference and CARMEN code solutions for є⁰ = 0.01, є = 0.05, є⁰ = 0.1, on x = 0.5 and y = 0.5.

Figure 4:
Resistive MHD cuts of variables ρ, _^uy and _^Bx (from top to bottom) obtained at t = 0.5 and L = 9, with reference and CARMEN code solutions for є⁰ = 0.01, є = 0.05, є⁰ = 0.1, on x = 0.5 and y = 0.5.

According to Figures 3 and 4, the solution obtained with the resistive MHD model approaches the reference solution. We can also notice that for sufficiently big values of є the convergence of the solution is not sufficient, which means that the number of cells used to represent the solution is not enough, thus it is important to find an acceptable threshold parameter in order to obtain a fair relation between CPU time, number of cells, and accuracy.

To illustrate this situation, we can observe the cuts for є = 0.05, which do not converge properly to the reference solution when compared to the other two cases. However, this behavior is expected since the percentage of cells needed for those simulations are 11% and 5% for the ideal and resistive cases, respectively. This reinforces the importance of the choice of parameter є.

The adaptive meshes for the ideal and resistive cases at time t = 0.5 are presented for є⁰ = 0.1 in Figure 5. In this simulation the use of memory (percentage of cells) over time is 32% for the ideal case and 24% for the resistive case, when compared to a uniform mesh. The percentages of CPU time needed are 67% and 72% smaller for the ideal and resistive MHD, respectively. To quantitatively compare the presented results, Table 2 and 3 show the L ₁, L ₂, and L _∞ error values for every variable of the MHD model when the solution is compared with the reference results. The errors are computed for MR and FV CARMEN solutions and the reference solution. The values for the FV approach are slighty smaller when compared to the MR results. This result is expected since we need a full mesh in order to obtain them. The errors for the resistive case are even smaller when compared to the ideal case. That happens due to the diffusive behavior of the solution, since the greater values occur where the sharpest structures are located. In Figure 6, the diffusive behavior of the resistive terms of the MHD model is illustrated. We can observe how the topology of the solution of the variable ρ changes according to the value of η. In this case we choose η = 0 (ideal case), η = 0.005, η = 0.02 and η = 0.05, and it shows that the solution becomes more diffusive as η is increased. As long as the local structures become more diffuse, the adaptive mesh for greater values of η will have less cells in the most refined level.

Figure 5:
Adaptive MR meshes at t = 0.5 obtained with CARMEN code with є⁰ = 0.1 for the ideal and resistive cases.

Table 2:
Errors for the ideal Kelvin-Helmholtz simulation with MR (є⁰ = 0.1) and FV compared with referece results.

CARMEN code	2*Variables	Errors
		ℒ1 (× 10-2 )	ℒ2 (× 10-4 )	ℒ∞ (× 10-1 )
FV	ρ	0.29	0.14	0.89
	p	18.25	10.53	67.87
	ux	6.57	4.32	22.86
	uy	3.34	1.74	8.38
	Bx	4.04	2.09	10.69
	By	1.65	0.76	3.22
MR	ρ	0.30	0.14	1.00
	p	18.38	9.10	60.11
	ux	6.72	4.14	21.98
	uy	3.82	1.83	8.43
	Bx	4.08	1.81	8.63
	By	1.91	0.80	3.37

Table 3:
Errors for the resistive Kelvin-Helmholtz simulation with MR (є⁰ = 0.1) and FV compared with reference results.

CARMEN code	Variables	Errors
		ℒ1 (× 10-3 )	ℒ2 (× 10-4 )	ℒ∞ (× 10-2 )
FV	ρ	0.13	0.06	0.37
	p	7.57	3.56	7.39
	ux	2.72	1.47	4.05
	uy	1.56	0.90	2.33
	Bx	2.26	1.27	1.94
	By	1.18	0.64	1.07
MR	ρ	0.81	0.38	0.66
	p	55.89	26.04	46.27
	ux	10.04	5.58	13.57
	uy	9.91	5.62	12.37
	Bx	7.27	4.06	5.82
	By	3.38	1.89	2.62

Figure 6:
Variable ρ obtained with CARMEN code at t = 0.5 for the ideal case and η = 0.005, η = 0.02, η = 0.05 (from left to right) and є⁰ = 0.1.

The error obtained for variable p is larger when compared to the other variables of the model. That happens because the pressure depends on the values of B, e and u to be computed, thus the approximation errors of these variables accumulate in the pressure. In Figure 7 we observe that the solution of p is close to the reference solution on ideal (top) and resistive cases (bottom), for є⁰ = 0.01, є = 0.05 and є⁰ = 0.1, at x = 0.5 and y = 0.5.

Figure 7:
Ideal and resistive MHD cuts of variable p (from top to bottom) obtained at t = 0.5 and L = 9, with reference and CARMEN code solutions for є⁰ = 0.01, є = 0.05, є⁰ = 0.1, on x = 0.5 and y = 0.5.

5 FINAL REMARKS

We have presented a wavelet based multiresolution approach to compute the Kelvin-Helmholtz instabilities for the ideal and resistive MHD system with parabolic-hyperbolic divergence cleaning approach and compared its results to a reference solution obtained with FLASH code. New implementations of the resistive MHD model are done in CARMEN code for the MR and FV approaches. It is shown that the topology of the solution of the problem achieved with the adaptive MR is very close to the reference solution and with the respective cuts in the domain and the error values we ensure that it leads to the expected solution quantitatively. The adaptive method has proven to be efficient to keep the accuracy of the solution even decreasing the number of cells in the simulation. In our case we needed an average of 30% of cells for the simulations, which implied in a reduction of 67% of the CPU time (compared to the FV CARMEN code). We conclude that the presented MR methodology for the ideal and resistive MHD models is relevant in this instability context and it would be interesting to extend it to more complex space physics problems in future works.

ACKNOWLEDGEMENTS

M.O.D. and O.M. thankfully acknowledge financial support from MCTI/FINEP/INFRINPE-1 (grant 01.12.0527.00), CNPq (grants 30603-8/2015-3, 31224-6/2013-7) and FAPESP (grants 2015/50403-0, 2015/25624-2, 2016/50016-9). A.G. thankfully acknowledges financial support of MCTI/INPE-PCI and PhD scholarships from CNPq (grants 132045/2010-9, 312479/2012-3, 141741/2013-9). The software used in this work was in part developed by the DOE NNSA-ASC OASCR Flash Center at the University of Chicago. We thank Dr. Olivier Roussel for developing the original MR CARMEN code for hydrodynamics and Dr. Kai Schneider fruitful for scientific discussions. We also want to thank to Dr. Edgard Evangelista for his help with FLASH code and Eng. Varlei Menconi for the computational assistance.

REFERENCES

• ¹
  A. Dedner, F. Kemm, D. Kröner, C.-D. Munz, T. Schnitzer & M. Wesenberg. Hyperbolic divergence cleaning for the MHD equations. Journal of Computational Physics, 175 (2002), 645-673.
• ²
  M.O. Domingues, S.M. Gomes, O. Roussel & K. Schneider. Adaptive multiresolution methods. ESAIM Proceedings, 34 (2011), 1-96.
• ³
  M.O. Domingues, A.K.F. Gomes, S.M. Gomes, O. Mendes, B. Di Pierro & K. Schneider. Extended generalized lagrangian multipliers for magnetohydrodynamics using adaptive multiresolution methods. ESAIM Proceedings, 43 (2013), 95-107.
• ⁴
  R. Deiterding, M.O. Domingues, S.M. Gomes & K. Schneider. Comparison of adaptive multiresolution and adaptive mesh refinement applied to simulations of the compressible euler equations. SIAM Journal on Scientific Computing, 38 (2016), S173-S193.
• ⁵
  E.F.D. Evangelista, M.O. Domingues, O. Mendes & O.D. Miranda. A brief study of instabilities in the context of space magnetohydrodynamic simulations. Revista Brasileira de Ensino de Fìsica, 38(1) (2016).
• ⁶
  A. Frank, T.W. Jones, D. Ryu & J.B. Gaalaas. The magnetohydrodynamic Kelvin-Helmholtz instability: A two-dimensional numerical study. The Astrophysical Journal, 460 (1996), 777.
• ⁷
  B. Fryxell, K. Olson, P. Ricker, F.X. Timmes, M. Zingale, D.Q. Lamb, P. MacNeice, R. Rosner, J.W. Truran & H. Tufo. FLASH: An adaptive mesh hydrodynamics code for modeling astrophysical thermonuclear flashes. The Astrophysical Journal Supplement Series, 131 (2000), 273-334.
• ⁸
  A.K.F. Gomes. Análise multirresolução adaptativa no contexto da resolução numérica de um modelo de magnetohidrodinâmica ideal. Master’s thesis, Instituto Nacional de Pesquisas Espaciais (INPE), São José dos Campos, (2012).
• ⁹
  A.K F. Gomes, M.O. Domingues, K. Schneider, O. Mendes & R. Deiterding. An adaptive multiresolution method for ideal magnetohydrodynamics using divergence cleaning with parabolic-hyperbolic correction. Applied Numerical Mathematics, 95 (2015), 199-213.
• ¹⁰
  A. Harten. Multiresolution representation of data: a general framework. SIAM Journal of Numerical Analysis, 33(3) (1996), 385-394.
• ¹¹
  R.J. Hosking & R.L. Dewar. “Fundamental Fluid Mechanics and Magnetohydrodynamics”. Springer Singapore, (2016).
• ¹²
  T. Miyoshi & K. Kusano. A multi-state HLL approximate Riemann solver for ideal magnetohydrodynamics. Journal of Computational Physics, 208 (2005), 315-344.
• ¹³
  A. Mignone & P. Tzeferacos. A second-order unsplit Godunov scheme for cell-centered MHD: The CTU-GLM scheme. Journal of Computational Physics, 229(6) (2010), 2117-2138.
• ¹⁴
  K.G. Powell, P.L. Roe, T.J. Linde, T.I. Gombosi & D.L. De Zeeuw. A Solution-Adaptative Upwind Scheme for Ideal Magnetohydrodynamics. Journal of Computational Physics, 154 (1999), 284-309.
• ¹⁵
  O. Roussel, K. Schneider, A. Tsigulin & H. Bockhorn. A conservative fully adaptive multiresolution algorithm for parabolic {PDEs}. Journal of Computational Physics, 188(2) (2003), 493-523.
• ¹⁶
  E.F. Toro. “Riemann Solvers and Numerical Methods for Fluid Dynamics: A Practical Introduction”. Springer-Verlag Berlin Heidelberg, (1999).

Publication Dates

• Publication in this collection
  Aug 2017

History

• Received
  28 Nov 2016
• Accepted
  10 May 2017