Modification of toy drift

examples/ToyGyrokinetic/diocotron_instability/params_diocotron.py

@@ -22,6 +22,11 @@
 # Import Struphy API
 # ------------------
 
+
+import logging
+from struphy import set_logging_level
+set_logging_level(logging.INFO)
+
 from struphy import (
     BaseUnits,
     DerhamOptions,
@@ -52,7 +57,14 @@
 # ---------------------
 
 from struphy.models import ToyDrift
-model = ToyDrift(epsilon=1.0)
+
+
+base_units = BaseUnits(kBT=1.0)
+model = ToyDrift(
+    epsilon=1.0,
+    alpha=1.0,
+    base_units=base_units,
+    )
 
 # List all variables and decide whether to save their data
 model.em_fields.phi.save_data = True
@@ -63,10 +75,10 @@
 # --------------------------
 
 # Environment options
-env = EnvironmentOptions(sim_folder="simdata")
+env = EnvironmentOptions(sim_folder="simdata",profiling_activated=True, profiling_trace=True)
 
 # Time stepping
-time_opts = Time(dt=0.05, Tend=20.0, split_algo="LieTrotter")
+time_opts = Time(dt=0.05, Tend=5.2, split_algo="LieTrotter")
 
 # Geometry
 domain = domains.HollowCylinder(a1=1.0, a2=10.0, Lz=10.0)
@@ -75,14 +87,15 @@
 equil = equils.HomogenSlab()
 
 # Grid
-grid = grids.TensorProductGrid(num_elements=(32,64,1), mpi_dims_mask=(False,True,False))
+grid = grids.TensorProductGrid(num_elements=(64,128,1), mpi_dims_mask=(False,True,False))
 
 # Derham options
 derham_opts = DerhamOptions(
     degree=(3,3,1), 
     bcs=(("dirichlet", "dirichlet"), None, None),
     )
 
+
 # Simulation object
 sim = Simulation(
     model=model,
@@ -101,15 +114,16 @@
 # Particle parameters
 # -------------------
 
-loading_params = LoadingParameters(ppc = 500, seed=1234)
-weights_params = WeightsParameters(control_variate=True)
+Np=200000
+loading_params = LoadingParameters(Np = Np, loading="sobol_standard", spatial="disc")
+weights_params = WeightsParameters(control_variate=True, reject_weights=True, threshold=0.00001)
 boundary_params = BoundaryParameters()
 model.kinetic_ions.set_markers(loading_params=loading_params,
                                weights_params=weights_params,
                                boundary_params=boundary_params,
                                bufsize=2.0,
                                )
-model.kinetic_ions.set_sorting_boxes(boxes_per_dim=(16,16,1), do_sort=True)
+model.kinetic_ions.set_sorting_boxes(boxes_per_dim=(24,24,1), do_sort=True)
 
 # density binning
 eta_bin = BinningPlot(slice='e1_e2', n_bins= (128,128), ranges= ((0.0, 1.0), (0.0,1.0)))
@@ -120,7 +134,7 @@
 # ------------------
 
 model.propagators.gc_poisson.options = model.propagators.gc_poisson.Options()
-model.propagators.push_gc_bxe.options = model.propagators.push_gc_bxe.Options(phi=model.em_fields.phi, evaluate_e_field=True)
+model.propagators.push_gc_bxe.options = model.propagators.push_gc_bxe.Options(algo="explicit", phi=model.em_fields.phi, evaluate_e_field=True)
 
 # ------------------
 # Initial conditions

examples/ToyGyrokinetic/diocotron_instability/pproc_diocotron.py

@@ -11,7 +11,8 @@
 # Post process simulation data
 # ------------------
 def main():
-    sim_path = os.path.join(os.getcwd(), "simdata")
+    sim_name = "simdata"
+    sim_path = os.path.join(os.getcwd(), sim_name)
 
     pp = PostProcessor(sim=params.sim)
     pp.process(physical=True)
@@ -49,10 +50,12 @@ def main():
     else:
         tf = 2*ti
     print(f"{ti = }, {tf = }")
+    #ti, tf = 2.5, 5.1
 
     xi = xp.abs(pdata.t_grid - ti).argmin() + 1 # index of time 100 [a.lu.] (observed end of growth rate)
     xf = xp.abs(pdata.t_grid - tf).argmin() + 1 # index of time 200 [a.lu.] (observed end of growth rate)
-
+    phi_init=en_phi[1]
+    en_phi = en_phi - phi_init
     fitting = xp.polyfit(time[xi:xf], xp.log10(en_phi[xi:xf]), deg=1)
 
     fig, ax = plt.subplots(1, figsize = (18, 12))
@@ -79,6 +82,7 @@ def main():
     # plt.savefig(os.path.join(save_path, "growth_rate.png"))
     # plt.close()
 
+    en_phi = en_phi + phi_init
 
     # ------------------
     # Show evolution of mass density distribution
@@ -208,6 +212,22 @@ def extract_images(bin_name, quantity, img_dir):
             plt.close(fig)
 
     # extract_images("e1_e2_density", "f_binned", os.path.join(save_path, "video"))
+    save_video_pngs = True
+    if save_video_pngs:
+        if not os.path.exists(sim_path+"/video"):
+            os.mkdir(sim_path+"/video")
+        # create .png for video
+        jump = 1
+        fig = plt.figure(figsize=(8, 8))
+        for n in range(ntime):
+            if n % jump == 0:
+                color_mapped = pdata.f.kinetic_ions.e1_e2_density.f_binned[n].T
+                plt.pcolor(phy_bin[0], phy_bin[1], pdata.f.kinetic_ions.e1_e2_density.f_binned[n])
+
+                plt.xlabel("x position")
+                plt.ylabel("y position")
+                plt.title(f"t = {pdata.t_grid[n]:4.2e}")
+                plt.savefig(sim_path+"/video"+f"/fig_{n:04.0f}.png", transparent=False, bbox_inches='tight', pad_inches=0)
 
 if __name__ == "__main__":
     main()

src/struphy/models/toy_drift.py

@@ -53,12 +53,14 @@ def __init__(
             charge_number: int = 1,
             mass_number: float = 1.0,
             epsilon: float = None,
+            alpha: float = None,
         ):
             self.var = PICVariable(space="Particles5D")
             self.init_variables(
                 charge_number=charge_number,
                 mass_number=mass_number,
                 epsilon=epsilon,
+                alpha=alpha,
             )
 
     ## propagators
@@ -76,6 +78,7 @@ def __init__(
         charge_number: int = 1,
         mass_number: float = 1.0,
         epsilon: float = None,
+        alpha: float = None,
     ):
 
         # 1. instantiate all species
@@ -84,6 +87,7 @@ def __init__(
             charge_number,
             mass_number,
             epsilon,
+            alpha,
         )
 
         # 2. derive units (must be done after instantiating species to access charge and mass numbers)

src/struphy/simulation/sim.py

@@ -816,24 +816,19 @@ def compute_plasma_params(self, verbose: bool = True):
         if verbose and MPI.COMM_WORLD.Get_rank() == 0:
             logger.info("\nPLASMA PARAMETERS:")
             logger.info(
-                "Plasma volume:".ljust(25),
-                "{:4.3e}".format(plasma_volume) + units_affix["plasma volume"],
+                "Plasma volume:".ljust(25) + "{:4.3e}".format(plasma_volume) + units_affix["plasma volume"],
             )
             logger.info(
-                "Transit length:".ljust(25),
-                "{:4.3e}".format(transit_length) + units_affix["transit length"],
+                "Transit length:".ljust(25) + "{:4.3e}".format(transit_length) + units_affix["transit length"],
             )
             logger.info(
-                "Avg. magnetic field:".ljust(25),
-                "{:4.3e}".format(magnetic_field) + units_affix["magnetic field"],
+                "Avg. magnetic field:".ljust(25) + "{:4.3e}".format(magnetic_field) + units_affix["magnetic field"],
             )
             logger.info(
-                "Max magnetic field:".ljust(25),
-                "{:4.3e}".format(B_max) + units_affix["magnetic field"],
+                "Max magnetic field:".ljust(25) + "{:4.3e}".format(B_max) + units_affix["magnetic field"],
             )
             logger.info(
-                "Min magnetic field:".ljust(25),
-                "{:4.3e}".format(B_min) + units_affix["magnetic field"],
+                "Min magnetic field:".ljust(25) + "{:4.3e}".format(B_min) + units_affix["magnetic field"],
             )
 
     def spawn_sister(