From b47602506d2bba67798f9bb57f68e8c480971f5c Mon Sep 17 00:00:00 2001 From: bknaadmin Date: Mon, 8 Dec 2025 13:26:02 +0100 Subject: [PATCH] upload parameter.py file for ITPA benchmark simulations --- ITPA_benchmark/gc.py | 145 ++++++++++++ ITPA_benchmark/linearmhd.py | 126 +++++++++++ ITPA_benchmark/tae_2n64_NC2_debug.py | 206 +++++++++++++++++ ITPA_benchmark/tae_2n64_debug.py | 205 +++++++++++++++++ ITPA_benchmark/tae_hybrid.py | 207 ++++++++++++++++++ ITPA_benchmark/tae_hybrid_nonlinear.py | 207 ++++++++++++++++++ .../tae_hybrid_nonlinear_sonic_dt0125.py | 206 +++++++++++++++++ ...ae_hybrid_nonlinear_sonic_gradp0_dt0125.py | 206 +++++++++++++++++ .../tae_hybrid_nonlinear_sonic_m2jboff.py | 206 +++++++++++++++++ ...e_hybrid_nonlinear_sonic_m2jboff_dt0125.py | 206 +++++++++++++++++ ITPA_benchmark/tae_hybrid_sonic.py | 204 +++++++++++++++++ ITPA_benchmark/tae_hybrid_sonic_dt0125.py | 206 +++++++++++++++++ ITPA_benchmark/tae_hybrid_sonic_gradp0.py | 205 +++++++++++++++++ .../tae_hybrid_sonic_gradp0_dt0125.py | 206 +++++++++++++++++ ITPA_benchmark/tae_hybrid_sonic_m2jboff.py | 206 +++++++++++++++++ .../tae_hybrid_sonic_m2jboff_dt0125.py | 206 +++++++++++++++++ ITPA_benchmark/tae_tp.py | 207 ++++++++++++++++++ ITPA_benchmark/tae_tp_nonlinear.py | 207 ++++++++++++++++++ .../tae_tp_nonlinear_sonic_dt0125.py | 206 +++++++++++++++++ .../tae_tp_nonlinear_sonic_gradp0_dt0125.py | 206 +++++++++++++++++ .../tae_tp_nonlinear_sonic_m2jboff.py | 206 +++++++++++++++++ .../tae_tp_nonlinear_sonic_m2jboff_dt0125.py | 206 +++++++++++++++++ ITPA_benchmark/tae_tp_sonic.py | 204 +++++++++++++++++ ITPA_benchmark/tae_tp_sonic_dt0125.py | 206 +++++++++++++++++ ITPA_benchmark/tae_tp_sonic_gradp0.py | 205 +++++++++++++++++ ITPA_benchmark/tae_tp_sonic_gradp0_dt0125.py | 206 +++++++++++++++++ ITPA_benchmark/tae_tp_sonic_m2jboff.py | 206 +++++++++++++++++ ITPA_benchmark/tae_tp_sonic_m2jboff_dt0125.py | 206 +++++++++++++++++ 28 files changed, 5624 insertions(+) create mode 100644 ITPA_benchmark/gc.py create mode 100644 ITPA_benchmark/linearmhd.py create mode 100644 ITPA_benchmark/tae_2n64_NC2_debug.py create mode 100644 ITPA_benchmark/tae_2n64_debug.py create mode 100644 ITPA_benchmark/tae_hybrid.py create mode 100644 ITPA_benchmark/tae_hybrid_nonlinear.py create mode 100644 ITPA_benchmark/tae_hybrid_nonlinear_sonic_dt0125.py create mode 100644 ITPA_benchmark/tae_hybrid_nonlinear_sonic_gradp0_dt0125.py create mode 100644 ITPA_benchmark/tae_hybrid_nonlinear_sonic_m2jboff.py create mode 100644 ITPA_benchmark/tae_hybrid_nonlinear_sonic_m2jboff_dt0125.py create mode 100644 ITPA_benchmark/tae_hybrid_sonic.py create mode 100644 ITPA_benchmark/tae_hybrid_sonic_dt0125.py create mode 100644 ITPA_benchmark/tae_hybrid_sonic_gradp0.py create mode 100644 ITPA_benchmark/tae_hybrid_sonic_gradp0_dt0125.py create mode 100644 ITPA_benchmark/tae_hybrid_sonic_m2jboff.py create mode 100644 ITPA_benchmark/tae_hybrid_sonic_m2jboff_dt0125.py create mode 100644 ITPA_benchmark/tae_tp.py create mode 100644 ITPA_benchmark/tae_tp_nonlinear.py create mode 100644 ITPA_benchmark/tae_tp_nonlinear_sonic_dt0125.py create mode 100644 ITPA_benchmark/tae_tp_nonlinear_sonic_gradp0_dt0125.py create mode 100644 ITPA_benchmark/tae_tp_nonlinear_sonic_m2jboff.py create mode 100644 ITPA_benchmark/tae_tp_nonlinear_sonic_m2jboff_dt0125.py create mode 100644 ITPA_benchmark/tae_tp_sonic.py create mode 100644 ITPA_benchmark/tae_tp_sonic_dt0125.py create mode 100644 ITPA_benchmark/tae_tp_sonic_gradp0.py create mode 100644 ITPA_benchmark/tae_tp_sonic_gradp0_dt0125.py create mode 100644 ITPA_benchmark/tae_tp_sonic_m2jboff.py create mode 100644 ITPA_benchmark/tae_tp_sonic_m2jboff_dt0125.py diff --git a/ITPA_benchmark/gc.py b/ITPA_benchmark/gc.py new file mode 100644 index 0000000..f80c408 --- /dev/null +++ b/ITPA_benchmark/gc.py @@ -0,0 +1,145 @@ +from struphy.io.options import EnvironmentOptions, BaseUnits, Time +from struphy.geometry import domains +from struphy.fields_background import equils +from struphy.topology import grids +from struphy.io.options import DerhamOptions +from struphy.io.options import FieldsBackground +from struphy.initial import perturbations +from struphy.kinetic_background import maxwellians +from struphy.pic.utilities import (LoadingParameters, + WeightsParameters, + BoundaryParameters, + BinningPlot, + KernelDensityPlot, + ) +from struphy import main + +# import model, set verbosity +from struphy.models.toy import GuidingCenter + +# environment options +env = EnvironmentOptions( + out_folders='/tokp/work/bkna/tae', + sim_folder='gc', + save_step=8, max_runtime=1440) + +# units +base_units = BaseUnits(n=0.2) + +# time stepping +time_opts = Time(dt=0.25, Tend=200., split_algo='Strang') + +# geometry +domain = domains.HollowTorus(a1=0.1, a2=1., R0=10., sfl=False, tor_period=6) + +# fluid equilibrium (can be used as part of initial conditions) +equil = equils.AdhocTorus( + a=1., + R0=10., + B0=3., + q_kind=0, + q0=1.71, + q1=1.87, + n1=0., + n2=0., + na=1., + p_kind=1., + p0=1., + p1=0.95, + p2=0.05, + beta=0.0018, +) + +# grid +grid = grids.TensorProductGrid( + Nel=(24,96,16), + mpi_dims_mask=(True, True, False) +) + +# derham options +derham_opts = DerhamOptions( + p=(3, 3, 3), + spl_kind=(False, True, True), + dirichlet_bc=((True, True), (False, False), (False, False)), + nq_pr=(6, 6, 6), + nquads=(6, 6, 6), + polar_ck=-1, +) + +# light-weight model instance +model = GuidingCenter() + +# species parameters +model.kinetic_ions.set_phys_params( + charge_number=2., +) + +loading_params = LoadingParameters( + ppc=100, + moments=(0., 0., 0.89740839, 0.89740839) + ) +weights_params = WeightsParameters( + control_variate=True, +) +boundary_params = BoundaryParameters( + bc=("remove", "periodic", "periodic"), + bc_refill=("outer", "inner"), +) +model.kinetic_ions.set_markers(loading_params=loading_params, + weights_params=weights_params, + boundary_params=boundary_params, + bufsize=1., + n_cols_diag=3, + n_cols_aux=5, + ) +model.kinetic_ions.set_sorting_boxes(boxes_per_dim=None) +binning_plots = ( + BinningPlot( + slice="e1", + n_bins=128, + ranges=(0., 1.), + divide_by_jac=True), + BinningPlot( + slice="v1", + n_bins=128, + ranges=(-2.5, 2.5), + divide_by_jac=True), + ) +model.kinetic_ions.set_save_data( + binning_plots=binning_plots, +) + +# propagator options +model.propagators.push_bxe.options = model.propagators.push_bxe.Options( + algo = 'explicit',) +model.propagators.push_parallel.options = model.propagators.push_parallel.Options( + algo = 'explicit',) +# background, perturbations and initial conditions +n_profile = perturbations.ITPA_density( + n0=0.00720652, + c=(0.47023, 0.20323, 0.13273, 0.521298) +) +maxwellian_1 = maxwellians.CanonicalMaxwellian2D( + equil=equil, + volume_form=True, + n=(0., n_profile), +) + +background = maxwellian_1 +model.kinetic_ions.var.add_background(background, n_as_volume_form=False) + +if __name__ == "__main__": + # start run + verbose = True + + main.run(model, + params_path=__file__, + env=env, + base_units=base_units, + time_opts=time_opts, + domain=domain, + equil=equil, + grid=grid, + derham_opts=derham_opts, + verbose=verbose, + ) \ No newline at end of file diff --git a/ITPA_benchmark/linearmhd.py b/ITPA_benchmark/linearmhd.py new file mode 100644 index 0000000..17d934e --- /dev/null +++ b/ITPA_benchmark/linearmhd.py @@ -0,0 +1,126 @@ +from struphy.io.options import EnvironmentOptions, BaseUnits, Time +from struphy.geometry import domains +from struphy.fields_background import equils +from struphy.topology import grids +from struphy.io.options import DerhamOptions +from struphy.io.options import FieldsBackground +from struphy.initial import perturbations +from struphy.kinetic_background import maxwellians +from struphy.pic.utilities import (LoadingParameters, + WeightsParameters, + BoundaryParameters, + BinningPlot, + KernelDensityPlot, + ) +from struphy import main + +# import model, set verbosity +from struphy.models.fluid import LinearMHD + +# environment options +env = EnvironmentOptions( + out_folders='/tokp/work/bkna/tae', + sim_folder='linearmhd', + save_step=8, max_runtime=1440) + + +# units +base_units = BaseUnits(n=0.2) + +# time stepping +time_opts = Time(dt=0.25, Tend=200., split_algo='Strang') + +# geometry +domain = domains.HollowTorus(a1=0.1, a2=1., R0=10., sfl=False, tor_period=6) + +# fluid equilibrium (can be used as part of initial conditions) +equil = equils.AdhocTorus( + a=1., + R0=10., + B0=3., + q_kind=0, + q0=1.71, + q1=1.87, + n1=0., + n2=0., + na=1., + p_kind=1., + p0=1., + p1=0.95, + p2=0.05, + beta=0.0018, +) + +# grid +grid = grids.TensorProductGrid( + Nel=(24,96,16), + mpi_dims_mask=(True, True, False) +) + +# derham options +derham_opts = DerhamOptions( + p=(3, 3, 3), + spl_kind=(False, True, True), + dirichlet_bc=((True, True), (False, False), (False, False)), + nq_pr=(6, 6, 6), + nquads=(6, 6, 6), + polar_ck=-1, +) + +# light-weight model instance +model = LinearMHD() + +# species parameters +model.mhd.set_phys_params() + +# propagator options +model.propagators.shear_alf.options = model.propagators.shear_alf.Options( + algo='implicit', + precond='MassMatrixDiagonalPreconditioner', +) +model.propagators.mag_sonic.options = model.propagators.mag_sonic.Options(b_field=model.em_fields.b_field) + +# background, perturbations and initial conditions +model.mhd.velocity.add_background(FieldsBackground( + type="LogicalConst", + values=(0., 0., 0.) +)) +model.mhd.velocity.add_perturbation( + perturbations.TorusModesSin( + given_in_basis="2", + amps=(0.01,0.01), + ms=(10,11), + ns=(-1,-1), + pfuns=('exp','exp'), + pfun_params=((0.44444444,0.1),(0.44444444,0.1)), + comp=0) +) +model.mhd.velocity.add_perturbation( + perturbations.TorusModesCos( + given_in_basis="2", + amps=(0.00015915,0.00015915), + ms=(10,11), + ns=(-1,-1), + pfuns=('d_exp','d_exp'), + pfun_params=((0.44444444,0.1),(0.44444444,0.1)), + comp=1), +) + +# optional: exclude variables from saving +# model.mhd.pressure.save_data = False + +if __name__ == "__main__": + # start run + verbose = True + + main.run(model, + params_path=__file__, + env=env, + base_units=base_units, + time_opts=time_opts, + domain=domain, + equil=equil, + grid=grid, + derham_opts=derham_opts, + verbose=verbose, + ) \ No newline at end of file diff --git a/ITPA_benchmark/tae_2n64_NC2_debug.py b/ITPA_benchmark/tae_2n64_NC2_debug.py new file mode 100644 index 0000000..4373488 --- /dev/null +++ b/ITPA_benchmark/tae_2n64_NC2_debug.py @@ -0,0 +1,206 @@ +from struphy.io.options import EnvironmentOptions, BaseUnits, Time +from struphy.geometry import domains +from struphy.fields_background import equils +from struphy.topology import grids +from struphy.io.options import DerhamOptions +from struphy.io.options import FieldsBackground +from struphy.initial import perturbations +from struphy.kinetic_background import maxwellians +from struphy.pic.utilities import (LoadingParameters, + WeightsParameters, + BoundaryParameters, + BinningPlot, + KernelDensityPlot, + ) +from struphy import main + +# import model, set verbosity +from struphy.models.hybrid import LinearMHDDriftkineticCC + +# environment options +env = EnvironmentOptions( + out_folders='/tokp/work/bkna/tae', + sim_folder='tae_2n64_debug', + save_step=1, max_runtime=110, + num_clones=2, + restart=False) + +# units +base_units = BaseUnits(n=0.2) + +# time stepping +time_opts = Time(dt=0.25, Tend=2000., split_algo='Strang') + +# geometry +domain = domains.HollowTorus(a1=0.1, a2=1., R0=10., sfl=False, tor_period=6) + +# fluid equilibrium (can be used as part of initial conditions) +equil = equils.AdhocTorus( + a=1., + R0=10., + B0=3., + q_kind=0, + q0=1.71, + q1=1.87, + n1=0., + n2=0., + na=1., + p_kind=1., + p0=1., + p1=0.95, + p2=0.05, + beta=0.0018, +) + +# grid +grid = grids.TensorProductGrid( + Nel=(24,96,16), + mpi_dims_mask=(True, True, False) +) + +# derham options +derham_opts = DerhamOptions( + p=(3, 3, 3), + spl_kind=(False, True, True), + dirichlet_bc=((True, True), (False, False), (False, False)), + nq_pr=(6, 6, 6), + nquads=(6, 6, 6), + polar_ck=-1, +) + +# light-weight model instance +model = LinearMHDDriftkineticCC( + turn_off = (#'PushGuidingCenterBxEstar', + #'PushGuidingCenterParallel', + #'ShearAlfvenCurrentCoupling5D', + 'Magnetosonic', + #'CurrentCoupling5DDensity', + #'CurrentCoupling5DGradB', + #'CurrentCoupling5DCurlb', + ) +) + +# species parameters +model.mhd.set_phys_params() +model.energetic_ions.set_phys_params( + mass_number=2., +) + +loading_params = LoadingParameters( + ppc=100, + moments=(0., 0., 0.89740839, 0.89740839), + seed=1234, + ) +weights_params = WeightsParameters( + control_variate=True, +) +boundary_params = BoundaryParameters( + bc=("remove", "periodic", "periodic"), + bc_refill=("outer", "inner"), +) +model.energetic_ions.set_markers(loading_params=loading_params, + weights_params=weights_params, + boundary_params=boundary_params, + bufsize=1., + n_cols_diag=3, + n_cols_aux=5, + ) +model.energetic_ions.set_sorting_boxes(boxes_per_dim=None) +model.energetic_ions.set_save_data() + +# calculate constant +ep_scale = model.energetic_ions.var.species.mass_number / model.mhd.mass_number + +# params else +from struphy.pic.accumulation.filter import FilterParameters +from struphy.linear_algebra.solver import SolverParameters, DiscreteGradientSolverParameters +filter_params = FilterParameters(use_filter="hybrid",) +solver_params = SolverParameters() + +# propagator options +model.propagators.push_bxe.options = model.propagators.push_bxe.Options( + algo = 'explicit', + b_tilde = model.em_fields.b_field,) +model.propagators.push_parallel.options = model.propagators.push_parallel.Options( + algo = 'explicit', + b_tilde = model.em_fields.b_field,) +model.propagators.shearalfen_cc5d.options = model.propagators.shearalfen_cc5d.Options( + filter_params=filter_params, + ep_scale=ep_scale, + energetic_ions = model.energetic_ions.var, + solver_params=solver_params, + nonlinear=False,) +# model.propagators.magnetosonic.options = model.propagators.magnetosonic.Options( +# b_field=model.em_fields.b_field,) +model.propagators.cc5d_density.options = model.propagators.cc5d_density.Options( + filter_params=filter_params, + ep_scale=ep_scale, + energetic_ions = model.energetic_ions.var, + b_tilde = model.em_fields.b_field,) +model.propagators.cc5d_gradb.options = model.propagators.cc5d_gradb.Options( + filter_params=filter_params, + ep_scale=ep_scale, + b_tilde = model.em_fields.b_field,) +model.propagators.cc5d_curlb.options = model.propagators.cc5d_curlb.Options( + filter_params=filter_params, + ep_scale=ep_scale, + b_tilde = model.em_fields.b_field,) + +# background, perturbations and initial conditions +model.mhd.velocity.add_background(FieldsBackground( + type="LogicalConst", + values=(0., 0., 0.) +)) +model.mhd.velocity.add_perturbation( + perturbations.TorusModesSin( + given_in_basis="2", + amps=(0.001,0.001), + ms=(10,11), + ns=(-1,-1), + pfuns=('exp','exp'), + pfun_params=((0.44444444,0.1),(0.44444444,0.1)), + comp=0) +) +model.mhd.velocity.add_perturbation( + perturbations.TorusModesCos( + given_in_basis="2", + amps=(0.000015915494309189535,0.000014468631190172303), + ms=(10,11), + ns=(-1,-1), + pfuns=('d_exp','d_exp'), + pfun_params=((0.44444444,0.1),(0.44444444,0.1)), + comp=1), +) +n_profile = perturbations.ITPA_density( + n0=0.00720652, + c=(0.47023, 0.20323, 0.13273, 0.521298) +) +maxwellian_1 = maxwellians.CanonicalMaxwellian2D( + vth=(0.89740839, None), + equil=equil, + volume_form=True, + n=(0., n_profile), +) +background = maxwellian_1 +model.energetic_ions.var.add_background( + background=background, + n_as_volume_form=False) + +# optional: exclude variables from saving +# model.energetic_ions.var.save_data = False + +if __name__ == "__main__": + # start run + verbose = True + + main.run(model, + params_path=__file__, + env=env, + base_units=base_units, + time_opts=time_opts, + domain=domain, + equil=equil, + grid=grid, + derham_opts=derham_opts, + verbose=verbose, + ) diff --git a/ITPA_benchmark/tae_2n64_debug.py b/ITPA_benchmark/tae_2n64_debug.py new file mode 100644 index 0000000..ee683eb --- /dev/null +++ b/ITPA_benchmark/tae_2n64_debug.py @@ -0,0 +1,205 @@ +from struphy.io.options import EnvironmentOptions, BaseUnits, Time +from struphy.geometry import domains +from struphy.fields_background import equils +from struphy.topology import grids +from struphy.io.options import DerhamOptions +from struphy.io.options import FieldsBackground +from struphy.initial import perturbations +from struphy.kinetic_background import maxwellians +from struphy.pic.utilities import (LoadingParameters, + WeightsParameters, + BoundaryParameters, + BinningPlot, + KernelDensityPlot, + ) +from struphy import main + +# import model, set verbosity +from struphy.models.hybrid import LinearMHDDriftkineticCC + +# environment options +env = EnvironmentOptions( + out_folders='/tokp/work/bkna/tae', + sim_folder='tae_2n64_debug', + save_step=1, max_runtime=110, + restart=True) + +# units +base_units = BaseUnits(n=0.2) + +# time stepping +time_opts = Time(dt=0.25, Tend=2000., split_algo='Strang') + +# geometry +domain = domains.HollowTorus(a1=0.1, a2=1., R0=10., sfl=False, tor_period=6) + +# fluid equilibrium (can be used as part of initial conditions) +equil = equils.AdhocTorus( + a=1., + R0=10., + B0=3., + q_kind=0, + q0=1.71, + q1=1.87, + n1=0., + n2=0., + na=1., + p_kind=1., + p0=1., + p1=0.95, + p2=0.05, + beta=0.0018, +) + +# grid +grid = grids.TensorProductGrid( + Nel=(24,96,16), + mpi_dims_mask=(True, True, False) +) + +# derham options +derham_opts = DerhamOptions( + p=(3, 3, 3), + spl_kind=(False, True, True), + dirichlet_bc=((True, True), (False, False), (False, False)), + nq_pr=(6, 6, 6), + nquads=(6, 6, 6), + polar_ck=-1, +) + +# light-weight model instance +model = LinearMHDDriftkineticCC( + turn_off = (#'PushGuidingCenterBxEstar', + #'PushGuidingCenterParallel', + #'ShearAlfvenCurrentCoupling5D', + 'Magnetosonic', + #'CurrentCoupling5DDensity', + #'CurrentCoupling5DGradB', + #'CurrentCoupling5DCurlb', + ) +) + +# species parameters +model.mhd.set_phys_params() +model.energetic_ions.set_phys_params( + mass_number=2., +) + +loading_params = LoadingParameters( + ppc=100, + moments=(0., 0., 0.89740839, 0.89740839), + seed=1234, + ) +weights_params = WeightsParameters( + control_variate=True, +) +boundary_params = BoundaryParameters( + bc=("remove", "periodic", "periodic"), + bc_refill=("outer", "inner"), +) +model.energetic_ions.set_markers(loading_params=loading_params, + weights_params=weights_params, + boundary_params=boundary_params, + bufsize=1., + n_cols_diag=3, + n_cols_aux=5, + ) +model.energetic_ions.set_sorting_boxes(boxes_per_dim=None) +model.energetic_ions.set_save_data() + +# calculate constant +ep_scale = model.energetic_ions.var.species.mass_number / model.mhd.mass_number + +# params else +from struphy.pic.accumulation.filter import FilterParameters +from struphy.linear_algebra.solver import SolverParameters, DiscreteGradientSolverParameters +filter_params = FilterParameters(use_filter="hybrid",) +solver_params = SolverParameters() + +# propagator options +model.propagators.push_bxe.options = model.propagators.push_bxe.Options( + algo = 'explicit', + b_tilde = model.em_fields.b_field,) +model.propagators.push_parallel.options = model.propagators.push_parallel.Options( + algo = 'explicit', + b_tilde = model.em_fields.b_field,) +model.propagators.shearalfen_cc5d.options = model.propagators.shearalfen_cc5d.Options( + filter_params=filter_params, + ep_scale=ep_scale, + energetic_ions = model.energetic_ions.var, + solver_params=solver_params, + nonlinear=False,) +# model.propagators.magnetosonic.options = model.propagators.magnetosonic.Options( +# b_field=model.em_fields.b_field,) +model.propagators.cc5d_density.options = model.propagators.cc5d_density.Options( + filter_params=filter_params, + ep_scale=ep_scale, + energetic_ions = model.energetic_ions.var, + b_tilde = model.em_fields.b_field,) +model.propagators.cc5d_gradb.options = model.propagators.cc5d_gradb.Options( + filter_params=filter_params, + ep_scale=ep_scale, + b_tilde = model.em_fields.b_field,) +model.propagators.cc5d_curlb.options = model.propagators.cc5d_curlb.Options( + filter_params=filter_params, + ep_scale=ep_scale, + b_tilde = model.em_fields.b_field,) + +# background, perturbations and initial conditions +model.mhd.velocity.add_background(FieldsBackground( + type="LogicalConst", + values=(0., 0., 0.) +)) +model.mhd.velocity.add_perturbation( + perturbations.TorusModesSin( + given_in_basis="2", + amps=(0.001,0.001), + ms=(10,11), + ns=(-1,-1), + pfuns=('exp','exp'), + pfun_params=((0.44444444,0.1),(0.44444444,0.1)), + comp=0) +) +model.mhd.velocity.add_perturbation( + perturbations.TorusModesCos( + given_in_basis="2", + amps=(0.000015915494309189535,0.000014468631190172303), + ms=(10,11), + ns=(-1,-1), + pfuns=('d_exp','d_exp'), + pfun_params=((0.44444444,0.1),(0.44444444,0.1)), + comp=1), +) +n_profile = perturbations.ITPA_density( + n0=0.00720652, + c=(0.47023, 0.20323, 0.13273, 0.521298) +) +maxwellian_1 = maxwellians.CanonicalMaxwellian2D( + vth=(0.89740839, None), + equil=equil, + volume_form=True, + n=(0., n_profile), +) +background = maxwellian_1 +model.energetic_ions.var.add_background( + background=background, + n_as_volume_form=False) + +# optional: exclude variables from saving +# model.energetic_ions.var.save_data = False + +if __name__ == "__main__": + # start run + verbose = True + + main.run(model, + params_path=__file__, + env=env, + base_units=base_units, + time_opts=time_opts, + domain=domain, + equil=equil, + grid=grid, + derham_opts=derham_opts, + verbose=verbose, + ) diff --git a/ITPA_benchmark/tae_hybrid.py b/ITPA_benchmark/tae_hybrid.py new file mode 100644 index 0000000..7ffc57e --- /dev/null +++ b/ITPA_benchmark/tae_hybrid.py @@ -0,0 +1,207 @@ +from struphy.io.options import EnvironmentOptions, BaseUnits, Time +from struphy.geometry import domains +from struphy.fields_background import equils +from struphy.topology import grids +from struphy.io.options import DerhamOptions +from struphy.io.options import FieldsBackground +from struphy.initial import perturbations +from struphy.kinetic_background import maxwellians +from struphy.pic.utilities import (LoadingParameters, + WeightsParameters, + BoundaryParameters, + BinningPlot, + KernelDensityPlot, + ) +from struphy import main + +# import model, set verbosity +from struphy.models.hybrid import LinearMHDDriftkineticCC + +# environment options +# 32 +env = EnvironmentOptions( + out_folders='/tokp/work/bkna/tae', + sim_folder='tae_hybrid', + save_step=8, + max_runtime=2840, + restart=True) + +# units +base_units = BaseUnits(n=0.2) + +# time stepping +time_opts = Time(dt=0.25, Tend=2000., split_algo='Strang') + +# geometry +domain = domains.HollowTorus(a1=0.1, a2=1., R0=10., sfl=False, tor_period=6) + +# fluid equilibrium (can be used as part of initial conditions) +equil = equils.AdhocTorus( + a=1., + R0=10., + B0=3., + q_kind=0, + q0=1.71, + q1=1.87, + n1=0., + n2=0., + na=1., + p_kind=1., + p0=1., + p1=0.95, + p2=0.05, + beta=0.0018, +) + +# grid +grid = grids.TensorProductGrid( + Nel=(24,96,16), + mpi_dims_mask=(True, True, False) +) + +# derham options +derham_opts = DerhamOptions( + p=(3, 3, 3), + spl_kind=(False, True, True), + dirichlet_bc=((True, True), (False, False), (False, False)), + nq_pr=(6, 6, 6), + nquads=(6, 6, 6), + polar_ck=-1, +) + +# light-weight model instance +model = LinearMHDDriftkineticCC( + turn_off = (#'PushGuidingCenterBxEstar', + #'PushGuidingCenterParallel', + #'ShearAlfvenCurrentCoupling5D', + 'Magnetosonic', + #'CurrentCoupling5DDensity', + #'CurrentCoupling5DGradB', + #'CurrentCoupling5DCurlb', + ) +) + +# species parameters +model.mhd.set_phys_params() +model.energetic_ions.set_phys_params( + mass_number=2., +) + +loading_params = LoadingParameters( + ppc=100, + moments=(0., 0., 0.89740839, 0.89740839), + seed=1234, + ) +weights_params = WeightsParameters( + control_variate=True, +) +boundary_params = BoundaryParameters( + bc=("remove", "periodic", "periodic"), + bc_refill=("outer", "inner"), +) +model.energetic_ions.set_markers(loading_params=loading_params, + weights_params=weights_params, + boundary_params=boundary_params, + bufsize=1., + n_cols_diag=3, + n_cols_aux=5, + ) +model.energetic_ions.set_sorting_boxes(boxes_per_dim=None) +model.energetic_ions.set_save_data() + +# calculate constant +ep_scale = model.energetic_ions.var.species.mass_number / model.mhd.mass_number + +# params else +from struphy.pic.accumulation.filter import FilterParameters +from struphy.linear_algebra.solver import SolverParameters, DiscreteGradientSolverParameters +filter_params = FilterParameters(use_filter="hybrid",) +solver_params = SolverParameters() + +# propagator options +model.propagators.push_bxe.options = model.propagators.push_bxe.Options( + algo = 'explicit', + b_tilde = model.em_fields.b_field,) +model.propagators.push_parallel.options = model.propagators.push_parallel.Options( + algo = 'explicit', + b_tilde = model.em_fields.b_field,) +model.propagators.shearalfen_cc5d.options = model.propagators.shearalfen_cc5d.Options( + filter_params=filter_params, + ep_scale=ep_scale, + energetic_ions = model.energetic_ions.var, + solver_params=solver_params, + nonlinear=False,) +# model.propagators.magnetosonic.options = model.propagators.magnetosonic.Options( +# b_field=model.em_fields.b_field,) +model.propagators.cc5d_density.options = model.propagators.cc5d_density.Options( + filter_params=filter_params, + ep_scale=ep_scale, + energetic_ions = model.energetic_ions.var, + b_tilde = model.em_fields.b_field,) +model.propagators.cc5d_gradb.options = model.propagators.cc5d_gradb.Options( + filter_params=filter_params, + ep_scale=ep_scale, + b_tilde = model.em_fields.b_field,) +model.propagators.cc5d_curlb.options = model.propagators.cc5d_curlb.Options( + filter_params=filter_params, + ep_scale=ep_scale, + b_tilde = model.em_fields.b_field,) + +# background, perturbations and initial conditions +model.mhd.velocity.add_background(FieldsBackground( + type="LogicalConst", + values=(0., 0., 0.) +)) +model.mhd.velocity.add_perturbation( + perturbations.TorusModesSin( + given_in_basis="2", + amps=(0.001,0.001), + ms=(10,11), + ns=(-1,-1), + pfuns=('exp','exp'), + pfun_params=((0.44444444,0.1),(0.44444444,0.1)), + comp=0) +) +model.mhd.velocity.add_perturbation( + perturbations.TorusModesCos( + given_in_basis="2", + amps=(0.000015915494309189535,0.000014468631190172303), + ms=(10,11), + ns=(-1,-1), + pfuns=('d_exp','d_exp'), + pfun_params=((0.44444444,0.1),(0.44444444,0.1)), + comp=1), +) +n_profile = perturbations.ITPA_density( + n0=0.00720652, + c=(0.47023, 0.20323, 0.13273, 0.521298) +) +maxwellian_1 = maxwellians.CanonicalMaxwellian2D( + vth=(0.89740839, None), + equil=equil, + volume_form=True, + n=(0., n_profile), +) +background = maxwellian_1 +model.energetic_ions.var.add_background( + background=background, + n_as_volume_form=False) + +# optional: exclude variables from saving +# model.energetic_ions.var.save_data = False + +if __name__ == "__main__": + # start run + verbose = True + + main.run(model, + params_path=__file__, + env=env, + base_units=base_units, + time_opts=time_opts, + domain=domain, + equil=equil, + grid=grid, + derham_opts=derham_opts, + verbose=verbose, + ) diff --git a/ITPA_benchmark/tae_hybrid_nonlinear.py b/ITPA_benchmark/tae_hybrid_nonlinear.py new file mode 100644 index 0000000..fd74bd7 --- /dev/null +++ b/ITPA_benchmark/tae_hybrid_nonlinear.py @@ -0,0 +1,207 @@ +from struphy.io.options import EnvironmentOptions, BaseUnits, Time +from struphy.geometry import domains +from struphy.fields_background import equils +from struphy.topology import grids +from struphy.io.options import DerhamOptions +from struphy.io.options import FieldsBackground +from struphy.initial import perturbations +from struphy.kinetic_background import maxwellians +from struphy.pic.utilities import (LoadingParameters, + WeightsParameters, + BoundaryParameters, + BinningPlot, + KernelDensityPlot, + ) +from struphy import main + +# import model, set verbosity +from struphy.models.hybrid import LinearMHDDriftkineticCC + +# environment options +# 32 +env = EnvironmentOptions( + out_folders='/tokp/work/bkna/tae', + sim_folder='tae_hybrid_nonlinear', + save_step=8, + max_runtime=2840, + restart=True) + +# units +base_units = BaseUnits(n=0.2) + +# time stepping +time_opts = Time(dt=0.25, Tend=2000., split_algo='Strang') + +# geometry +domain = domains.HollowTorus(a1=0.1, a2=1., R0=10., sfl=False, tor_period=6) + +# fluid equilibrium (can be used as part of initial conditions) +equil = equils.AdhocTorus( + a=1., + R0=10., + B0=3., + q_kind=0, + q0=1.71, + q1=1.87, + n1=0., + n2=0., + na=1., + p_kind=1., + p0=1., + p1=0.95, + p2=0.05, + beta=0.0018, +) + +# grid +grid = grids.TensorProductGrid( + Nel=(24,96,16), + mpi_dims_mask=(True, True, False) +) + +# derham options +derham_opts = DerhamOptions( + p=(3, 3, 3), + spl_kind=(False, True, True), + dirichlet_bc=((True, True), (False, False), (False, False)), + nq_pr=(6, 6, 6), + nquads=(6, 6, 6), + polar_ck=-1, +) + +# light-weight model instance +model = LinearMHDDriftkineticCC( + turn_off = (#'PushGuidingCenterBxEstar', + #'PushGuidingCenterParallel', + #'ShearAlfvenCurrentCoupling5D', + 'Magnetosonic', + #'CurrentCoupling5DDensity', + #'CurrentCoupling5DGradB', + #'CurrentCoupling5DCurlb', + ) +) + +# species parameters +model.mhd.set_phys_params() +model.energetic_ions.set_phys_params( + mass_number=2., +) + +loading_params = LoadingParameters( + ppc=100, + moments=(0., 0., 0.89740839, 0.89740839), + seed=1234, + ) +weights_params = WeightsParameters( + control_variate=True, +) +boundary_params = BoundaryParameters( + bc=("remove", "periodic", "periodic"), + bc_refill=("outer", "inner"), +) +model.energetic_ions.set_markers(loading_params=loading_params, + weights_params=weights_params, + boundary_params=boundary_params, + bufsize=1., + n_cols_diag=3, + n_cols_aux=5, + ) +model.energetic_ions.set_sorting_boxes(boxes_per_dim=None) +model.energetic_ions.set_save_data() + +# calculate constant +ep_scale = model.energetic_ions.var.species.mass_number / model.mhd.mass_number + +# params else +from struphy.pic.accumulation.filter import FilterParameters +from struphy.linear_algebra.solver import SolverParameters, DiscreteGradientSolverParameters +filter_params = FilterParameters(use_filter="hybrid",) +solver_params = SolverParameters() + +# propagator options +model.propagators.push_bxe.options = model.propagators.push_bxe.Options( + algo = 'explicit', + b_tilde = model.em_fields.b_field,) +model.propagators.push_parallel.options = model.propagators.push_parallel.Options( + algo = 'explicit', + b_tilde = model.em_fields.b_field,) +model.propagators.shearalfen_cc5d.options = model.propagators.shearalfen_cc5d.Options( + filter_params=filter_params, + ep_scale=ep_scale, + energetic_ions = model.energetic_ions.var, + solver_params=solver_params, + nonlinear=True,) +# model.propagators.magnetosonic.options = model.propagators.magnetosonic.Options( +# b_field=model.em_fields.b_field,) +model.propagators.cc5d_density.options = model.propagators.cc5d_density.Options( + filter_params=filter_params, + ep_scale=ep_scale, + energetic_ions = model.energetic_ions.var, + b_tilde = model.em_fields.b_field,) +model.propagators.cc5d_gradb.options = model.propagators.cc5d_gradb.Options( + filter_params=filter_params, + ep_scale=ep_scale, + b_tilde = model.em_fields.b_field,) +model.propagators.cc5d_curlb.options = model.propagators.cc5d_curlb.Options( + filter_params=filter_params, + ep_scale=ep_scale, + b_tilde = model.em_fields.b_field,) + +# background, perturbations and initial conditions +model.mhd.velocity.add_background(FieldsBackground( + type="LogicalConst", + values=(0., 0., 0.) +)) +model.mhd.velocity.add_perturbation( + perturbations.TorusModesSin( + given_in_basis="2", + amps=(0.001,0.001), + ms=(10,11), + ns=(-1,-1), + pfuns=('exp','exp'), + pfun_params=((0.44444444,0.1),(0.44444444,0.1)), + comp=0) +) +model.mhd.velocity.add_perturbation( + perturbations.TorusModesCos( + given_in_basis="2", + amps=(0.000015915494309189535,0.000014468631190172303), + ms=(10,11), + ns=(-1,-1), + pfuns=('d_exp','d_exp'), + pfun_params=((0.44444444,0.1),(0.44444444,0.1)), + comp=1), +) +n_profile = perturbations.ITPA_density( + n0=0.00720652, + c=(0.47023, 0.20323, 0.13273, 0.521298) +) +maxwellian_1 = maxwellians.CanonicalMaxwellian2D( + vth=(0.89740839, None), + equil=equil, + volume_form=True, + n=(0., n_profile), +) +background = maxwellian_1 +model.energetic_ions.var.add_background( + background=background, + n_as_volume_form=False) + +# optional: exclude variables from saving +# model.energetic_ions.var.save_data = False + +if __name__ == "__main__": + # start run + verbose = True + + main.run(model, + params_path=__file__, + env=env, + base_units=base_units, + time_opts=time_opts, + domain=domain, + equil=equil, + grid=grid, + derham_opts=derham_opts, + verbose=verbose, + ) diff --git a/ITPA_benchmark/tae_hybrid_nonlinear_sonic_dt0125.py b/ITPA_benchmark/tae_hybrid_nonlinear_sonic_dt0125.py new file mode 100644 index 0000000..ba8c02a --- /dev/null +++ b/ITPA_benchmark/tae_hybrid_nonlinear_sonic_dt0125.py @@ -0,0 +1,206 @@ +from struphy.io.options import EnvironmentOptions, BaseUnits, Time +from struphy.geometry import domains +from struphy.fields_background import equils +from struphy.topology import grids +from struphy.io.options import DerhamOptions +from struphy.io.options import FieldsBackground +from struphy.initial import perturbations +from struphy.kinetic_background import maxwellians +from struphy.pic.utilities import (LoadingParameters, + WeightsParameters, + BoundaryParameters, + BinningPlot, + KernelDensityPlot, + ) +from struphy import main + +# import model, set verbosity +from struphy.models.hybrid import LinearMHDDriftkineticCC + +# environment options +env = EnvironmentOptions( + out_folders='/tokp/work/bkna/tae', + sim_folder='tae_hybrid_nonlinear_sonic_dt0125', + save_step=8, + max_runtime=2840, + restart=False) + +# units +base_units = BaseUnits(n=0.2) + +# time stepping +time_opts = Time(dt=0.125, Tend=2000., split_algo='Strang') + +# geometry +domain = domains.HollowTorus(a1=0.1, a2=1., R0=10., sfl=False, tor_period=6) + +# fluid equilibrium (can be used as part of initial conditions) +equil = equils.AdhocTorus( + a=1., + R0=10., + B0=3., + q_kind=0, + q0=1.71, + q1=1.87, + n1=0., + n2=0., + na=1., + p_kind=1., + p0=1., + p1=0.95, + p2=0.05, + beta=0.0018, +) + +# grid +grid = grids.TensorProductGrid( + Nel=(24,96,16), + mpi_dims_mask=(True, True, False) +) + +# derham options +derham_opts = DerhamOptions( + p=(3, 3, 3), + spl_kind=(False, True, True), + dirichlet_bc=((True, True), (False, False), (False, False)), + nq_pr=(6, 6, 6), + nquads=(6, 6, 6), + polar_ck=-1, +) + +# light-weight model instance +model = LinearMHDDriftkineticCC( + turn_off = (#'PushGuidingCenterBxEstar', + #'PushGuidingCenterParallel', + #'ShearAlfvenCurrentCoupling5D', + #'Magnetosonic', + #'CurrentCoupling5DDensity', + #'CurrentCoupling5DGradB', + #'CurrentCoupling5DCurlb', + ) +) + +# species parameters +model.mhd.set_phys_params() +model.energetic_ions.set_phys_params( + mass_number=2., +) + +loading_params = LoadingParameters( + ppc=100, + moments=(0., 0., 0.89740839, 0.89740839), + seed=1234, + ) +weights_params = WeightsParameters( + control_variate=True, +) +boundary_params = BoundaryParameters( + bc=("remove", "periodic", "periodic"), + bc_refill=("outer", "inner"), +) +model.energetic_ions.set_markers(loading_params=loading_params, + weights_params=weights_params, + boundary_params=boundary_params, + bufsize=1., + n_cols_diag=3, + n_cols_aux=5, + ) +model.energetic_ions.set_sorting_boxes(boxes_per_dim=None) +model.energetic_ions.set_save_data() + +# calculate constant +ep_scale = model.energetic_ions.var.species.mass_number / model.mhd.mass_number + +# params else +from struphy.pic.accumulation.filter import FilterParameters +from struphy.linear_algebra.solver import SolverParameters, DiscreteGradientSolverParameters +filter_params = FilterParameters(use_filter="hybrid",) +solver_params = SolverParameters() + +# propagator options +model.propagators.push_bxe.options = model.propagators.push_bxe.Options( + algo = 'explicit', + b_tilde = model.em_fields.b_field,) +model.propagators.push_parallel.options = model.propagators.push_parallel.Options( + algo = 'explicit', + b_tilde = model.em_fields.b_field,) +model.propagators.shearalfen_cc5d.options = model.propagators.shearalfen_cc5d.Options( + filter_params=filter_params, + ep_scale=ep_scale, + energetic_ions = model.energetic_ions.var, + solver_params=solver_params, + nonlinear=True,) +model.propagators.magnetosonic.options = model.propagators.magnetosonic.Options( + b_field=model.em_fields.b_field,) +model.propagators.cc5d_density.options = model.propagators.cc5d_density.Options( + filter_params=filter_params, + ep_scale=ep_scale, + energetic_ions = model.energetic_ions.var, + b_tilde = model.em_fields.b_field,) +model.propagators.cc5d_gradb.options = model.propagators.cc5d_gradb.Options( + filter_params=filter_params, + ep_scale=ep_scale, + b_tilde = model.em_fields.b_field,) +model.propagators.cc5d_curlb.options = model.propagators.cc5d_curlb.Options( + filter_params=filter_params, + ep_scale=ep_scale, + b_tilde = model.em_fields.b_field,) + +# background, perturbations and initial conditions +model.mhd.velocity.add_background(FieldsBackground( + type="LogicalConst", + values=(0., 0., 0.) +)) +model.mhd.velocity.add_perturbation( + perturbations.TorusModesSin( + given_in_basis="2", + amps=(0.001,0.001), + ms=(10,11), + ns=(-1,-1), + pfuns=('exp','exp'), + pfun_params=((0.44444444,0.1),(0.44444444,0.1)), + comp=0) +) +model.mhd.velocity.add_perturbation( + perturbations.TorusModesCos( + given_in_basis="2", + amps=(0.000015915494309189535,0.000014468631190172303), + ms=(10,11), + ns=(-1,-1), + pfuns=('d_exp','d_exp'), + pfun_params=((0.44444444,0.1),(0.44444444,0.1)), + comp=1), +) +n_profile = perturbations.ITPA_density( + n0=0.00720652, + c=(0.47023, 0.20323, 0.13273, 0.521298) +) +maxwellian_1 = maxwellians.CanonicalMaxwellian2D( + vth=(0.89740839, None), + equil=equil, + volume_form=True, + n=(0., n_profile), +) +background = maxwellian_1 +model.energetic_ions.var.add_background( + background=background, + n_as_volume_form=False) + +# optional: exclude variables from saving +# model.energetic_ions.var.save_data = False + +if __name__ == "__main__": + # start run + verbose = True + + main.run(model, + params_path=__file__, + env=env, + base_units=base_units, + time_opts=time_opts, + domain=domain, + equil=equil, + grid=grid, + derham_opts=derham_opts, + verbose=verbose, + ) diff --git a/ITPA_benchmark/tae_hybrid_nonlinear_sonic_gradp0_dt0125.py b/ITPA_benchmark/tae_hybrid_nonlinear_sonic_gradp0_dt0125.py new file mode 100644 index 0000000..a9a5e24 --- /dev/null +++ b/ITPA_benchmark/tae_hybrid_nonlinear_sonic_gradp0_dt0125.py @@ -0,0 +1,206 @@ +from struphy.io.options import EnvironmentOptions, BaseUnits, Time +from struphy.geometry import domains +from struphy.fields_background import equils +from struphy.topology import grids +from struphy.io.options import DerhamOptions +from struphy.io.options import FieldsBackground +from struphy.initial import perturbations +from struphy.kinetic_background import maxwellians +from struphy.pic.utilities import (LoadingParameters, + WeightsParameters, + BoundaryParameters, + BinningPlot, + KernelDensityPlot, + ) +from struphy import main + +# import model, set verbosity +from struphy.models.hybrid import LinearMHDDriftkineticCC + +# environment options +env = EnvironmentOptions( + out_folders='/tokp/work/bkna/tae', + sim_folder='tae_hybrid_nonlinear_sonic_gradp0_dt0125', + save_step=8, + max_runtime=2840, + restart=False) + +# units +base_units = BaseUnits(n=0.2) + +# time stepping +time_opts = Time(dt=0.125, Tend=2000., split_algo='Strang') + +# geometry +domain = domains.HollowTorus(a1=0.1, a2=1., R0=10., sfl=False, tor_period=6) + +# fluid equilibrium (can be used as part of initial conditions) +equil = equils.AdhocTorus( + a=1., + R0=10., + B0=3., + q_kind=0, + q0=1.71, + q1=1.87, + n1=0., + n2=0., + na=1., + p_kind=1., + p0=1., + p1=0., + p2=0., + beta=0.0018, +) + +# grid +grid = grids.TensorProductGrid( + Nel=(24,96,16), + mpi_dims_mask=(True, True, False) +) + +# derham options +derham_opts = DerhamOptions( + p=(3, 3, 3), + spl_kind=(False, True, True), + dirichlet_bc=((True, True), (False, False), (False, False)), + nq_pr=(6, 6, 6), + nquads=(6, 6, 6), + polar_ck=-1, +) + +# light-weight model instance +model = LinearMHDDriftkineticCC( + turn_off = (#'PushGuidingCenterBxEstar', + #'PushGuidingCenterParallel', + #'ShearAlfvenCurrentCoupling5D', + #'Magnetosonic', + #'CurrentCoupling5DDensity', + #'CurrentCoupling5DGradB', + #'CurrentCoupling5DCurlb', + ) +) + +# species parameters +model.mhd.set_phys_params() +model.energetic_ions.set_phys_params( + mass_number=2., +) + +loading_params = LoadingParameters( + ppc=100, + moments=(0., 0., 0.89740839, 0.89740839), + seed=1234, + ) +weights_params = WeightsParameters( + control_variate=True, +) +boundary_params = BoundaryParameters( + bc=("remove", "periodic", "periodic"), + bc_refill=("outer", "inner"), +) +model.energetic_ions.set_markers(loading_params=loading_params, + weights_params=weights_params, + boundary_params=boundary_params, + bufsize=1., + n_cols_diag=3, + n_cols_aux=5, + ) +model.energetic_ions.set_sorting_boxes(boxes_per_dim=None) +model.energetic_ions.set_save_data() + +# calculate constant +ep_scale = model.energetic_ions.var.species.mass_number / model.mhd.mass_number + +# params else +from struphy.pic.accumulation.filter import FilterParameters +from struphy.linear_algebra.solver import SolverParameters, DiscreteGradientSolverParameters +filter_params = FilterParameters(use_filter="hybrid",) +solver_params = SolverParameters() + +# propagator options +model.propagators.push_bxe.options = model.propagators.push_bxe.Options( + algo = 'explicit', + b_tilde = model.em_fields.b_field,) +model.propagators.push_parallel.options = model.propagators.push_parallel.Options( + algo = 'explicit', + b_tilde = model.em_fields.b_field,) +model.propagators.shearalfen_cc5d.options = model.propagators.shearalfen_cc5d.Options( + filter_params=filter_params, + ep_scale=ep_scale, + energetic_ions = model.energetic_ions.var, + solver_params=solver_params, + nonlinear=True,) +model.propagators.magnetosonic.options = model.propagators.magnetosonic.Options( + b_field=model.em_fields.b_field,) +model.propagators.cc5d_density.options = model.propagators.cc5d_density.Options( + filter_params=filter_params, + ep_scale=ep_scale, + energetic_ions = model.energetic_ions.var, + b_tilde = model.em_fields.b_field,) +model.propagators.cc5d_gradb.options = model.propagators.cc5d_gradb.Options( + filter_params=filter_params, + ep_scale=ep_scale, + b_tilde = model.em_fields.b_field,) +model.propagators.cc5d_curlb.options = model.propagators.cc5d_curlb.Options( + filter_params=filter_params, + ep_scale=ep_scale, + b_tilde = model.em_fields.b_field,) + +# background, perturbations and initial conditions +model.mhd.velocity.add_background(FieldsBackground( + type="LogicalConst", + values=(0., 0., 0.) +)) +model.mhd.velocity.add_perturbation( + perturbations.TorusModesSin( + given_in_basis="2", + amps=(0.001,0.001), + ms=(10,11), + ns=(-1,-1), + pfuns=('exp','exp'), + pfun_params=((0.44444444,0.1),(0.44444444,0.1)), + comp=0) +) +model.mhd.velocity.add_perturbation( + perturbations.TorusModesCos( + given_in_basis="2", + amps=(0.000015915494309189535,0.000014468631190172303), + ms=(10,11), + ns=(-1,-1), + pfuns=('d_exp','d_exp'), + pfun_params=((0.44444444,0.1),(0.44444444,0.1)), + comp=1), +) +n_profile = perturbations.ITPA_density( + n0=0.00720652, + c=(0.47023, 0.20323, 0.13273, 0.521298) +) +maxwellian_1 = maxwellians.CanonicalMaxwellian2D( + vth=(0.89740839, None), + equil=equil, + volume_form=True, + n=(0., n_profile), +) +background = maxwellian_1 +model.energetic_ions.var.add_background( + background=background, + n_as_volume_form=False) + +# optional: exclude variables from saving +# model.energetic_ions.var.save_data = False + +if __name__ == "__main__": + # start run + verbose = True + + main.run(model, + params_path=__file__, + env=env, + base_units=base_units, + time_opts=time_opts, + domain=domain, + equil=equil, + grid=grid, + derham_opts=derham_opts, + verbose=verbose, + ) diff --git a/ITPA_benchmark/tae_hybrid_nonlinear_sonic_m2jboff.py b/ITPA_benchmark/tae_hybrid_nonlinear_sonic_m2jboff.py new file mode 100644 index 0000000..6624745 --- /dev/null +++ b/ITPA_benchmark/tae_hybrid_nonlinear_sonic_m2jboff.py @@ -0,0 +1,206 @@ +from struphy.io.options import EnvironmentOptions, BaseUnits, Time +from struphy.geometry import domains +from struphy.fields_background import equils +from struphy.topology import grids +from struphy.io.options import DerhamOptions +from struphy.io.options import FieldsBackground +from struphy.initial import perturbations +from struphy.kinetic_background import maxwellians +from struphy.pic.utilities import (LoadingParameters, + WeightsParameters, + BoundaryParameters, + BinningPlot, + KernelDensityPlot, + ) +from struphy import main + +# import model, set verbosity +from struphy.models.hybrid import LinearMHDDriftkineticCC + +# environment options +# 64 +env = EnvironmentOptions( + out_folders='/tokp/work/bkna/tae', + sim_folder='tae_hybrid_nonlinear_sonic_m2jboff', + save_step=8, + max_runtime=2840, + restart=True) + +# units +base_units = BaseUnits(n=0.2) + +# time stepping +time_opts = Time(dt=0.25, Tend=2000., split_algo='Strang') + +# geometry +domain = domains.HollowTorus(a1=0.1, a2=1., R0=10., sfl=False, tor_period=6) + +# fluid equilibrium (can be used as part of initial conditions) +equil = equils.AdhocTorus( + a=1., + R0=10., + B0=3., + q_kind=0, + q0=1.71, + q1=1.87, + n1=0., + n2=0., + na=1., + p_kind=1., + p0=1., + p1=0.95, + p2=0.05, + beta=0.0018, +) + +# grid +grid = grids.TensorProductGrid( + Nel=(24,96,16), + mpi_dims_mask=(True, True, False) +) + +# derham options +derham_opts = DerhamOptions( + p=(3, 3, 3), + spl_kind=(False, True, True), + dirichlet_bc=((True, True), (False, False), (False, False)), + nq_pr=(6, 6, 6), + nquads=(6, 6, 6), + polar_ck=-1, +) + +# light-weight model instance +model = LinearMHDDriftkineticCC( + turn_off = (#'PushGuidingCenterBxEstar', + #'PushGuidingCenterParallel', + #'ShearAlfvenCurrentCoupling5D', + #'Magnetosonic', + #'CurrentCoupling5DDensity', + #'CurrentCoupling5DGradB', + #'CurrentCoupling5DCurlb', + ) +) + +# species parameters +model.mhd.set_phys_params() +model.energetic_ions.set_phys_params( + mass_number=2., +) + +loading_params = LoadingParameters( + ppc=100, + moments=(0., 0., 0.89740839, 0.89740839), + seed=1234, + ) +weights_params = WeightsParameters( + control_variate=True, +) +boundary_params = BoundaryParameters( + bc=("remove", "periodic", "periodic"), + bc_refill=("outer", "inner"), +) +model.energetic_ions.set_markers(loading_params=loading_params, + weights_params=weights_params, + boundary_params=boundary_params, + bufsize=1., + n_cols_diag=3, + n_cols_aux=5, + ) +model.energetic_ions.set_sorting_boxes(boxes_per_dim=None) +model.energetic_ions.set_save_data() + +# calculate constant +ep_scale = model.energetic_ions.var.species.mass_number / model.mhd.mass_number + +# params else +from struphy.pic.accumulation.filter import FilterParameters +from struphy.linear_algebra.solver import SolverParameters, DiscreteGradientSolverParameters +filter_params = FilterParameters(use_filter="hybrid",) +solver_params = SolverParameters() + +# propagator options +model.propagators.push_bxe.options = model.propagators.push_bxe.Options( + algo = 'explicit', + b_tilde = model.em_fields.b_field,) +model.propagators.push_parallel.options = model.propagators.push_parallel.Options( + algo = 'explicit', + b_tilde = model.em_fields.b_field,) +model.propagators.shearalfen_cc5d.options = model.propagators.shearalfen_cc5d.Options( + filter_params=filter_params, + ep_scale=ep_scale, + energetic_ions = model.energetic_ions.var, + solver_params=solver_params, + nonlinear=True,) +model.propagators.magnetosonic.options = model.propagators.magnetosonic.Options() +model.propagators.cc5d_density.options = model.propagators.cc5d_density.Options( + filter_params=filter_params, + ep_scale=ep_scale, + energetic_ions = model.energetic_ions.var, + b_tilde = model.em_fields.b_field,) +model.propagators.cc5d_gradb.options = model.propagators.cc5d_gradb.Options( + filter_params=filter_params, + ep_scale=ep_scale, + b_tilde = model.em_fields.b_field,) +model.propagators.cc5d_curlb.options = model.propagators.cc5d_curlb.Options( + filter_params=filter_params, + ep_scale=ep_scale, + b_tilde = model.em_fields.b_field,) + +# background, perturbations and initial conditions +model.mhd.velocity.add_background(FieldsBackground( + type="LogicalConst", + values=(0., 0., 0.) +)) +model.mhd.velocity.add_perturbation( + perturbations.TorusModesSin( + given_in_basis="2", + amps=(0.001,0.001), + ms=(10,11), + ns=(-1,-1), + pfuns=('exp','exp'), + pfun_params=((0.44444444,0.1),(0.44444444,0.1)), + comp=0) +) +model.mhd.velocity.add_perturbation( + perturbations.TorusModesCos( + given_in_basis="2", + amps=(0.000015915494309189535,0.000014468631190172303), + ms=(10,11), + ns=(-1,-1), + pfuns=('d_exp','d_exp'), + pfun_params=((0.44444444,0.1),(0.44444444,0.1)), + comp=1), +) +n_profile = perturbations.ITPA_density( + n0=0.00720652, + c=(0.47023, 0.20323, 0.13273, 0.521298) +) +maxwellian_1 = maxwellians.CanonicalMaxwellian2D( + vth=(0.89740839, None), + equil=equil, + volume_form=True, + n=(0., n_profile), +) +background = maxwellian_1 +model.energetic_ions.var.add_background( + background=background, + n_as_volume_form=False) + +# optional: exclude variables from saving +# model.energetic_ions.var.save_data = False + +if __name__ == "__main__": + # start run + verbose = True + + main.run(model, + params_path=__file__, + env=env, + base_units=base_units, + time_opts=time_opts, + domain=domain, + equil=equil, + grid=grid, + derham_opts=derham_opts, + verbose=verbose, + ) diff --git a/ITPA_benchmark/tae_hybrid_nonlinear_sonic_m2jboff_dt0125.py b/ITPA_benchmark/tae_hybrid_nonlinear_sonic_m2jboff_dt0125.py new file mode 100644 index 0000000..163d505 --- /dev/null +++ b/ITPA_benchmark/tae_hybrid_nonlinear_sonic_m2jboff_dt0125.py @@ -0,0 +1,206 @@ +from struphy.io.options import EnvironmentOptions, BaseUnits, Time +from struphy.geometry import domains +from struphy.fields_background import equils +from struphy.topology import grids +from struphy.io.options import DerhamOptions +from struphy.io.options import FieldsBackground +from struphy.initial import perturbations +from struphy.kinetic_background import maxwellians +from struphy.pic.utilities import (LoadingParameters, + WeightsParameters, + BoundaryParameters, + BinningPlot, + KernelDensityPlot, + ) +from struphy import main + +# import model, set verbosity +from struphy.models.hybrid import LinearMHDDriftkineticCC + +# environment options +# 64 +env = EnvironmentOptions( + out_folders='/tokp/work/bkna/tae', + sim_folder='tae_hybrid_nonlinear_sonic_m2jboff_dt0125', + save_step=8, + max_runtime=2840, + restart=False) + +# units +base_units = BaseUnits(n=0.2) + +# time stepping +time_opts = Time(dt=0.125, Tend=2000., split_algo='Strang') + +# geometry +domain = domains.HollowTorus(a1=0.1, a2=1., R0=10., sfl=False, tor_period=6) + +# fluid equilibrium (can be used as part of initial conditions) +equil = equils.AdhocTorus( + a=1., + R0=10., + B0=3., + q_kind=0, + q0=1.71, + q1=1.87, + n1=0., + n2=0., + na=1., + p_kind=1., + p0=1., + p1=0.95, + p2=0.05, + beta=0.0018, +) + +# grid +grid = grids.TensorProductGrid( + Nel=(24,96,16), + mpi_dims_mask=(True, True, False) +) + +# derham options +derham_opts = DerhamOptions( + p=(3, 3, 3), + spl_kind=(False, True, True), + dirichlet_bc=((True, True), (False, False), (False, False)), + nq_pr=(6, 6, 6), + nquads=(6, 6, 6), + polar_ck=-1, +) + +# light-weight model instance +model = LinearMHDDriftkineticCC( + turn_off = (#'PushGuidingCenterBxEstar', + #'PushGuidingCenterParallel', + #'ShearAlfvenCurrentCoupling5D', + #'Magnetosonic', + #'CurrentCoupling5DDensity', + #'CurrentCoupling5DGradB', + #'CurrentCoupling5DCurlb', + ) +) + +# species parameters +model.mhd.set_phys_params() +model.energetic_ions.set_phys_params( + mass_number=2., +) + +loading_params = LoadingParameters( + ppc=100, + moments=(0., 0., 0.89740839, 0.89740839), + seed=1234, + ) +weights_params = WeightsParameters( + control_variate=True, +) +boundary_params = BoundaryParameters( + bc=("remove", "periodic", "periodic"), + bc_refill=("outer", "inner"), +) +model.energetic_ions.set_markers(loading_params=loading_params, + weights_params=weights_params, + boundary_params=boundary_params, + bufsize=1., + n_cols_diag=3, + n_cols_aux=5, + ) +model.energetic_ions.set_sorting_boxes(boxes_per_dim=None) +model.energetic_ions.set_save_data() + +# calculate constant +ep_scale = model.energetic_ions.var.species.mass_number / model.mhd.mass_number + +# params else +from struphy.pic.accumulation.filter import FilterParameters +from struphy.linear_algebra.solver import SolverParameters, DiscreteGradientSolverParameters +filter_params = FilterParameters(use_filter="hybrid",) +solver_params = SolverParameters() + +# propagator options +model.propagators.push_bxe.options = model.propagators.push_bxe.Options( + algo = 'explicit', + b_tilde = model.em_fields.b_field,) +model.propagators.push_parallel.options = model.propagators.push_parallel.Options( + algo = 'explicit', + b_tilde = model.em_fields.b_field,) +model.propagators.shearalfen_cc5d.options = model.propagators.shearalfen_cc5d.Options( + filter_params=filter_params, + ep_scale=ep_scale, + energetic_ions = model.energetic_ions.var, + solver_params=solver_params, + nonlinear=True,) +model.propagators.magnetosonic.options = model.propagators.magnetosonic.Options() +model.propagators.cc5d_density.options = model.propagators.cc5d_density.Options( + filter_params=filter_params, + ep_scale=ep_scale, + energetic_ions = model.energetic_ions.var, + b_tilde = model.em_fields.b_field,) +model.propagators.cc5d_gradb.options = model.propagators.cc5d_gradb.Options( + filter_params=filter_params, + ep_scale=ep_scale, + b_tilde = model.em_fields.b_field,) +model.propagators.cc5d_curlb.options = model.propagators.cc5d_curlb.Options( + filter_params=filter_params, + ep_scale=ep_scale, + b_tilde = model.em_fields.b_field,) + +# background, perturbations and initial conditions +model.mhd.velocity.add_background(FieldsBackground( + type="LogicalConst", + values=(0., 0., 0.) +)) +model.mhd.velocity.add_perturbation( + perturbations.TorusModesSin( + given_in_basis="2", + amps=(0.001,0.001), + ms=(10,11), + ns=(-1,-1), + pfuns=('exp','exp'), + pfun_params=((0.44444444,0.1),(0.44444444,0.1)), + comp=0) +) +model.mhd.velocity.add_perturbation( + perturbations.TorusModesCos( + given_in_basis="2", + amps=(0.000015915494309189535,0.000014468631190172303), + ms=(10,11), + ns=(-1,-1), + pfuns=('d_exp','d_exp'), + pfun_params=((0.44444444,0.1),(0.44444444,0.1)), + comp=1), +) +n_profile = perturbations.ITPA_density( + n0=0.00720652, + c=(0.47023, 0.20323, 0.13273, 0.521298) +) +maxwellian_1 = maxwellians.CanonicalMaxwellian2D( + vth=(0.89740839, None), + equil=equil, + volume_form=True, + n=(0., n_profile), +) +background = maxwellian_1 +model.energetic_ions.var.add_background( + background=background, + n_as_volume_form=False) + +# optional: exclude variables from saving +# model.energetic_ions.var.save_data = False + +if __name__ == "__main__": + # start run + verbose = True + + main.run(model, + params_path=__file__, + env=env, + base_units=base_units, + time_opts=time_opts, + domain=domain, + equil=equil, + grid=grid, + derham_opts=derham_opts, + verbose=verbose, + ) diff --git a/ITPA_benchmark/tae_hybrid_sonic.py b/ITPA_benchmark/tae_hybrid_sonic.py new file mode 100644 index 0000000..c66ba7a --- /dev/null +++ b/ITPA_benchmark/tae_hybrid_sonic.py @@ -0,0 +1,204 @@ +from struphy.io.options import EnvironmentOptions, BaseUnits, Time +from struphy.geometry import domains +from struphy.fields_background import equils +from struphy.topology import grids +from struphy.io.options import DerhamOptions +from struphy.io.options import FieldsBackground +from struphy.initial import perturbations +from struphy.kinetic_background import maxwellians +from struphy.pic.utilities import (LoadingParameters, + WeightsParameters, + BoundaryParameters, + BinningPlot, + KernelDensityPlot, + ) +from struphy import main + +# import model, set verbosity +from struphy.models.hybrid import LinearMHDDriftkineticCC + +# environment options +env = EnvironmentOptions( + out_folders='/tokp/work/bkna/tae', + sim_folder='tae_hybrid_sonic', + save_step=8, max_runtime=2840) + +# units +base_units = BaseUnits(n=0.2) + +# time stepping +time_opts = Time(dt=0.25, Tend=2000., split_algo='Strang') + +# geometry +domain = domains.HollowTorus(a1=0.1, a2=1., R0=10., sfl=False, tor_period=6) + +# fluid equilibrium (can be used as part of initial conditions) +equil = equils.AdhocTorus( + a=1., + R0=10., + B0=3., + q_kind=0, + q0=1.71, + q1=1.87, + n1=0., + n2=0., + na=1., + p_kind=1., + p0=1., + p1=0.95, + p2=0.05, + beta=0.0018, +) + +# grid +grid = grids.TensorProductGrid( + Nel=(24,96,16), + mpi_dims_mask=(True, True, False) +) + +# derham options +derham_opts = DerhamOptions( + p=(3, 3, 3), + spl_kind=(False, True, True), + dirichlet_bc=((True, True), (False, False), (False, False)), + nq_pr=(6, 6, 6), + nquads=(6, 6, 6), + polar_ck=-1, +) + +# light-weight model instance +model = LinearMHDDriftkineticCC( + turn_off = (#'PushGuidingCenterBxEstar', + #'PushGuidingCenterParallel', + #'ShearAlfvenCurrentCoupling5D', + #'Magnetosonic', + #'CurrentCoupling5DDensity', + #'CurrentCoupling5DGradB', + #'CurrentCoupling5DCurlb', + ) +) + +# species parameters +model.mhd.set_phys_params() +model.energetic_ions.set_phys_params( + mass_number=2., +) + +loading_params = LoadingParameters( + ppc=100, + moments=(0., 0., 0.89740839, 0.89740839), + seed=1234, + ) +weights_params = WeightsParameters( + control_variate=True, +) +boundary_params = BoundaryParameters( + bc=("remove", "periodic", "periodic"), + bc_refill=("outer", "inner"), +) +model.energetic_ions.set_markers(loading_params=loading_params, + weights_params=weights_params, + boundary_params=boundary_params, + bufsize=1., + n_cols_diag=3, + n_cols_aux=5, + ) +model.energetic_ions.set_sorting_boxes(boxes_per_dim=None) +model.energetic_ions.set_save_data() + +# calculate constant +ep_scale = model.energetic_ions.var.species.mass_number / model.mhd.mass_number + +# params else +from struphy.pic.accumulation.filter import FilterParameters +from struphy.linear_algebra.solver import SolverParameters, DiscreteGradientSolverParameters +filter_params = FilterParameters(use_filter="hybrid",) +solver_params = SolverParameters() + +# propagator options +model.propagators.push_bxe.options = model.propagators.push_bxe.Options( + algo = 'explicit', + b_tilde = model.em_fields.b_field,) +model.propagators.push_parallel.options = model.propagators.push_parallel.Options( + algo = 'explicit', + b_tilde = model.em_fields.b_field,) +model.propagators.shearalfen_cc5d.options = model.propagators.shearalfen_cc5d.Options( + filter_params=filter_params, + ep_scale=ep_scale, + energetic_ions = model.energetic_ions.var, + solver_params=solver_params, + nonlinear=False,) +model.propagators.magnetosonic.options = model.propagators.magnetosonic.Options( + b_field=model.em_fields.b_field,) +model.propagators.cc5d_density.options = model.propagators.cc5d_density.Options( + filter_params=filter_params, + ep_scale=ep_scale, + energetic_ions = model.energetic_ions.var, + b_tilde = model.em_fields.b_field,) +model.propagators.cc5d_gradb.options = model.propagators.cc5d_gradb.Options( + filter_params=filter_params, + ep_scale=ep_scale, + b_tilde = model.em_fields.b_field,) +model.propagators.cc5d_curlb.options = model.propagators.cc5d_curlb.Options( + filter_params=filter_params, + ep_scale=ep_scale, + b_tilde = model.em_fields.b_field,) + +# background, perturbations and initial conditions +model.mhd.velocity.add_background(FieldsBackground( + type="LogicalConst", + values=(0., 0., 0.) +)) +model.mhd.velocity.add_perturbation( + perturbations.TorusModesSin( + given_in_basis="2", + amps=(0.001,0.001), + ms=(10,11), + ns=(-1,-1), + pfuns=('exp','exp'), + pfun_params=((0.44444444,0.1),(0.44444444,0.1)), + comp=0) +) +model.mhd.velocity.add_perturbation( + perturbations.TorusModesCos( + given_in_basis="2", + amps=(0.000015915494309189535,0.000014468631190172303), + ms=(10,11), + ns=(-1,-1), + pfuns=('d_exp','d_exp'), + pfun_params=((0.44444444,0.1),(0.44444444,0.1)), + comp=1), +) +n_profile = perturbations.ITPA_density( + n0=0.00720652, + c=(0.47023, 0.20323, 0.13273, 0.521298) +) +maxwellian_1 = maxwellians.CanonicalMaxwellian2D( + vth=(0.89740839, None), + equil=equil, + volume_form=True, + n=(0., n_profile), +) +background = maxwellian_1 +model.energetic_ions.var.add_background( + background=background, + n_as_volume_form=False) + +# optional: exclude variables from saving +# model.energetic_ions.var.save_data = False + +if __name__ == "__main__": + # start run + verbose = True + + main.run(model, + params_path=__file__, + env=env, + base_units=base_units, + time_opts=time_opts, + domain=domain, + equil=equil, + grid=grid, + derham_opts=derham_opts, + verbose=verbose, + ) diff --git a/ITPA_benchmark/tae_hybrid_sonic_dt0125.py b/ITPA_benchmark/tae_hybrid_sonic_dt0125.py new file mode 100644 index 0000000..8a06bd1 --- /dev/null +++ b/ITPA_benchmark/tae_hybrid_sonic_dt0125.py @@ -0,0 +1,206 @@ +from struphy.io.options import EnvironmentOptions, BaseUnits, Time +from struphy.geometry import domains +from struphy.fields_background import equils +from struphy.topology import grids +from struphy.io.options import DerhamOptions +from struphy.io.options import FieldsBackground +from struphy.initial import perturbations +from struphy.kinetic_background import maxwellians +from struphy.pic.utilities import (LoadingParameters, + WeightsParameters, + BoundaryParameters, + BinningPlot, + KernelDensityPlot, + ) +from struphy import main + +# import model, set verbosity +from struphy.models.hybrid import LinearMHDDriftkineticCC + +# environment options +env = EnvironmentOptions( + out_folders='/tokp/work/bkna/tae', + sim_folder='tae_hybrid_sonic_dt0125', + save_step=8, + max_runtime=2840, + restart=False) + +# units +base_units = BaseUnits(n=0.2) + +# time stepping +time_opts = Time(dt=0.125, Tend=2000., split_algo='Strang') + +# geometry +domain = domains.HollowTorus(a1=0.1, a2=1., R0=10., sfl=False, tor_period=6) + +# fluid equilibrium (can be used as part of initial conditions) +equil = equils.AdhocTorus( + a=1., + R0=10., + B0=3., + q_kind=0, + q0=1.71, + q1=1.87, + n1=0., + n2=0., + na=1., + p_kind=1., + p0=1., + p1=0.95, + p2=0.05, + beta=0.0018, +) + +# grid +grid = grids.TensorProductGrid( + Nel=(24,96,16), + mpi_dims_mask=(True, True, False) +) + +# derham options +derham_opts = DerhamOptions( + p=(3, 3, 3), + spl_kind=(False, True, True), + dirichlet_bc=((True, True), (False, False), (False, False)), + nq_pr=(6, 6, 6), + nquads=(6, 6, 6), + polar_ck=-1, +) + +# light-weight model instance +model = LinearMHDDriftkineticCC( + turn_off = (#'PushGuidingCenterBxEstar', + #'PushGuidingCenterParallel', + #'ShearAlfvenCurrentCoupling5D', + #'Magnetosonic', + #'CurrentCoupling5DDensity', + #'CurrentCoupling5DGradB', + #'CurrentCoupling5DCurlb', + ) +) + +# species parameters +model.mhd.set_phys_params() +model.energetic_ions.set_phys_params( + mass_number=2., +) + +loading_params = LoadingParameters( + ppc=100, + moments=(0., 0., 0.89740839, 0.89740839), + seed=1234, + ) +weights_params = WeightsParameters( + control_variate=True, +) +boundary_params = BoundaryParameters( + bc=("remove", "periodic", "periodic"), + bc_refill=("outer", "inner"), +) +model.energetic_ions.set_markers(loading_params=loading_params, + weights_params=weights_params, + boundary_params=boundary_params, + bufsize=1., + n_cols_diag=3, + n_cols_aux=5, + ) +model.energetic_ions.set_sorting_boxes(boxes_per_dim=None) +model.energetic_ions.set_save_data() + +# calculate constant +ep_scale = model.energetic_ions.var.species.mass_number / model.mhd.mass_number + +# params else +from struphy.pic.accumulation.filter import FilterParameters +from struphy.linear_algebra.solver import SolverParameters, DiscreteGradientSolverParameters +filter_params = FilterParameters(use_filter="hybrid",) +solver_params = SolverParameters() + +# propagator options +model.propagators.push_bxe.options = model.propagators.push_bxe.Options( + algo = 'explicit', + b_tilde = model.em_fields.b_field,) +model.propagators.push_parallel.options = model.propagators.push_parallel.Options( + algo = 'explicit', + b_tilde = model.em_fields.b_field,) +model.propagators.shearalfen_cc5d.options = model.propagators.shearalfen_cc5d.Options( + filter_params=filter_params, + ep_scale=ep_scale, + energetic_ions = model.energetic_ions.var, + solver_params=solver_params, + nonlinear=False,) +model.propagators.magnetosonic.options = model.propagators.magnetosonic.Options( + b_field=model.em_fields.b_field,) +model.propagators.cc5d_density.options = model.propagators.cc5d_density.Options( + filter_params=filter_params, + ep_scale=ep_scale, + energetic_ions = model.energetic_ions.var, + b_tilde = model.em_fields.b_field,) +model.propagators.cc5d_gradb.options = model.propagators.cc5d_gradb.Options( + filter_params=filter_params, + ep_scale=ep_scale, + b_tilde = model.em_fields.b_field,) +model.propagators.cc5d_curlb.options = model.propagators.cc5d_curlb.Options( + filter_params=filter_params, + ep_scale=ep_scale, + b_tilde = model.em_fields.b_field,) + +# background, perturbations and initial conditions +model.mhd.velocity.add_background(FieldsBackground( + type="LogicalConst", + values=(0., 0., 0.) +)) +model.mhd.velocity.add_perturbation( + perturbations.TorusModesSin( + given_in_basis="2", + amps=(0.001,0.001), + ms=(10,11), + ns=(-1,-1), + pfuns=('exp','exp'), + pfun_params=((0.44444444,0.1),(0.44444444,0.1)), + comp=0) +) +model.mhd.velocity.add_perturbation( + perturbations.TorusModesCos( + given_in_basis="2", + amps=(0.000015915494309189535,0.000014468631190172303), + ms=(10,11), + ns=(-1,-1), + pfuns=('d_exp','d_exp'), + pfun_params=((0.44444444,0.1),(0.44444444,0.1)), + comp=1), +) +n_profile = perturbations.ITPA_density( + n0=0.00720652, + c=(0.47023, 0.20323, 0.13273, 0.521298) +) +maxwellian_1 = maxwellians.CanonicalMaxwellian2D( + vth=(0.89740839, None), + equil=equil, + volume_form=True, + n=(0., n_profile), +) +background = maxwellian_1 +model.energetic_ions.var.add_background( + background=background, + n_as_volume_form=False) + +# optional: exclude variables from saving +# model.energetic_ions.var.save_data = False + +if __name__ == "__main__": + # start run + verbose = True + + main.run(model, + params_path=__file__, + env=env, + base_units=base_units, + time_opts=time_opts, + domain=domain, + equil=equil, + grid=grid, + derham_opts=derham_opts, + verbose=verbose, + ) diff --git a/ITPA_benchmark/tae_hybrid_sonic_gradp0.py b/ITPA_benchmark/tae_hybrid_sonic_gradp0.py new file mode 100644 index 0000000..aa231d8 --- /dev/null +++ b/ITPA_benchmark/tae_hybrid_sonic_gradp0.py @@ -0,0 +1,205 @@ +from struphy.io.options import EnvironmentOptions, BaseUnits, Time +from struphy.geometry import domains +from struphy.fields_background import equils +from struphy.topology import grids +from struphy.io.options import DerhamOptions +from struphy.io.options import FieldsBackground +from struphy.initial import perturbations +from struphy.kinetic_background import maxwellians +from struphy.pic.utilities import (LoadingParameters, + WeightsParameters, + BoundaryParameters, + BinningPlot, + KernelDensityPlot, + ) +from struphy import main + +# import model, set verbosity +from struphy.models.hybrid import LinearMHDDriftkineticCC + +# environment options +env = EnvironmentOptions( + out_folders='/tokp/work/bkna/tae', + sim_folder='tae_hybrid_sonic_gradp0', + save_step=8, + max_runtime=2840) + +# units +base_units = BaseUnits(n=0.2) + +# time stepping +time_opts = Time(dt=0.25, Tend=2000., split_algo='Strang') + +# geometry +domain = domains.HollowTorus(a1=0.1, a2=1., R0=10., sfl=False, tor_period=6) + +# fluid equilibrium (can be used as part of initial conditions) +equil = equils.AdhocTorus( + a=1., + R0=10., + B0=3., + q_kind=0, + q0=1.71, + q1=1.87, + n1=0., + n2=0., + na=1., + p_kind=1., + p0=1., + p1=0., + p2=0., + beta=0.0018, +) + +# grid +grid = grids.TensorProductGrid( + Nel=(24,96,16), + mpi_dims_mask=(True, True, False) +) + +# derham options +derham_opts = DerhamOptions( + p=(3, 3, 3), + spl_kind=(False, True, True), + dirichlet_bc=((True, True), (False, False), (False, False)), + nq_pr=(6, 6, 6), + nquads=(6, 6, 6), + polar_ck=-1, +) + +# light-weight model instance +model = LinearMHDDriftkineticCC( + turn_off = (#'PushGuidingCenterBxEstar', + #'PushGuidingCenterParallel', + #'ShearAlfvenCurrentCoupling5D', + #'Magnetosonic', + #'CurrentCoupling5DDensity', + #'CurrentCoupling5DGradB', + #'CurrentCoupling5DCurlb', + ) +) + +# species parameters +model.mhd.set_phys_params() +model.energetic_ions.set_phys_params( + mass_number=2., +) + +loading_params = LoadingParameters( + ppc=100, + moments=(0., 0., 0.89740839, 0.89740839), + seed=1234, + ) +weights_params = WeightsParameters( + control_variate=True, +) +boundary_params = BoundaryParameters( + bc=("remove", "periodic", "periodic"), + bc_refill=("outer", "inner"), +) +model.energetic_ions.set_markers(loading_params=loading_params, + weights_params=weights_params, + boundary_params=boundary_params, + bufsize=1., + n_cols_diag=3, + n_cols_aux=5, + ) +model.energetic_ions.set_sorting_boxes(boxes_per_dim=None) +model.energetic_ions.set_save_data() + +# calculate constant +ep_scale = model.energetic_ions.var.species.mass_number / model.mhd.mass_number + +# params else +from struphy.pic.accumulation.filter import FilterParameters +from struphy.linear_algebra.solver import SolverParameters, DiscreteGradientSolverParameters +filter_params = FilterParameters(use_filter="hybrid",) +solver_params = SolverParameters() + +# propagator options +model.propagators.push_bxe.options = model.propagators.push_bxe.Options( + algo = 'explicit', + b_tilde = model.em_fields.b_field,) +model.propagators.push_parallel.options = model.propagators.push_parallel.Options( + algo = 'explicit', + b_tilde = model.em_fields.b_field,) +model.propagators.shearalfen_cc5d.options = model.propagators.shearalfen_cc5d.Options( + filter_params=filter_params, + ep_scale=ep_scale, + energetic_ions = model.energetic_ions.var, + solver_params=solver_params, + nonlinear=False,) +model.propagators.magnetosonic.options = model.propagators.magnetosonic.Options( + b_field=model.em_fields.b_field,) +model.propagators.cc5d_density.options = model.propagators.cc5d_density.Options( + filter_params=filter_params, + ep_scale=ep_scale, + energetic_ions = model.energetic_ions.var, + b_tilde = model.em_fields.b_field,) +model.propagators.cc5d_gradb.options = model.propagators.cc5d_gradb.Options( + filter_params=filter_params, + ep_scale=ep_scale, + b_tilde = model.em_fields.b_field,) +model.propagators.cc5d_curlb.options = model.propagators.cc5d_curlb.Options( + filter_params=filter_params, + ep_scale=ep_scale, + b_tilde = model.em_fields.b_field,) + +# background, perturbations and initial conditions +model.mhd.velocity.add_background(FieldsBackground( + type="LogicalConst", + values=(0., 0., 0.) +)) +model.mhd.velocity.add_perturbation( + perturbations.TorusModesSin( + given_in_basis="2", + amps=(0.001,0.001), + ms=(10,11), + ns=(-1,-1), + pfuns=('exp','exp'), + pfun_params=((0.44444444,0.1),(0.44444444,0.1)), + comp=0) +) +model.mhd.velocity.add_perturbation( + perturbations.TorusModesCos( + given_in_basis="2", + amps=(0.000015915494309189535,0.000014468631190172303), + ms=(10,11), + ns=(-1,-1), + pfuns=('d_exp','d_exp'), + pfun_params=((0.44444444,0.1),(0.44444444,0.1)), + comp=1), +) +n_profile = perturbations.ITPA_density( + n0=0.00720652, + c=(0.47023, 0.20323, 0.13273, 0.521298) +) +maxwellian_1 = maxwellians.CanonicalMaxwellian2D( + vth=(0.89740839, None), + equil=equil, + volume_form=True, + n=(0., n_profile), +) +background = maxwellian_1 +model.energetic_ions.var.add_background( + background=background, + n_as_volume_form=False) + +# optional: exclude variables from saving +# model.energetic_ions.var.save_data = False + +if __name__ == "__main__": + # start run + verbose = True + + main.run(model, + params_path=__file__, + env=env, + base_units=base_units, + time_opts=time_opts, + domain=domain, + equil=equil, + grid=grid, + derham_opts=derham_opts, + verbose=verbose, + ) diff --git a/ITPA_benchmark/tae_hybrid_sonic_gradp0_dt0125.py b/ITPA_benchmark/tae_hybrid_sonic_gradp0_dt0125.py new file mode 100644 index 0000000..1647fc8 --- /dev/null +++ b/ITPA_benchmark/tae_hybrid_sonic_gradp0_dt0125.py @@ -0,0 +1,206 @@ +from struphy.io.options import EnvironmentOptions, BaseUnits, Time +from struphy.geometry import domains +from struphy.fields_background import equils +from struphy.topology import grids +from struphy.io.options import DerhamOptions +from struphy.io.options import FieldsBackground +from struphy.initial import perturbations +from struphy.kinetic_background import maxwellians +from struphy.pic.utilities import (LoadingParameters, + WeightsParameters, + BoundaryParameters, + BinningPlot, + KernelDensityPlot, + ) +from struphy import main + +# import model, set verbosity +from struphy.models.hybrid import LinearMHDDriftkineticCC + +# environment options +env = EnvironmentOptions( + out_folders='/tokp/work/bkna/tae', + sim_folder='tae_hybrid_sonic_gradp0_dt0125', + save_step=8, + max_runtime=2840, + restart=False) + +# units +base_units = BaseUnits(n=0.2) + +# time stepping +time_opts = Time(dt=0.125, Tend=2000., split_algo='Strang') + +# geometry +domain = domains.HollowTorus(a1=0.1, a2=1., R0=10., sfl=False, tor_period=6) + +# fluid equilibrium (can be used as part of initial conditions) +equil = equils.AdhocTorus( + a=1., + R0=10., + B0=3., + q_kind=0, + q0=1.71, + q1=1.87, + n1=0., + n2=0., + na=1., + p_kind=1., + p0=1., + p1=0., + p2=0., + beta=0.0018, +) + +# grid +grid = grids.TensorProductGrid( + Nel=(24,96,16), + mpi_dims_mask=(True, True, False) +) + +# derham options +derham_opts = DerhamOptions( + p=(3, 3, 3), + spl_kind=(False, True, True), + dirichlet_bc=((True, True), (False, False), (False, False)), + nq_pr=(6, 6, 6), + nquads=(6, 6, 6), + polar_ck=-1, +) + +# light-weight model instance +model = LinearMHDDriftkineticCC( + turn_off = (#'PushGuidingCenterBxEstar', + #'PushGuidingCenterParallel', + #'ShearAlfvenCurrentCoupling5D', + #'Magnetosonic', + #'CurrentCoupling5DDensity', + #'CurrentCoupling5DGradB', + #'CurrentCoupling5DCurlb', + ) +) + +# species parameters +model.mhd.set_phys_params() +model.energetic_ions.set_phys_params( + mass_number=2., +) + +loading_params = LoadingParameters( + ppc=100, + moments=(0., 0., 0.89740839, 0.89740839), + seed=1234, + ) +weights_params = WeightsParameters( + control_variate=True, +) +boundary_params = BoundaryParameters( + bc=("remove", "periodic", "periodic"), + bc_refill=("outer", "inner"), +) +model.energetic_ions.set_markers(loading_params=loading_params, + weights_params=weights_params, + boundary_params=boundary_params, + bufsize=1., + n_cols_diag=3, + n_cols_aux=5, + ) +model.energetic_ions.set_sorting_boxes(boxes_per_dim=None) +model.energetic_ions.set_save_data() + +# calculate constant +ep_scale = model.energetic_ions.var.species.mass_number / model.mhd.mass_number + +# params else +from struphy.pic.accumulation.filter import FilterParameters +from struphy.linear_algebra.solver import SolverParameters, DiscreteGradientSolverParameters +filter_params = FilterParameters(use_filter="hybrid",) +solver_params = SolverParameters() + +# propagator options +model.propagators.push_bxe.options = model.propagators.push_bxe.Options( + algo = 'explicit', + b_tilde = model.em_fields.b_field,) +model.propagators.push_parallel.options = model.propagators.push_parallel.Options( + algo = 'explicit', + b_tilde = model.em_fields.b_field,) +model.propagators.shearalfen_cc5d.options = model.propagators.shearalfen_cc5d.Options( + filter_params=filter_params, + ep_scale=ep_scale, + energetic_ions = model.energetic_ions.var, + solver_params=solver_params, + nonlinear=False,) +model.propagators.magnetosonic.options = model.propagators.magnetosonic.Options( + b_field=model.em_fields.b_field,) +model.propagators.cc5d_density.options = model.propagators.cc5d_density.Options( + filter_params=filter_params, + ep_scale=ep_scale, + energetic_ions = model.energetic_ions.var, + b_tilde = model.em_fields.b_field,) +model.propagators.cc5d_gradb.options = model.propagators.cc5d_gradb.Options( + filter_params=filter_params, + ep_scale=ep_scale, + b_tilde = model.em_fields.b_field,) +model.propagators.cc5d_curlb.options = model.propagators.cc5d_curlb.Options( + filter_params=filter_params, + ep_scale=ep_scale, + b_tilde = model.em_fields.b_field,) + +# background, perturbations and initial conditions +model.mhd.velocity.add_background(FieldsBackground( + type="LogicalConst", + values=(0., 0., 0.) +)) +model.mhd.velocity.add_perturbation( + perturbations.TorusModesSin( + given_in_basis="2", + amps=(0.001,0.001), + ms=(10,11), + ns=(-1,-1), + pfuns=('exp','exp'), + pfun_params=((0.44444444,0.1),(0.44444444,0.1)), + comp=0) +) +model.mhd.velocity.add_perturbation( + perturbations.TorusModesCos( + given_in_basis="2", + amps=(0.000015915494309189535,0.000014468631190172303), + ms=(10,11), + ns=(-1,-1), + pfuns=('d_exp','d_exp'), + pfun_params=((0.44444444,0.1),(0.44444444,0.1)), + comp=1), +) +n_profile = perturbations.ITPA_density( + n0=0.00720652, + c=(0.47023, 0.20323, 0.13273, 0.521298) +) +maxwellian_1 = maxwellians.CanonicalMaxwellian2D( + vth=(0.89740839, None), + equil=equil, + volume_form=True, + n=(0., n_profile), +) +background = maxwellian_1 +model.energetic_ions.var.add_background( + background=background, + n_as_volume_form=False) + +# optional: exclude variables from saving +# model.energetic_ions.var.save_data = False + +if __name__ == "__main__": + # start run + verbose = True + + main.run(model, + params_path=__file__, + env=env, + base_units=base_units, + time_opts=time_opts, + domain=domain, + equil=equil, + grid=grid, + derham_opts=derham_opts, + verbose=verbose, + ) diff --git a/ITPA_benchmark/tae_hybrid_sonic_m2jboff.py b/ITPA_benchmark/tae_hybrid_sonic_m2jboff.py new file mode 100644 index 0000000..6a058e3 --- /dev/null +++ b/ITPA_benchmark/tae_hybrid_sonic_m2jboff.py @@ -0,0 +1,206 @@ +from struphy.io.options import EnvironmentOptions, BaseUnits, Time +from struphy.geometry import domains +from struphy.fields_background import equils +from struphy.topology import grids +from struphy.io.options import DerhamOptions +from struphy.io.options import FieldsBackground +from struphy.initial import perturbations +from struphy.kinetic_background import maxwellians +from struphy.pic.utilities import (LoadingParameters, + WeightsParameters, + BoundaryParameters, + BinningPlot, + KernelDensityPlot, + ) +from struphy import main + +# import model, set verbosity +from struphy.models.hybrid import LinearMHDDriftkineticCC + +# environment options +# 64 +env = EnvironmentOptions( + out_folders='/tokp/work/bkna/tae', + sim_folder='tae_hybrid_sonic_m2jboff', + save_step=8, + max_runtime=2840, + restart=True) + +# units +base_units = BaseUnits(n=0.2) + +# time stepping +time_opts = Time(dt=0.25, Tend=2000., split_algo='Strang') + +# geometry +domain = domains.HollowTorus(a1=0.1, a2=1., R0=10., sfl=False, tor_period=6) + +# fluid equilibrium (can be used as part of initial conditions) +equil = equils.AdhocTorus( + a=1., + R0=10., + B0=3., + q_kind=0, + q0=1.71, + q1=1.87, + n1=0., + n2=0., + na=1., + p_kind=1., + p0=1., + p1=0.95, + p2=0.05, + beta=0.0018, +) + +# grid +grid = grids.TensorProductGrid( + Nel=(24,96,16), + mpi_dims_mask=(True, True, False) +) + +# derham options +derham_opts = DerhamOptions( + p=(3, 3, 3), + spl_kind=(False, True, True), + dirichlet_bc=((True, True), (False, False), (False, False)), + nq_pr=(6, 6, 6), + nquads=(6, 6, 6), + polar_ck=-1, +) + +# light-weight model instance +model = LinearMHDDriftkineticCC( + turn_off = (#'PushGuidingCenterBxEstar', + #'PushGuidingCenterParallel', + #'ShearAlfvenCurrentCoupling5D', + #'Magnetosonic', + #'CurrentCoupling5DDensity', + #'CurrentCoupling5DGradB', + #'CurrentCoupling5DCurlb', + ) +) + +# species parameters +model.mhd.set_phys_params() +model.energetic_ions.set_phys_params( + mass_number=2., +) + +loading_params = LoadingParameters( + ppc=100, + moments=(0., 0., 0.89740839, 0.89740839), + seed=1234, + ) +weights_params = WeightsParameters( + control_variate=True, +) +boundary_params = BoundaryParameters( + bc=("remove", "periodic", "periodic"), + bc_refill=("outer", "inner"), +) +model.energetic_ions.set_markers(loading_params=loading_params, + weights_params=weights_params, + boundary_params=boundary_params, + bufsize=1., + n_cols_diag=3, + n_cols_aux=5, + ) +model.energetic_ions.set_sorting_boxes(boxes_per_dim=None) +model.energetic_ions.set_save_data() + +# calculate constant +ep_scale = model.energetic_ions.var.species.mass_number / model.mhd.mass_number + +# params else +from struphy.pic.accumulation.filter import FilterParameters +from struphy.linear_algebra.solver import SolverParameters, DiscreteGradientSolverParameters +filter_params = FilterParameters(use_filter="hybrid",) +solver_params = SolverParameters() + +# propagator options +model.propagators.push_bxe.options = model.propagators.push_bxe.Options( + algo = 'explicit', + b_tilde = model.em_fields.b_field,) +model.propagators.push_parallel.options = model.propagators.push_parallel.Options( + algo = 'explicit', + b_tilde = model.em_fields.b_field,) +model.propagators.shearalfen_cc5d.options = model.propagators.shearalfen_cc5d.Options( + filter_params=filter_params, + ep_scale=ep_scale, + energetic_ions = model.energetic_ions.var, + solver_params=solver_params, + nonlinear=False,) +model.propagators.magnetosonic.options = model.propagators.magnetosonic.Options() +model.propagators.cc5d_density.options = model.propagators.cc5d_density.Options( + filter_params=filter_params, + ep_scale=ep_scale, + energetic_ions = model.energetic_ions.var, + b_tilde = model.em_fields.b_field,) +model.propagators.cc5d_gradb.options = model.propagators.cc5d_gradb.Options( + filter_params=filter_params, + ep_scale=ep_scale, + b_tilde = model.em_fields.b_field,) +model.propagators.cc5d_curlb.options = model.propagators.cc5d_curlb.Options( + filter_params=filter_params, + ep_scale=ep_scale, + b_tilde = model.em_fields.b_field,) + +# background, perturbations and initial conditions +model.mhd.velocity.add_background(FieldsBackground( + type="LogicalConst", + values=(0., 0., 0.) +)) +model.mhd.velocity.add_perturbation( + perturbations.TorusModesSin( + given_in_basis="2", + amps=(0.001,0.001), + ms=(10,11), + ns=(-1,-1), + pfuns=('exp','exp'), + pfun_params=((0.44444444,0.1),(0.44444444,0.1)), + comp=0) +) +model.mhd.velocity.add_perturbation( + perturbations.TorusModesCos( + given_in_basis="2", + amps=(0.000015915494309189535,0.000014468631190172303), + ms=(10,11), + ns=(-1,-1), + pfuns=('d_exp','d_exp'), + pfun_params=((0.44444444,0.1),(0.44444444,0.1)), + comp=1), +) +n_profile = perturbations.ITPA_density( + n0=0.00720652, + c=(0.47023, 0.20323, 0.13273, 0.521298) +) +maxwellian_1 = maxwellians.CanonicalMaxwellian2D( + vth=(0.89740839, None), + equil=equil, + volume_form=True, + n=(0., n_profile), +) +background = maxwellian_1 +model.energetic_ions.var.add_background( + background=background, + n_as_volume_form=False) + +# optional: exclude variables from saving +# model.energetic_ions.var.save_data = False + +if __name__ == "__main__": + # start run + verbose = True + + main.run(model, + params_path=__file__, + env=env, + base_units=base_units, + time_opts=time_opts, + domain=domain, + equil=equil, + grid=grid, + derham_opts=derham_opts, + verbose=verbose, + ) diff --git a/ITPA_benchmark/tae_hybrid_sonic_m2jboff_dt0125.py b/ITPA_benchmark/tae_hybrid_sonic_m2jboff_dt0125.py new file mode 100644 index 0000000..6f13515 --- /dev/null +++ b/ITPA_benchmark/tae_hybrid_sonic_m2jboff_dt0125.py @@ -0,0 +1,206 @@ +from struphy.io.options import EnvironmentOptions, BaseUnits, Time +from struphy.geometry import domains +from struphy.fields_background import equils +from struphy.topology import grids +from struphy.io.options import DerhamOptions +from struphy.io.options import FieldsBackground +from struphy.initial import perturbations +from struphy.kinetic_background import maxwellians +from struphy.pic.utilities import (LoadingParameters, + WeightsParameters, + BoundaryParameters, + BinningPlot, + KernelDensityPlot, + ) +from struphy import main + +# import model, set verbosity +from struphy.models.hybrid import LinearMHDDriftkineticCC + +# environment options +# 64 +env = EnvironmentOptions( + out_folders='/tokp/work/bkna/tae', + sim_folder='tae_hybrid_sonic_m2jboff_dt0125', + save_step=8, + max_runtime=2840, + restart=False) + +# units +base_units = BaseUnits(n=0.2) + +# time stepping +time_opts = Time(dt=0.125, Tend=2000., split_algo='Strang') + +# geometry +domain = domains.HollowTorus(a1=0.1, a2=1., R0=10., sfl=False, tor_period=6) + +# fluid equilibrium (can be used as part of initial conditions) +equil = equils.AdhocTorus( + a=1., + R0=10., + B0=3., + q_kind=0, + q0=1.71, + q1=1.87, + n1=0., + n2=0., + na=1., + p_kind=1., + p0=1., + p1=0.95, + p2=0.05, + beta=0.0018, +) + +# grid +grid = grids.TensorProductGrid( + Nel=(24,96,16), + mpi_dims_mask=(True, True, False) +) + +# derham options +derham_opts = DerhamOptions( + p=(3, 3, 3), + spl_kind=(False, True, True), + dirichlet_bc=((True, True), (False, False), (False, False)), + nq_pr=(6, 6, 6), + nquads=(6, 6, 6), + polar_ck=-1, +) + +# light-weight model instance +model = LinearMHDDriftkineticCC( + turn_off = (#'PushGuidingCenterBxEstar', + #'PushGuidingCenterParallel', + #'ShearAlfvenCurrentCoupling5D', + #'Magnetosonic', + #'CurrentCoupling5DDensity', + #'CurrentCoupling5DGradB', + #'CurrentCoupling5DCurlb', + ) +) + +# species parameters +model.mhd.set_phys_params() +model.energetic_ions.set_phys_params( + mass_number=2., +) + +loading_params = LoadingParameters( + ppc=100, + moments=(0., 0., 0.89740839, 0.89740839), + seed=1234, + ) +weights_params = WeightsParameters( + control_variate=True, +) +boundary_params = BoundaryParameters( + bc=("remove", "periodic", "periodic"), + bc_refill=("outer", "inner"), +) +model.energetic_ions.set_markers(loading_params=loading_params, + weights_params=weights_params, + boundary_params=boundary_params, + bufsize=1., + n_cols_diag=3, + n_cols_aux=5, + ) +model.energetic_ions.set_sorting_boxes(boxes_per_dim=None) +model.energetic_ions.set_save_data() + +# calculate constant +ep_scale = model.energetic_ions.var.species.mass_number / model.mhd.mass_number + +# params else +from struphy.pic.accumulation.filter import FilterParameters +from struphy.linear_algebra.solver import SolverParameters, DiscreteGradientSolverParameters +filter_params = FilterParameters(use_filter="hybrid",) +solver_params = SolverParameters() + +# propagator options +model.propagators.push_bxe.options = model.propagators.push_bxe.Options( + algo = 'explicit', + b_tilde = model.em_fields.b_field,) +model.propagators.push_parallel.options = model.propagators.push_parallel.Options( + algo = 'explicit', + b_tilde = model.em_fields.b_field,) +model.propagators.shearalfen_cc5d.options = model.propagators.shearalfen_cc5d.Options( + filter_params=filter_params, + ep_scale=ep_scale, + energetic_ions = model.energetic_ions.var, + solver_params=solver_params, + nonlinear=False,) +model.propagators.magnetosonic.options = model.propagators.magnetosonic.Options() +model.propagators.cc5d_density.options = model.propagators.cc5d_density.Options( + filter_params=filter_params, + ep_scale=ep_scale, + energetic_ions = model.energetic_ions.var, + b_tilde = model.em_fields.b_field,) +model.propagators.cc5d_gradb.options = model.propagators.cc5d_gradb.Options( + filter_params=filter_params, + ep_scale=ep_scale, + b_tilde = model.em_fields.b_field,) +model.propagators.cc5d_curlb.options = model.propagators.cc5d_curlb.Options( + filter_params=filter_params, + ep_scale=ep_scale, + b_tilde = model.em_fields.b_field,) + +# background, perturbations and initial conditions +model.mhd.velocity.add_background(FieldsBackground( + type="LogicalConst", + values=(0., 0., 0.) +)) +model.mhd.velocity.add_perturbation( + perturbations.TorusModesSin( + given_in_basis="2", + amps=(0.001,0.001), + ms=(10,11), + ns=(-1,-1), + pfuns=('exp','exp'), + pfun_params=((0.44444444,0.1),(0.44444444,0.1)), + comp=0) +) +model.mhd.velocity.add_perturbation( + perturbations.TorusModesCos( + given_in_basis="2", + amps=(0.000015915494309189535,0.000014468631190172303), + ms=(10,11), + ns=(-1,-1), + pfuns=('d_exp','d_exp'), + pfun_params=((0.44444444,0.1),(0.44444444,0.1)), + comp=1), +) +n_profile = perturbations.ITPA_density( + n0=0.00720652, + c=(0.47023, 0.20323, 0.13273, 0.521298) +) +maxwellian_1 = maxwellians.CanonicalMaxwellian2D( + vth=(0.89740839, None), + equil=equil, + volume_form=True, + n=(0., n_profile), +) +background = maxwellian_1 +model.energetic_ions.var.add_background( + background=background, + n_as_volume_form=False) + +# optional: exclude variables from saving +# model.energetic_ions.var.save_data = False + +if __name__ == "__main__": + # start run + verbose = True + + main.run(model, + params_path=__file__, + env=env, + base_units=base_units, + time_opts=time_opts, + domain=domain, + equil=equil, + grid=grid, + derham_opts=derham_opts, + verbose=verbose, + ) diff --git a/ITPA_benchmark/tae_tp.py b/ITPA_benchmark/tae_tp.py new file mode 100644 index 0000000..ba8cdda --- /dev/null +++ b/ITPA_benchmark/tae_tp.py @@ -0,0 +1,207 @@ +from struphy.io.options import EnvironmentOptions, BaseUnits, Time +from struphy.geometry import domains +from struphy.fields_background import equils +from struphy.topology import grids +from struphy.io.options import DerhamOptions +from struphy.io.options import FieldsBackground +from struphy.initial import perturbations +from struphy.kinetic_background import maxwellians +from struphy.pic.utilities import (LoadingParameters, + WeightsParameters, + BoundaryParameters, + BinningPlot, + KernelDensityPlot, + ) +from struphy import main + +# import model, set verbosity +# 32 +from struphy.models.hybrid import LinearMHDDriftkineticCC + +# environment options +env = EnvironmentOptions( + out_folders='/tokp/work/bkna/tae', + sim_folder='tae_tp', + save_step=8, + max_runtime=2840, + restart=True) + +# units +base_units = BaseUnits(n=0.2) + +# time stepping +time_opts = Time(dt=0.25, Tend=2000., split_algo='Strang') + +# geometry +domain = domains.HollowTorus(a1=0.1, a2=1., R0=10., sfl=False, tor_period=6) + +# fluid equilibrium (can be used as part of initial conditions) +equil = equils.AdhocTorus( + a=1., + R0=10., + B0=3., + q_kind=0, + q0=1.71, + q1=1.87, + n1=0., + n2=0., + na=1., + p_kind=1., + p0=1., + p1=0.95, + p2=0.05, + beta=0.0018, +) + +# grid +grid = grids.TensorProductGrid( + Nel=(24,96,16), + mpi_dims_mask=(True, True, False) +) + +# derham options +derham_opts = DerhamOptions( + p=(3, 3, 3), + spl_kind=(False, True, True), + dirichlet_bc=((True, True), (False, False), (False, False)), + nq_pr=(6, 6, 6), + nquads=(6, 6, 6), + polar_ck=-1, +) + +# light-weight model instance +model = LinearMHDDriftkineticCC( + turn_off = (#'PushGuidingCenterBxEstar', + #'PushGuidingCenterParallel', + #'ShearAlfvenCurrentCoupling5D', + 'Magnetosonic', + #'CurrentCoupling5DDensity', + #'CurrentCoupling5DGradB', + #'CurrentCoupling5DCurlb', + ) +) + +# species parameters +model.mhd.set_phys_params() +model.energetic_ions.set_phys_params( + mass_number=2., +) + +loading_params = LoadingParameters( + ppc=100, + moments=(0., 0., 0.89740839, 0.89740839), + seed=1234, + ) +weights_params = WeightsParameters( + control_variate=True, +) +boundary_params = BoundaryParameters( + bc=("remove", "periodic", "periodic"), + bc_refill=("outer", "inner"), +) +model.energetic_ions.set_markers(loading_params=loading_params, + weights_params=weights_params, + boundary_params=boundary_params, + bufsize=1., + n_cols_diag=3, + n_cols_aux=5, + ) +model.energetic_ions.set_sorting_boxes(boxes_per_dim=None) +model.energetic_ions.set_save_data() + +# calculate constant +ep_scale = model.energetic_ions.var.species.mass_number / model.mhd.mass_number + +# params else +from struphy.pic.accumulation.filter import FilterParameters +from struphy.linear_algebra.solver import SolverParameters, DiscreteGradientSolverParameters +filter_params = FilterParameters(use_filter="three_point",) +solver_params = SolverParameters() + +# propagator options +model.propagators.push_bxe.options = model.propagators.push_bxe.Options( + algo = 'explicit', + b_tilde = model.em_fields.b_field,) +model.propagators.push_parallel.options = model.propagators.push_parallel.Options( + algo = 'explicit', + b_tilde = model.em_fields.b_field,) +model.propagators.shearalfen_cc5d.options = model.propagators.shearalfen_cc5d.Options( + filter_params=filter_params, + ep_scale=ep_scale, + energetic_ions = model.energetic_ions.var, + solver_params=solver_params, + nonlinear=False,) +# model.propagators.magnetosonic.options = model.propagators.magnetosonic.Options( +# b_field=model.em_fields.b_field,) +model.propagators.cc5d_density.options = model.propagators.cc5d_density.Options( + filter_params=filter_params, + ep_scale=ep_scale, + energetic_ions = model.energetic_ions.var, + b_tilde = model.em_fields.b_field,) +model.propagators.cc5d_gradb.options = model.propagators.cc5d_gradb.Options( + filter_params=filter_params, + ep_scale=ep_scale, + b_tilde = model.em_fields.b_field,) +model.propagators.cc5d_curlb.options = model.propagators.cc5d_curlb.Options( + filter_params=filter_params, + ep_scale=ep_scale, + b_tilde = model.em_fields.b_field,) + +# background, perturbations and initial conditions +model.mhd.velocity.add_background(FieldsBackground( + type="LogicalConst", + values=(0., 0., 0.) +)) +model.mhd.velocity.add_perturbation( + perturbations.TorusModesSin( + given_in_basis="2", + amps=(0.001,0.001), + ms=(10,11), + ns=(-1,-1), + pfuns=('exp','exp'), + pfun_params=((0.44444444,0.1),(0.44444444,0.1)), + comp=0) +) +model.mhd.velocity.add_perturbation( + perturbations.TorusModesCos( + given_in_basis="2", + amps=(0.000015915494309189535,0.000014468631190172303), + ms=(10,11), + ns=(-1,-1), + pfuns=('d_exp','d_exp'), + pfun_params=((0.44444444,0.1),(0.44444444,0.1)), + comp=1), +) +n_profile = perturbations.ITPA_density( + n0=0.00720652, + c=(0.47023, 0.20323, 0.13273, 0.521298) +) +maxwellian_1 = maxwellians.CanonicalMaxwellian2D( + vth=(0.89740839, None), + equil=equil, + volume_form=True, + n=(0., n_profile), +) +background = maxwellian_1 +model.energetic_ions.var.add_background( + background=background, + n_as_volume_form=False) + +# optional: exclude variables from saving +# model.energetic_ions.var.save_data = False + +if __name__ == "__main__": + # start run + verbose = True + + main.run(model, + params_path=__file__, + env=env, + base_units=base_units, + time_opts=time_opts, + domain=domain, + equil=equil, + grid=grid, + derham_opts=derham_opts, + verbose=verbose, + ) diff --git a/ITPA_benchmark/tae_tp_nonlinear.py b/ITPA_benchmark/tae_tp_nonlinear.py new file mode 100644 index 0000000..52926a5 --- /dev/null +++ b/ITPA_benchmark/tae_tp_nonlinear.py @@ -0,0 +1,207 @@ +from struphy.io.options import EnvironmentOptions, BaseUnits, Time +from struphy.geometry import domains +from struphy.fields_background import equils +from struphy.topology import grids +from struphy.io.options import DerhamOptions +from struphy.io.options import FieldsBackground +from struphy.initial import perturbations +from struphy.kinetic_background import maxwellians +from struphy.pic.utilities import (LoadingParameters, + WeightsParameters, + BoundaryParameters, + BinningPlot, + KernelDensityPlot, + ) +from struphy import main + +# import model, set verbosity +# 32 +from struphy.models.hybrid import LinearMHDDriftkineticCC + +# environment options +env = EnvironmentOptions( + out_folders='/tokp/work/bkna/tae', + sim_folder='tae_tp_nonlinear', + save_step=8, + max_runtime=2840, + restart=True) + +# units +base_units = BaseUnits(n=0.2) + +# time stepping +time_opts = Time(dt=0.25, Tend=2000., split_algo='Strang') + +# geometry +domain = domains.HollowTorus(a1=0.1, a2=1., R0=10., sfl=False, tor_period=6) + +# fluid equilibrium (can be used as part of initial conditions) +equil = equils.AdhocTorus( + a=1., + R0=10., + B0=3., + q_kind=0, + q0=1.71, + q1=1.87, + n1=0., + n2=0., + na=1., + p_kind=1., + p0=1., + p1=0.95, + p2=0.05, + beta=0.0018, +) + +# grid +grid = grids.TensorProductGrid( + Nel=(24,96,16), + mpi_dims_mask=(True, True, False) +) + +# derham options +derham_opts = DerhamOptions( + p=(3, 3, 3), + spl_kind=(False, True, True), + dirichlet_bc=((True, True), (False, False), (False, False)), + nq_pr=(6, 6, 6), + nquads=(6, 6, 6), + polar_ck=-1, +) + +# light-weight model instance +model = LinearMHDDriftkineticCC( + turn_off = (#'PushGuidingCenterBxEstar', + #'PushGuidingCenterParallel', + #'ShearAlfvenCurrentCoupling5D', + 'Magnetosonic', + #'CurrentCoupling5DDensity', + #'CurrentCoupling5DGradB', + #'CurrentCoupling5DCurlb', + ) +) + +# species parameters +model.mhd.set_phys_params() +model.energetic_ions.set_phys_params( + mass_number=2., +) + +loading_params = LoadingParameters( + ppc=100, + moments=(0., 0., 0.89740839, 0.89740839), + seed=1234, + ) +weights_params = WeightsParameters( + control_variate=True, +) +boundary_params = BoundaryParameters( + bc=("remove", "periodic", "periodic"), + bc_refill=("outer", "inner"), +) +model.energetic_ions.set_markers(loading_params=loading_params, + weights_params=weights_params, + boundary_params=boundary_params, + bufsize=1., + n_cols_diag=3, + n_cols_aux=5, + ) +model.energetic_ions.set_sorting_boxes(boxes_per_dim=None) +model.energetic_ions.set_save_data() + +# calculate constant +ep_scale = model.energetic_ions.var.species.mass_number / model.mhd.mass_number + +# params else +from struphy.pic.accumulation.filter import FilterParameters +from struphy.linear_algebra.solver import SolverParameters, DiscreteGradientSolverParameters +filter_params = FilterParameters(use_filter="three_point",) +solver_params = SolverParameters() + +# propagator options +model.propagators.push_bxe.options = model.propagators.push_bxe.Options( + algo = 'explicit', + b_tilde = model.em_fields.b_field,) +model.propagators.push_parallel.options = model.propagators.push_parallel.Options( + algo = 'explicit', + b_tilde = model.em_fields.b_field,) +model.propagators.shearalfen_cc5d.options = model.propagators.shearalfen_cc5d.Options( + filter_params=filter_params, + ep_scale=ep_scale, + energetic_ions = model.energetic_ions.var, + solver_params=solver_params, + nonlinear=True,) +# model.propagators.magnetosonic.options = model.propagators.magnetosonic.Options( +# b_field=model.em_fields.b_field,) +model.propagators.cc5d_density.options = model.propagators.cc5d_density.Options( + filter_params=filter_params, + ep_scale=ep_scale, + energetic_ions = model.energetic_ions.var, + b_tilde = model.em_fields.b_field,) +model.propagators.cc5d_gradb.options = model.propagators.cc5d_gradb.Options( + filter_params=filter_params, + ep_scale=ep_scale, + b_tilde = model.em_fields.b_field,) +model.propagators.cc5d_curlb.options = model.propagators.cc5d_curlb.Options( + filter_params=filter_params, + ep_scale=ep_scale, + b_tilde = model.em_fields.b_field,) + +# background, perturbations and initial conditions +model.mhd.velocity.add_background(FieldsBackground( + type="LogicalConst", + values=(0., 0., 0.) +)) +model.mhd.velocity.add_perturbation( + perturbations.TorusModesSin( + given_in_basis="2", + amps=(0.001,0.001), + ms=(10,11), + ns=(-1,-1), + pfuns=('exp','exp'), + pfun_params=((0.44444444,0.1),(0.44444444,0.1)), + comp=0) +) +model.mhd.velocity.add_perturbation( + perturbations.TorusModesCos( + given_in_basis="2", + amps=(0.000015915494309189535,0.000014468631190172303), + ms=(10,11), + ns=(-1,-1), + pfuns=('d_exp','d_exp'), + pfun_params=((0.44444444,0.1),(0.44444444,0.1)), + comp=1), +) +n_profile = perturbations.ITPA_density( + n0=0.00720652, + c=(0.47023, 0.20323, 0.13273, 0.521298) +) +maxwellian_1 = maxwellians.CanonicalMaxwellian2D( + vth=(0.89740839, None), + equil=equil, + volume_form=True, + n=(0., n_profile), +) +background = maxwellian_1 +model.energetic_ions.var.add_background( + background=background, + n_as_volume_form=False) + +# optional: exclude variables from saving +# model.energetic_ions.var.save_data = False + +if __name__ == "__main__": + # start run + verbose = True + + main.run(model, + params_path=__file__, + env=env, + base_units=base_units, + time_opts=time_opts, + domain=domain, + equil=equil, + grid=grid, + derham_opts=derham_opts, + verbose=verbose, + ) diff --git a/ITPA_benchmark/tae_tp_nonlinear_sonic_dt0125.py b/ITPA_benchmark/tae_tp_nonlinear_sonic_dt0125.py new file mode 100644 index 0000000..6a3751c --- /dev/null +++ b/ITPA_benchmark/tae_tp_nonlinear_sonic_dt0125.py @@ -0,0 +1,206 @@ +from struphy.io.options import EnvironmentOptions, BaseUnits, Time +from struphy.geometry import domains +from struphy.fields_background import equils +from struphy.topology import grids +from struphy.io.options import DerhamOptions +from struphy.io.options import FieldsBackground +from struphy.initial import perturbations +from struphy.kinetic_background import maxwellians +from struphy.pic.utilities import (LoadingParameters, + WeightsParameters, + BoundaryParameters, + BinningPlot, + KernelDensityPlot, + ) +from struphy import main + +# import model, set verbosity +from struphy.models.hybrid import LinearMHDDriftkineticCC + +# environment options +env = EnvironmentOptions( + out_folders='/tokp/work/bkna/tae', + sim_folder='tae_tp_nonlinear_sonic_dt0125', + save_step=8, + max_runtime=2840, + restart=False) + +# units +base_units = BaseUnits(n=0.2) + +# time stepping +time_opts = Time(dt=0.125, Tend=1000., split_algo='Strang') + +# geometry +domain = domains.HollowTorus(a1=0.1, a2=1., R0=10., sfl=False, tor_period=6) + +# fluid equilibrium (can be used as part of initial conditions) +equil = equils.AdhocTorus( + a=1., + R0=10., + B0=3., + q_kind=0, + q0=1.71, + q1=1.87, + n1=0., + n2=0., + na=1., + p_kind=1., + p0=1., + p1=0.95, + p2=0.05, + beta=0.0018, +) + +# grid +grid = grids.TensorProductGrid( + Nel=(24,96,16), + mpi_dims_mask=(True, True, False) +) + +# derham options +derham_opts = DerhamOptions( + p=(3, 3, 3), + spl_kind=(False, True, True), + dirichlet_bc=((True, True), (False, False), (False, False)), + nq_pr=(6, 6, 6), + nquads=(6, 6, 6), + polar_ck=-1, +) + +# light-weight model instance +model = LinearMHDDriftkineticCC( + turn_off = (#'PushGuidingCenterBxEstar', + #'PushGuidingCenterParallel', + #'ShearAlfvenCurrentCoupling5D', + #'Magnetosonic', + #'CurrentCoupling5DDensity', + #'CurrentCoupling5DGradB', + #'CurrentCoupling5DCurlb', + ) +) + +# species parameters +model.mhd.set_phys_params() +model.energetic_ions.set_phys_params( + mass_number=2., +) + +loading_params = LoadingParameters( + ppc=100, + moments=(0., 0., 0.89740839, 0.89740839), + seed=1234, + ) +weights_params = WeightsParameters( + control_variate=True, +) +boundary_params = BoundaryParameters( + bc=("remove", "periodic", "periodic"), + bc_refill=("outer", "inner"), +) +model.energetic_ions.set_markers(loading_params=loading_params, + weights_params=weights_params, + boundary_params=boundary_params, + bufsize=1., + n_cols_diag=3, + n_cols_aux=5, + ) +model.energetic_ions.set_sorting_boxes(boxes_per_dim=None) +model.energetic_ions.set_save_data() + +# calculate constant +ep_scale = model.energetic_ions.var.species.mass_number / model.mhd.mass_number + +# params else +from struphy.pic.accumulation.filter import FilterParameters +from struphy.linear_algebra.solver import SolverParameters, DiscreteGradientSolverParameters +filter_params = FilterParameters(use_filter="three_point",) +solver_params = SolverParameters() + +# propagator options +model.propagators.push_bxe.options = model.propagators.push_bxe.Options( + algo = 'explicit', + b_tilde = model.em_fields.b_field,) +model.propagators.push_parallel.options = model.propagators.push_parallel.Options( + algo = 'explicit', + b_tilde = model.em_fields.b_field,) +model.propagators.shearalfen_cc5d.options = model.propagators.shearalfen_cc5d.Options( + filter_params=filter_params, + ep_scale=ep_scale, + energetic_ions = model.energetic_ions.var, + solver_params=solver_params, + nonlinear=True,) +model.propagators.magnetosonic.options = model.propagators.magnetosonic.Options( + b_field=model.em_fields.b_field,) +model.propagators.cc5d_density.options = model.propagators.cc5d_density.Options( + filter_params=filter_params, + ep_scale=ep_scale, + energetic_ions = model.energetic_ions.var, + b_tilde = model.em_fields.b_field,) +model.propagators.cc5d_gradb.options = model.propagators.cc5d_gradb.Options( + filter_params=filter_params, + ep_scale=ep_scale, + b_tilde = model.em_fields.b_field,) +model.propagators.cc5d_curlb.options = model.propagators.cc5d_curlb.Options( + filter_params=filter_params, + ep_scale=ep_scale, + b_tilde = model.em_fields.b_field,) + +# background, perturbations and initial conditions +model.mhd.velocity.add_background(FieldsBackground( + type="LogicalConst", + values=(0., 0., 0.) +)) +model.mhd.velocity.add_perturbation( + perturbations.TorusModesSin( + given_in_basis="2", + amps=(0.001,0.001), + ms=(10,11), + ns=(-1,-1), + pfuns=('exp','exp'), + pfun_params=((0.44444444,0.1),(0.44444444,0.1)), + comp=0) +) +model.mhd.velocity.add_perturbation( + perturbations.TorusModesCos( + given_in_basis="2", + amps=(0.000015915494309189535,0.000014468631190172303), + ms=(10,11), + ns=(-1,-1), + pfuns=('d_exp','d_exp'), + pfun_params=((0.44444444,0.1),(0.44444444,0.1)), + comp=1), +) +n_profile = perturbations.ITPA_density( + n0=0.00720652, + c=(0.47023, 0.20323, 0.13273, 0.521298) +) +maxwellian_1 = maxwellians.CanonicalMaxwellian2D( + vth=(0.89740839, None), + equil=equil, + volume_form=True, + n=(0., n_profile), +) +background = maxwellian_1 +model.energetic_ions.var.add_background( + background=background, + n_as_volume_form=False) + +# optional: exclude variables from saving +# model.energetic_ions.var.save_data = False + +if __name__ == "__main__": + # start run + verbose = True + + main.run(model, + params_path=__file__, + env=env, + base_units=base_units, + time_opts=time_opts, + domain=domain, + equil=equil, + grid=grid, + derham_opts=derham_opts, + verbose=verbose, + ) diff --git a/ITPA_benchmark/tae_tp_nonlinear_sonic_gradp0_dt0125.py b/ITPA_benchmark/tae_tp_nonlinear_sonic_gradp0_dt0125.py new file mode 100644 index 0000000..f34ff79 --- /dev/null +++ b/ITPA_benchmark/tae_tp_nonlinear_sonic_gradp0_dt0125.py @@ -0,0 +1,206 @@ +from struphy.io.options import EnvironmentOptions, BaseUnits, Time +from struphy.geometry import domains +from struphy.fields_background import equils +from struphy.topology import grids +from struphy.io.options import DerhamOptions +from struphy.io.options import FieldsBackground +from struphy.initial import perturbations +from struphy.kinetic_background import maxwellians +from struphy.pic.utilities import (LoadingParameters, + WeightsParameters, + BoundaryParameters, + BinningPlot, + KernelDensityPlot, + ) +from struphy import main + +# import model, set verbosity +from struphy.models.hybrid import LinearMHDDriftkineticCC + +# environment options +env = EnvironmentOptions( + out_folders='/tokp/work/bkna/tae', + sim_folder='tae_tp_nonlinear_sonic_gradp0_dt0125', + save_step=8, + max_runtime=2840, + restart=False) + +# units +base_units = BaseUnits(n=0.2) + +# time stepping +time_opts = Time(dt=0.125, Tend=1000., split_algo='Strang') + +# geometry +domain = domains.HollowTorus(a1=0.1, a2=1., R0=10., sfl=False, tor_period=6) + +# fluid equilibrium (can be used as part of initial conditions) +equil = equils.AdhocTorus( + a=1., + R0=10., + B0=3., + q_kind=0, + q0=1.71, + q1=1.87, + n1=0., + n2=0., + na=1., + p_kind=1., + p0=1., + p1=0., + p2=0., + beta=0.0018, +) + +# grid +grid = grids.TensorProductGrid( + Nel=(24,96,16), + mpi_dims_mask=(True, True, False) +) + +# derham options +derham_opts = DerhamOptions( + p=(3, 3, 3), + spl_kind=(False, True, True), + dirichlet_bc=((True, True), (False, False), (False, False)), + nq_pr=(6, 6, 6), + nquads=(6, 6, 6), + polar_ck=-1, +) + +# light-weight model instance +model = LinearMHDDriftkineticCC( + turn_off = (#'PushGuidingCenterBxEstar', + #'PushGuidingCenterParallel', + #'ShearAlfvenCurrentCoupling5D', + #'Magnetosonic', + #'CurrentCoupling5DDensity', + #'CurrentCoupling5DGradB', + #'CurrentCoupling5DCurlb', + ) +) + +# species parameters +model.mhd.set_phys_params() +model.energetic_ions.set_phys_params( + mass_number=2., +) + +loading_params = LoadingParameters( + ppc=100, + moments=(0., 0., 0.89740839, 0.89740839), + seed=1234, + ) +weights_params = WeightsParameters( + control_variate=True, +) +boundary_params = BoundaryParameters( + bc=("remove", "periodic", "periodic"), + bc_refill=("outer", "inner"), +) +model.energetic_ions.set_markers(loading_params=loading_params, + weights_params=weights_params, + boundary_params=boundary_params, + bufsize=1., + n_cols_diag=3, + n_cols_aux=5, + ) +model.energetic_ions.set_sorting_boxes(boxes_per_dim=None) +model.energetic_ions.set_save_data() + +# calculate constant +ep_scale = model.energetic_ions.var.species.mass_number / model.mhd.mass_number + +# params else +from struphy.pic.accumulation.filter import FilterParameters +from struphy.linear_algebra.solver import SolverParameters, DiscreteGradientSolverParameters +filter_params = FilterParameters(use_filter="three_point",) +solver_params = SolverParameters() + +# propagator options +model.propagators.push_bxe.options = model.propagators.push_bxe.Options( + algo = 'explicit', + b_tilde = model.em_fields.b_field,) +model.propagators.push_parallel.options = model.propagators.push_parallel.Options( + algo = 'explicit', + b_tilde = model.em_fields.b_field,) +model.propagators.shearalfen_cc5d.options = model.propagators.shearalfen_cc5d.Options( + filter_params=filter_params, + ep_scale=ep_scale, + energetic_ions = model.energetic_ions.var, + solver_params=solver_params, + nonlinear=True,) +model.propagators.magnetosonic.options = model.propagators.magnetosonic.Options( + b_field=model.em_fields.b_field,) +model.propagators.cc5d_density.options = model.propagators.cc5d_density.Options( + filter_params=filter_params, + ep_scale=ep_scale, + energetic_ions = model.energetic_ions.var, + b_tilde = model.em_fields.b_field,) +model.propagators.cc5d_gradb.options = model.propagators.cc5d_gradb.Options( + filter_params=filter_params, + ep_scale=ep_scale, + b_tilde = model.em_fields.b_field,) +model.propagators.cc5d_curlb.options = model.propagators.cc5d_curlb.Options( + filter_params=filter_params, + ep_scale=ep_scale, + b_tilde = model.em_fields.b_field,) + +# background, perturbations and initial conditions +model.mhd.velocity.add_background(FieldsBackground( + type="LogicalConst", + values=(0., 0., 0.) +)) +model.mhd.velocity.add_perturbation( + perturbations.TorusModesSin( + given_in_basis="2", + amps=(0.001,0.001), + ms=(10,11), + ns=(-1,-1), + pfuns=('exp','exp'), + pfun_params=((0.44444444,0.1),(0.44444444,0.1)), + comp=0) +) +model.mhd.velocity.add_perturbation( + perturbations.TorusModesCos( + given_in_basis="2", + amps=(0.000015915494309189535,0.000014468631190172303), + ms=(10,11), + ns=(-1,-1), + pfuns=('d_exp','d_exp'), + pfun_params=((0.44444444,0.1),(0.44444444,0.1)), + comp=1), +) +n_profile = perturbations.ITPA_density( + n0=0.00720652, + c=(0.47023, 0.20323, 0.13273, 0.521298) +) +maxwellian_1 = maxwellians.CanonicalMaxwellian2D( + vth=(0.89740839, None), + equil=equil, + volume_form=True, + n=(0., n_profile), +) +background = maxwellian_1 +model.energetic_ions.var.add_background( + background=background, + n_as_volume_form=False) + +# optional: exclude variables from saving +# model.energetic_ions.var.save_data = False + +if __name__ == "__main__": + # start run + verbose = True + + main.run(model, + params_path=__file__, + env=env, + base_units=base_units, + time_opts=time_opts, + domain=domain, + equil=equil, + grid=grid, + derham_opts=derham_opts, + verbose=verbose, + ) diff --git a/ITPA_benchmark/tae_tp_nonlinear_sonic_m2jboff.py b/ITPA_benchmark/tae_tp_nonlinear_sonic_m2jboff.py new file mode 100644 index 0000000..b656797 --- /dev/null +++ b/ITPA_benchmark/tae_tp_nonlinear_sonic_m2jboff.py @@ -0,0 +1,206 @@ +from struphy.io.options import EnvironmentOptions, BaseUnits, Time +from struphy.geometry import domains +from struphy.fields_background import equils +from struphy.topology import grids +from struphy.io.options import DerhamOptions +from struphy.io.options import FieldsBackground +from struphy.initial import perturbations +from struphy.kinetic_background import maxwellians +from struphy.pic.utilities import (LoadingParameters, + WeightsParameters, + BoundaryParameters, + BinningPlot, + KernelDensityPlot, + ) +from struphy import main + +# import model, set verbosity +# 64 +from struphy.models.hybrid import LinearMHDDriftkineticCC + +# environment options +env = EnvironmentOptions( + out_folders='/tokp/work/bkna/tae', + sim_folder='tae_tp_nonlinear_sonic_m2jboff', + save_step=8, + max_runtime=2840, + restart=True) + +# units +base_units = BaseUnits(n=0.2) + +# time stepping +time_opts = Time(dt=0.25, Tend=2000., split_algo='Strang') + +# geometry +domain = domains.HollowTorus(a1=0.1, a2=1., R0=10., sfl=False, tor_period=6) + +# fluid equilibrium (can be used as part of initial conditions) +equil = equils.AdhocTorus( + a=1., + R0=10., + B0=3., + q_kind=0, + q0=1.71, + q1=1.87, + n1=0., + n2=0., + na=1., + p_kind=1., + p0=1., + p1=0.95, + p2=0.05, + beta=0.0018, +) + +# grid +grid = grids.TensorProductGrid( + Nel=(24,96,16), + mpi_dims_mask=(True, True, False) +) + +# derham options +derham_opts = DerhamOptions( + p=(3, 3, 3), + spl_kind=(False, True, True), + dirichlet_bc=((True, True), (False, False), (False, False)), + nq_pr=(6, 6, 6), + nquads=(6, 6, 6), + polar_ck=-1, +) + +# light-weight model instance +model = LinearMHDDriftkineticCC( + turn_off = (#'PushGuidingCenterBxEstar', + #'PushGuidingCenterParallel', + #'ShearAlfvenCurrentCoupling5D', + #'Magnetosonic', + #'CurrentCoupling5DDensity', + #'CurrentCoupling5DGradB', + #'CurrentCoupling5DCurlb', + ) +) + +# species parameters +model.mhd.set_phys_params() +model.energetic_ions.set_phys_params( + mass_number=2., +) + +loading_params = LoadingParameters( + ppc=100, + moments=(0., 0., 0.89740839, 0.89740839), + seed=1234, + ) +weights_params = WeightsParameters( + control_variate=True, +) +boundary_params = BoundaryParameters( + bc=("remove", "periodic", "periodic"), + bc_refill=("outer", "inner"), +) +model.energetic_ions.set_markers(loading_params=loading_params, + weights_params=weights_params, + boundary_params=boundary_params, + bufsize=1., + n_cols_diag=3, + n_cols_aux=5, + ) +model.energetic_ions.set_sorting_boxes(boxes_per_dim=None) +model.energetic_ions.set_save_data() + +# calculate constant +ep_scale = model.energetic_ions.var.species.mass_number / model.mhd.mass_number + +# params else +from struphy.pic.accumulation.filter import FilterParameters +from struphy.linear_algebra.solver import SolverParameters, DiscreteGradientSolverParameters +filter_params = FilterParameters(use_filter="three_point",) +solver_params = SolverParameters() + +# propagator options +model.propagators.push_bxe.options = model.propagators.push_bxe.Options( + algo = 'explicit', + b_tilde = model.em_fields.b_field,) +model.propagators.push_parallel.options = model.propagators.push_parallel.Options( + algo = 'explicit', + b_tilde = model.em_fields.b_field,) +model.propagators.shearalfen_cc5d.options = model.propagators.shearalfen_cc5d.Options( + filter_params=filter_params, + ep_scale=ep_scale, + energetic_ions = model.energetic_ions.var, + solver_params=solver_params, + nonlinear=True,) +model.propagators.magnetosonic.options = model.propagators.magnetosonic.Options() +model.propagators.cc5d_density.options = model.propagators.cc5d_density.Options( + filter_params=filter_params, + ep_scale=ep_scale, + energetic_ions = model.energetic_ions.var, + b_tilde = model.em_fields.b_field,) +model.propagators.cc5d_gradb.options = model.propagators.cc5d_gradb.Options( + filter_params=filter_params, + ep_scale=ep_scale, + b_tilde = model.em_fields.b_field,) +model.propagators.cc5d_curlb.options = model.propagators.cc5d_curlb.Options( + filter_params=filter_params, + ep_scale=ep_scale, + b_tilde = model.em_fields.b_field,) + +# background, perturbations and initial conditions +model.mhd.velocity.add_background(FieldsBackground( + type="LogicalConst", + values=(0., 0., 0.) +)) +model.mhd.velocity.add_perturbation( + perturbations.TorusModesSin( + given_in_basis="2", + amps=(0.001,0.001), + ms=(10,11), + ns=(-1,-1), + pfuns=('exp','exp'), + pfun_params=((0.44444444,0.1),(0.44444444,0.1)), + comp=0) +) +model.mhd.velocity.add_perturbation( + perturbations.TorusModesCos( + given_in_basis="2", + amps=(0.000015915494309189535,0.000014468631190172303), + ms=(10,11), + ns=(-1,-1), + pfuns=('d_exp','d_exp'), + pfun_params=((0.44444444,0.1),(0.44444444,0.1)), + comp=1), +) +n_profile = perturbations.ITPA_density( + n0=0.00720652, + c=(0.47023, 0.20323, 0.13273, 0.521298) +) +maxwellian_1 = maxwellians.CanonicalMaxwellian2D( + vth=(0.89740839, None), + equil=equil, + volume_form=True, + n=(0., n_profile), +) +background = maxwellian_1 +model.energetic_ions.var.add_background( + background=background, + n_as_volume_form=False) + +# optional: exclude variables from saving +# model.energetic_ions.var.save_data = False + +if __name__ == "__main__": + # start run + verbose = True + + main.run(model, + params_path=__file__, + env=env, + base_units=base_units, + time_opts=time_opts, + domain=domain, + equil=equil, + grid=grid, + derham_opts=derham_opts, + verbose=verbose, + ) diff --git a/ITPA_benchmark/tae_tp_nonlinear_sonic_m2jboff_dt0125.py b/ITPA_benchmark/tae_tp_nonlinear_sonic_m2jboff_dt0125.py new file mode 100644 index 0000000..5c45483 --- /dev/null +++ b/ITPA_benchmark/tae_tp_nonlinear_sonic_m2jboff_dt0125.py @@ -0,0 +1,206 @@ +from struphy.io.options import EnvironmentOptions, BaseUnits, Time +from struphy.geometry import domains +from struphy.fields_background import equils +from struphy.topology import grids +from struphy.io.options import DerhamOptions +from struphy.io.options import FieldsBackground +from struphy.initial import perturbations +from struphy.kinetic_background import maxwellians +from struphy.pic.utilities import (LoadingParameters, + WeightsParameters, + BoundaryParameters, + BinningPlot, + KernelDensityPlot, + ) +from struphy import main + +# import model, set verbosity +# 64 +from struphy.models.hybrid import LinearMHDDriftkineticCC + +# environment options +env = EnvironmentOptions( + out_folders='/tokp/work/bkna/tae', + sim_folder='tae_tp_nonlinear_sonic_m2jboff_dt0125', + save_step=8, + max_runtime=2840, + restart=False) + +# units +base_units = BaseUnits(n=0.2) + +# time stepping +time_opts = Time(dt=0.125, Tend=2000., split_algo='Strang') + +# geometry +domain = domains.HollowTorus(a1=0.1, a2=1., R0=10., sfl=False, tor_period=6) + +# fluid equilibrium (can be used as part of initial conditions) +equil = equils.AdhocTorus( + a=1., + R0=10., + B0=3., + q_kind=0, + q0=1.71, + q1=1.87, + n1=0., + n2=0., + na=1., + p_kind=1., + p0=1., + p1=0.95, + p2=0.05, + beta=0.0018, +) + +# grid +grid = grids.TensorProductGrid( + Nel=(24,96,16), + mpi_dims_mask=(True, True, False) +) + +# derham options +derham_opts = DerhamOptions( + p=(3, 3, 3), + spl_kind=(False, True, True), + dirichlet_bc=((True, True), (False, False), (False, False)), + nq_pr=(6, 6, 6), + nquads=(6, 6, 6), + polar_ck=-1, +) + +# light-weight model instance +model = LinearMHDDriftkineticCC( + turn_off = (#'PushGuidingCenterBxEstar', + #'PushGuidingCenterParallel', + #'ShearAlfvenCurrentCoupling5D', + #'Magnetosonic', + #'CurrentCoupling5DDensity', + #'CurrentCoupling5DGradB', + #'CurrentCoupling5DCurlb', + ) +) + +# species parameters +model.mhd.set_phys_params() +model.energetic_ions.set_phys_params( + mass_number=2., +) + +loading_params = LoadingParameters( + ppc=100, + moments=(0., 0., 0.89740839, 0.89740839), + seed=1234, + ) +weights_params = WeightsParameters( + control_variate=True, +) +boundary_params = BoundaryParameters( + bc=("remove", "periodic", "periodic"), + bc_refill=("outer", "inner"), +) +model.energetic_ions.set_markers(loading_params=loading_params, + weights_params=weights_params, + boundary_params=boundary_params, + bufsize=1., + n_cols_diag=3, + n_cols_aux=5, + ) +model.energetic_ions.set_sorting_boxes(boxes_per_dim=None) +model.energetic_ions.set_save_data() + +# calculate constant +ep_scale = model.energetic_ions.var.species.mass_number / model.mhd.mass_number + +# params else +from struphy.pic.accumulation.filter import FilterParameters +from struphy.linear_algebra.solver import SolverParameters, DiscreteGradientSolverParameters +filter_params = FilterParameters(use_filter="three_point",) +solver_params = SolverParameters() + +# propagator options +model.propagators.push_bxe.options = model.propagators.push_bxe.Options( + algo = 'explicit', + b_tilde = model.em_fields.b_field,) +model.propagators.push_parallel.options = model.propagators.push_parallel.Options( + algo = 'explicit', + b_tilde = model.em_fields.b_field,) +model.propagators.shearalfen_cc5d.options = model.propagators.shearalfen_cc5d.Options( + filter_params=filter_params, + ep_scale=ep_scale, + energetic_ions = model.energetic_ions.var, + solver_params=solver_params, + nonlinear=True,) +model.propagators.magnetosonic.options = model.propagators.magnetosonic.Options() +model.propagators.cc5d_density.options = model.propagators.cc5d_density.Options( + filter_params=filter_params, + ep_scale=ep_scale, + energetic_ions = model.energetic_ions.var, + b_tilde = model.em_fields.b_field,) +model.propagators.cc5d_gradb.options = model.propagators.cc5d_gradb.Options( + filter_params=filter_params, + ep_scale=ep_scale, + b_tilde = model.em_fields.b_field,) +model.propagators.cc5d_curlb.options = model.propagators.cc5d_curlb.Options( + filter_params=filter_params, + ep_scale=ep_scale, + b_tilde = model.em_fields.b_field,) + +# background, perturbations and initial conditions +model.mhd.velocity.add_background(FieldsBackground( + type="LogicalConst", + values=(0., 0., 0.) +)) +model.mhd.velocity.add_perturbation( + perturbations.TorusModesSin( + given_in_basis="2", + amps=(0.001,0.001), + ms=(10,11), + ns=(-1,-1), + pfuns=('exp','exp'), + pfun_params=((0.44444444,0.1),(0.44444444,0.1)), + comp=0) +) +model.mhd.velocity.add_perturbation( + perturbations.TorusModesCos( + given_in_basis="2", + amps=(0.000015915494309189535,0.000014468631190172303), + ms=(10,11), + ns=(-1,-1), + pfuns=('d_exp','d_exp'), + pfun_params=((0.44444444,0.1),(0.44444444,0.1)), + comp=1), +) +n_profile = perturbations.ITPA_density( + n0=0.00720652, + c=(0.47023, 0.20323, 0.13273, 0.521298) +) +maxwellian_1 = maxwellians.CanonicalMaxwellian2D( + vth=(0.89740839, None), + equil=equil, + volume_form=True, + n=(0., n_profile), +) +background = maxwellian_1 +model.energetic_ions.var.add_background( + background=background, + n_as_volume_form=False) + +# optional: exclude variables from saving +# model.energetic_ions.var.save_data = False + +if __name__ == "__main__": + # start run + verbose = True + + main.run(model, + params_path=__file__, + env=env, + base_units=base_units, + time_opts=time_opts, + domain=domain, + equil=equil, + grid=grid, + derham_opts=derham_opts, + verbose=verbose, + ) diff --git a/ITPA_benchmark/tae_tp_sonic.py b/ITPA_benchmark/tae_tp_sonic.py new file mode 100644 index 0000000..31c2d07 --- /dev/null +++ b/ITPA_benchmark/tae_tp_sonic.py @@ -0,0 +1,204 @@ +from struphy.io.options import EnvironmentOptions, BaseUnits, Time +from struphy.geometry import domains +from struphy.fields_background import equils +from struphy.topology import grids +from struphy.io.options import DerhamOptions +from struphy.io.options import FieldsBackground +from struphy.initial import perturbations +from struphy.kinetic_background import maxwellians +from struphy.pic.utilities import (LoadingParameters, + WeightsParameters, + BoundaryParameters, + BinningPlot, + KernelDensityPlot, + ) +from struphy import main + +# import model, set verbosity +from struphy.models.hybrid import LinearMHDDriftkineticCC + +# environment options +env = EnvironmentOptions( + out_folders='/tokp/work/bkna/tae', + sim_folder='tae_tp_sonic', + save_step=8, max_runtime=2840) + +# units +base_units = BaseUnits(n=0.2) + +# time stepping +time_opts = Time(dt=0.25, Tend=1000., split_algo='Strang') + +# geometry +domain = domains.HollowTorus(a1=0.1, a2=1., R0=10., sfl=False, tor_period=6) + +# fluid equilibrium (can be used as part of initial conditions) +equil = equils.AdhocTorus( + a=1., + R0=10., + B0=3., + q_kind=0, + q0=1.71, + q1=1.87, + n1=0., + n2=0., + na=1., + p_kind=1., + p0=1., + p1=0.95, + p2=0.05, + beta=0.0018, +) + +# grid +grid = grids.TensorProductGrid( + Nel=(24,96,16), + mpi_dims_mask=(True, True, False) +) + +# derham options +derham_opts = DerhamOptions( + p=(3, 3, 3), + spl_kind=(False, True, True), + dirichlet_bc=((True, True), (False, False), (False, False)), + nq_pr=(6, 6, 6), + nquads=(6, 6, 6), + polar_ck=-1, +) + +# light-weight model instance +model = LinearMHDDriftkineticCC( + turn_off = (#'PushGuidingCenterBxEstar', + #'PushGuidingCenterParallel', + #'ShearAlfvenCurrentCoupling5D', + #'Magnetosonic', + #'CurrentCoupling5DDensity', + #'CurrentCoupling5DGradB', + #'CurrentCoupling5DCurlb', + ) +) + +# species parameters +model.mhd.set_phys_params() +model.energetic_ions.set_phys_params( + mass_number=2., +) + +loading_params = LoadingParameters( + ppc=100, + moments=(0., 0., 0.89740839, 0.89740839), + seed=1234, + ) +weights_params = WeightsParameters( + control_variate=True, +) +boundary_params = BoundaryParameters( + bc=("remove", "periodic", "periodic"), + bc_refill=("outer", "inner"), +) +model.energetic_ions.set_markers(loading_params=loading_params, + weights_params=weights_params, + boundary_params=boundary_params, + bufsize=1., + n_cols_diag=3, + n_cols_aux=5, + ) +model.energetic_ions.set_sorting_boxes(boxes_per_dim=None) +model.energetic_ions.set_save_data() + +# calculate constant +ep_scale = model.energetic_ions.var.species.mass_number / model.mhd.mass_number + +# params else +from struphy.pic.accumulation.filter import FilterParameters +from struphy.linear_algebra.solver import SolverParameters, DiscreteGradientSolverParameters +filter_params = FilterParameters(use_filter="three_point",) +solver_params = SolverParameters() + +# propagator options +model.propagators.push_bxe.options = model.propagators.push_bxe.Options( + algo = 'explicit', + b_tilde = model.em_fields.b_field,) +model.propagators.push_parallel.options = model.propagators.push_parallel.Options( + algo = 'explicit', + b_tilde = model.em_fields.b_field,) +model.propagators.shearalfen_cc5d.options = model.propagators.shearalfen_cc5d.Options( + filter_params=filter_params, + ep_scale=ep_scale, + energetic_ions = model.energetic_ions.var, + solver_params=solver_params, + nonlinear=False,) +model.propagators.magnetosonic.options = model.propagators.magnetosonic.Options( + b_field=model.em_fields.b_field,) +model.propagators.cc5d_density.options = model.propagators.cc5d_density.Options( + filter_params=filter_params, + ep_scale=ep_scale, + energetic_ions = model.energetic_ions.var, + b_tilde = model.em_fields.b_field,) +model.propagators.cc5d_gradb.options = model.propagators.cc5d_gradb.Options( + filter_params=filter_params, + ep_scale=ep_scale, + b_tilde = model.em_fields.b_field,) +model.propagators.cc5d_curlb.options = model.propagators.cc5d_curlb.Options( + filter_params=filter_params, + ep_scale=ep_scale, + b_tilde = model.em_fields.b_field,) + +# background, perturbations and initial conditions +model.mhd.velocity.add_background(FieldsBackground( + type="LogicalConst", + values=(0., 0., 0.) +)) +model.mhd.velocity.add_perturbation( + perturbations.TorusModesSin( + given_in_basis="2", + amps=(0.001,0.001), + ms=(10,11), + ns=(-1,-1), + pfuns=('exp','exp'), + pfun_params=((0.44444444,0.1),(0.44444444,0.1)), + comp=0) +) +model.mhd.velocity.add_perturbation( + perturbations.TorusModesCos( + given_in_basis="2", + amps=(0.000015915494309189535,0.000014468631190172303), + ms=(10,11), + ns=(-1,-1), + pfuns=('d_exp','d_exp'), + pfun_params=((0.44444444,0.1),(0.44444444,0.1)), + comp=1), +) +n_profile = perturbations.ITPA_density( + n0=0.00720652, + c=(0.47023, 0.20323, 0.13273, 0.521298) +) +maxwellian_1 = maxwellians.CanonicalMaxwellian2D( + vth=(0.89740839, None), + equil=equil, + volume_form=True, + n=(0., n_profile), +) +background = maxwellian_1 +model.energetic_ions.var.add_background( + background=background, + n_as_volume_form=False) + +# optional: exclude variables from saving +# model.energetic_ions.var.save_data = False + +if __name__ == "__main__": + # start run + verbose = True + + main.run(model, + params_path=__file__, + env=env, + base_units=base_units, + time_opts=time_opts, + domain=domain, + equil=equil, + grid=grid, + derham_opts=derham_opts, + verbose=verbose, + ) diff --git a/ITPA_benchmark/tae_tp_sonic_dt0125.py b/ITPA_benchmark/tae_tp_sonic_dt0125.py new file mode 100644 index 0000000..2142ed6 --- /dev/null +++ b/ITPA_benchmark/tae_tp_sonic_dt0125.py @@ -0,0 +1,206 @@ +from struphy.io.options import EnvironmentOptions, BaseUnits, Time +from struphy.geometry import domains +from struphy.fields_background import equils +from struphy.topology import grids +from struphy.io.options import DerhamOptions +from struphy.io.options import FieldsBackground +from struphy.initial import perturbations +from struphy.kinetic_background import maxwellians +from struphy.pic.utilities import (LoadingParameters, + WeightsParameters, + BoundaryParameters, + BinningPlot, + KernelDensityPlot, + ) +from struphy import main + +# import model, set verbosity +from struphy.models.hybrid import LinearMHDDriftkineticCC + +# environment options +env = EnvironmentOptions( + out_folders='/tokp/work/bkna/tae', + sim_folder='tae_tp_sonic_dt0125', + save_step=8, + max_runtime=2840, + restart=False) + +# units +base_units = BaseUnits(n=0.2) + +# time stepping +time_opts = Time(dt=0.125, Tend=1000., split_algo='Strang') + +# geometry +domain = domains.HollowTorus(a1=0.1, a2=1., R0=10., sfl=False, tor_period=6) + +# fluid equilibrium (can be used as part of initial conditions) +equil = equils.AdhocTorus( + a=1., + R0=10., + B0=3., + q_kind=0, + q0=1.71, + q1=1.87, + n1=0., + n2=0., + na=1., + p_kind=1., + p0=1., + p1=0.95, + p2=0.05, + beta=0.0018, +) + +# grid +grid = grids.TensorProductGrid( + Nel=(24,96,16), + mpi_dims_mask=(True, True, False) +) + +# derham options +derham_opts = DerhamOptions( + p=(3, 3, 3), + spl_kind=(False, True, True), + dirichlet_bc=((True, True), (False, False), (False, False)), + nq_pr=(6, 6, 6), + nquads=(6, 6, 6), + polar_ck=-1, +) + +# light-weight model instance +model = LinearMHDDriftkineticCC( + turn_off = (#'PushGuidingCenterBxEstar', + #'PushGuidingCenterParallel', + #'ShearAlfvenCurrentCoupling5D', + #'Magnetosonic', + #'CurrentCoupling5DDensity', + #'CurrentCoupling5DGradB', + #'CurrentCoupling5DCurlb', + ) +) + +# species parameters +model.mhd.set_phys_params() +model.energetic_ions.set_phys_params( + mass_number=2., +) + +loading_params = LoadingParameters( + ppc=100, + moments=(0., 0., 0.89740839, 0.89740839), + seed=1234, + ) +weights_params = WeightsParameters( + control_variate=True, +) +boundary_params = BoundaryParameters( + bc=("remove", "periodic", "periodic"), + bc_refill=("outer", "inner"), +) +model.energetic_ions.set_markers(loading_params=loading_params, + weights_params=weights_params, + boundary_params=boundary_params, + bufsize=1., + n_cols_diag=3, + n_cols_aux=5, + ) +model.energetic_ions.set_sorting_boxes(boxes_per_dim=None) +model.energetic_ions.set_save_data() + +# calculate constant +ep_scale = model.energetic_ions.var.species.mass_number / model.mhd.mass_number + +# params else +from struphy.pic.accumulation.filter import FilterParameters +from struphy.linear_algebra.solver import SolverParameters, DiscreteGradientSolverParameters +filter_params = FilterParameters(use_filter="three_point",) +solver_params = SolverParameters() + +# propagator options +model.propagators.push_bxe.options = model.propagators.push_bxe.Options( + algo = 'explicit', + b_tilde = model.em_fields.b_field,) +model.propagators.push_parallel.options = model.propagators.push_parallel.Options( + algo = 'explicit', + b_tilde = model.em_fields.b_field,) +model.propagators.shearalfen_cc5d.options = model.propagators.shearalfen_cc5d.Options( + filter_params=filter_params, + ep_scale=ep_scale, + energetic_ions = model.energetic_ions.var, + solver_params=solver_params, + nonlinear=False,) +model.propagators.magnetosonic.options = model.propagators.magnetosonic.Options( + b_field=model.em_fields.b_field,) +model.propagators.cc5d_density.options = model.propagators.cc5d_density.Options( + filter_params=filter_params, + ep_scale=ep_scale, + energetic_ions = model.energetic_ions.var, + b_tilde = model.em_fields.b_field,) +model.propagators.cc5d_gradb.options = model.propagators.cc5d_gradb.Options( + filter_params=filter_params, + ep_scale=ep_scale, + b_tilde = model.em_fields.b_field,) +model.propagators.cc5d_curlb.options = model.propagators.cc5d_curlb.Options( + filter_params=filter_params, + ep_scale=ep_scale, + b_tilde = model.em_fields.b_field,) + +# background, perturbations and initial conditions +model.mhd.velocity.add_background(FieldsBackground( + type="LogicalConst", + values=(0., 0., 0.) +)) +model.mhd.velocity.add_perturbation( + perturbations.TorusModesSin( + given_in_basis="2", + amps=(0.001,0.001), + ms=(10,11), + ns=(-1,-1), + pfuns=('exp','exp'), + pfun_params=((0.44444444,0.1),(0.44444444,0.1)), + comp=0) +) +model.mhd.velocity.add_perturbation( + perturbations.TorusModesCos( + given_in_basis="2", + amps=(0.000015915494309189535,0.000014468631190172303), + ms=(10,11), + ns=(-1,-1), + pfuns=('d_exp','d_exp'), + pfun_params=((0.44444444,0.1),(0.44444444,0.1)), + comp=1), +) +n_profile = perturbations.ITPA_density( + n0=0.00720652, + c=(0.47023, 0.20323, 0.13273, 0.521298) +) +maxwellian_1 = maxwellians.CanonicalMaxwellian2D( + vth=(0.89740839, None), + equil=equil, + volume_form=True, + n=(0., n_profile), +) +background = maxwellian_1 +model.energetic_ions.var.add_background( + background=background, + n_as_volume_form=False) + +# optional: exclude variables from saving +# model.energetic_ions.var.save_data = False + +if __name__ == "__main__": + # start run + verbose = True + + main.run(model, + params_path=__file__, + env=env, + base_units=base_units, + time_opts=time_opts, + domain=domain, + equil=equil, + grid=grid, + derham_opts=derham_opts, + verbose=verbose, + ) diff --git a/ITPA_benchmark/tae_tp_sonic_gradp0.py b/ITPA_benchmark/tae_tp_sonic_gradp0.py new file mode 100644 index 0000000..5a942fe --- /dev/null +++ b/ITPA_benchmark/tae_tp_sonic_gradp0.py @@ -0,0 +1,205 @@ +from struphy.io.options import EnvironmentOptions, BaseUnits, Time +from struphy.geometry import domains +from struphy.fields_background import equils +from struphy.topology import grids +from struphy.io.options import DerhamOptions +from struphy.io.options import FieldsBackground +from struphy.initial import perturbations +from struphy.kinetic_background import maxwellians +from struphy.pic.utilities import (LoadingParameters, + WeightsParameters, + BoundaryParameters, + BinningPlot, + KernelDensityPlot, + ) +from struphy import main + +# import model, set verbosity +from struphy.models.hybrid import LinearMHDDriftkineticCC + +# environment options +env = EnvironmentOptions( + out_folders='/tokp/work/bkna/tae', + sim_folder='tae_tp_sonic_gradp0', + save_step=8, + max_runtime=2840) + +# units +base_units = BaseUnits(n=0.2) + +# time stepping +time_opts = Time(dt=0.25, Tend=1000., split_algo='Strang') + +# geometry +domain = domains.HollowTorus(a1=0.1, a2=1., R0=10., sfl=False, tor_period=6) + +# fluid equilibrium (can be used as part of initial conditions) +equil = equils.AdhocTorus( + a=1., + R0=10., + B0=3., + q_kind=0, + q0=1.71, + q1=1.87, + n1=0., + n2=0., + na=1., + p_kind=1., + p0=1., + p1=0., + p2=0., + beta=0.0018, +) + +# grid +grid = grids.TensorProductGrid( + Nel=(24,96,16), + mpi_dims_mask=(True, True, False) +) + +# derham options +derham_opts = DerhamOptions( + p=(3, 3, 3), + spl_kind=(False, True, True), + dirichlet_bc=((True, True), (False, False), (False, False)), + nq_pr=(6, 6, 6), + nquads=(6, 6, 6), + polar_ck=-1, +) + +# light-weight model instance +model = LinearMHDDriftkineticCC( + turn_off = (#'PushGuidingCenterBxEstar', + #'PushGuidingCenterParallel', + #'ShearAlfvenCurrentCoupling5D', + #'Magnetosonic', + #'CurrentCoupling5DDensity', + #'CurrentCoupling5DGradB', + #'CurrentCoupling5DCurlb', + ) +) + +# species parameters +model.mhd.set_phys_params() +model.energetic_ions.set_phys_params( + mass_number=2., +) + +loading_params = LoadingParameters( + ppc=100, + moments=(0., 0., 0.89740839, 0.89740839), + seed=1234, + ) +weights_params = WeightsParameters( + control_variate=True, +) +boundary_params = BoundaryParameters( + bc=("remove", "periodic", "periodic"), + bc_refill=("outer", "inner"), +) +model.energetic_ions.set_markers(loading_params=loading_params, + weights_params=weights_params, + boundary_params=boundary_params, + bufsize=1., + n_cols_diag=3, + n_cols_aux=5, + ) +model.energetic_ions.set_sorting_boxes(boxes_per_dim=None) +model.energetic_ions.set_save_data() + +# calculate constant +ep_scale = model.energetic_ions.var.species.mass_number / model.mhd.mass_number + +# params else +from struphy.pic.accumulation.filter import FilterParameters +from struphy.linear_algebra.solver import SolverParameters, DiscreteGradientSolverParameters +filter_params = FilterParameters(use_filter="three_point",) +solver_params = SolverParameters() + +# propagator options +model.propagators.push_bxe.options = model.propagators.push_bxe.Options( + algo = 'explicit', + b_tilde = model.em_fields.b_field,) +model.propagators.push_parallel.options = model.propagators.push_parallel.Options( + algo = 'explicit', + b_tilde = model.em_fields.b_field,) +model.propagators.shearalfen_cc5d.options = model.propagators.shearalfen_cc5d.Options( + filter_params=filter_params, + ep_scale=ep_scale, + energetic_ions = model.energetic_ions.var, + solver_params=solver_params, + nonlinear=False,) +model.propagators.magnetosonic.options = model.propagators.magnetosonic.Options( + b_field=model.em_fields.b_field,) +model.propagators.cc5d_density.options = model.propagators.cc5d_density.Options( + filter_params=filter_params, + ep_scale=ep_scale, + energetic_ions = model.energetic_ions.var, + b_tilde = model.em_fields.b_field,) +model.propagators.cc5d_gradb.options = model.propagators.cc5d_gradb.Options( + filter_params=filter_params, + ep_scale=ep_scale, + b_tilde = model.em_fields.b_field,) +model.propagators.cc5d_curlb.options = model.propagators.cc5d_curlb.Options( + filter_params=filter_params, + ep_scale=ep_scale, + b_tilde = model.em_fields.b_field,) + +# background, perturbations and initial conditions +model.mhd.velocity.add_background(FieldsBackground( + type="LogicalConst", + values=(0., 0., 0.) +)) +model.mhd.velocity.add_perturbation( + perturbations.TorusModesSin( + given_in_basis="2", + amps=(0.001,0.001), + ms=(10,11), + ns=(-1,-1), + pfuns=('exp','exp'), + pfun_params=((0.44444444,0.1),(0.44444444,0.1)), + comp=0) +) +model.mhd.velocity.add_perturbation( + perturbations.TorusModesCos( + given_in_basis="2", + amps=(0.000015915494309189535,0.000014468631190172303), + ms=(10,11), + ns=(-1,-1), + pfuns=('d_exp','d_exp'), + pfun_params=((0.44444444,0.1),(0.44444444,0.1)), + comp=1), +) +n_profile = perturbations.ITPA_density( + n0=0.00720652, + c=(0.47023, 0.20323, 0.13273, 0.521298) +) +maxwellian_1 = maxwellians.CanonicalMaxwellian2D( + vth=(0.89740839, None), + equil=equil, + volume_form=True, + n=(0., n_profile), +) +background = maxwellian_1 +model.energetic_ions.var.add_background( + background=background, + n_as_volume_form=False) + +# optional: exclude variables from saving +# model.energetic_ions.var.save_data = False + +if __name__ == "__main__": + # start run + verbose = True + + main.run(model, + params_path=__file__, + env=env, + base_units=base_units, + time_opts=time_opts, + domain=domain, + equil=equil, + grid=grid, + derham_opts=derham_opts, + verbose=verbose, + ) diff --git a/ITPA_benchmark/tae_tp_sonic_gradp0_dt0125.py b/ITPA_benchmark/tae_tp_sonic_gradp0_dt0125.py new file mode 100644 index 0000000..ce9eacb --- /dev/null +++ b/ITPA_benchmark/tae_tp_sonic_gradp0_dt0125.py @@ -0,0 +1,206 @@ +from struphy.io.options import EnvironmentOptions, BaseUnits, Time +from struphy.geometry import domains +from struphy.fields_background import equils +from struphy.topology import grids +from struphy.io.options import DerhamOptions +from struphy.io.options import FieldsBackground +from struphy.initial import perturbations +from struphy.kinetic_background import maxwellians +from struphy.pic.utilities import (LoadingParameters, + WeightsParameters, + BoundaryParameters, + BinningPlot, + KernelDensityPlot, + ) +from struphy import main + +# import model, set verbosity +from struphy.models.hybrid import LinearMHDDriftkineticCC + +# environment options +env = EnvironmentOptions( + out_folders='/tokp/work/bkna/tae', + sim_folder='tae_tp_sonic_gradp0_dt0125', + save_step=8, + max_runtime=2840, + restart=False) + +# units +base_units = BaseUnits(n=0.2) + +# time stepping +time_opts = Time(dt=0.125, Tend=1000., split_algo='Strang') + +# geometry +domain = domains.HollowTorus(a1=0.1, a2=1., R0=10., sfl=False, tor_period=6) + +# fluid equilibrium (can be used as part of initial conditions) +equil = equils.AdhocTorus( + a=1., + R0=10., + B0=3., + q_kind=0, + q0=1.71, + q1=1.87, + n1=0., + n2=0., + na=1., + p_kind=1., + p0=1., + p1=0., + p2=0., + beta=0.0018, +) + +# grid +grid = grids.TensorProductGrid( + Nel=(24,96,16), + mpi_dims_mask=(True, True, False) +) + +# derham options +derham_opts = DerhamOptions( + p=(3, 3, 3), + spl_kind=(False, True, True), + dirichlet_bc=((True, True), (False, False), (False, False)), + nq_pr=(6, 6, 6), + nquads=(6, 6, 6), + polar_ck=-1, +) + +# light-weight model instance +model = LinearMHDDriftkineticCC( + turn_off = (#'PushGuidingCenterBxEstar', + #'PushGuidingCenterParallel', + #'ShearAlfvenCurrentCoupling5D', + #'Magnetosonic', + #'CurrentCoupling5DDensity', + #'CurrentCoupling5DGradB', + #'CurrentCoupling5DCurlb', + ) +) + +# species parameters +model.mhd.set_phys_params() +model.energetic_ions.set_phys_params( + mass_number=2., +) + +loading_params = LoadingParameters( + ppc=100, + moments=(0., 0., 0.89740839, 0.89740839), + seed=1234, + ) +weights_params = WeightsParameters( + control_variate=True, +) +boundary_params = BoundaryParameters( + bc=("remove", "periodic", "periodic"), + bc_refill=("outer", "inner"), +) +model.energetic_ions.set_markers(loading_params=loading_params, + weights_params=weights_params, + boundary_params=boundary_params, + bufsize=1., + n_cols_diag=3, + n_cols_aux=5, + ) +model.energetic_ions.set_sorting_boxes(boxes_per_dim=None) +model.energetic_ions.set_save_data() + +# calculate constant +ep_scale = model.energetic_ions.var.species.mass_number / model.mhd.mass_number + +# params else +from struphy.pic.accumulation.filter import FilterParameters +from struphy.linear_algebra.solver import SolverParameters, DiscreteGradientSolverParameters +filter_params = FilterParameters(use_filter="three_point",) +solver_params = SolverParameters() + +# propagator options +model.propagators.push_bxe.options = model.propagators.push_bxe.Options( + algo = 'explicit', + b_tilde = model.em_fields.b_field,) +model.propagators.push_parallel.options = model.propagators.push_parallel.Options( + algo = 'explicit', + b_tilde = model.em_fields.b_field,) +model.propagators.shearalfen_cc5d.options = model.propagators.shearalfen_cc5d.Options( + filter_params=filter_params, + ep_scale=ep_scale, + energetic_ions = model.energetic_ions.var, + solver_params=solver_params, + nonlinear=False,) +model.propagators.magnetosonic.options = model.propagators.magnetosonic.Options( + b_field=model.em_fields.b_field,) +model.propagators.cc5d_density.options = model.propagators.cc5d_density.Options( + filter_params=filter_params, + ep_scale=ep_scale, + energetic_ions = model.energetic_ions.var, + b_tilde = model.em_fields.b_field,) +model.propagators.cc5d_gradb.options = model.propagators.cc5d_gradb.Options( + filter_params=filter_params, + ep_scale=ep_scale, + b_tilde = model.em_fields.b_field,) +model.propagators.cc5d_curlb.options = model.propagators.cc5d_curlb.Options( + filter_params=filter_params, + ep_scale=ep_scale, + b_tilde = model.em_fields.b_field,) + +# background, perturbations and initial conditions +model.mhd.velocity.add_background(FieldsBackground( + type="LogicalConst", + values=(0., 0., 0.) +)) +model.mhd.velocity.add_perturbation( + perturbations.TorusModesSin( + given_in_basis="2", + amps=(0.001,0.001), + ms=(10,11), + ns=(-1,-1), + pfuns=('exp','exp'), + pfun_params=((0.44444444,0.1),(0.44444444,0.1)), + comp=0) +) +model.mhd.velocity.add_perturbation( + perturbations.TorusModesCos( + given_in_basis="2", + amps=(0.000015915494309189535,0.000014468631190172303), + ms=(10,11), + ns=(-1,-1), + pfuns=('d_exp','d_exp'), + pfun_params=((0.44444444,0.1),(0.44444444,0.1)), + comp=1), +) +n_profile = perturbations.ITPA_density( + n0=0.00720652, + c=(0.47023, 0.20323, 0.13273, 0.521298) +) +maxwellian_1 = maxwellians.CanonicalMaxwellian2D( + vth=(0.89740839, None), + equil=equil, + volume_form=True, + n=(0., n_profile), +) +background = maxwellian_1 +model.energetic_ions.var.add_background( + background=background, + n_as_volume_form=False) + +# optional: exclude variables from saving +# model.energetic_ions.var.save_data = False + +if __name__ == "__main__": + # start run + verbose = True + + main.run(model, + params_path=__file__, + env=env, + base_units=base_units, + time_opts=time_opts, + domain=domain, + equil=equil, + grid=grid, + derham_opts=derham_opts, + verbose=verbose, + ) diff --git a/ITPA_benchmark/tae_tp_sonic_m2jboff.py b/ITPA_benchmark/tae_tp_sonic_m2jboff.py new file mode 100644 index 0000000..dec0b86 --- /dev/null +++ b/ITPA_benchmark/tae_tp_sonic_m2jboff.py @@ -0,0 +1,206 @@ +from struphy.io.options import EnvironmentOptions, BaseUnits, Time +from struphy.geometry import domains +from struphy.fields_background import equils +from struphy.topology import grids +from struphy.io.options import DerhamOptions +from struphy.io.options import FieldsBackground +from struphy.initial import perturbations +from struphy.kinetic_background import maxwellians +from struphy.pic.utilities import (LoadingParameters, + WeightsParameters, + BoundaryParameters, + BinningPlot, + KernelDensityPlot, + ) +from struphy import main + +# import model, set verbosity +# 64 +from struphy.models.hybrid import LinearMHDDriftkineticCC + +# environment options +env = EnvironmentOptions( + out_folders='/tokp/work/bkna/tae', + sim_folder='tae_tp_sonic_m2jboff', + save_step=8, + max_runtime=2840, + restart=True) + +# units +base_units = BaseUnits(n=0.2) + +# time stepping +time_opts = Time(dt=0.25, Tend=2000., split_algo='Strang') + +# geometry +domain = domains.HollowTorus(a1=0.1, a2=1., R0=10., sfl=False, tor_period=6) + +# fluid equilibrium (can be used as part of initial conditions) +equil = equils.AdhocTorus( + a=1., + R0=10., + B0=3., + q_kind=0, + q0=1.71, + q1=1.87, + n1=0., + n2=0., + na=1., + p_kind=1., + p0=1., + p1=0.95, + p2=0.05, + beta=0.0018, +) + +# grid +grid = grids.TensorProductGrid( + Nel=(24,96,16), + mpi_dims_mask=(True, True, False) +) + +# derham options +derham_opts = DerhamOptions( + p=(3, 3, 3), + spl_kind=(False, True, True), + dirichlet_bc=((True, True), (False, False), (False, False)), + nq_pr=(6, 6, 6), + nquads=(6, 6, 6), + polar_ck=-1, +) + +# light-weight model instance +model = LinearMHDDriftkineticCC( + turn_off = (#'PushGuidingCenterBxEstar', + #'PushGuidingCenterParallel', + #'ShearAlfvenCurrentCoupling5D', + #'Magnetosonic', + #'CurrentCoupling5DDensity', + #'CurrentCoupling5DGradB', + #'CurrentCoupling5DCurlb', + ) +) + +# species parameters +model.mhd.set_phys_params() +model.energetic_ions.set_phys_params( + mass_number=2., +) + +loading_params = LoadingParameters( + ppc=100, + moments=(0., 0., 0.89740839, 0.89740839), + seed=1234, + ) +weights_params = WeightsParameters( + control_variate=True, +) +boundary_params = BoundaryParameters( + bc=("remove", "periodic", "periodic"), + bc_refill=("outer", "inner"), +) +model.energetic_ions.set_markers(loading_params=loading_params, + weights_params=weights_params, + boundary_params=boundary_params, + bufsize=1., + n_cols_diag=3, + n_cols_aux=5, + ) +model.energetic_ions.set_sorting_boxes(boxes_per_dim=None) +model.energetic_ions.set_save_data() + +# calculate constant +ep_scale = model.energetic_ions.var.species.mass_number / model.mhd.mass_number + +# params else +from struphy.pic.accumulation.filter import FilterParameters +from struphy.linear_algebra.solver import SolverParameters, DiscreteGradientSolverParameters +filter_params = FilterParameters(use_filter="three_point",) +solver_params = SolverParameters() + +# propagator options +model.propagators.push_bxe.options = model.propagators.push_bxe.Options( + algo = 'explicit', + b_tilde = model.em_fields.b_field,) +model.propagators.push_parallel.options = model.propagators.push_parallel.Options( + algo = 'explicit', + b_tilde = model.em_fields.b_field,) +model.propagators.shearalfen_cc5d.options = model.propagators.shearalfen_cc5d.Options( + filter_params=filter_params, + ep_scale=ep_scale, + energetic_ions = model.energetic_ions.var, + solver_params=solver_params, + nonlinear=False,) +model.propagators.magnetosonic.options = model.propagators.magnetosonic.Options() +model.propagators.cc5d_density.options = model.propagators.cc5d_density.Options( + filter_params=filter_params, + ep_scale=ep_scale, + energetic_ions = model.energetic_ions.var, + b_tilde = model.em_fields.b_field,) +model.propagators.cc5d_gradb.options = model.propagators.cc5d_gradb.Options( + filter_params=filter_params, + ep_scale=ep_scale, + b_tilde = model.em_fields.b_field,) +model.propagators.cc5d_curlb.options = model.propagators.cc5d_curlb.Options( + filter_params=filter_params, + ep_scale=ep_scale, + b_tilde = model.em_fields.b_field,) + +# background, perturbations and initial conditions +model.mhd.velocity.add_background(FieldsBackground( + type="LogicalConst", + values=(0., 0., 0.) +)) +model.mhd.velocity.add_perturbation( + perturbations.TorusModesSin( + given_in_basis="2", + amps=(0.001,0.001), + ms=(10,11), + ns=(-1,-1), + pfuns=('exp','exp'), + pfun_params=((0.44444444,0.1),(0.44444444,0.1)), + comp=0) +) +model.mhd.velocity.add_perturbation( + perturbations.TorusModesCos( + given_in_basis="2", + amps=(0.000015915494309189535,0.000014468631190172303), + ms=(10,11), + ns=(-1,-1), + pfuns=('d_exp','d_exp'), + pfun_params=((0.44444444,0.1),(0.44444444,0.1)), + comp=1), +) +n_profile = perturbations.ITPA_density( + n0=0.00720652, + c=(0.47023, 0.20323, 0.13273, 0.521298) +) +maxwellian_1 = maxwellians.CanonicalMaxwellian2D( + vth=(0.89740839, None), + equil=equil, + volume_form=True, + n=(0., n_profile), +) +background = maxwellian_1 +model.energetic_ions.var.add_background( + background=background, + n_as_volume_form=False) + +# optional: exclude variables from saving +# model.energetic_ions.var.save_data = False + +if __name__ == "__main__": + # start run + verbose = True + + main.run(model, + params_path=__file__, + env=env, + base_units=base_units, + time_opts=time_opts, + domain=domain, + equil=equil, + grid=grid, + derham_opts=derham_opts, + verbose=verbose, + ) diff --git a/ITPA_benchmark/tae_tp_sonic_m2jboff_dt0125.py b/ITPA_benchmark/tae_tp_sonic_m2jboff_dt0125.py new file mode 100644 index 0000000..e0a6cfb --- /dev/null +++ b/ITPA_benchmark/tae_tp_sonic_m2jboff_dt0125.py @@ -0,0 +1,206 @@ +from struphy.io.options import EnvironmentOptions, BaseUnits, Time +from struphy.geometry import domains +from struphy.fields_background import equils +from struphy.topology import grids +from struphy.io.options import DerhamOptions +from struphy.io.options import FieldsBackground +from struphy.initial import perturbations +from struphy.kinetic_background import maxwellians +from struphy.pic.utilities import (LoadingParameters, + WeightsParameters, + BoundaryParameters, + BinningPlot, + KernelDensityPlot, + ) +from struphy import main + +# import model, set verbosity +# 64 +from struphy.models.hybrid import LinearMHDDriftkineticCC + +# environment options +env = EnvironmentOptions( + out_folders='/tokp/work/bkna/tae', + sim_folder='tae_tp_sonic_m2jboff_dt0125', + save_step=8, + max_runtime=2840, + restart=False) + +# units +base_units = BaseUnits(n=0.2) + +# time stepping +time_opts = Time(dt=0.125, Tend=2000., split_algo='Strang') + +# geometry +domain = domains.HollowTorus(a1=0.1, a2=1., R0=10., sfl=False, tor_period=6) + +# fluid equilibrium (can be used as part of initial conditions) +equil = equils.AdhocTorus( + a=1., + R0=10., + B0=3., + q_kind=0, + q0=1.71, + q1=1.87, + n1=0., + n2=0., + na=1., + p_kind=1., + p0=1., + p1=0.95, + p2=0.05, + beta=0.0018, +) + +# grid +grid = grids.TensorProductGrid( + Nel=(24,96,16), + mpi_dims_mask=(True, True, False) +) + +# derham options +derham_opts = DerhamOptions( + p=(3, 3, 3), + spl_kind=(False, True, True), + dirichlet_bc=((True, True), (False, False), (False, False)), + nq_pr=(6, 6, 6), + nquads=(6, 6, 6), + polar_ck=-1, +) + +# light-weight model instance +model = LinearMHDDriftkineticCC( + turn_off = (#'PushGuidingCenterBxEstar', + #'PushGuidingCenterParallel', + #'ShearAlfvenCurrentCoupling5D', + #'Magnetosonic', + #'CurrentCoupling5DDensity', + #'CurrentCoupling5DGradB', + #'CurrentCoupling5DCurlb', + ) +) + +# species parameters +model.mhd.set_phys_params() +model.energetic_ions.set_phys_params( + mass_number=2., +) + +loading_params = LoadingParameters( + ppc=100, + moments=(0., 0., 0.89740839, 0.89740839), + seed=1234, + ) +weights_params = WeightsParameters( + control_variate=True, +) +boundary_params = BoundaryParameters( + bc=("remove", "periodic", "periodic"), + bc_refill=("outer", "inner"), +) +model.energetic_ions.set_markers(loading_params=loading_params, + weights_params=weights_params, + boundary_params=boundary_params, + bufsize=1., + n_cols_diag=3, + n_cols_aux=5, + ) +model.energetic_ions.set_sorting_boxes(boxes_per_dim=None) +model.energetic_ions.set_save_data() + +# calculate constant +ep_scale = model.energetic_ions.var.species.mass_number / model.mhd.mass_number + +# params else +from struphy.pic.accumulation.filter import FilterParameters +from struphy.linear_algebra.solver import SolverParameters, DiscreteGradientSolverParameters +filter_params = FilterParameters(use_filter="three_point",) +solver_params = SolverParameters() + +# propagator options +model.propagators.push_bxe.options = model.propagators.push_bxe.Options( + algo = 'explicit', + b_tilde = model.em_fields.b_field,) +model.propagators.push_parallel.options = model.propagators.push_parallel.Options( + algo = 'explicit', + b_tilde = model.em_fields.b_field,) +model.propagators.shearalfen_cc5d.options = model.propagators.shearalfen_cc5d.Options( + filter_params=filter_params, + ep_scale=ep_scale, + energetic_ions = model.energetic_ions.var, + solver_params=solver_params, + nonlinear=False,) +model.propagators.magnetosonic.options = model.propagators.magnetosonic.Options() +model.propagators.cc5d_density.options = model.propagators.cc5d_density.Options( + filter_params=filter_params, + ep_scale=ep_scale, + energetic_ions = model.energetic_ions.var, + b_tilde = model.em_fields.b_field,) +model.propagators.cc5d_gradb.options = model.propagators.cc5d_gradb.Options( + filter_params=filter_params, + ep_scale=ep_scale, + b_tilde = model.em_fields.b_field,) +model.propagators.cc5d_curlb.options = model.propagators.cc5d_curlb.Options( + filter_params=filter_params, + ep_scale=ep_scale, + b_tilde = model.em_fields.b_field,) + +# background, perturbations and initial conditions +model.mhd.velocity.add_background(FieldsBackground( + type="LogicalConst", + values=(0., 0., 0.) +)) +model.mhd.velocity.add_perturbation( + perturbations.TorusModesSin( + given_in_basis="2", + amps=(0.001,0.001), + ms=(10,11), + ns=(-1,-1), + pfuns=('exp','exp'), + pfun_params=((0.44444444,0.1),(0.44444444,0.1)), + comp=0) +) +model.mhd.velocity.add_perturbation( + perturbations.TorusModesCos( + given_in_basis="2", + amps=(0.000015915494309189535,0.000014468631190172303), + ms=(10,11), + ns=(-1,-1), + pfuns=('d_exp','d_exp'), + pfun_params=((0.44444444,0.1),(0.44444444,0.1)), + comp=1), +) +n_profile = perturbations.ITPA_density( + n0=0.00720652, + c=(0.47023, 0.20323, 0.13273, 0.521298) +) +maxwellian_1 = maxwellians.CanonicalMaxwellian2D( + vth=(0.89740839, None), + equil=equil, + volume_form=True, + n=(0., n_profile), +) +background = maxwellian_1 +model.energetic_ions.var.add_background( + background=background, + n_as_volume_form=False) + +# optional: exclude variables from saving +# model.energetic_ions.var.save_data = False + +if __name__ == "__main__": + # start run + verbose = True + + main.run(model, + params_path=__file__, + env=env, + base_units=base_units, + time_opts=time_opts, + domain=domain, + equil=equil, + grid=grid, + derham_opts=derham_opts, + verbose=verbose, + )