Simulation inputs

Python initialization blocks

PHARE is initialized from Python objects. These objects are declared by users in a script to define simulation parameters, and to describe the initial conditionn.

The following shows the various objects to declare in a pseudo-script:

# ------------ MANDATORY BLOCKS

Simulation(
    # some parameters
    # configuring numerical and AMR
    # parameters
)


MawellianFluidModel(
   # some parameters
   # configuring the magnetic field profile
   # and ion initial condition
   # as fluid moment profiles
   # ion particles are assumed to follow
   # locally Maxwellian distributions with these moments
   )

# ------------ END OF MANDATORY BLOCKS


# ------------ OPTIONAL BLOCKS

ElectronModel(
    # configures electron fluid properties
)



ElectromagDiagnostics(
   # parameters configuring outputs
   # of E and B
)


FluidDiagnostics(
    # parameters configuring ion moment outputs
)


ParticleDiagnostics(
    # some parameters configuring particle outputs
)


# ------------ END OF OPTIONAL BLOCKS

The Simulation block

The Simulation block is used to set general parameters of the simulation like the integration time, domain size, interpolation order, adaptive meshing, restart and diagnostcs options. The Simulation must be the first block defined in an input script

class pyphare.pharein.Simulation(**kwargs_in)[source]

Usage example:

This declares a 2D simulation of 100x100 cells with an isotropic mesh size of 0.2. The simulation will run 1000 time steps of dt=0.001 and so will stop at t=1

The simulation will run with adaptive mesh refinement and evolve up to 3 levels, i.e. a base coarsest mesh with up to 2 additional refinement levels.

Diagnostics, if declared, will be saved as native PHARE HDF5 (default format) in the directory diag_outputs and will overwrite any existing h5 files there.

The resistivity and hyper-resistivity are set to constant custom values.

from pyphare.pharein import Simulation



Simulation(
           # time evolution parameters
           time_step_nbr=1000,
           time_step=0.001,

           # Domain parameters
           cells=(100,100),
           dl=(0.2, 0.2),

           # AMR parameters
           max_nbr_levels=3,

            # general physics parameters
           hyper_resistivity=0.001,
           resistivity=0.001,

           # diagnostics and restart parameters
           diag_options={"format": "phareh5",
                         "options": {"dir": diag_outputs,
                                     "mode":"overwrite"}},
          )

Time parameters:

final_time (float), final simulation time. Use with time_step OR time_step_nbr

time_step (float), simulation time step. Use with time_step_nbr OR final_time

time_step_nbr (int), number of time step to perform. Use with final_time OR time_step
Simulation(
           time_step_nbr=1000,
           time_step=0.001,

        # ....
        )
is equivalent to:
Simulation(
           time_step_nbr=1000,
           final_time=1.,

        # ....
        )
or even to:
Simulation(
           time_step=0.001,
           final_time=1.,

        # ....
        )

Domain/grid parameters:

The number of dimensions of the simulation is deduced from the length of these parameters.

dl (float, tuple) grid spacing dx, (dx,dy) or (dx,dy,dz) in 1D, 2D or 3D. A single float value in 2D or 3D are assumed equal for each dimension. Use this parameter with either domain_size OR cells

domain_size (float or tuple) size of the physical domain Lx, (Lx,Ly), (Lx,Ly,Lz) in 1D, 2D or 3D Single float is assumed equal for each direction. Use this parameter with either cells or dl

cells (int or tuple) number of cells nx or (nx, ny) or (nx, ny, nz) in 1, 2 and 3D. Single int is assumed equal for each direction Use this parameter with either domain_size or dl

Expert parameters:

These parameters are more advanced, modify them at your own risk

layout (str), layout of the physical quantities on the mesh (default = “yee”)

For instance:

Simulation(
           dl=(0.1, 0.1),
           domain_size=(100,100),

        # ....
        )

is equivalent to:

Simulation(
           cells=(1000, 1000),
           domain_size=(100,100),

        # ....
        )

or even to:

Simulation(
           cells=(1000, 1000),
           dl=(0.1,0.1),

        # ....
        )

Macro-particle parameters:

interp_order (int), 1, 2 or 3 (default=1) particle b-spline order

particle_pusher (str), algo to push particles (default = “modifiedBoris”)

Diagnostics output parameters:

diag_options (dict)

path (str) path for outputs (default : ‘./’)

diag_export_format (str) format of the output diagnostics (default= “phareh5”)

mode (str) mode of the output diagnostics (default= “overwrite” will write over existing files)

Adaptive Mesh Refinement (AMR) parameters:

max_nbr_levels (int), default=1, max number of levels in the hierarchy. Used if no refinement_boxes are set

tag_buffer (int), default=1, value representing the number of cells by which tagged cells are buffered before clustering into boxes. The larger tag_buffer, the wider refined regions will be around tagged cells.

clustering (str), {“berger”, “tile” (default)}, type of clustering to use for AMR. tile results in wider patches, less artifacts and better scalability

Expert parameters:

These parameters are more advanced, modify them at your own risk

refined_particle_nbr (ìnt), number of refined particle per coarse particle.

nesting_buffer (ìnt), default=0 minimum gap in coarse cells from the border of a level and any refined patch border

refinement_boxes, default=None, {“L0”:{“B0”:[(lox,loy,loz),(upx,upy,upz)],…,”Bi”:[(),()]},…”Li”:{B0:[(),()]}}

smallest_patch_size (int or tuple), minimum number of cells in a patch in each direction This parameter cannot be smaller than the number of field ghost nodes

largest_patch_size (int or tuple), maximum size of a patch in each direction

Restart parameters:

These parameters are used to restart a simulation from a previous state and dump chekpoints.

restart_options (dict)

dir (str) path for restart files (default : ‘./’)

mode (str) mode of the restart files

“conserve” - (default), will conserve existing files

“overwrite” - will overwrite existing files

restart_time (float) time at which to restart the simulation (default=0)

timestamps (list) list of timestamps at which to restart the simulation

Misc:

description (string), [default=None] arbitrary string for per simulation context - injected in output files when feasible

strict (bool), turns warnings into errors (default False)

resistivity (float), resistivity value (default=0.0)

hyper-resistivity (float), hyper-resistivity value (default=0.0)

boundary_types (str or tuple) type of boundary conditions (default is “periodic” for each direction)

Magnetic field and ions

The typical way to initialize a Hybrid-PIC model is to parametrize the different ion populations with fluid moments, assuming an underlying velocity distributotion function. The MaxwellianFluidModel block allows just that, assuming a Maxwellian distribution for each population. It also allows to set the magnetic field profile.

Below is an example of a MaxwellianFluidModel block for which a single population of protons is initialized.

ph.MaxwellianFluidModel(
    bx=bx,
    by=by,
    bz=bz,
    protons={"charge": 1,
             "density": density,
             "vbulkx": vx,
             "vbulky": vy,
             "vbulkz": vz,
             "vthx": vthx,
             "vthy": vthy,
             "vthz": vthz,
             "nbr_part_per_cell": 100}
)

In the example above, bx, by, bz are the components of the magnetic field, density is the density of the protons, vx, vy, vz are the bulk velocities and vthx, vthy, vthz are the thermal velocities of the protons. All these parameters are (previously defined) functions of the spatial coordinates.

As an example, the function below defines a density profile. This function will be called by the C++ code to load the population to which this density is assigned. In this example, we create two current sheets at y = 0.3 and y = 0.7 of the simulation domain along the y direction, of half-width 0.5.

# assume the simulation is created and accessed by the variable `sim`

import numpy as np

def density(x, y):

    # L[1] is the length of the simulation domain in the y direction
    L = sim.simulation_domain()[1]
    return (
        0.2
        + 1.0 / np.cosh((y - L * 0.3) / 0.5) ** 2
        + 1.0 / np.cosh((y - L * 0.7) / 0.5) ** 2
    )

The magnetic field function could be defined as follows:

def S(y, y0, l):
    return 0.5 * (1.0 + np.tanh((y - y0) / l))

def bx(x, y):
    Ly = sim.simulation_domain()[1]
    v1 = -1.0
    v2 =  1.0
    return (
        v1
        + (v2 - v1) * (S(y, Ly * 0.3, 0.5) - S(y, Ly * 0.7, 0.5))
    )

Note here that the function bx uses the other function S. This underlines the great power that comes with initializing the simulation with Python since the initialization script can be as complex as needed.

Adding a population is simple. In the example below a beam proton poulation is added with more particles per cells and different moments.:

ph.MaxwellianFluidModel(
    bx=bx,
    by=by,
    bz=bz,
    protons={"charge": 1,
             "density": density,
             "vbulkx": vx,
             "vbulky": vy,
             "vbulkz": vz,
             "vthx": vthx,
             "vthy": vthy,
             "vthz": vthz,
             "nbr_part_per_cell": 100},

    beam={"charge": 1,
             "density": beam_density,
             "vbulkx": vx_beam,
             "vbulky": vy_beam,
             "vbulkz": vz_beam,
             "vthx": vthx_beam,
             "vthy": vthy_beam,
             "vthz": vthz_beam,
             "nbr_part_per_cell": 500}
)

Examples

See tests/functional/harris/harris_2d.py for a complete example with one population.
See tests/functional/ionIon/beam_ions.py for a complete example with two populations.

Details

See the Maxwellian Fluid Model class for more details.

Electron model

The ElectronModel block is used to set the electron fluid properties. Below is an example of an ElectronModel block:

from pyphare.pharein import ElectronModel

ElectronModel(closure="isothermal", Te=0.2)

For now, the only closure available is the isothermal closure, the given temperatureis thus constantin space and time.

Diagnostics

Diagnostics blocks are used to set parameters of the different diagnostics output PHARE produces. There are three types of diagnostics blocks:

ElectromagDiagnostics for electromagnetic field outputs
FluidDiagnostics for ion moment outputs
ParticleDiagnostics for particle outputs

Electromagnetic Diagnostics

Electromagnetic diagnostics are used to write either the electric or magnetic field to disk. Below is an example of an ElectromagDiagnostics block.

from pyphare.pharein import ElectromagDiagnostics

time_step_nbr = 1000
time_step = 0.001
final_time = time_step * time_step_nbr
dt = 10 * time_step
nt = final_time / dt + 1
timestamps = dt * np.arange(nt)

ElectromagDiagnostics(quantity="E", write_timestamps=timestamps)
ElectromagDiagnostics(quantity="B", write_timestamps=timestamps)

In this example, the electric and magnetic field will be written to disk every 10 time steps. Note that although here the timestamps are equally spaced, they do not have to be. The write_timestamps parameter can be any list of times. Times are provided in simulation units.

Fluid Diagnostics

Fluid quantities represent the moments of the ion populations. Quantities available are of two kinds:

Total ion quantities:

These diagnostics are properties of the ion populations as a whole.

density: represents the total particle density of the ions
mass density: represents the total mass density of the ions
bulk velocity: represents the total bulk velocity of the ions
pressure_tensory: represent the total pressure tensor of the ions

Ion population quantities:

These diagnostics are properties of each ion population. The name of the population must be provided.

flux: represents the flux of a given population
momentum_tensor: represents the momentum tensor of a given population

Example

from pyphare.pharein import FluidDiagnostics
from pyphare.pharein import MaxwellianFluidModel

# assume bx, by, ... are defined
# ...
MaxwellianFluidModel(
    bx=bx,
    by=by,
    bz=bz,
    protons={"charge": 1,
             "density": density,
             "vbulkx": vx,
             "vbulky": vy,
             "vbulkz": vz,
             "vthx": vthx,
             "vthy": vthy,
             "vthz": vthz,
             "nbr_part_per_cell": 100}
)

time_step_nbr = 1000
time_step = 0.001
final_time = time_step * time_step_nbr
dt = 10 * time_step
nt = final_time / dt + 1
timestamps = dt * np.arange(nt)

FluidDiagnostics(
    quantity="density",
    write_timestamps=timestamps
)

FluidDiagnostics(
    quantity="flux",
    write_timestamps=timestamps,
    population_name="protons"
)

Particle Diagnostics

These diagnostics are used to write particle data to disk. They are typically much heavier than any other diagnostics.

The block below declares a particle diagnostics so that “protons” are written to disk at the given timestamps. Note the quantity=”domain” parameter. This is used to write all the particles living within the interior of the simulation patches.

from pyphare.pharein import ParticleDiagnostics

time_step_nbr = 1000
time_step = 0.001
final_time = time_step * time_step_nbr
dt = 10 * time_step
nt = final_time / dt + 1
timestamps = dt * np.arange(nt)

ph.ParticleDiagnostics(
        quantity="domain",
        population_name="protons",
        write_timestamps=timestamps
    )