population fluid diagnostics not written

[nicolasaunai]

The following script does not produce any pop fluid diag h5 file whereas it should.

#!/usr/bin/env python3

import pyphare.pharein as ph #lgtm [py/import-and-import-from]
from pyphare.pharein import Simulation
from pyphare.pharein import MaxwellianFluidModel
from pyphare.pharein import ElectromagDiagnostics,FluidDiagnostics, ParticleDiagnostics, InfoDiagnostics
from pyphare.pharein import MetaDiagnostics
from pyphare.pharein import ElectronModel
from pyphare.simulator.simulator import Simulator, startMPI
from pyphare.pharein import global_vars as gv
from pyphare.pharein import LoadBalancer
import sys
import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
mpl.use('Agg')
from pyphare.cpp import cpp_lib
cpp = cpp_lib()
startMPI()



def config():

    start_time = 0.
    L=1.
    Simulation(
        time_step=0.001,
        final_time=1.,
        #boundary_types="periodic",
        cells=(320,128),
        dl=(0.40, 0.40),
        refinement="tagging",
        max_nbr_levels = 3,
        nesting_buffer=1,
        clustering="tile",
        tag_buffer="10",
	tagging_threshold=0.4,
        hyper_resistivity=0.0001,
	hyper_mode="spatial",
        resistivity=0.001,
        diag_options={"format": "phareh5",
                       "options": {"dir": "run079a",
                                  "mode":"overwrite"}},
        restart_options={"dir":"checkpoints", "mode":"overwrite", "elapsed_timestamps":[36000, 79000]}# ,"restart_time":start_time }
    )

    theta=0.
    nsh0 = 1
    nsp0 = 0.1
    nr = nsh0/nsp0
    theta=0.
    Bsp = 2
    Bsh = 1
    Br = Bsh/Bsp
    K = 1/Br**2

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


    def density_msh(x, y):
        from pyphare.pharein.global_vars import sim
        Ly = sim.simulation_domain()[1]
        n = S(y, 0.3*Ly, L) - S(y, 0.7*Ly, L)
        return n



    def by(x, y):
        from pyphare.pharein.global_vars import sim
        Lx = sim.simulation_domain()[0]
        Ly = sim.simulation_domain()[1]
        sigma = 1.
        dB = 0.1

        x0 = (x - 0.5 * Lx)
        y1 = (y - 0.3 * Ly)
        y2 = (y - 0.7 * Ly)

        dBy1 = 2*dB*x0 * np.exp(-(x0**2 + y1**2)/(sigma)**2)
        dBy2 = -2*dB*x0 * np.exp(-(x0**2 + y2**2)/(sigma)**2)

        return dBy1 + dBy2


    def bx(x, y):
        from pyphare.pharein.global_vars import sim
        Lx = sim.simulation_domain()[0]
        Ly = sim.simulation_domain()[1]
        sigma = 1.
        dB = 0.1

        x0 = (x - 0.5 * Lx)
        y1 = (y - 0.3 * Ly)
        y2 = (y - 0.7 * Ly)

        dBx1 = -2*dB*y1 * np.exp(-(x0**2 + y1**2)/(sigma)**2)
        dBx2 =  2*dB*y2 * np.exp(-(x0**2 + y2**2)/(sigma)**2)

        v1 = 2
        v2 = -1
        return v1 + (v2-v1)*S(y, Ly*0.3, L)   - (v2-v1)*S(y, Ly*0.7, L)+ dBx1 + dBx2


    def bz(x, y):
        return 0.


    def b2(x, y):
        return bx(x,y)**2 + by(x, y)**2 + bz(x, y)**2

    
    def Tish():
        K = 4
        return (K-0.5)/(1+theta)
    
    
    def Tisp():
        K = 4
        return nr/(1+theta)* (K - 1/(2*Br**2))
    
    
    def density_msp(x, y):
        from pyphare.pharein.global_vars import sim
        Ly = sim.simulation_domain()[1]
        n = 1/Tisp() * ( (K - b2(x,y) /2 )/(1+theta) - density_msh(x,y)*Tish())
        return n
    
    
    def density_total(x, y):
        return density_msh(x,y) + density_msp(x,y)
    

    def vx(x, y):
        return 0.


    def vy(x, y):
        return 0.


    def vz(x, y):
        return 0.


    def vth_sp_x(x, y):
        return np.sqrt(Tisp())

    def vth_sp_y(x, y):
        return np.sqrt(Tisp())

    def vth_sp_z(x, y):
        return np.sqrt(Tisp())

    def vth_sh_x(x, y):
        return np.sqrt(Tish())

    def vth_sh_y(x, y):
        return np.sqrt(Tish())

    def vth_sh_z(x, y):
        return np.sqrt(Tish())

    vvv_sp = {
        "vbulkx": vx, "vbulky": vy, "vbulkz": vz,
        "vthx": vth_sp_x, "vthy": vth_sp_y, "vthz": vth_sp_z,
        "nbr_part_per_cell":100
    }
    vvv_sh = {
        "vbulkx": vx, "vbulky": vy, "vbulkz": vz,
        "vthx": vth_sh_x, "vthy": vth_sh_y, "vthz": vth_sh_z,
        "nbr_part_per_cell":100
    }

    MaxwellianFluidModel(
        bx=bx, by=by, bz=bz,
        msp={"charge": 1, "density": density_msp,  **vvv_sp},
        msh={"charge": 1, "density": density_msh,  **vvv_sh}
    )

    ElectronModel(closure="isothermal", Te=0.0)
 
    LoadBalancer(active=True, mode="nppc", tol=0.05, every=1000)


    sim = ph.global_vars.sim
    dt =   500.*sim.time_step
    nt = (sim.final_time-start_time)/dt
    timestamps = start_time +dt * np.arange(nt)
    print(timestamps)
    print(f"Ion temperatures are Tish = {Tish()} and Tisp = {Tisp()}")


    for quantity in ["E", "B"]:
        ElectromagDiagnostics(
            quantity=quantity,
            write_timestamps=timestamps
        )

    for pop in ("msh", "msp"):
        FluidDiagnostics(
            quantity="density",
            write_timestamps=timestamps
            )



    for quantity in ["density", "bulkVelocity"]:
        FluidDiagnostics(
            quantity=quantity,
            write_timestamps=timestamps
            )

    InfoDiagnostics(quantity="particle_count", write_timestamps=timestamps)

def main():

    config()
    simulator = Simulator(gv.sim, print_one_line=False)
    simulator.initialize()
    simulator.run()


if __name__=="__main__":
    main()

[PhilipDeegan]

Something to note, this is I guess on your hyper branch #893