Cyclone test case

doc/sections/subsections/propagators-avail.rst

@@ -11,6 +11,11 @@ Field solvers
     :exclude-members: options, allocate
     :show-inheritance:
 
+.. automodule:: struphy.propagators.poisson_adiabatic_gyrokinetic
+    :members:
+    :exclude-members: options, allocate
+    :show-inheritance:
+
 .. automodule:: struphy.propagators.curl_curl_solve
     :members:
     :exclude-members: options, allocate

examples/ToyGyrokinetic/diocotron_instability/params_diocotron.py

@@ -70,17 +70,17 @@
 
 # List all variables and decide whether to save their data
 model.em_fields.phi.save_data = True
-model.kinetic_ions.var.save_data = True
+model.kinetic_ions.var.save_data = False
 
 # --------------------------
 # Instance of the simulation
 # --------------------------
 
 # Environment options
-env = EnvironmentOptions(sim_folder="simdata",profiling_activated=True, profiling_trace=True)
+env = EnvironmentOptions(sim_folder="sim_1",profiling_activated=True, profiling_trace=True, restart=False)
 
 # Time stepping
-time_opts = Time(dt=0.05, Tend=5.2, split_algo="LieTrotter")
+time_opts = Time(dt=0.01, Tend=51.0, split_algo="LieTrotter")
 
 # Geometry
 domain = domains.HollowCylinder(a1=1.0, a2=10.0, Lz=10.0)
@@ -116,15 +116,15 @@
 # Particle parameters
 # -------------------
 
-Np=200000
+Np=500000
 loading_params = LoadingParameters(Np = Np, loading="sobol_standard", spatial="disc")
-weights_params = WeightsParameters(control_variate=True, reject_weights=True, threshold=0.00001)
+weights_params = WeightsParameters(control_variate=True, reject_weights=True, threshold=0.0001)
 boundary_params = BoundaryParameters()
-sorting_params = SortingParameters(boxes_per_dim=(24,24,1), do_sort=True)
+sorting_params = SortingParameters(boxes_per_dim=(12,12,1), do_sort=True, sorting_frequency=5)
 
 # density binning
 eta_bin = BinningPlot(slice='e1_e2', n_bins= (128,128), ranges= ((0.0, 1.0), (0.0,1.0)))
-saving_params = SavingParameters(binning_plots=(eta_bin, ))
+saving_params = SavingParameters(binning_plots=(eta_bin,)) # (binning_plots=(eta_bin,)) if you want to save the density binning data
 
 model.kinetic_ions.set_markers(loading_params=loading_params,
                                weights_params=weights_params,
@@ -139,7 +139,7 @@
 # ------------------
 
 model.propagators.gc_poisson.options = model.propagators.gc_poisson.Options()
-model.propagators.push_gc_bxe.options = model.propagators.push_gc_bxe.Options(algo="explicit", phi=model.em_fields.phi, evaluate_e_field=True)
+model.propagators.push_gc_bxe.options = model.propagators.push_gc_bxe.Options(algo="discrete_gradient_1st_order_newton", phi=model.em_fields.phi, evaluate_e_field=True)
 
 # ------------------
 # Initial conditions
@@ -172,8 +172,8 @@ def n_init(etas,r_minus=r_minus,r_plus=r_plus):
 
 # Perturbations for (some) kinetic species
 
-# for linear case amps = (1e-6,)
-perturbation = perturbations.ModesCos(amps=(0.5,), ms=(ms,), perb_domain=((eta_minus,eta_plus), None, None))
+# for non linear case amps = (0.5,)
+perturbation = perturbations.ModesCos(amps=(1e-6,), ms=(ms,), perb_domain=((eta_minus,eta_plus), None, None))
 init = maxwellians.GyroMaxwellian2D(n=(n_init, perturbation), equil=equil)
 model.kinetic_ions.var.add_initial_condition(init)
 

examples/ToyGyrokinetic/diocotron_instability/pproc_diocotron.py

@@ -1,76 +1,100 @@
-import params_diocotron as params
+import importlib.util
 from struphy import PlottingData, PostProcessor
 
 import os
+import sys
 import cunumpy as xp
+import scipy.optimize as sc
 from matplotlib import pyplot as plt
 import h5py
 
+import logging
+from struphy import set_logging_level
+set_logging_level(logging.INFO)
+
 
 # ------------------
 # Post process simulation data
+# In order to compare different simulations, execute this file as `python pproc_diocotron.py sim_1 sim_2 ...` 
+# where `sim_1`, `sim_2`, etc. are the names of the simulation folders to be post-processed and plotted together.
+# If only one argument, the 2D plots will be shown. If multiple arguments, only the growth rate plot will be shown.
 # ------------------
 def main():
-    sim_name = "simdata"
-    sim_path = os.path.join(os.getcwd(), sim_name)
-
-    pp = PostProcessor(sim=params.sim)
-    pp.process(physical=True)
-
-    pdata = PlottingData(sim=params.sim)
-    pdata.load()
-
-    # path to save plots
-    # save_path = os.path.join(os.getcwd(), "images", "sim")
-    # os.makedirs(save_path, exist_ok=True)
-
-    # ------------------
-    # Check simulation domain
-    # ------------------
-
-    params.domain.show()
-
-    # ------------------
-    # Determine electrical potentail growth rate
-    # ------------------
-
-    # get scalar data (post processing not needed for scalar data)
-    pa_data = os.path.join(sim_path, "data")
-    with h5py.File(os.path.join(pa_data, "data_proc0.hdf5"), "r") as f:
-        time = f["time"]["value"][()]
-        en_phi = f["scalar"]["en_phi"][()]
-
-    # determine growth rate
-    exp_func = lambda x,m,b: 10**(m*x+b)
-
-    # time interval to determine growth rate
-    ti = pdata.t_grid[-1]//4 
-    if ti == 0.0:
-        tf = pdata.t_grid[-1]
+    if __name__ == "main":
+        if len(sys.argv)>1:
+            sim_names = sys.argv[1:]
+        else:
+            sim_names = ["sim_1"]
     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)
+        sim_names = ["sim_1"]
+    en_phis = []
+    times = []
+    sls = []
+    params_opts = []
+    fitting = []
+    for i, sim_name in enumerate(sim_names):
+        sim_path = os.path.join(os.getcwd(), sim_name)
+
+        spec = importlib.util.spec_from_file_location("params", os.path.join(sim_path, "parameters.py"))
+        params = importlib.util.module_from_spec(spec)
+        spec.loader.exec_module(params)
+
+        if not os.path.isdir(os.path.join(sim_path, "post_processing")):
+            pp = PostProcessor(sim=params.sim)
+            pp.process(physical=True)
+
+        pdata = PlottingData(sim=params.sim)
+        pdata.load()
+
+        logging.warning(f"arriving at post processing with these sim_names : {sim_names}, knowing that {__name__=}")
+        # ------------------
+        # Determine electrical potentail growth rate
+        # ------------------
+
+        # get scalar data (post processing not needed for scalar data)
+        pa_data = os.path.join(sim_path, "data")
+        with h5py.File(os.path.join(pa_data, "data_proc0.hdf5"), "r") as f:
+            times.append(f["time"]["value"][()])
+            en_phis.append(xp.power(f["scalar"]["en_phi"][()], 1.0))
+
+        # time interval to determine growth rate
+        ti, tf = 0.0, 42.0
+        if tf>times[i][-1]: tf = times[i][-1]
+        if ti>tf:
+            ti = tf/2
+        xi = xp.abs(pdata.t_grid - ti).argmin() # 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)
+        if xi==0:
+            xi=1 # avoid including t=0 in fit
+
+        sls.append(tuple([slice(xi, xf)]))
+
+        if len(times) > 3:
+            fitting.append(True)
+            # determine growth rate
+            fitting_func = lambda x,m,b,c0: xp.exp(m*x+b)+c0
+            jac_func = lambda x,m,b,c0: xp.array([x*xp.exp(m*x+b), xp.exp(m*x+b), xp.ones_like(x)]).transpose()
+
+            params_opt, _ = sc.curve_fit(fitting_func, times[i][sls[i]], en_phis[i][sls[i]], p0=(1e-3, -5, en_phis[i][1]), jac=jac_func, maxfev=10000)#3.07e2
+            params_opts.append(params_opt)
+
+            logging.info(f"Fitted growth rate for {sim_name}: {params_opt[0]:.4e}")
+        else:
+            fitting.append(False)
 
     fig, ax = plt.subplots(1, figsize = (18, 12))
-
-    # plot
-    ax.plot(time, en_phi, label=r"$\phi$")
-    ax.plot(
-        pdata.t_grid, 
-        exp_func(pdata.t_grid, *fitting), 
-        label=f"fitted growth rate {ti=}, {tf=}, growth_rate={fitting[0]:.4e}"
-    )
+    for i in range(len(sim_names)):
+        ax.scatter(times[i][1:], en_phis[i][1:], label=r"$\phi_{"+sim_names[i][4:]+r"}$", marker='x', s=0.05)
+        if fitting[i]:
+            ax.plot(
+                times[i][sls[i]], 
+                fitting_func(times[i][sls[i]], *params_opts[i]), 
+                label=f"fitted growth rate {ti=}, {tf=}, growth_rate={params_opts[i][0]:.4e}, b={params_opts[i][1]:.4e}, c0={params_opts[i][2]:.4e}"
+            )
     ax.axvline(ti, color="gray", linestyle="--", alpha=0.5)
     ax.axvline(tf, color="gray", linestyle="--", alpha=0.5)
 
-    ax.set_yscale('log')
+    #ax.set_yscale('log')
     ax.legend()
 
     ax.set_title(f"{params.time_opts.dt=}, {params.time_opts.split_algo=}, {params.grid.num_elements=}, {params.derham_opts.degree=}, {params.loading_params.ppc=}")
@@ -79,11 +103,8 @@ def main():
 
     plt.tight_layout()
     plt.show()
-    # plt.savefig(os.path.join(save_path, "growth_rate.png"))
-    # plt.close()
-
-    en_phi = en_phi + phi_init
-
+    if len(sim_names)>1:
+        exit()
     # ------------------
     # Show evolution of mass density distribution
     # ------------------
@@ -212,7 +233,7 @@ 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
+    save_video_pngs = False
     if save_video_pngs:
         if not os.path.exists(sim_path+"/video"):
             os.mkdir(sim_path+"/video")

src/struphy/fields_background/equils.py

@@ -871,6 +871,9 @@ class AdhocTorus(AxisymmMHDequilibrium):
         &q_0 + ( q_1 - q_0 )\frac{r^2}{a^2} \quad &&\textnormal{if} \quad q_\textnormal{kind}=0\,,
 
         &\frac{q_0}{1-\left(1-\frac{r^2}{a^2}\right)^{\frac{q_1}{q_0}}}\frac{r^2}{a^2} \quad &&\textnormal{if} \quad q_\textnormal{kind}=1\,.
+
+
+        &q_0 + l\frac{r}{a} ( q_1 - q_0 )\frac{r^2}{a^2} \quad &&\textnormal{if} \quad q_\textnormal{kind}=2\,,
         \end{aligned}\right.
 
     The pressure profile
@@ -903,6 +906,8 @@ class AdhocTorus(AxisymmMHDequilibrium):
         Safety factor at r=0 (default: 1.71).
     q1 : float
         Safety factor at r=a (default: 1.87).
+    l : float
+        Linear term factor for q profile if q_kind=2 (default: 0.).
     n1 : float
         1st shape factor for ion number density profile (default: 0.).
     n2 : float
@@ -932,9 +937,10 @@ class AdhocTorus(AxisymmMHDequilibrium):
             a       : 1.   # minor radius
             R0      : 3.   # major radius
             B0      : 2.   # on-axis toroidal magnetic field
-            q_kind  : 0    # which profile (0 : parabolic, 1 : other, see documentation)
+            q_kind  : 0    # which profile (0 : parabolic, 1 : other, 2 : parabolic with linear term, see documentation)
             q0      : 1.05 # safety factor at r=0
             q1      : 1.80 # safety factor at r=a
+            l       : 0.   # linear term factor for q profile if q_kind=2
             n1      : .5   # 1st shape factor for number density profile
             n2      : 1.   # 2nd shape factor for number density profile
             na      : .2   # number density at r=a
@@ -955,6 +961,7 @@ def __init__(
         q_kind: int = 0,
         q0: float = 1.71,
         q1: float = 1.87,
+        l: float = 0.0,
         n1: float = 2.0,
         n2: float = 1.0,
         na: float = 0.2,
@@ -983,7 +990,7 @@ def __init__(
             self._psi_i = None
             self._p_i = None
 
-        else:
+        elif self.params["q_kind"] == 1 or self.params["q_kind"] == 2:
             r_i = xp.linspace(0.0, self.params["a"], self.params["psi_nel"] + 1)
 
             def dpsi_dr(r):
@@ -1088,7 +1095,7 @@ def psi_r(self, r, der=0):
                 out = self.params["B0"] * (q_bar_0 - r * q_bar_1) / q_bar_0**2
 
         # alternative profile (interpolated)
-        else:
+        elif self.params["q_kind"] == 1 or self.params["q_kind"] == 2:
             out = self._psi_i(r, nu=der)
 
             # remove all "dimensions" for point-wise evaluation
@@ -1105,6 +1112,7 @@ def q_r(self, r, der=0):
 
         q0 = self.params["q0"]
         q1 = self.params["q1"]
+        l = self.params["l"]
 
         a = self.params["a"]
 
@@ -1115,8 +1123,14 @@ def q_r(self, r, der=0):
             else:
                 qout = 2 * (q1 - q0) * r / a**2
 
+        elif self.params["q_kind"] == 2:
+            if der == 0:
+                qout = q0 + (q1 - q0) * (r / a) ** 2 + l * r / a
+            else:
+                qout = 2 * (q1 - q0) * r / a**2 + l / a
+
         # alternative profile
-        else:
+        elif self.params["q_kind"] == 1:
             # int/float input
             if isinstance(r, (int, float)):
                 if r == 0:
@@ -1210,7 +1224,7 @@ def p_r(self, r):
                     )
 
             # alternative profile
-            else:
+            elif self.params["q_kind"] == 1:
                 pout = self._p_i(r)
 
                 # remove all "dimensions" for point-wise evaluation
@@ -1219,7 +1233,7 @@ def p_r(self, r):
                     pout = pout.item()
 
         # ad-hoc profile
-        else:
+        elif self.params["p_kind"] == 1:
             pout = (
                 self.params["B0"] ** 2
                 * self.params["beta"]