diff --git a/examples/ToyGyrokinetic/diocotron_instability/params_diocotron.py b/examples/ToyGyrokinetic/diocotron_instability/params_diocotron.py index a2d253951..bbd0fc2e2 100644 --- a/examples/ToyGyrokinetic/diocotron_instability/params_diocotron.py +++ b/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,7 +87,7 @@ 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( @@ -83,6 +95,7 @@ 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 diff --git a/examples/ToyGyrokinetic/diocotron_instability/pproc_diocotron.py b/examples/ToyGyrokinetic/diocotron_instability/pproc_diocotron.py index e99b17726..96c031330 100644 --- a/examples/ToyGyrokinetic/diocotron_instability/pproc_diocotron.py +++ b/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() \ No newline at end of file diff --git a/src/struphy/models/toy_drift.py b/src/struphy/models/toy_drift.py index f67862799..417a63d6d 100644 --- a/src/struphy/models/toy_drift.py +++ b/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) diff --git a/src/struphy/simulation/sim.py b/src/struphy/simulation/sim.py index 5e2ff11ec..3b5601f46 100644 --- a/src/struphy/simulation/sim.py +++ b/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( diff --git a/tutorials/tutorial_03_test_particles.ipynb b/tutorials/tutorial_03_test_particles.ipynb index 77c7e0d7c..ed3ab9c89 100644 --- a/tutorials/tutorial_03_test_particles.ipynb +++ b/tutorials/tutorial_03_test_particles.ipynb @@ -359,13 +359,90 @@ "ax.set_aspect(\"equal\")\n", "plt.xlabel(\"x\")\n", "plt.ylabel(\"y\")\n", - "plt.title(\"sim_2: draw uniform on disc\");" + "plt.title(\"sim_2: draw uniform on disc\")" ] }, { "cell_type": "markdown", "id": "27", "metadata": {}, + "source": [ + "Sometimes you need to load particles more uniformly. It is usefull to reach an homogeneous description of the plasma with a limited number of particles.\n", + "\n", + "It is performed using `loading=\"sobol_standard\"` or `loading=\"sobol_antithetic\"` in the `LoadingParameters` function:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "28", + "metadata": {}, + "outputs": [], + "source": [ + "env_3 = EnvironmentOptions(sim_folder=\"sim_3\")\n", + "model_3 = Vlasov()\n", + "sim_3 = Simulation.spawn_sister(sim, env=env_3, model=model_3)\n", + "\n", + "# species parameters\n", + "loading_params_3 = LoadingParameters(Np=1000, spatial=\"disc\", loading=\"sobol_standard\")\n", + "\n", + "model_3.kinetic_ions.set_markers(\n", + " loading_params=loading_params_3, weights_params=weights_params, boundary_params=boundary_params\n", + ")\n", + "model_3.kinetic_ions.set_sorting_boxes()\n", + "\n", + "model_3.kinetic_ions.set_save_data(n_markers=1.0)\n", + "\n", + "# propagator options\n", + "model_3.propagators.push_vxb.options = model_3.propagators.push_vxb.Options()\n", + "model_3.propagators.push_eta.options = model_3.propagators.push_eta.Options()\n", + "\n", + "# initial conditions (background + perturbation)\n", + "model_3.kinetic_ions.var.add_background(background)\n", + "\n", + "sim_3.run()\n", + "sim_3.pproc()\n", + "sim_3.load_plotting_data()\n", + "\n", + "fig = plt.figure(figsize=(15, 6))\n", + "\n", + "orbits_standard = sim_3.orbits.kinetic_ions\n", + "\n", + "plt.subplot(1, 3, 1)\n", + "plt.scatter(orbits[0, :, 0], orbits[0, :, 1], s=2.0)\n", + "circle1 = plt.Circle((0, 0), a2, color=\"k\", fill=False)\n", + "ax = plt.gca()\n", + "ax.add_patch(circle1)\n", + "ax.set_aspect(\"equal\")\n", + "plt.xlabel(\"x\")\n", + "plt.ylabel(\"y\")\n", + "plt.title(\"sim_1: draw uniform in logical space\")\n", + "\n", + "plt.subplot(1, 3, 2)\n", + "plt.scatter(orbits_uni[0, :, 0], orbits_uni[0, :, 1], s=2.0)\n", + "circle2 = plt.Circle((0, 0), a2, color=\"k\", fill=False)\n", + "ax = plt.gca()\n", + "ax.add_patch(circle2)\n", + "ax.set_aspect(\"equal\")\n", + "plt.xlabel(\"x\")\n", + "plt.ylabel(\"y\")\n", + "plt.title(\"sim_2: draw uniform on disc\")\n", + "\n", + "plt.subplot(1, 3, 3)\n", + "plt.scatter(orbits_standard[0, :, 0], orbits_standard[0, :, 1], s=2.0)\n", + "circle3 = plt.Circle((0, 0), a2, color=\"k\", fill=False)\n", + "ax = plt.gca()\n", + "ax.add_patch(circle3)\n", + "ax.set_aspect(\"equal\")\n", + "plt.xlabel(\"x\")\n", + "plt.ylabel(\"y\")\n", + "plt.title(\"sim_3: draw using sobol_standard\")" + ] + }, + { + "cell_type": "markdown", + "id": "29", + "metadata": {}, "source": [ "## Reflecting boundary conditions \n", "\n", @@ -376,7 +453,7 @@ { "cell_type": "code", "execution_count": null, - "id": "28", + "id": "30", "metadata": {}, "outputs": [], "source": [ @@ -388,7 +465,7 @@ { "cell_type": "code", "execution_count": null, - "id": "29", + "id": "31", "metadata": {}, "outputs": [], "source": [ @@ -407,7 +484,7 @@ { "cell_type": "code", "execution_count": null, - "id": "30", + "id": "32", "metadata": {}, "outputs": [], "source": [ @@ -420,7 +497,7 @@ }, { "cell_type": "markdown", - "id": "31", + "id": "33", "metadata": {}, "source": [ "We still have to set the propagator options and the initial conditions:" @@ -429,7 +506,7 @@ { "cell_type": "code", "execution_count": null, - "id": "32", + "id": "34", "metadata": {}, "outputs": [], "source": [ @@ -446,7 +523,7 @@ }, { "cell_type": "markdown", - "id": "33", + "id": "35", "metadata": {}, "source": [ "We can now run the simulation, then post-process the data and plot the resulting orbits:" @@ -455,7 +532,7 @@ { "cell_type": "code", "execution_count": null, - "id": "34", + "id": "36", "metadata": {}, "outputs": [], "source": [ @@ -465,7 +542,7 @@ { "cell_type": "code", "execution_count": null, - "id": "35", + "id": "37", "metadata": {}, "outputs": [], "source": [ @@ -475,7 +552,7 @@ { "cell_type": "code", "execution_count": null, - "id": "36", + "id": "38", "metadata": {}, "outputs": [], "source": [ @@ -484,7 +561,7 @@ }, { "cell_type": "markdown", - "id": "37", + "id": "39", "metadata": {}, "source": [ "Under `sim.orbits[]` one finds a three-dimensional numpy array; the first index refers to the time step, the second index to the particle and the third index to the particle attribute. The first three attributes are the particle positions, followed by the velocities and the (initial and time-dependent) weights." @@ -493,7 +570,7 @@ { "cell_type": "code", "execution_count": null, - "id": "38", + "id": "40", "metadata": {}, "outputs": [], "source": [ @@ -507,7 +584,7 @@ { "cell_type": "code", "execution_count": null, - "id": "39", + "id": "41", "metadata": {}, "outputs": [], "source": [ @@ -537,7 +614,7 @@ }, { "cell_type": "markdown", - "id": "40", + "id": "42", "metadata": {}, "source": [ "## Particles in a cylinder with a magnetic field\n", @@ -561,7 +638,7 @@ { "cell_type": "code", "execution_count": null, - "id": "41", + "id": "43", "metadata": {}, "outputs": [], "source": [ @@ -573,7 +650,7 @@ }, { "cell_type": "markdown", - "id": "42", + "id": "44", "metadata": {}, "source": [ "In order to project the equilibrium on the spline basis for fast evaluation in the particle kernels, we need a Derham complex:" @@ -582,7 +659,7 @@ { "cell_type": "code", "execution_count": null, - "id": "43", + "id": "45", "metadata": {}, "outputs": [], "source": [ @@ -592,7 +669,7 @@ }, { "cell_type": "markdown", - "id": "44", + "id": "46", "metadata": {}, "source": [ "Now we create the light-weight instance of the model and set the species options. We shall `remove` particles that hit the boundary in $\\eta_1$ (radial) direction:" @@ -601,7 +678,7 @@ { "cell_type": "code", "execution_count": null, - "id": "45", + "id": "47", "metadata": {}, "outputs": [], "source": [ @@ -616,7 +693,7 @@ { "cell_type": "code", "execution_count": null, - "id": "46", + "id": "48", "metadata": {}, "outputs": [], "source": [ @@ -641,7 +718,7 @@ }, { "cell_type": "markdown", - "id": "47", + "id": "49", "metadata": {}, "source": [ "Now the usual procedure: run, post-process, load data and finally plot the orbits:" @@ -650,7 +727,7 @@ { "cell_type": "code", "execution_count": null, - "id": "48", + "id": "50", "metadata": {}, "outputs": [], "source": [ @@ -660,7 +737,7 @@ { "cell_type": "code", "execution_count": null, - "id": "49", + "id": "51", "metadata": {}, "outputs": [], "source": [ @@ -670,7 +747,7 @@ { "cell_type": "code", "execution_count": null, - "id": "50", + "id": "52", "metadata": {}, "outputs": [], "source": [ @@ -680,7 +757,7 @@ { "cell_type": "code", "execution_count": null, - "id": "51", + "id": "53", "metadata": {}, "outputs": [], "source": [ @@ -693,7 +770,7 @@ { "cell_type": "code", "execution_count": null, - "id": "52", + "id": "54", "metadata": {}, "outputs": [], "source": [ @@ -721,7 +798,7 @@ }, { "cell_type": "markdown", - "id": "53", + "id": "55", "metadata": {}, "source": [ "## Particles in a Tokamak equilibrium\n", @@ -732,7 +809,7 @@ { "cell_type": "code", "execution_count": null, - "id": "54", + "id": "56", "metadata": {}, "outputs": [], "source": [ @@ -744,7 +821,7 @@ }, { "cell_type": "markdown", - "id": "55", + "id": "57", "metadata": {}, "source": [ "Since `EQDSKequilibrium` is an `AxisymmMHDequilibrium`, which is a subclass of `CartesianMHDequilibrium`, we are free to choose any mapping for the simulation (e.g. a Cuboid for Cartesian coordinates). In order to be conforming to the boundary of the equilibrium, we shall choose the `Tokamak` mapping:" @@ -753,7 +830,7 @@ { "cell_type": "code", "execution_count": null, - "id": "56", + "id": "58", "metadata": {}, "outputs": [], "source": [ @@ -766,7 +843,7 @@ }, { "cell_type": "markdown", - "id": "57", + "id": "59", "metadata": {}, "source": [ "The `Tokamak` domain is a `PoloidalSplineTorus`, hence\n", @@ -817,7 +894,7 @@ { "cell_type": "code", "execution_count": null, - "id": "58", + "id": "60", "metadata": {}, "outputs": [], "source": [ @@ -835,7 +912,7 @@ { "cell_type": "code", "execution_count": null, - "id": "59", + "id": "61", "metadata": {}, "outputs": [], "source": [ @@ -850,7 +927,7 @@ { "cell_type": "code", "execution_count": null, - "id": "60", + "id": "62", "metadata": {}, "outputs": [], "source": [ @@ -867,7 +944,7 @@ { "cell_type": "code", "execution_count": null, - "id": "61", + "id": "63", "metadata": {}, "outputs": [], "source": [ @@ -924,7 +1001,7 @@ }, { "cell_type": "markdown", - "id": "62", + "id": "64", "metadata": {}, "source": [ "We need a Derham complex for the projection of the equilibrium onto the spline basis:" @@ -933,7 +1010,7 @@ { "cell_type": "code", "execution_count": null, - "id": "63", + "id": "65", "metadata": {}, "outputs": [], "source": [ @@ -947,7 +1024,7 @@ }, { "cell_type": "markdown", - "id": "64", + "id": "66", "metadata": {}, "source": [ "We aim to simulate 15000 time steps with a second-order splitting algorithm:" @@ -956,7 +1033,7 @@ { "cell_type": "code", "execution_count": null, - "id": "65", + "id": "67", "metadata": {}, "outputs": [], "source": [ @@ -965,7 +1042,7 @@ }, { "cell_type": "markdown", - "id": "66", + "id": "68", "metadata": {}, "source": [ "We now set up a simulation of 4 specific particle orbits in this equilibrium:" @@ -974,7 +1051,7 @@ { "cell_type": "code", "execution_count": null, - "id": "67", + "id": "69", "metadata": {}, "outputs": [], "source": [ @@ -993,7 +1070,7 @@ { "cell_type": "code", "execution_count": null, - "id": "68", + "id": "70", "metadata": {}, "outputs": [], "source": [ @@ -1027,7 +1104,7 @@ { "cell_type": "code", "execution_count": null, - "id": "69", + "id": "71", "metadata": {}, "outputs": [], "source": [ @@ -1037,7 +1114,7 @@ { "cell_type": "code", "execution_count": null, - "id": "70", + "id": "72", "metadata": {}, "outputs": [], "source": [ @@ -1047,7 +1124,7 @@ { "cell_type": "code", "execution_count": null, - "id": "71", + "id": "73", "metadata": {}, "outputs": [], "source": [ @@ -1057,7 +1134,7 @@ { "cell_type": "code", "execution_count": null, - "id": "72", + "id": "74", "metadata": {}, "outputs": [], "source": [ @@ -1070,7 +1147,7 @@ { "cell_type": "code", "execution_count": null, - "id": "73", + "id": "75", "metadata": {}, "outputs": [], "source": [ @@ -1095,7 +1172,7 @@ }, { "cell_type": "markdown", - "id": "74", + "id": "76", "metadata": {}, "source": [ "## Guiding-centers in a Tokamak equilibrium\n", @@ -1106,7 +1183,7 @@ { "cell_type": "code", "execution_count": null, - "id": "75", + "id": "77", "metadata": {}, "outputs": [], "source": [ @@ -1119,7 +1196,7 @@ { "cell_type": "code", "execution_count": null, - "id": "76", + "id": "78", "metadata": {}, "outputs": [], "source": [ @@ -1131,7 +1208,7 @@ { "cell_type": "code", "execution_count": null, - "id": "77", + "id": "79", "metadata": {}, "outputs": [], "source": [ @@ -1165,7 +1242,7 @@ { "cell_type": "code", "execution_count": null, - "id": "78", + "id": "80", "metadata": {}, "outputs": [], "source": [ @@ -1223,7 +1300,7 @@ { "cell_type": "code", "execution_count": null, - "id": "79", + "id": "81", "metadata": {}, "outputs": [], "source": [ @@ -1233,7 +1310,7 @@ { "cell_type": "code", "execution_count": null, - "id": "80", + "id": "82", "metadata": {}, "outputs": [], "source": [ @@ -1243,7 +1320,7 @@ { "cell_type": "code", "execution_count": null, - "id": "81", + "id": "83", "metadata": {}, "outputs": [], "source": [ @@ -1253,7 +1330,7 @@ { "cell_type": "code", "execution_count": null, - "id": "82", + "id": "84", "metadata": {}, "outputs": [], "source": [ @@ -1266,7 +1343,7 @@ { "cell_type": "code", "execution_count": null, - "id": "83", + "id": "85", "metadata": {}, "outputs": [], "source": [