diff --git a/.gitignore b/.gitignore index deec3fc9b..69cd8144b 100644 --- a/.gitignore +++ b/.gitignore @@ -24,4 +24,4 @@ perf.* PHARE_REPORT.zip .cache .gdbinit -.phlop \ No newline at end of file +.phlop diff --git a/tests/simulator/refinement/CMakeLists.txt b/tests/simulator/refinement/CMakeLists.txt index 31772ce5d..b7f48be8f 100644 --- a/tests/simulator/refinement/CMakeLists.txt +++ b/tests/simulator/refinement/CMakeLists.txt @@ -17,5 +17,6 @@ if(HighFive) endif(testMPI) phare_python3_exec(9 simple_2d_refinement test_2d_2_core.py ${CMAKE_CURRENT_BINARY_DIR}) + phare_python3_exec(9 test_2d_overlaps test_2d_overlaps.py ${CMAKE_CURRENT_BINARY_DIR}) endif() diff --git a/tests/simulator/refinement/test_2d_overlaps.py b/tests/simulator/refinement/test_2d_overlaps.py new file mode 100644 index 000000000..bffe70654 --- /dev/null +++ b/tests/simulator/refinement/test_2d_overlaps.py @@ -0,0 +1,188 @@ +#!/usr/bin/env python3 + +import numpy as np +import matplotlib as mpl +from pathlib import Path + +import pyphare.pharein as ph +from pyphare.cpp import cpp_lib +from pyphare.simulator.simulator import Simulator, startMPI + +from tests.simulator import SimulatorTest +from tests.simulator.test_advance import AdvanceTestBase +from tests.diagnostic import dump_all_diags + +mpl.use("Agg") + +cpp = cpp_lib() + + +cells = (40, 80) +time_step = 0.005 +final_time = 1 +timestamps = np.arange(0, final_time + time_step, time_step) +diag_dir = "phare_outputs/test/refinenment/2d/overlaps" +test = AdvanceTestBase(rethrow=True) + + +def config(): + L = 0.5 + + sim = ph.Simulation( + time_step=time_step, + final_time=final_time, + largest_patch_size=10, + cells=cells, + dl=(0.40, 0.40), + refinement="tagging", + max_nbr_levels=2, + hyper_resistivity=0.002, + resistivity=0.001, + diag_options={ + "format": "phareh5", + "options": {"dir": diag_dir, "mode": "overwrite"}, + }, + strict=False, + ) + + def density(x, y): + Ly = sim.simulation_domain()[1] + return ( + 0.4 + + 1.0 / np.cosh((y - Ly * 0.3) / L) ** 2 + + 1.0 / np.cosh((y - Ly * 0.7) / L) ** 2 + ) + + def S(y, y0, l): + return 0.5 * (1.0 + np.tanh((y - y0) / l)) + + def by(x, y): + Lx = sim.simulation_domain()[0] + Ly = sim.simulation_domain()[1] + sigma = 1.0 + 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 0.0 # dBy1 + dBy2 + + def bx(x, y): + Lx = sim.simulation_domain()[0] + Ly = sim.simulation_domain()[1] + sigma = 1.0 + 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 = -1 + v2 = 1.0 + return v1 + (v2 - v1) * (S(y, Ly * 0.3, L) - S(y, Ly * 0.7, L)) # + dBx1 + dBx2 + + def bz(x, y): + return 0.0 + + def b2(x, y): + return bx(x, y) ** 2 + by(x, y) ** 2 + bz(x, y) ** 2 + + def T(x, y): + K = 0.7 + temp = 1.0 / density(x, y) * (K - b2(x, y) * 0.5) + assert np.all(temp > 0) + return temp + + def vx(x, y): + return 0.0 + + def vy(x, y): + return 0.0 + + def vz(x, y): + return 0.0 + + def vthx(x, y): + return np.sqrt(T(x, y)) + + def vthy(x, y): + return np.sqrt(T(x, y)) + + def vthz(x, y): + return np.sqrt(T(x, y)) + + vvv = { + "vbulkx": vx, + "vbulky": vy, + "vbulkz": vz, + "vthx": vthx, + "vthy": vthy, + "vthz": vthz, + "nbr_part_per_cell": 20, + } + + model = ph.MaxwellianFluidModel( + bx=bx, + by=by, + bz=bz, + protons={"charge": 1, "density": density, **vvv, "init": {"seed": 12334}}, + ) + ph.ElectronModel(closure="isothermal", Te=0.0) + + dump_all_diags(model.populations, flush_every=1) + + return sim + + +def get_time(path, time=None, datahier=None): + if time is not None: + time = "{:.10f}".format(time) + from pyphare.pharesee.hierarchy import hierarchy_from + + datahier = hierarchy_from(h5_filename=path + "/EM_E.h5", times=time, hier=datahier) + datahier = hierarchy_from(h5_filename=path + "/EM_B.h5", times=time, hier=datahier) + + return datahier + + +def get_hier(path): + return get_time(path) + + +class HarrisTest(SimulatorTest): + def __init__(self, *args, **kwargs): + super(HarrisTest, self).__init__(*args, **kwargs) + self.simulator = None + self.plot_dir = Path(f"{diag_dir}_plots") / str(cpp.mpi_size()) + self.plot_dir.mkdir(parents=True, exist_ok=True) + + def tearDown(self): + super(HarrisTest, self).tearDown() + if self.simulator is not None: + self.simulator.reset() + self.simulator = None + ph.global_vars.sim = None + + def test_run(self): + sim = Simulator(config(), post_advance=self.post_advance) + sim.initialize() + self.post_advance(0) + sim.run().reset() + return self + + def post_advance(self, new_time): + if cpp.mpi_rank() == 0: + L0L1_datahier = get_hier(diag_dir) + test.base_test_overlaped_fields_are_equal(L0L1_datahier, new_time) + + +if __name__ == "__main__": + startMPI() + HarrisTest().test_run().tearDown()