PHARE Parallel hybrid particle-in-cell code with patch-based adaptive mesh refinement

2025-05-02 0 0 2.13MB 21 页 10玖币

PHARE : Parallel hybrid particle-in-cell code with patch-based adaptive

mesh reﬁnement

Nicolas Aunaia,∗, Roch Smetsa, Andrea Ciardib, Philip Deegana, Alexis Jeandeta, Thibault Payetc, Nathan

Guyota, Loic Darrieumerloua

aLaboratoire de Physique des Plasmas (LPP), CNRS, Observatoire de Paris, Sorbonne Universit´e, Universit´e Paris-Saclay,

Ecole polytechnique, Institut Polytechnique de Paris, 91120 Palaiseau, France

bSorbonne Universit´e, Observatoire de Paris, Universit´e´e PSL, CNRS, LERMA, F-75005 Paris, France

cAix Marseille Univ, CNRS, LIS, Marseille, France

Abstract

Modeling multi-scale collisionless magnetized processes constitutes an important numerical challenge. By

treating electrons as a ﬂuid and ions kinetically, the so-called hybrid Particle-In-Cell (PIC) codes represent

a promising intermediary between fully kinetic codes, limited to model small scales and short durations,

and magnetohydrodynamic codes used large scale. However, simulating processes at scales signiﬁcantly

larger than typical ion particle dynamics while resolving sub-ion dissipative current sheets remain extremely

diﬃcult. This paper presents a new hybrid PIC code with patch-based adaptive mesh reﬁnement. Here,

hybrid PIC equations are solved on a hierarchy of an arbitrary number of Cartesian meshes of incrementally

ﬁner resolution dynamically mapping regions of interest, and with a reﬁned time stepping. This paper

presents how the hybrid PIC algorithm is adapted to evolve such mesh hierarchy and the validation of the

code on a uniform mesh, ﬁxed reﬁned mesh and dynamically reﬁned mesh.

1. Introduction

The numerical modeling of magnetized space and

laboratory plasmas represents an important com-

putational challenge that is generally linked to the

need to capture the multi-scale nature of plasma

dynamics. While kinetic equations such as the

Vlasov-Maxwell system for collisionless plasmas can

in principle describe such dynamics, from electron

scales to global system scales, their use in multi-

dimensions is in practice restricted to relatively

small scales. For example, fully kinetic simulations

of magnetic reconnection usually handle scales of

a few tens of ion inertial lengths δi, quite rarely

in three dimensions, and very often with artiﬁ-

cially reduced scale separations through modiﬁed

ion to electron mass ratio and/or reduced speed

of light [1, 2]. Modeling larger scales is typically

done with more approximate physical formalisms.

Among them is the hybrid formalism, in which ions

∗Corresponding author

Email address: nicolas.aunai@lpp.polytechnique.fr

(Nicolas Aunai)

are treated kinetically but electrons are modeled as

a ﬂuid [3, 4]. This allows, in principle, to transfer

the computational power required to resolve elec-

tron scales into solving larger domains for longer

times.

Hybrid codes have been used for modeling large

systems such as (small or artiﬁcially reduced) plan-

etary magnetospheres. However such large scale

models typically poorly resolve the Hall scale [5,

6, 7, 8]. Not resolving such sub-ion scale cur-

rent sheets may result in errors leading to spuri-

ous macroscopic behavior [9]. The main advantage

such models have compared to a well resolved yet

much lighter magnetohydrodynamics (MHD) coun-

terparts is that they account for collisionless popu-

lation mixing.

Modeling scales much larger than typical ion

scales and in three dimensional systems yet having

sub-ion scale space/time resolution remains a very

challenging task for hybrid codes. The dispersive

nature of the whistler waves supported by the colli-

sionless Ohm’s law imposes very strong constraints

on the time step, which scales as the square of the

Preprint submitted to Computer Physics Communications October 17, 2023

arXiv:2210.14580v2 [physics.comp-ph] 14 Oct 2023

mesh size. Therefore, high resolution large scale

models (∆ ∼0.1δi) remain computationally inten-

sive. In practice, such high resolution hybrid simu-

lations have typically been done in domains that are

not vastly larger than those accessible with full-PIC

codes using reduced physical normalized constants.

The small gain in domain size or evolution time

oﬀered by hybrid codes is thus often quickly bal-

anced by the drawback of losing the electron kinetic

physics oﬀered by full PIC ones. In practice, hybrid

codes have not been so competitive once high reso-

lution is needed.

The hybrid equations can in practice either be

solved in a Eulerian way, with so-called Vlasov hy-

brid codes [10], or with a Lagrangian perspective,

with the so-called hybrid Particle-In-Cell (PIC)

codes [11]. PIC codes present the advantage over

Vlasov ones to be inherently adaptive in their de-

scription of the particle distribution. This is more

eﬃcient computationally than resolving ﬁne struc-

tures developing for diﬀerent reasons in position

space and velocity space and for each of the pop-

ulations composing the plasma. Some features are

also easier and/or lighter to handle in PIC than in

Vlasov models, such as the possible interactions be-

tween the populations. On the other hand, Vlasov

codes present the advantage to be noise-free con-

trary to PIC results which are inherently based on

statistical sampling of the continuous distribution

function. Recent advances have been achieved in

using the Vlasov hybrid formalism to model large

scales in the context of the Earth’s magnetosphere

[10]. Yet, simulations remain heavier than their

PIC counterpart for which noise can moreover be

reduced through the use of higher order interpola-

tion kernels and increased number of macroparti-

cles per cell. PIC codes thus constitute a promising

approach for developing multi-scale models.

There has been several attempts to make hybrid

PIC codes more adapted to model multi-scale sys-

tems. Time scale disparities have been for instance

addressed by decomposing the domain into diﬀer-

ent time zones, each evolving with a proper time

step [12]. Evolving quantities based on local event

triggers rather than at regular and global time step

intervals [13, 14, 15] has also been considered.

Spatial scale disparities, on the other hand, are

generally handled by adopting a non-uniform mesh-

ing. This can be achieved either by reﬁning the grid

in the so-called adaptive mesh reﬁnement (AMR)

methods [16, 17, 18, 19, 20, 21, 22, 23], or adapting

the position of the nodes of a single mesh in the

so-called moving mesh adaptation (MMA) method

[24, 25]. In both cases the goal is to focus the mesh-

ing in regions developing small scale features where

more computational power will be dedicated, and

for which one needs to deﬁne ad-hoc detection cri-

teria.

Adaptive mesh reﬁnement (AMR) has now be-

come a mainstream feature of ﬂuid codes [26, 27,

28]. It is however much less used in PIC codes. The

main reason probably stems from the algorithmic

complexity in handling macroparticles across multi-

ple grid levels. Only one AMR fully kinetic electro-

magnetic code is currently in use in our knowledge

[18, 19, 22, 23].

An AMR Hybrid PIC code following a hybrid

block AMR method [29] has also been proposed

[21]. It is worth noticing however that the code

has been mostly used with multiple grids that were

ﬁxed in time [30, 31, 32, 33].

In existing AMR implementations for PIC codes,

macroparticles interact only with the ﬁnest mesh

at their location. This method has several draw-

backs. The ﬁrst one, shared with MMA, is that

macroparticles see multiple mesh spacings as they

evolve. Either because they cross a grid boundary

along their trajectory or because the region in which

they currently are is being reﬁned or coarsened in

time (or streched/contracted in MMA). Because

the compact support of the macroparticle interpola-

tion kernel is usually constant in index-space, these

macroparticles will see their physical extent change

in time, potentially leading to a spurious evolution

of their momentum [11]. A more important draw-

back of having particles exploring multiple grids is

that the number of particles per cell strongly de-

creases in reﬁned cells if nothing is done. In prac-

tice this number is usually kept roughly constant

by splitting macroparticles when they enter in re-

ﬁned regions [18, 19, 22, 21]. If the splitting pro-

cedure can trivially conserve the velocity distribu-

tion of macroparticles, it may not conserve their

spatial distribution, and thus their contribution to

moments on the mesh. Several splitting procedures

have been proposed [34, 23, 35, 21]. They all dis-

patch reﬁned macroparticles within some ad hoc

distance from the original coarse one, usually in

the same cell, randomly or not. If macroparti-

cles only interact with the ﬁnest mesh at their lo-

cation, and splitting is used not to decrease their

number per cell, the opposite operation, i.e. the

merging of macroparticles, needs to be done as well

not to stall the simulation. Merging macroparti-

cles is more hazardous than splitting as it generally

does not preserve the local structure of the distri-

bution function [34], the very reason why kinetic

simulations are done in the ﬁrst place.

The Multi-Level-Multi-Domain (MLMD)

method has later been proposed as a way to pre-

vent previous issues [35, 36, 37, 38, 39]. Contrary

to other AMR codes, the MLMD method solves

all equations on all grid levels, as in a patch-based

AMR approach. This means that not only elec-

tromagnetic ﬁelds and moments are deﬁned on

all nodes of all levels, but also macroparticles.

While this may appear as an overhead to deal with

macroparticles in coarse regions where there is

also a ﬁne mesh and its associated macroparticles,

it comes with several advantages. First, the

macroparticles only see one mesh resolution, their

shape is thus perfectly constant in time. Then,

merging macroparticles is not required anymore

since macroparticles can simply be deleted as they

exit a reﬁned level because their is a self-consistent

kinetic ﬂux in the overlapped region of the next

coarser level. As in other AMR PIC codes, reﬁned

levels are fed with macroparticles from splitting

those from coarser levels at level boundaries. And

once the ﬁne solution is obtained, the electromag-

netic ﬁeld is coarsened onto the next coarser level,

and in practice it either overwrites the coarser

solution in the overlapped region, or is averaged

with it.

To our knowledge, existing AMR PIC codes use

in-house developed code for the adaptive meshing

mechanism and evolve the system with a time step

uniform across all mesh levels [19, 21]. In such a

case, the time step is thus constrained by the ﬁnest

grid of the model, which leads to much heavier sim-

ulations than necessary. The MLMD method was

also originally proposed with uniform time stepping

across grid levels. It was then updated to consider

a proper stepping per level [39]. Coarser levels can

evolve much fewer cycles than reﬁned ones, which,

considering the dispersive nature of kinetic plasma

waves, is much more advantageous than codes based

on a uniform and ﬁxed time stepping. Published

results[35, 36, 37, 38, 39] however only demonstrate

the method with only one reﬁned level, consisting

in a single reﬁned patch with a predeﬁned position

which is ﬁxed in time. Such a code is thus useful

when the region in which to enhance the resolution

is known in advance and does not evolve with time.

In complex systems, where critical small scale re-

gions are moving, appear and disappear, the lack of

adaptivity imposes the reﬁnement of a substantial

part of the domain, which may become prohibitive.

Contrary to MHD codes, AMR kinetic codes can-

not have arbitrarily large mesh spacing. Kinetic

plasmas include intrinsic particle scales that need to

be correctly resolved even in regions where ”noth-

ing” happens. Solving explicit fully kinetic equa-

tions on a mesh much coarser than the electron De-

bye length is unstable. Solving hybrid equations

with a mesh and time step much coarser than the

local ion scales is irrelevant, if not wrong. Such an

upper bound to the coarsest mesh resolution make

modeling very large domains expensive even with

AMR. It thus appears interesting to not only con-

sider reﬁning the mesh and the time step, but also

the physical formalism that is resolved. The MLMD

method also appears promising in that regard since

each reﬁnement level, having its own macroparti-

cles, could in principle be coupled to levels evolving

diﬀerent equations, given that one knows how to

transfer information from one to the other.

Coupling a kinetic solver, operating on critical

regions, with a ﬂuid solver, evolving less impor-

tant regions of the domain, has been an important

goal over the last decade. The ﬁrst coupling be-

tween MHD and fully kinetic PIC equations was

achieved for local simulations of magnetic reconnec-

tion [40, 41, 42] along the inﬂow direction. A 2D

coupling was later achieved using anisotropic MHD

and Hall MHD and implicit fully kinetic equations

using the codes BATS-R-US and iPIC 3D [43]. This

method was then extended to 3D [44]. Simula-

tions embedding one or several rectangular full-

PIC regions in a global MHD domain were then

performed for modeling the Earth’s magnetosphere

[45], Ganymede’s magnetosphere [46, 47], the Mars’

magnetosphere [48], or Mercury’s magnetosphere

[49]. The method has later been implemented be-

tween iPIC3D and MPI-AMRVAC [50, 51].

These works provided for the ﬁrst time a large

scale context to otherwise and until then isolated

kinetic boxes. However, the choice of using fully

kinetic equations, even through an implicit scheme,

forces the resolution of some electron particle scales

which prevents a kinetic treatment of ions over large

scales. Coupling ﬂuid equations to hybrid ones

would instead allow to cover much broader regions

with kinetic ions.

In this paper we present the design and imple-

mentation of a new hybrid kinetic PIC code with

adaptive mesh reﬁnement inspired from the MLMD

method. AMR will make high resolution more af-

fordable for large domains and set the code architec-

ture for future multi-formalisms simulations. Sec-

tion 2 presents the hybrid formalism, its physical

assumptions and associated equations. Section 3

details the AMR method we employ. The proposed

code named PHARE, handles an arbitrary number of

reﬁnement levels and patches which position and

shape is adapted dynamically according to the cur-

rent state of the solution and reﬁnement criteria.

A reﬁned time stepping is used to accommodate

the local Courant–Friedrichs–Lewy condition [52]

required by our explicit scheme, and use the largest

possible time step on each reﬁnement level. Sec-

tion 4 presents the validation of the hybrid solver

with and without AMR. Section 5 explains the main

points of the code implementation, in particular re-

lated to the AMR mechanism. Finally, section 6

concludes the paper.

2. Hybrid Particle-in-Cell formalism

2.1. Physical assumptions and model equations

In this section we present the physical assump-

tions concerning the plasma dynamics and the rele-

vant mathematical equations of the hybrid formal-

ism. The magnetic ﬁeld is evolved from the electric

ﬁeld via Faraday’s equation:

∂B

∂t =−∇ × E(1)

The total current density jis obtained from Am-

pere’s equation where the displacement current is

neglected:

µ0j=∇ × B(2)

where µ0is the vacuum permittivity. The total ion

density, ni, is the sum of the density npof each

ion population p, obtained from integration of their

distribution function fp:

ni(r, t) = X

pZfp(r,v, t)dv(3)

The total ion bulk velocity uis computed from

the particle ﬂux of all populations and the total ion

density:

u(r, t) = 1

niX

pZvfp(r,v, t)dv(4)

The hybrid formalism assumes that the plasma

dynamics occurs at spatial and temporal scales for

which quasi-neutrality holds, i.e. the total ion

charge density ρi=Ppnpqp=PpeZp, where np,

qpare the particle density and the electric charge

of the population p, respectively, is equal to the

electron charge density ene. With this assumption

the electron bulk velocity can be calculated from

the known electrical current density j, and ion mo-

ments:

ve=u−1

ene

j(5)

The major physical assumption of the hybrid for-

malism, that is the less justiﬁable in collisionless

systems, is the closure of the electron moment hi-

erarchy at the pressure level. Currently PHARE uses

the isothermal closure:

Pe=nekBTe(6)

where kBis the Boltzmann constant. This closure is

the simplest and most commonly used. It assumes

that the electron pressure tensor reduces to a sim-

ple scalar, Pe, proportional to the electron density.

The coeﬃcient Tedeﬁnes the temperature which

is a constant in both space and time. Since the

ﬁrst three electron moments are known, one can use

and rewrite the electron ﬂuid momentum equation

to calculate an electric ﬁeld that is consistent with

them. We further assume time and spatial scales

at which the electron ﬂuid is accelerated and thus

neglect the bulk inertia in the momentum equation.

The electric ﬁeld is then given by:

E=−ve×B−1

ene∇Pe(7)

In practice, eq. 7 can be problematic because

it cannot prevent the possibly inﬁnite thinning of

current sheets. In reality, the thickness of any cur-

rent sheet would be limited by intrinsic electron

particle orbit scales, leading to non-negligible non-

gyrotropic pressure tensor components and bulk in-

ertia. Since we do not include these scales, we

need to add terms in eq. 7 to include a dissipative

scale controlling the minimal current sheet thick-

ness in our systems. Two terms are added, the

classical Joule resistive term ηj, where ηis constant,

which diﬀuses large scale structures, and the hyper-

resistive term −ν∇2j, where νis constant, which on

the contrary is able to dominate over the Hall term

and deﬁnes a dissipation-dominated small scale [53].

The inclusion of these two terms leads to :

E=−ve×B−1

ene∇Pe+ηj−ν∇2j(8)

Finally, in the hybrid formalism, the evolution

of the distribution function of each ion population

p, used in equations 3 and 4, is described by the

Vlasov equation:

∂fp

∂t +v·∂fp

∂r+qp

(E+v×B)·∂fp

∂v= 0 (9)

2.2. Macroparticles and their interaction with the

mesh

In practice, equation 9 is not directly solved

in PIC codes, but its solution is equivalently ob-

tained by the evolution of a collection of so-

called macroparticles, which, according to the PIC

method, follow the same dynamical equations as

real particles:

dt =qp(E+v×B) (10)

dt =v(11)

The moments of the ion distribution used in 8

are obtained from macroparticles. The distribution

function fpof an ion population pat a given point

in phase (r,v) space writes:

fp(r,v) =

l=1

S(r−rl)δ(v−vl) (12)

where the sum is performed over the Np

macroparticles of the population p,Sis a symmet-

ric kernel normalized to 1 over its compact support,

and δis the Dirac function. Macroparticles thus

have a ﬁnite spatial extent and a unique velocity.

Each of the Ncmacroparticles in a cell is assigned

a statistical weight W=Vcnc/Ncrepresenting its

contribution to the cell’s density ncof volume Vc.

The function Sis chosen to be a B-spline Sνof or-

der ν. The particle number density at the mesh

point rijk is obtained from the Ncmacroparticles

with position rl:

nijk =1

l=1

WlSν(rijk −rl) (13)

Similarly, the average particle density ﬂux is ob-

tained from particles velocity vl:

Φijk =

l=1

WlvlSν(rijk −rl) (14)

PHARE supports ﬁrst, second and third order B-

splines S1,S2and S3:

S1(x) = 1− |x|if |x| ≤ 1

0 elsewhere (15)

S2(x) = 





4−x2if |x| ≤ 1

23

2+|x|2if 1

2≤ |x| ≤ 3

0 elsewhere

(16)

S3(x) = 









|x|3

2−x2+ 2/3 if |x| ≤ 1

31−|x|

23if 1 ≤ |x| ≤ 2

0 elsewhere

(17)

Increasing the B-spline order uses more mesh

nodes. It thus decreases the noise level, at the price

of heavier computations and an increased diﬀusion

of gradients that may exist at a scale close to that

of the mesh.

2.3. Normalization of physical quantities

The code evolves dimensionless quantities that

are obtained with the following normalization. The

magnetic ﬁeld and the particle density are normal-

ized by arbitrary constants B0and n0, respectively.

Masses are normalized to the proton mass mpand

charges to the elementary charge e. It follows that

the time is normalized to the inverse proton gy-

rofrequency Ω0= (eB0/mp). Velocities are nor-

malized to the proton Alfv´en speed in the reference

magnetic ﬁeld and density VA0=B0/√mpn0µ0.

The distances are thus normalized to the ion iner-

tial length δi0=VA0/Ωci0

2.4. Spatial discretization

The spatial discretization of equations 1 to 8 in

PHARE follow the Yee layout [54], which preserves

the divergence-free character of the magnetic ﬁeld.

In each dimension, the components of the electric

and magnetic ﬁeld are positioned at integer or half

integer multiples of the mesh size, so-called primal

and dual positions, respectively. The 3D version of

the Yee lattice is represented on Fig. 1. By choice,

plasma moments are deﬁned on primal positions.

文档加载中……请稍候！
如果长时间未打开，您也可以点击刷新试试。

下载文档到电脑，查找使用更方便

10 玖币 0人已下载

立即下载

摘要：

PHARE:Parallelhybridparticle-in-cellcodewithpatch-basedadaptivemeshrefinementNicolasAunaia,∗,RochSmetsa,AndreaCiardib,PhilipDeegana,AlexisJeandeta,ThibaultPayetc,NathanGuyota,LoicDarrieumerlouaaLaboratoiredePhysiquedesPlasmas(LPP),CNRS,ObservatoiredeParis,SorbonneUniversit´e,Universit´eParis-Saclay,...

展开>> 收起<<

PHARE Parallel hybrid particle-in-cell code with patch-based adaptive mesh refinement.pdf

共21页,预览5页

还剩页未读，继续阅读

声明：本站为文档C2C交易模式，即用户上传的文档直接被用户下载，本站只是中间服务平台，本站所有文档下载所得的收益归上传人(含作者)所有。玖贝云文库仅提供信息存储空间，仅对用户上传内容的表现方式做保护处理，对上载内容本身不做任何修改或编辑。若文档所含内容侵犯了您的版权或隐私，请立即通知玖贝云文库，我们立即给予删除！

PHARE Parallel hybrid particle-in-cell code with patch-based adaptive mesh refinement

相关推荐

开通VIP享超值会员特权

作者详情

相关内容

热门标签

举报选择: