diff --git a/tests/functional/harris/harris_2d.py b/tests/functional/harris/harris_2d.py index ae54b4c9e..bcfd5559c 100644 --- a/tests/functional/harris/harris_2d.py +++ b/tests/functional/harris/harris_2d.py @@ -1,64 +1,54 @@ #!/usr/bin/env python3 -import os +import os 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 -from pyphare.simulator.simulator import startMPI - -os.environ["PHARE_SCOPE_TIMING"] = "1" # turn on scope timing -""" - For scope timings to work - The env var PHARE_SCOPE_TIMING must be == "1" (or "true") - See src/phare/phare.hpp - CMake must be configured with: -DwithPhlop=ON - And a LOG_LEVEL must be defined via compile args: -DPHARE_LOG_LEVEL=1 - Or change the default value in src/core/logger.hpp - And phlop must be available on PYTHONPATH either from subprojects - or install phlop via pip -""" - - -ph.NO_GUI() +from pyphare.pharesee.run import Run +from pyphare.simulator.simulator import Simulator, startMPI + +from tests.simulator import SimulatorTest + +mpl.use("Agg") + cpp = cpp_lib() -startMPI() -diag_outputs = "phare_outputs/test/harris/2d" -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) + +cells = (200, 100) +time_step = 0.005 +final_time = 50 +timestamps = np.arange(0, final_time + time_step, final_time / 5) +diag_dir = "phare_outputs/harris" def config(): + L = 0.5 + sim = ph.Simulation( - smallest_patch_size=15, - largest_patch_size=25, - time_step_nbr=time_step_nbr, time_step=time_step, - # boundary_types="periodic", - cells=(200, 400), - dl=(0.2, 0.2), + final_time=final_time, + cells=cells, + dl=(0.40, 0.40), refinement="tagging", - max_nbr_levels=1, - hyper_resistivity=0.001, + max_nbr_levels=2, + hyper_resistivity=0.002, resistivity=0.001, diag_options={ "format": "phareh5", - "options": {"dir": diag_outputs, "mode": "overwrite"}, + "options": {"dir": diag_dir, "mode": "overwrite"}, }, strict=True, ) def density(x, y): - L = sim.simulation_domain()[1] + Ly = 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 + 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): @@ -67,35 +57,34 @@ def S(y, y0, l): def by(x, y): Lx = sim.simulation_domain()[0] Ly = sim.simulation_domain()[1] - w1 = 0.2 - w2 = 1.0 + sigma = 1.0 + dB = 0.1 + x0 = x - 0.5 * Lx y1 = y - 0.3 * Ly y2 = y - 0.7 * Ly - w3 = np.exp(-(x0 * x0 + y1 * y1) / (w2 * w2)) - w4 = np.exp(-(x0 * x0 + y2 * y2) / (w2 * w2)) - w5 = 2.0 * w1 / w2 - return (w5 * x0 * w3) + (-w5 * x0 * w4) + + 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): Lx = sim.simulation_domain()[0] Ly = sim.simulation_domain()[1] - w1 = 0.2 - w2 = 1.0 + sigma = 1.0 + dB = 0.1 + x0 = x - 0.5 * Lx y1 = y - 0.3 * Ly y2 = y - 0.7 * Ly - w3 = np.exp(-(x0 * x0 + y1 * y1) / (w2 * w2)) - w4 = np.exp(-(x0 * x0 + y2 * y2) / (w2 * w2)) - w5 = 2.0 * w1 / w2 + + 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, 0.5) - S(y, Ly * 0.7, 0.5)) - + (-w5 * y1 * w3) - + (+w5 * y2 * w4) - ) + 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 @@ -104,7 +93,7 @@ def b2(x, y): return bx(x, y) ** 2 + by(x, y) ** 2 + bz(x, y) ** 2 def T(x, y): - K = 1 + K = 0.7 temp = 1.0 / density(x, y) * (K - b2(x, y) * 0.5) assert np.all(temp > 0) return temp @@ -143,26 +132,80 @@ def vthz(x, y): bz=bz, protons={"charge": 1, "density": density, **vvv, "init": {"seed": 12334}}, ) - ph.ElectronModel(closure="isothermal", Te=0.0) for quantity in ["E", "B"]: ph.ElectromagDiagnostics(quantity=quantity, write_timestamps=timestamps) - ph.InfoDiagnostics(quantity="particle_count") # defaults all coarse time steps + for quantity in ["density", "bulkVelocity"]: + ph.FluidDiagnostics(quantity=quantity, write_timestamps=timestamps) + + ph.FluidDiagnostics( + quantity="density", write_timestamps=timestamps, population_name="protons" + ) + ph.InfoDiagnostics(quantity="particle_count") + + ph.LoadBalancer(active=True, auto=True, mode="nppc", tol=0.05) return sim -def main(): - Simulator(config()).run() - try: - from tools.python3 import plotting as m_plotting +def plot_file_for_qty(plot_dir, qty, time): + return f"{plot_dir}/harris_{qty}_t{time}.png" + + +def plot(diag_dir, plot_dir): + run = Run(diag_dir) + for time in timestamps: + run.GetDivB(time).plot( + filename=plot_file_for_qty(plot_dir, "divb", time), + plot_patches=True, + vmin=1e-11, + vmax=2e-10, + ) + run.GetRanks(time).plot( + filename=plot_file_for_qty(plot_dir, "Ranks", time), plot_patches=True + ) + run.GetN(time, pop_name="protons").plot( + filename=plot_file_for_qty(plot_dir, "N", time), plot_patches=True + ) + for c in ["x", "y", "z"]: + run.GetB(time).plot( + filename=plot_file_for_qty(plot_dir, f"b{c}", time), + plot_patches=True, + qty=f"{c}", + ) + run.GetJ(time).plot( + filename=plot_file_for_qty(plot_dir, "jz", time), + qty="z", + plot_patches=True, + vmin=-2, + vmax=2, + ) + + +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 - m_plotting.plot_run_timer_data(diag_outputs, cpp.mpi_rank()) - except ImportError: - print("Phlop not found - install with: `pip install phlop`") - cpp.mpi_barrier() + def test_run(self): + self.register_diag_dir_for_cleanup(diag_dir) + Simulator(config()).run().reset() + if cpp.mpi_rank() == 0: + plot(diag_dir, self.plot_dir) + cpp.mpi_barrier() + return self if __name__ == "__main__": - main() + startMPI() + HarrisTest().test_run().tearDown() diff --git a/tests/functional/harris/harris_2d_lb.py b/tests/functional/harris/harris_2d_lb.py deleted file mode 100644 index 1153b6e8b..000000000 --- a/tests/functional/harris/harris_2d_lb.py +++ /dev/null @@ -1,211 +0,0 @@ -#!/usr/bin/env python3 - -import os -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.pharesee.run import Run -from pyphare.simulator.simulator import Simulator, startMPI - -from tests.simulator import SimulatorTest -from tools.python3 import plotting as m_plotting - -mpl.use("Agg") - -SCOPE_TIMING = os.getenv("PHARE_SCOPE_TIMING", "True").lower() in ("true", "1", "t") -LOAD_BALANCE = os.getenv("LOAD_BALANCE", "True").lower() in ("true", "1", "t") - -cpp = cpp_lib() -startMPI() - -cells = (800, 800) -time_step = 0.005 -final_time = 50 -timestamps = np.arange(0, final_time + time_step, final_time / 5) - -if cpp.mpi_rank() == 0: - print(LOAD_BALANCE, "diag timestamps:", timestamps) - -diag_dir = "phare_outputs/harris_lb" -if not LOAD_BALANCE: - diag_dir = "phare_outputs/harris" - -plot_dir = Path(f"{diag_dir}_plots") -plot_dir.mkdir(parents=True, exist_ok=True) - - -def config(): - L = 0.5 - - sim = ph.Simulation( - time_step=time_step, - final_time=final_time, - cells=cells, - dl=(0.40, 0.40), - refinement="tagging", - max_nbr_levels=2, - nesting_buffer=1, - clustering="tile", - tag_buffer="1", - hyper_resistivity=0.002, - resistivity=0.001, - diag_options={ - "format": "phareh5", - "options": {"dir": diag_dir, "mode": "overwrite"}, - }, - restart_options={ - "dir": "checkpoints", - "mode": "overwrite", - "timestamps": timestamps, - # "restart_time": 0.0, - }, - ) - - 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 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 vxyz(x, y): - return 0.0 - - def vthxyz(x, y): - return np.sqrt(T(x, y)) - - vvv = {**{f"vbulk{c}": vxyz for c in "xyz"}, **{f"vth{c}": vthxyz for c in "xyz"}} - - ph.MaxwellianFluidModel( - bx=bx, by=by, bz=bz, protons={"charge": 1, "density": density, **vvv} - ) - ph.ElectronModel(closure="isothermal", Te=0.0) - - for quantity in ["E", "B"]: - ph.ElectromagDiagnostics(quantity=quantity, write_timestamps=timestamps) - for quantity in ["density", "bulkVelocity"]: - ph.FluidDiagnostics(quantity=quantity, write_timestamps=timestamps) - - ph.FluidDiagnostics( - quantity="density", write_timestamps=timestamps, population_name="protons" - ) - ph.InfoDiagnostics(quantity="particle_count") - - if LOAD_BALANCE: - ph.LoadBalancer(active=True, auto=True, mode="nppc", tol=0.05) - - return sim - - -def plot_file_for_qty(qty, time): - return f"{plot_dir}/harris_{qty}_t{time}.png" - - -def plot(diag_dir): - run = Run(diag_dir) - for time in timestamps: - run.GetDivB(time).plot( - filename=plot_file_for_qty("divb", time), - plot_patches=True, - vmin=1e-11, - vmax=2e-10, - ) - run.GetRanks(time).plot( - filename=plot_file_for_qty("Ranks", time), - plot_patches=True, - ) - run.GetN(time, pop_name="protons").plot( - filename=plot_file_for_qty("N", time), - plot_patches=True, - ) - for c in ["x", "y", "z"]: - run.GetB(time).plot( - filename=plot_file_for_qty(f"b{c}", time), - qty=f"{c}", - plot_patches=True, - ) - run.GetJ(time).plot( - filename=plot_file_for_qty("jz", time), - qty="z", - plot_patches=True, - vmin=-2, - vmax=2, - ) - - -class HarrisTest(SimulatorTest): - def __init__(self, *args, **kwargs): - super(HarrisTest, self).__init__(*args, **kwargs) - self.simulator = None - - 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): - self.register_diag_dir_for_cleanup(diag_dir) - Simulator(config()).run().reset() - if cpp.mpi_rank() == 0: - plot(diag_dir) - if SCOPE_TIMING: - m_plotting.plot_run_timer_data(diag_dir, cpp.mpi_rank()) - cpp.mpi_barrier() - return self - - -if __name__ == "__main__": - HarrisTest().test_run().tearDown()