make harris make plots

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()