diff --git a/TestAllMHD.sh b/TestAllMHD.sh index 68d8fa91ebbcb06324e9be76214c754924b33674..ce3b543ec96c1a916229fe1c71689c3a571d83bf 100755 --- a/TestAllMHD.sh +++ b/TestAllMHD.sh @@ -12,7 +12,8 @@ SCHEME_NAME=("Vector Potential" "Direct Induction (Orestis)" "Direct Induction ( # DO ~ Direct Induction Orestis # DF ~ Direct Induction Fede -BASE_CONF="--with-kernel=wendland-C4 --disable-hand-vec" +#BASE_CONF="--with-kernel=wendland-C4 --disable-hand-vec" +BASE_CONF="" case $1 in all) diff --git a/configure.ac b/configure.ac index b04ff0d46e15e9eca21fcbe50bb270a1efc7c1fc..c6ec362cdabb1b28608c204e8044d07abeebe0c6 100644 --- a/configure.ac +++ b/configure.ac @@ -146,7 +146,7 @@ AC_ARG_ENABLE([mpi], [enable_mpi="$enableval"], [enable_mpi="yes"] ) -# Use extra flags set by AC_PROG_CC as part of $CC. Currently undocumented as ac_cv_prog_cc_stdc, +# Use extra flags set by AC_PROG_CC as part of $CC. Currently undocumented as ac_cv_prog_cc_stdc, # so could change. good_mpi="yes" if test "$enable_mpi" = "yes"; then @@ -814,6 +814,21 @@ if test "x$with_metis" != "xno"; then fi AC_CHECK_LIB([metis],[METIS_PartGraphKway], [have_metis="yes"], [have_metis="no"], $METIS_LIBS) + + # Recent METIS releases have an external GKlib, test for that before + # giving up. Assume sane and in the same directories. + if test "x$have_metis" = "xno"; then + if test "x$with_metis" != "xyes" -a "x$with_metis" != "x"; then + METIS_LIBS="-L$with_metis/lib -lmetis -lGKlib" + METIS_INCS="-I$with_metis/include" + else + METIS_LIBS="-lmetis -lGKlib" + METIS_INCS="" + fi + AC_CHECK_LIB([metis],[METIS_PartGraphRecursive], [have_metis="yes"], + [have_metis="no"], $METIS_LIBS) + fi + if test "$have_metis" = "yes"; then AC_DEFINE([HAVE_METIS],1,[The METIS library is present.]) else @@ -847,10 +862,9 @@ if test "x$with_parmetis" != "xno"; then fi AC_CHECK_LIB([parmetis],[ParMETIS_V3_RefineKway], [have_parmetis="yes"], [have_parmetis="no"], $PARMETIS_LIBS) - if test "$have_parmetis" = "no"; then # A build may use an external METIS library, check for that. - + if test "$have_parmetis" = "no"; then if test "x$with_parmetis" != "xyes" -a "x$with_parmetis" != "x"; then PARMETIS_LIBS="-L$with_parmetis/lib -lparmetis -lmetis" PARMETIS_INCS="-I$with_parmetis/include" @@ -861,6 +875,20 @@ if test "x$with_parmetis" != "xno"; then # Note use different function to avoid caching of first check. AC_CHECK_LIB([parmetis],[ParMETIS_V3_PartKway], [have_parmetis="yes"], [have_parmetis="no"], [$METIS_LIBS $PARMETIS_LIBS]) + fi + +# A build may use an external GKlib in later releases... + if test "$have_parmetis" = "no"; then + if test "x$with_parmetis" != "xyes" -a "x$with_parmetis" != "x"; then + PARMETIS_LIBS="-L$with_parmetis/lib -lparmetis -lGKlib" + PARMETIS_INCS="-I$with_parmetis/include" + else + PARMETIS_LIBS="-lparmetis -lGKlib" + PARMETIS_INCS="" + fi + # Note use different function to avoid caching of first check. + AC_CHECK_LIB([parmetis],[ParMETIS_V3_PartGeom], [have_parmetis="yes"], + [have_parmetis="no"], [$METIS_LIBS $PARMETIS_LIBS]) fi if test "$have_parmetis" = "yes"; then @@ -1500,7 +1528,7 @@ if test "$enable_lightcone" = "yes"; then # Also need to make sure we have GSL if we're making lightcones if test "$have_gsl" != "yes"; then AC_MSG_ERROR([Lightcone output requires GSL. Please configure with --with-gsl.]) - fi + fi else have_chealpix="no" fi @@ -1568,10 +1596,10 @@ AC_SUBST([NUMA_LIBS]) # Check for Sundials (required for the SPHM1RT library). -# There is a problems with the headers of this library -# as they do not pass the strict prototypes check when -# installed outside of the system directories. So we -# need to do this check in two phases. +# There is a problems with the headers of this library +# as they do not pass the strict prototypes check when +# installed outside of the system directories. So we +# need to do this check in two phases. have_sundials="no" SUNDIALS_LIBS="" SUNDIALS_INCS="" @@ -1600,22 +1628,22 @@ if test "x$with_sundials" != "xno"; then AC_DEFINE([HAVE_SUNDIALS],1,[The SUNDIALS library is present.]) else if test "x$with_sundials" != "xyes" -a "x$with_sundials" != "x"; then - # It might be that the libraries are in - # /lib64 rather than /lib + # It might be that the libraries are in + # /lib64 rather than /lib SUNDIALS_LIBS="-L$with_sundials/lib64 -lsundials_cvode -lsundials_nvecserial -lsundials_sunlinsoldense -lsundials_sunmatrixdense" # unset cached result of previous AC_CHECK_LIB unset ac_cv_lib_sundials_cvode_CVode - AC_CHECK_LIB([sundials_cvode], [CVode], [have_sundials="yes"], [have_sundials="no"], $SUNDIALS_LIBS) - + AC_CHECK_LIB([sundials_cvode], [CVode], [have_sundials="yes"], [have_sundials="no"], $SUNDIALS_LIBS) + if test "$have_sundials" == "yes"; then AC_DEFINE([HAVE_SUNDIALS],1,[The SUNDIALS library is present.]) - else + else AC_MSG_ERROR("Failed to find a SUNDIALS library") - fi + fi - else + else AC_MSG_ERROR("Failed to find a SUNDIALS library") fi fi @@ -1784,9 +1812,9 @@ fi AC_SUBST([NUMA_INCS]) -# Second part of the Sundials library checks. -# We now decide if we need to use -isystem to -# get around the strict-prototypes problem. Assumes +# Second part of the Sundials library checks. +# We now decide if we need to use -isystem to +# get around the strict-prototypes problem. Assumes # isystem is available when strict-prototypes is. if test "x$with_sundials" != "xno"; then if test "x$with_sundials" != "xyes" -a "x$with_sundials" != "x"; then @@ -1799,7 +1827,7 @@ if test "x$with_sundials" != "xno"; then ;; esac fi -fi +fi AC_SUBST([SUNDIALS_INCS]) @@ -1858,7 +1886,7 @@ case "$with_subgrid" in with_subgrid_sink=none with_subgrid_extra_io=none enable_fof=no - ;; + ;; QLA) with_subgrid_cooling=QLA with_subgrid_chemistry=QLA @@ -1993,7 +2021,7 @@ fi # Hydro scheme. AC_ARG_WITH([hydro], [AS_HELP_STRING([--with-hydro=<scheme>], - [Hydro dynamics to use @<:@gadget2, minimal, pressure-entropy, pressure-energy, pressure-energy-monaghan, phantom, gizmo-mfv, gizmo-mfm, shadowfax, planetary, sphenix, gasoline, anarchy-pu default: sphenix@:>@] + [Hydro dynamics to use @<:@gadget2, minimal, pressure-entropy, pressure-energy, pressure-energy-monaghan, phantom, gizmo-mfv, gizmo-mfm, shadowfax, planetary, sphenix, gasoline, anarchy-pu, magma2 default: sphenix@:>@] )], [with_hydro="$withval"], [with_hydro="sphenix"] @@ -2039,6 +2067,9 @@ case "$with_hydro" in sphenix) AC_DEFINE([SPHENIX_SPH], [1], [SPHENIX SPH]) ;; + magma2) + AC_DEFINE([MAGMA_SPH], [1], [MAGMA SPH]) + ;; gasoline) AC_DEFINE([GASOLINE_SPH], [1], [Gasoline SPH]) ;; @@ -2322,7 +2353,7 @@ case "$with_chemistry" in AGORA) AC_DEFINE([CHEMISTRY_AGORA], [1], [Chemistry taken from the AGORA model]) with_chemistry_name="AGORA" - ;; + ;; QLA) AC_DEFINE([CHEMISTRY_QLA], [1], [Chemistry taken from the Quick-Lyman-alpha model]) with_chemistry_name="QLA" @@ -2504,7 +2535,7 @@ case "$with_stars" in ;; GEAR) AC_DEFINE([STARS_GEAR], [1], [GEAR stellar model]) - ;; + ;; basic) AC_DEFINE([STARS_BASIC], [1], [Basic stellar model]) ;; @@ -2554,7 +2585,7 @@ case "$with_feedback" in AGORA) AC_DEFINE([FEEDBACK_AGORA], [1], [AGORA stellar feedback and evolution model]) with_feedback_name="AGORA" - ;; + ;; none) AC_DEFINE([FEEDBACK_NONE], [1], [No feedback]) ;; @@ -2720,9 +2751,9 @@ case "$with_potential" in esac # Entropy floor -AC_ARG_WITH([entropy-floor], +AC_ARG_WITH([entropy-floor], [AS_HELP_STRING([--with-entropy-floor=<floor>], - [entropy floor @<:@none, QLA, EAGLE, default: none@:>@] + [entropy floor @<:@none, QLA, EAGLE, default: none@:>@] )], [with_entropy_floor="$withval"], [with_entropy_floor="none"] @@ -2748,10 +2779,10 @@ case "$with_entropy_floor" in *) AC_MSG_ERROR([Unknown entropy floor model]) ;; -esac +esac # Pressure floor -AC_ARG_WITH([pressure-floor], +AC_ARG_WITH([pressure-floor], [AS_HELP_STRING([--with-pressure-floor=<floor>], [pressure floor @<:@none, GEAR, default: none@:>@ The hydro model needs to be compatible.] @@ -2777,12 +2808,12 @@ case "$with_pressure_floor" in *) AC_MSG_ERROR([Unknown pressure floor model]) ;; -esac +esac # Star formation -AC_ARG_WITH([star-formation], +AC_ARG_WITH([star-formation], [AS_HELP_STRING([--with-star-formation=<sfm>], - [star formation @<:@none, QLA, EAGLE, GEAR, default: none@:>@] + [star formation @<:@none, QLA, EAGLE, GEAR, default: none@:>@] )], [with_star_formation="$withval"], [with_star_formation="none"] @@ -2811,7 +2842,7 @@ case "$with_star_formation" in *) AC_MSG_ERROR([Unknown star formation model]) ;; -esac +esac AC_ARG_WITH([gadget2-physical-constants], [AS_HELP_STRING([--with-gadget2-physical-constants], @@ -2837,7 +2868,7 @@ AC_DEFINE_UNQUOTED([SELF_GRAVITY_MULTIPOLE_ORDER], [$with_multipole_order], [Mul # Radiative transfer scheme AC_ARG_WITH([rt], [AS_HELP_STRING([--with-rt=<scheme>], - [Radiative transfer scheme to use @<:@none, GEAR_*, SPHM1RT_*, debug default: none@:>@. + [Radiative transfer scheme to use @<:@none, GEAR_*, SPHM1RT_*, debug default: none@:>@. For GEAR and SPHM1RT, the number of photon groups (e.g. GEAR_4) needs to be provided.] )], [with_rt="$withval"], @@ -2847,7 +2878,7 @@ AC_ARG_WITH([rt], # For GEAR-RT scheme: Select a RT Riemann solver AC_ARG_WITH([rt-riemann-solver], [AS_HELP_STRING([--with-rt-riemann-solver=<scheme>], - [Riemann solver for the moments of the ratiadiative transfer equation with the M1 closure to use @<:@none, HLL, GLF, default: none@:>@. + [Riemann solver for the moments of the ratiadiative transfer equation with the M1 closure to use @<:@none, HLL, GLF, default: none@:>@. For the GEAR RT scheme, you need to select one Riemann solver.] )], [with_rt_riemann_solver="$withval"], @@ -2926,8 +2957,8 @@ case "$with_rt" in fi AC_MSG_CHECKING([for Sundials libraries]) AC_MSG_RESULT($have_sundials) - if test "$have_sundials" != "yes"; then - AC_MSG_ERROR([The Sundials library is not present. Sundials is required for the SPHM1RT module.]) + if test "$have_sundials" != "yes"; then + AC_MSG_ERROR([The Sundials library is not present. Sundials is required for the SPHM1RT module.]) fi ;; *) @@ -3059,8 +3090,8 @@ AC_MSG_RESULT([ HDF5 enabled : $with_hdf5 - parallel : $have_parallel_hdf5 METIS/ParMETIS : $have_metis / $have_parmetis - FFTW3 enabled : $have_fftw - - threaded/openmp : $have_threaded_fftw / $have_openmp_fftw + FFTW3 enabled : $have_fftw + - threaded/openmp : $have_threaded_fftw / $have_openmp_fftw - MPI : $have_mpi_fftw - ARM : $have_arm_fftw GSL enabled : $have_gsl @@ -3103,7 +3134,7 @@ AC_MSG_RESULT([ Black holes model : $with_black_holes Radiative transfer : $with_rt Extra i/o : $with_extra_io - + Atomic operations in tasks : $enable_atomics_within_tasks Individual timers : $enable_timers Task debugging : $enable_task_debugging diff --git a/doc/RTD/source/SubgridModels/AGNSpinJets/params.rst b/doc/RTD/source/SubgridModels/AGNSpinJets/params.rst index eef4286a1a52d0914cc5432873787f54ca8039ad..a4ac82e8de3e71c8dbc9d98c247a5ead87cd3cd3 100644 --- a/doc/RTD/source/SubgridModels/AGNSpinJets/params.rst +++ b/doc/RTD/source/SubgridModels/AGNSpinJets/params.rst @@ -50,7 +50,7 @@ Below we give an example of parameter choices applicable for e.g. a 50 Mpc box. include_jets: 1 # Global switch whether to include jet feedback [1] or not [0]. turn_off_radiative_feedback: 0 # Global switch whether to turn off radiative (thermal) feedback [1] or not [0]. This should only be used if 'include_jets' is set to 1, since we want feedback in some form or another. alpha_acc: 0.2 # Viscous alpha of the subgrid accretion disks. Likely to be within the 0.1-0.3 range. The main effect is that it sets the transition accretion rate between thin and thick disk, as dot(m) = 0.2 * alpha^2. - mdot_crit_ADAF: 0.01 # The transition normalized accretion rate (Eddington ratio) at which the disc goes from thick (low accretion rates) to thin (high accretion rates). The feedback also changes from kinetic jets to thermal isotropic, respectively. + mdot_crit_ADAF: 0.03 # The transition normalized accretion rate (Eddington ratio) at which the disc goes from thick (low accretion rates) to thin (high accretion rates). The feedback also changes from kinetic jets to thermal isotropic, respectively. seed_spin: 0.01 # The (randomly-directed) black hole spin assigned to BHs when they are seeded. Should be strictly between 0 and 1. AGN_jet_velocity_model: Constant # How AGN jet velocities are calculated. If 'Constant', a single value is used. If 'BlackHoleMass', then an empirical relation between halo mass and black hole mass is used to calculate jet velocities. 'HaloMass' is currently not supported. v_jet_km_p_s: 3160. # Jet velocity to use if 'AGN_jet_velocity_model' is 'Constant'. Units are km/s. @@ -69,7 +69,7 @@ Below we give an example of parameter choices applicable for e.g. a 50 Mpc box. fix_jet_efficiency: 0 # Global switch whether to fix jet efficiency to a particular value [1], or use a spin-dependant formula [0]. If used, jets will be launched exclusively along the z axis. Should be set to 1 only for tests. jet_efficiency: 0.1 # The constant jet efficiency used if 'fix_jet_efficiency' is set to 1. fix_jet_direction: 0 # Global switch whether to fix the jet direction to be along the z-axis, instead of along the spin vector. - accretion_efficiency: 1. # The accretion efficiency (suppression factor of the accretion rate) to use in the thick disc (ADAF), to represent the effects of subgrid winds that take away most of the mass flowing through the accretion disc. + accretion_efficiency: 0.1 # The accretion efficiency (suppression factor of the accretion rate) to use in the thick disc (ADAF), to represent the effects of subgrid winds that take away most of the mass flowing through the accretion disc. fix_radiative_efficiency: 0 # Global switch whether to fix the radiative efficiency to a particular value [1], or use a spin-dependant formula [0]. radiative_efficiency: 0.1 # The constant jet efficiency used if 'fix_radiative_efficiency' is set to 1. Otherwise, this value is used to define the Eddington accretion rate. TD_region: B # How to treat the subgrid accretion disk if it is thin, according to the Shakura & Sunyaev (1973) model. If set to B, region b will be used. If set to C, region c will be used. @@ -78,7 +78,7 @@ Below we give an example of parameter choices applicable for e.g. a 50 Mpc box. turn_off_secondary_feedback: 1 # If set to 1, there will be only radiative (thermal) feedback in the thin disk mode, and only jets in the thick disk mode. jet_h_r_slope: 1. # The slope of the dependence of jet efficiency on aspect ratio of the subgrid accretion disk, H/R. Default value is 1, and another reasonable value is 0 (same jet efficiency for all disks). Reality could be anything in between. This parameter is only used if turn_off_secondary_feedback is set to 0. delta_ADAF: 0.2 # Electron heating parameter, which controls the strength of radiative feedback in thick disks. Should be between 0.1 and 0.5. This parameter is only used if turn_off_secondary_feedback is set to 0. - include_slim_disk: 1 # Global switch whether to include super-Eddington accretion, modeled as the slim disk. If set to 0, disks will be considered thin even at very large accretion rates. + include_slim_disk: 0 # Global switch whether to include super-Eddington accretion, modeled as the slim disk. If set to 0, disks will be considered thin even at very large accretion rates. Most of these parameters should work well generally, and should not be changed except for tests. We will discuss only some of the more important ones. You can choose whether to have only the thick and thin disk (low and high BH accretion rates, respectively), or you can also include the slim disk at super-Eddington rates with ``include_slim_disk``. You can control what type of feedback you (do not) want with ``include_jets`` and ``turn_off_radiative_feedback``. If you choose to turn off jets, everything will be modeled as a thin disk (regardless of accretion rate), since jets go hand-in-hand with the thick and the slim disk. Similarly, if you turn off radiation, everything will be treated as a thick disk. diff --git a/doc/RTD/source/SubgridModels/GEAR/index.rst b/doc/RTD/source/SubgridModels/GEAR/index.rst index c1ac89f70be3bc2a75a99e2fb02ef515d7d67d5d..a3ada8807f5b85bd86783fd4034040d0ae98e6cd 100644 --- a/doc/RTD/source/SubgridModels/GEAR/index.rst +++ b/doc/RTD/source/SubgridModels/GEAR/index.rst @@ -91,18 +91,30 @@ The self shielding method is defined by ``GrackleCooling:self_shielding_method`` .. code:: YAML GrackleCooling: - cloudy_table: CloudyData_UVB=HM2012.h5 # Name of the Cloudy Table (available on the grackle bitbucket repository) - with_UV_background: 1 # Enable or not the UV background - redshift: 0 # Redshift to use (-1 means time based redshift) - with_metal_cooling: 1 # Enable or not the metal cooling - provide_volumetric_heating_rates: 0 # (optional) User provide volumetric heating rates - provide_specific_heating_rates: 0 # (optional) User provide specific heating rates - max_steps: 10000 # (optional) Max number of step when computing the initial composition - convergence_limit: 1e-2 # (optional) Convergence threshold (relative) for initial composition - thermal_time_myr: 5 # (optional) Time (in Myr) for adiabatic cooling after a feedback event. - self_shielding_method: -1 # (optional) Grackle (1->3 for Grackle's ones, 0 for none and -1 for GEAR) - self_shielding_threshold_atom_per_cm3: 0.007 # Required only with GEAR's self shielding. Density threshold of the self shielding - + cloudy_table: CloudyData_UVB=HM2012.h5 # Name of the Cloudy Table (available on the grackle bitbucket repository) + with_UV_background: 1 # Enable or not the UV background + redshift: 0 # Redshift to use (-1 means time based redshift) + with_metal_cooling: 1 # Enable or not the metal cooling + provide_volumetric_heating_rates: 0 # (optional) User provide volumetric heating rates + provide_specific_heating_rates: 0 # (optional) User provide specific heating rates + max_steps: 10000 # (optional) Max number of step when computing the initial composition + convergence_limit: 1e-2 # (optional) Convergence threshold (relative) for initial composition + thermal_time_myr: 5 # (optional) Time (in Myr) for adiabatic cooling after a feedback event. + self_shielding_method: -1 # (optional) Grackle (1->3 for Grackle's ones, 0 for none and -1 for GEAR) + self_shielding_threshold_atom_per_cm3: 0.007 # Required only with GEAR's self shielding. Density threshold of the self shielding + HydrogenFractionByMass : 1. # Hydrogen fraction by mass (default is 0.76) + use_radiative_transfer : 1 # Arrays of ionization and heating rates are provided + RT_heating_rate_cgs : 0 # heating rate in units of / nHI_cgs + RT_HI_ionization_rate_cgs : 0 # HI ionization rate in cgs [1/s] + RT_HeI_ionization_rate_cgs : 0 # HeI ionization rate in cgs [1/s] + RT_HeII_ionization_rate_cgs: 0 # HeII ionization rate in cgs [1/s] + RT_H2_dissociation_rate_cgs: 0 # H2 dissociation rate in cgs [1/s] + volumetric_heating_rates_cgs: 0 # Volumetric heating rate in cgs [erg/s/cm3] + specific_heating_rates_cgs: 0 # Specific heating rate in cgs [erg/s/g] + H2_three_body_rate : 1 # Specific the H2 formation three body rate (0->5,see Grackle documentation) + H2_cie_cooling : 0 # Enable/disable H2 collision-induced emission cooling from Ripamonti & Abel (2004) + cmb_temperature_floor : 1 # Enable/disable an effective CMB temperature floor + .. note:: A simple example running SWIFT with Grackle can be find in ``examples/Cooling/CoolingBox``. A more advanced example combining heating and cooling (with heating and ionization sources) is given in ``examples/Cooling/CoolingHeatingBox``. diff --git a/examples/Cooling/CoolingBox/README b/examples/Cooling/CoolingBox/README new file mode 100644 index 0000000000000000000000000000000000000000..49ab9afad7e5b893018c1d77b4690cd75ffcb16b --- /dev/null +++ b/examples/Cooling/CoolingBox/README @@ -0,0 +1,14 @@ +Runs a uniform box exposed to constant cooling rates computed with +the Grackle library. + +To run with Grackle mode 0, configure swift with: + ./configure --with-cooling=grackle_0 --with-grackle=$GRACKLE_ROOT + +To run with Grackle mode 1, configure swift with: + ./configure --with-cooling=grackle_1 --with-grackle=$GRACKLE_ROOT + +To run with Grackle mode 2, configure swift with: + ./configure --with-cooling=grackle_2 --with-grackle=$GRACKLE_ROOT + +To run with Grackle mode 3, configure swift with: + ./configure --with-cooling=grackle_3 --with-grackle=$GRACKLE_ROOT diff --git a/examples/Cooling/CoolingBox/run.sh b/examples/Cooling/CoolingBox/run.sh index d83974dd34c71310d1eca7c7376de83fd20d3210..ed66cfe12d469e6f7b756b7df7bfaa3fbaf07932 100755 --- a/examples/Cooling/CoolingBox/run.sh +++ b/examples/Cooling/CoolingBox/run.sh @@ -17,7 +17,7 @@ fi if [ ! -e CloudyData_UVB=HM2012.h5 ] then echo "Fetching the Cloudy tables required by Grackle..." - ../getCoolingTable.sh + ../getGrackleCoolingTable.sh fi # Run SWIFT diff --git a/examples/HydroTests/NFW_Hydrostatic/NFW_Hydrostatic.yml b/examples/HydroTests/NFW_Hydrostatic/NFW_Hydrostatic.yml new file mode 100644 index 0000000000000000000000000000000000000000..5e31dd03d6bc21c8d8429ab375c83e0a2df38546 --- /dev/null +++ b/examples/HydroTests/NFW_Hydrostatic/NFW_Hydrostatic.yml @@ -0,0 +1,58 @@ +# Define the system of units to use internally. +InternalUnitSystem: + UnitMass_in_cgs: 1.98848e43 # 10^10 M_sun in grams + UnitLength_in_cgs: 3.08567758e21 # kpc in centimeters + UnitVelocity_in_cgs: 1e5 # km/s in centimeters per second + UnitCurrent_in_cgs: 1 # Amperes + UnitTemp_in_cgs: 1 # Kelvin + + +# Parameters governing the time integration +TimeIntegration: + time_begin: 0. # The starting time of the simulation (in internal units). + time_end: 1.0 # The end time of the simulation (in internal units). + dt_min: 1e-10 # The minimal time-step size of the simulation (in internal units). + dt_max: 0.1 # The maximal time-step size of the simulation (in internal units). + +# Parameters governing the snapshots +Snapshots: + basename: snapshot # Common part of the name of output files + time_first: 0. # Time of the first output (in internal units) + delta_time: 2e-2 # Time difference between consecutive outputs (in internal units) + +# Parameters governing the conserved quantities statistics +Statistics: + delta_time: 1e-3 # Time between statistics output + +# Parameters for the self-gravity scheme +Gravity: + eta: 0.05 # Constant dimensionless multiplier for time integration. + MAC: geometric + theta_cr: 0.7 + comoving_DM_softening: 0.01 # Comoving softening length (in internal units). + max_physical_DM_softening: 0.01 # Physical softening length (in internal units). + comoving_baryon_softening: 0.01 # Comoving softening length (in internal units). + max_physical_baryon_softening: 0.01 # Physical softening length (in internal units). + +# Parameters for the hydrodynamics scheme +SPH: + resolution_eta: 1.2348 # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel). + CFL_condition: 0.1 # Courant-Friedrich-Levy condition for time integration. + minimal_temperature: 10. # Kelvin + +# Parameters related to the initial conditions +InitialConditions: + file_name: ./nfw.hdf5 # The file to read + periodic: 0 # Non-periodic BCs + shift: [0,0,0] # Centre the box + +NFWPotential: + useabspos: 0 # 0 -> positions based on centre, 1 -> absolute positions + position: [0.,0.,0.] # Location of centre of isothermal potential with respect to centre of the box (if 0) otherwise absolute (if 1) (internal units) + M_200: 9.5 # M200 of the galaxy disk + h: 0.72 # reduced Hubble constant (value does not specify the used units!) + concentration: 17.0 # concentration of the Halo + diskfraction: 0.15 # Disk mass fraction (here this is the gas) + bulgefraction: 0.0 # Bulge mass fraction + timestep_mult: 0.01 # Dimensionless pre-factor for the time-step condition, basically determines the fraction of the orbital time we use to do the time integration + epsilon: 0.001 # Softening size (internal units) diff --git a/examples/HydroTests/NFW_Hydrostatic/README b/examples/HydroTests/NFW_Hydrostatic/README new file mode 100644 index 0000000000000000000000000000000000000000..17655f5c5c83a3bd0894caa79c6d11e36451d67e --- /dev/null +++ b/examples/HydroTests/NFW_Hydrostatic/README @@ -0,0 +1,13 @@ + +Evolve gas at the static equilibrium in an NFW potential for several dynamical times. +This test allows to compare the impact of different hydro solvers on the initial +hydrostatic equilibrium. + +To run SWIFT configured with the sphenix hydro solver: + ./configure --with-ext-potential=nfw --with-hydro=sphenix + +To run SWIFT configured with the gadget2 hydro solver: + ./configure --with-ext-potential=nfw --with-hydro=gadget2 --disable-hand-vec + +To run SWIFT configured with the gizmo (MFV) hydro solver: + ./configure --with-ext-potential=nfw --with-hydro=gizmo-mfv --with-riemann-solver=exact diff --git a/examples/HydroTests/NFW_Hydrostatic/getICs.sh b/examples/HydroTests/NFW_Hydrostatic/getICs.sh new file mode 100755 index 0000000000000000000000000000000000000000..4f132dc2c82b01e80a0e051137d4bf52c3b91fe5 --- /dev/null +++ b/examples/HydroTests/NFW_Hydrostatic/getICs.sh @@ -0,0 +1,3 @@ +#!/bin/bash +wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/NFW_Hydrostatic/nfw.hdf5 + diff --git a/examples/HydroTests/NFW_Hydrostatic/plotGasDensityProfile.py b/examples/HydroTests/NFW_Hydrostatic/plotGasDensityProfile.py new file mode 100755 index 0000000000000000000000000000000000000000..beb677f5127f98626b3dbd2e8f544bbb78c8de56 --- /dev/null +++ b/examples/HydroTests/NFW_Hydrostatic/plotGasDensityProfile.py @@ -0,0 +1,126 @@ +############################################################################### +# This file is part of SWIFT. +# Copyright (c) 2023 Yves Revaz (yves.revaz@epfl.ch) +# +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU Lesser General Public License as published +# by the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU Lesser General Public License +# along with this program. If not, see <http://www.gnu.org/licenses/>. +# +############################################################################## +import matplotlib + +matplotlib.use("Agg") +import matplotlib.pyplot as plt +import numpy as np +from glob import glob +import h5py + + +plt.style.use("../../../tools/stylesheets/mnras.mplstyle") + +MyrInSec = 31557600000000.0 +gcm3InAcc = 1 / 5.978637406556783e23 + + +def ComputeDensity(snap): + + # Read the initial state of the gas + f = h5py.File(snap, "r") + + # Read the units parameters from the snapshot + units = f["InternalCodeUnits"] + unit_mass = units.attrs["Unit mass in cgs (U_M)"] + unit_length = units.attrs["Unit length in cgs (U_L)"] + unit_time = units.attrs["Unit time in cgs (U_t)"] + unit_density = unit_mass / unit_length ** 3 + + # Header + header = f["Header"] + BoxSize = header.attrs["BoxSize"] + Time = header.attrs["Time"] + + # Read data + pos = f["/PartType0/Coordinates"][:] + rho = f["/PartType0/Densities"][:] + mass = f["/PartType0/Masses"][:] + ids = f["/PartType0/ParticleIDs"][:] + + # Center the model and compute particle radius + pos = pos - BoxSize / 2 + x = pos[:, 0] + y = pos[:, 1] + z = pos[:, 2] + r = np.sqrt(x * x + y * y + z * z) + n = len(r) + + # sort particles according to their distance + idx = np.argsort(r) + r = r[idx] + mass = mass[idx] + ids = ids[idx] + + nparts_per_bin = 100 + + # loop over bins containing nparts_per_bin particles + idx = np.arange(n) + + # bins radius and mass + rs_beg = np.array([]) + rs_end = np.array([]) + ms = np.array([]) + + i = 0 + + while i + nparts_per_bin < n: + + rs_beg = np.concatenate((rs_beg, [r[i]])) + rs_end = np.concatenate((rs_end, [r[i + nparts_per_bin]])) + m_this_bin = np.sum(mass[i : i + nparts_per_bin]) + ms = np.concatenate((ms, [np.sum(m_this_bin)])) + + # shift + i = i + nparts_per_bin + + # compute density + vol = 4 / 3 * np.pi * (rs_end ** 3 - rs_beg ** 3) + rho = ms / vol + + # compute radius, we use the mean + rs = 0.5 * (rs_beg + rs_end) + + # convert rho to acc + rho = rho * unit_density / gcm3InAcc + + # convert time to Myr + Time = Time * unit_time / MyrInSec + + return rs, rho, Time + + +# Do the plot + +plt.figure() + +rs, rho, Time = ComputeDensity("snapshot_0000.hdf5") +plt.plot(rs, rho, c="b", label=r"$t=%5.1f\,\rm{[Myr]}$" % Time, lw=1) + +rs, rho, Time = ComputeDensity("snapshot_0050.hdf5") +plt.plot(rs, rho, c="r", label=r"$t=%5.1f\,\rm{[Myr]}$" % Time, lw=1) + +plt.loglog() + +plt.legend() +plt.xlabel("${\\rm{Radius~[kpc]}}$", labelpad=0) +plt.ylabel("${\\rm{Density~[atom/cm^3]}}$", labelpad=0) + + +plt.savefig("GasDensityProfile.png", dpi=200) diff --git a/examples/HydroTests/NFW_Hydrostatic/run.sh b/examples/HydroTests/NFW_Hydrostatic/run.sh new file mode 100755 index 0000000000000000000000000000000000000000..713957f30caaa90782902fcccab4d44b2ba366e4 --- /dev/null +++ b/examples/HydroTests/NFW_Hydrostatic/run.sh @@ -0,0 +1,15 @@ +#!/bin/bash + +# Generate the initial conditions if they are not present. +if [ ! -e nfw.hdf5 ] +then + echo "Fetching initial conditions file for the example..." + ./getICs.sh +fi + + +# Run SWIFT +../../../swift --hydro --external-gravity --self-gravity --threads=14 NFW_Hydrostatic.yml + +# Compute gas density profile +python3 plotGasDensityProfile.py diff --git a/examples/HydroTests/SedovBlast_3D/plotSolution.py b/examples/HydroTests/SedovBlast_3D/plotSolution.py index e6f5381df0397fc78e277faf5a676035e42ec92c..2dd850cb415d48fc80abd0dcd50f6f88bc447e84 100644 --- a/examples/HydroTests/SedovBlast_3D/plotSolution.py +++ b/examples/HydroTests/SedovBlast_3D/plotSolution.py @@ -391,3 +391,15 @@ yticks([]) tight_layout() savefig("Sedov.png") + +data = np.loadtxt("statistics.txt") +t = data[:,1] +E_kin = data[:,13] +E_int = data[:,14] +E_tot = E_kin + E_int +print(E_tot) +figure() +plot(t, E_tot) +plot(t, E_int, '--') +plot(t, E_kin, ':') +savefig("Energy.png") diff --git a/examples/IdealisedCluster/IdealisedCluster_M13/idealised_cluster_M13.yml b/examples/IdealisedCluster/IdealisedCluster_M13/idealised_cluster_M13.yml index 1e4987e51809807150234f11fdd58ecb1b35ea6b..e94b56a79e262b0836224f3ec7251185d09b60e7 100644 --- a/examples/IdealisedCluster/IdealisedCluster_M13/idealised_cluster_M13.yml +++ b/examples/IdealisedCluster/IdealisedCluster_M13/idealised_cluster_M13.yml @@ -284,6 +284,7 @@ SPINJETAGN: include_jets: 1 # Global switch whether to include jet feedback [1] or not [0]. Requires the 'include_spin' parameter to also be set to 1. turn_off_radiative_feedback: 1 # Global switch whether to turn off radiative feedback [1] or not [0]. This should only be used if 'include_jets' is set to 1. alpha_acc: 0.2 # Viscous alpha of the subgrid accretion disks. Should be within the 0.1-0.3 range. + mdot_crit_ADAF: 0.03 # The transition normalized accretion rate (Eddington ratio) at which the disc goes from thick (low accretion rates) to thin (high accretion rates). The feedback also changes from kinetic jets to thermal isotropic, respectively. delta_ADAF: 0.2 # Electron heating parameter, which controls the strength of radiative feedback in thick disks. Should be between 0.1 and 0.5. TD_region: B # How to treat the subgrid accretion disk if it is thin, according to the Shakura & Sunyaev (1973) model. If set to B, region b will be used. If set to C, region c will be used. N_jet: 2 # Target number of particles to kick as part of a single jet feedback event. Should be a multiple of 2 to ensure approximate momentum conservation. @@ -298,4 +299,6 @@ SPINJETAGN: turn_off_secondary_feedback: 1 # Whether to turn off radiative feedback at very low accretion rates and jets at moderate accretion rates. fix_jet_direction: 1 # Whether to fix the jet direction to be along the z-axis. fix_radiative_efficiency: 0 # Whether to fix the radiative efficiency to a constant value. - include_GRMHD_spindown: 0 # Use GRMHD simulation-derived formulas for spindown of BHs, or an analytical formula. \ No newline at end of file + include_GRMHD_spindown: 0 # Use GRMHD simulation-derived formulas for spindown of BHs, or an analytical formula. + accretion_efficiency: 1. # The accretion efficiency (suppression factor of the accretion rate) to use in the thick disc (ADAF), to represent the effects of subgrid winds that take away most of the mass flowing through the accretion disc. + diff --git a/examples/IdealisedCluster/IdealisedCluster_M135/idealised_cluster_M135.yml b/examples/IdealisedCluster/IdealisedCluster_M135/idealised_cluster_M135.yml index c0a8d256ed011b222ad848ad42ec97d636598e2f..ff3fbff68c733a702641229ed99717da452232dd 100644 --- a/examples/IdealisedCluster/IdealisedCluster_M135/idealised_cluster_M135.yml +++ b/examples/IdealisedCluster/IdealisedCluster_M135/idealised_cluster_M135.yml @@ -284,6 +284,7 @@ SPINJETAGN: include_jets: 1 # Global switch whether to include jet feedback [1] or not [0]. Requires the 'include_spin' parameter to also be set to 1. turn_off_radiative_feedback: 1 # Global switch whether to turn off radiative feedback [1] or not [0]. This should only be used if 'include_jets' is set to 1. alpha_acc: 0.2 # Viscous alpha of the subgrid accretion disks. Should be within the 0.1-0.3 range. + mdot_crit_ADAF: 0.03 # The transition normalized accretion rate (Eddington ratio) at which the disc goes from thick (low accretion rates) to thin (high accretion rates). The feedback also changes from kinetic jets to thermal isotropic, respectively. delta_ADAF: 0.2 # Electron heating parameter, which controls the strength of radiative feedback in thick disks. Should be between 0.1 and 0.5. TD_region: B # How to treat the subgrid accretion disk if it is thin, according to the Shakura & Sunyaev (1973) model. If set to B, region b will be used. If set to C, region c will be used. N_jet: 2 # Target number of particles to kick as part of a single jet feedback event. Should be a multiple of 2 to ensure approximate momentum conservation. @@ -299,3 +300,5 @@ SPINJETAGN: fix_jet_direction: 1 # Whether to fix the jet direction to be along the z-axis. fix_radiative_efficiency: 0 # Whether to fix the radiative efficiency to a constant value. include_GRMHD_spindown: 0 # Use GRMHD simulation-derived formulas for spindown of BHs, or an analytical formula. + accretion_efficiency: 1. # The accretion efficiency (suppression factor of the accretion rate) to use in the thick disc (ADAF), to represent the effects of subgrid winds that take away most of the mass flowing through the accretion disc. + diff --git a/examples/IdealisedCluster/IdealisedCluster_M14/idealised_cluster_M14.yml b/examples/IdealisedCluster/IdealisedCluster_M14/idealised_cluster_M14.yml index b29df144d0c3f5e98c918321dd64bfc7ef5a63f6..fd8534d25f42e8bbdf34d7abc4b0b08069c6fe1d 100644 --- a/examples/IdealisedCluster/IdealisedCluster_M14/idealised_cluster_M14.yml +++ b/examples/IdealisedCluster/IdealisedCluster_M14/idealised_cluster_M14.yml @@ -284,6 +284,7 @@ SPINJETAGN: include_jets: 1 # Global switch whether to include jet feedback [1] or not [0]. Requires the 'include_spin' parameter to also be set to 1. turn_off_radiative_feedback: 1 # Global switch whether to turn off radiative feedback [1] or not [0]. This should only be used if 'include_jets' is set to 1. alpha_acc: 0.2 # Viscous alpha of the subgrid accretion disks. Should be within the 0.1-0.3 range. + mdot_crit_ADAF: 0.03 # The transition normalized accretion rate (Eddington ratio) at which the disc goes from thick (low accretion rates) to thin (high accretion rates). The feedback also changes from kinetic jets to thermal isotropic, respectively. delta_ADAF: 0.2 # Electron heating parameter, which controls the strength of radiative feedback in thick disks. Should be between 0.1 and 0.5. TD_region: B # How to treat the subgrid accretion disk if it is thin, according to the Shakura & Sunyaev (1973) model. If set to B, region b will be used. If set to C, region c will be used. N_jet: 2 # Target number of particles to kick as part of a single jet feedback event. Should be a multiple of 2 to ensure approximate momentum conservation. @@ -299,3 +300,4 @@ SPINJETAGN: fix_jet_direction: 1 # Whether to fix the jet direction to be along the z-axis. fix_radiative_efficiency: 0 # Whether to fix the radiative efficiency to a constant value. include_GRMHD_spindown: 0 # Use GRMHD simulation-derived formulas for spindown of BHs, or an analytical formula. + accretion_efficiency: 1. # The accretion efficiency (suppression factor of the accretion rate) to use in the thick disc (ADAF), to represent the effects of subgrid winds that take away most of the mass flowing through the accretion disc. diff --git a/examples/IdealisedCluster/IdealisedCluster_M15/idealised_cluster_M15.yml b/examples/IdealisedCluster/IdealisedCluster_M15/idealised_cluster_M15.yml index f029c7c0194d025f233a176c9edc3dd05d30cdcb..edb92fa35b9507484c5b2e4bf9de34d190fbe7b8 100644 --- a/examples/IdealisedCluster/IdealisedCluster_M15/idealised_cluster_M15.yml +++ b/examples/IdealisedCluster/IdealisedCluster_M15/idealised_cluster_M15.yml @@ -284,6 +284,7 @@ SPINJETAGN: include_jets: 1 # Global switch whether to include jet feedback [1] or not [0]. Requires the 'include_spin' parameter to also be set to 1. turn_off_radiative_feedback: 1 # Global switch whether to turn off radiative feedback [1] or not [0]. This should only be used if 'include_jets' is set to 1. alpha_acc: 0.2 # Viscous alpha of the subgrid accretion disks. Should be within the 0.1-0.3 range. + mdot_crit_ADAF: 0.03 # The transition normalized accretion rate (Eddington ratio) at which the disc goes from thick (low accretion rates) to thin (high accretion rates). The feedback also changes from kinetic jets to thermal isotropic, respectively. delta_ADAF: 0.2 # Electron heating parameter, which controls the strength of radiative feedback in thick disks. Should be between 0.1 and 0.5. TD_region: B # How to treat the subgrid accretion disk if it is thin, according to the Shakura & Sunyaev (1973) model. If set to B, region b will be used. If set to C, region c will be used. N_jet: 2 # Target number of particles to kick as part of a single jet feedback event. Should be a multiple of 2 to ensure approximate momentum conservation. @@ -298,4 +299,5 @@ SPINJETAGN: turn_off_secondary_feedback: 1 # Whether to turn off radiative feedback at very low accretion rates and jets at moderate accretion rates. fix_jet_direction: 1 # Whether to fix the jet direction to be along the z-axis. fix_radiative_efficiency: 0 # Whether to fix the radiative efficiency to a constant value. - include_GRMHD_spindown: 0 # Use GRMHD simulation-derived formulas for spindown of BHs, or an analytical formula. \ No newline at end of file + include_GRMHD_spindown: 0 # Use GRMHD simulation-derived formulas for spindown of BHs, or an analytical formula. + accretion_efficiency: 1. # The accretion efficiency (suppression factor of the accretion rate) to use in the thick disc (ADAF), to represent the effects of subgrid winds that take away most of the mass flowing through the accretion disc. diff --git a/examples/IsolatedGalaxy/IsolatedGalaxy_feedback/getIC.sh b/examples/IsolatedGalaxy/IsolatedGalaxy_feedback/getIC.sh index 32195a97b154e849eacd781ddbe98f59ecf48311..227fcebe73fcf67aab9838f9eb332d32abc1e953 100755 --- a/examples/IsolatedGalaxy/IsolatedGalaxy_feedback/getIC.sh +++ b/examples/IsolatedGalaxy/IsolatedGalaxy_feedback/getIC.sh @@ -1,3 +1,7 @@ #!/bin/bash wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/IsolatedGalaxies/lowres8.hdf5 + +# Initial conditions from Nobels et al. arXiv:2309.13750 +#wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/IsolatedGalaxies_COLIBRE/M5_disk.hdf5 +#wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/IsolatedGalaxies_COLIBRE/M6_disk.hdf5 diff --git a/examples/MHDTests/CosmoVolume/L0150N0064/getIC.sh b/examples/MHDTests/CosmoVolume/L0150N0064/getIC.sh new file mode 100755 index 0000000000000000000000000000000000000000..c69a9cb103bffc86a2c053ac23d789e7a846d6c5 --- /dev/null +++ b/examples/MHDTests/CosmoVolume/L0150N0064/getIC.sh @@ -0,0 +1,2 @@ +#!/bin/bash +wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/MHD_cosmo_box/SWIFT_ICs_L0150N0064.hdf5 diff --git a/examples/MHDTests/CosmoVolume/L0150N0064/run.sh b/examples/MHDTests/CosmoVolume/L0150N0064/run.sh new file mode 100755 index 0000000000000000000000000000000000000000..93adeb99f3911f572b3ad569fdaa5eb96a49ef7a --- /dev/null +++ b/examples/MHDTests/CosmoVolume/L0150N0064/run.sh @@ -0,0 +1,16 @@ +#!/bin/bash + + # Generate the initial conditions if they are not present. +if [ ! -e SWIFT_ICs_L0150N0064.hdf5 ] +then + echo "Fetching initial conditions for the small cosmological volume example..." + ./getIC.sh + echo "Generating the Bfield" + python ../makeIC.py SWIFT_ICs_L0150N0064.hdf5 SWIFT_MHD_ICs_L0150N0064.hdf5 +fi + +# Run SWIFT +../../../../swift --cosmology --hydro --self-gravity --threads=8 swift_param_L0150N0064.yml 2>&1 | tee output.log + +# Plot the temperature evolution +# python3 plotTempEvolution.py diff --git a/examples/MHDTests/CosmoVolume/L0150N0064/swift_param_L0150N0064.yml b/examples/MHDTests/CosmoVolume/L0150N0064/swift_param_L0150N0064.yml new file mode 100644 index 0000000000000000000000000000000000000000..b99d65d132bcdac1e5a318b54f38d8f3c81645ae --- /dev/null +++ b/examples/MHDTests/CosmoVolume/L0150N0064/swift_param_L0150N0064.yml @@ -0,0 +1,89 @@ +# Define some meta-data about the simulation +MetaData: + run_name: SWIFT-L150N0064 + +InternalUnitSystem: + UnitMass_in_cgs: 1.98841e43 # 10^10 M_sun + UnitLength_in_cgs: 3.08567758e24 # 1 Mpc + UnitVelocity_in_cgs: 1e5 # 1 km/s + UnitCurrent_in_cgs: 1 # Amperes + UnitTemp_in_cgs: 1 # Kelvin + +# FLAMINGO cosmology (DES 3x2pt + all external) +Cosmology: +# Omega_cdm: 0.256011 + Omega_cdm: 0.257400 + Omega_b: 0.048600 + Omega_lambda: 0.693922 + h: 0.681 + a_begin: 0.015625 # z=63 + a_end: 1.00000 + +# Parameters governing the time integration +TimeIntegration: + dt_min: 1e-7 + dt_max: 1e-2 + +# Parameters for the self-gravity scheme +Gravity: + eta: 0.025 + MAC: adaptive + theta_cr: 0.7 + epsilon_fmm: 0.001 + comoving_DM_softening: 0.11 + max_physical_DM_softening: 0.035 + comoving_baryon_softening: 0.11 + max_physical_baryon_softening: 0.035 + mesh_side_length: 64 + +# Parameters of the hydro scheme +SPH: + resolution_eta: 1.2348 # "48 Ngb" with the cubic spline kernel + CFL_condition: 0.075 +# minimal_temperature: 100. + +MHD: + parabolic_dedner: 2.0 # + hyperbolic_dedner: 0.1 + resistive_eta: 0.00 + hyperbolic_dedner_divv: 0.0 + monopole_subtraction: 1.0 + artificial_diffusion: 1.0 + +# Parameters governing the snapshots +Snapshots: + subdir: snapshots + basename: snap + delta_time: 1.02 + scale_factor_first: 0.05 + compression: 4 + +# Parameters governing the conserved quantities statistics +Statistics: + delta_time: 1.01 + scale_factor_first: 0.016 + +Scheduler: + max_top_level_cells: 8 + +# Parameters related to the initial conditions +InitialConditions: + file_name: SWIFT_MHD_ICs_L0150N0064.hdf5 + periodic: 1 + cleanup_h_factors: 0 + cleanup_velocity_factors: 0 + generate_gas_in_ics: 0 + cleanup_smoothing_lengths: 0 + remap_ids: 0 + +# Parameters for the EAGLE "equation of state" +#EAGLEEntropyFloor: +# Jeans_density_threshold_H_p_cm3: 0.1 # Physical density above which the EAGLE Jeans limiter entropy floor kicks in expressed in Hydrogen atoms per cm^3. +# Jeans_over_density_threshold: 10. # Overdensity above which the EAGLE Jeans limiter entropy floor can kick in. +# Jeans_temperature_norm_K: 8000 # Temperature of the EAGLE Jeans limiter entropy floor at the density threshold expressed in Kelvin. +# Jeans_gamma_effective: 1.3333333 # Slope the of the EAGLE Jeans limiter entropy floor +# Cool_density_threshold_H_p_cm3: 1e-5 # Physical density above which the EAGLE Cool limiter entropy floor kicks in expressed in Hydrogen atoms per cm^3. +# Cool_over_density_threshold: 10. # Overdensity above which the EAGLE Cool limiter entropy floor can kick in. +# Cool_temperature_norm_K: 8000 # Temperature of the EAGLE Cool limiter entropy floor at the density threshold expressed in Kelvin. +# Cool_gamma_effective: 1. # Slope the of the EAGLE Cool limiter entropy floor + diff --git a/examples/MHDTests/CosmoVolume/L0150N0128/getIC.sh b/examples/MHDTests/CosmoVolume/L0150N0128/getIC.sh new file mode 100755 index 0000000000000000000000000000000000000000..ba48b9c98f27df51ef880b9eac4578e9117a210d --- /dev/null +++ b/examples/MHDTests/CosmoVolume/L0150N0128/getIC.sh @@ -0,0 +1,2 @@ +#!/bin/bash +wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/MHD_cosmo_box/SWIFT_ICs_L0150N0128.hdf5 diff --git a/examples/MHDTests/CosmoVolume/L0150N0128/run.sh b/examples/MHDTests/CosmoVolume/L0150N0128/run.sh new file mode 100755 index 0000000000000000000000000000000000000000..2028d97e3b7f802365e9ae307fb42a1ba4c1fc3f --- /dev/null +++ b/examples/MHDTests/CosmoVolume/L0150N0128/run.sh @@ -0,0 +1,16 @@ +#!/bin/bash + + # Generate the initial conditions if they are not present. +if [ ! -e SWIFT_ICs_L0150N0128.hdf5 ] +then + echo "Fetching initial conditions for the small cosmological volume example..." + ./getIC.sh + echo "Generating the Bfield" + python ../makeIC.py SWIFT_ICs_L0150N0128.hdf5 SWIFT_MHD_ICs_L0150N0128.hdf5 +fi + +# Run SWIFT +../../../../swift --cosmology --hydro --self-gravity --threads=8 swift_param_L0150N0128.yml 2>&1 | tee output.log + +# Plot the temperature evolution +# python3 plotTempEvolution.py diff --git a/examples/MHDTests/CosmoVolume/L0150N0128/swift_param_L0150N0128.yml b/examples/MHDTests/CosmoVolume/L0150N0128/swift_param_L0150N0128.yml new file mode 100644 index 0000000000000000000000000000000000000000..5c45f6090d4868117cd961e7ea2a03632eb784d5 --- /dev/null +++ b/examples/MHDTests/CosmoVolume/L0150N0128/swift_param_L0150N0128.yml @@ -0,0 +1,76 @@ +# Define some meta-data about the simulation +MetaData: + run_name: SWIFT-L150N0128 + +InternalUnitSystem: + UnitMass_in_cgs: 1.98841e43 # 10^10 M_sun + UnitLength_in_cgs: 3.08567758e24 # 1 Mpc + UnitVelocity_in_cgs: 1e5 # 1 km/s + UnitCurrent_in_cgs: 1 # Amperes + UnitTemp_in_cgs: 1 # Kelvin + +# FLAMINGO cosmology (DES 3x2pt + all external) +Cosmology: +# Omega_cdm: 0.256011 + Omega_cdm: 0.257400 + Omega_b: 0.048600 + Omega_lambda: 0.693922 + h: 0.681 + a_begin: 0.015625 # z=63 + a_end: 1.00000 + +# Parameters governing the time integration +TimeIntegration: + dt_min: 1e-7 + dt_max: 1e-2 + +# Parameters for the self-gravity scheme +Gravity: + eta: 0.025 + MAC: adaptive + theta_cr: 0.7 + epsilon_fmm: 0.001 + comoving_DM_softening: 0.058 + max_physical_DM_softening: 0.02 + comoving_baryon_softening: 0.058 + max_physical_baryon_softening: 0.02 + mesh_side_length: 128 + +# Parameters of the hydro scheme +SPH: + resolution_eta: 1.2348 # "48 Ngb" with the cubic spline kernel + CFL_condition: 0.075 + +MHD: + parabolic_dedner: 2.0 # + hyperbolic_dedner: 0.1 + resistive_eta: 0.00 + hyperbolic_dedner_divv: 0.0 + monopole_subtraction: 1.0 + artificial_diffusion: 1.0 + +# Parameters governing the snapshots +Snapshots: + subdir: snapshots + basename: snap + delta_time: 1.02 + scale_factor_first: 0.05 + compression: 4 + +# Parameters governing the conserved quantities statistics +Statistics: + delta_time: 1.01 + scale_factor_first: 0.016 + +Scheduler: + max_top_level_cells: 8 + +# Parameters related to the initial conditions +InitialConditions: + file_name: SWIFT_ICs_L0150N0128.hdf5 + periodic: 1 + cleanup_h_factors: 0 + cleanup_velocity_factors: 0 + generate_gas_in_ics: 0 + cleanup_smoothing_lengths: 0 + remap_ids: 0 diff --git a/examples/MHDTests/CosmoVolume/L0150N0256/getIC.sh b/examples/MHDTests/CosmoVolume/L0150N0256/getIC.sh new file mode 100755 index 0000000000000000000000000000000000000000..b83f323209b0516d10746eb017936a8bd4473902 --- /dev/null +++ b/examples/MHDTests/CosmoVolume/L0150N0256/getIC.sh @@ -0,0 +1,2 @@ +#!/bin/bash +wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/MHD_cosmo_box/SWIFT_ICs_L0150N0256.hdf5 diff --git a/examples/MHDTests/CosmoVolume/L0150N0256/run.sh b/examples/MHDTests/CosmoVolume/L0150N0256/run.sh new file mode 100755 index 0000000000000000000000000000000000000000..b48b7bb6394f5ac92e30e1f61ef3fc5831c84df1 --- /dev/null +++ b/examples/MHDTests/CosmoVolume/L0150N0256/run.sh @@ -0,0 +1,16 @@ +#!/bin/bash + + # Generate the initial conditions if they are not present. +if [ ! -e SWIFT_ICs_L0150N0256.hdf5 ] +then + echo "Fetching initial conditions for the small cosmological volume example..." + ./getIC.sh + echo "Generating the Bfield" + python ../makeIC.py SWIFT_ICs_L0150N0256.hdf5 SWIFT_MHD_ICs_L0150N0256.hdf5 +fi + +# Run SWIFT +../../../../swift --cosmology --hydro --self-gravity --threads=16 swift_param_L0150N0256.yml 2>&1 | tee output.log + +# Plot the temperature evolution +# python3 plotTempEvolution.py diff --git a/examples/MHDTests/CosmoVolume/L0150N0256/swift_param_L0150N0256.yml b/examples/MHDTests/CosmoVolume/L0150N0256/swift_param_L0150N0256.yml new file mode 100644 index 0000000000000000000000000000000000000000..601dcf4329a6090adff6fc6afb0f02446a95c95b --- /dev/null +++ b/examples/MHDTests/CosmoVolume/L0150N0256/swift_param_L0150N0256.yml @@ -0,0 +1,76 @@ +# Define some meta-data about the simulation +MetaData: + run_name: SWIFT-L150N0256 + +InternalUnitSystem: + UnitMass_in_cgs: 1.98841e43 # 10^10 M_sun + UnitLength_in_cgs: 3.08567758e24 # 1 Mpc + UnitVelocity_in_cgs: 1e5 # 1 km/s + UnitCurrent_in_cgs: 1 # Amperes + UnitTemp_in_cgs: 1 # Kelvin + +# FLAMINGO cosmology (DES 3x2pt + all external) +Cosmology: +# Omega_cdm: 0.256011 + Omega_cdm: 0.257400 + Omega_b: 0.048600 + Omega_lambda: 0.693922 + h: 0.681 + a_begin: 0.015625 # z=63 + a_end: 1.00000 + +# Parameters governing the time integration +TimeIntegration: + dt_min: 1e-7 + dt_max: 1e-2 + +# Parameters for the self-gravity scheme +Gravity: + eta: 0.025 + MAC: adaptive + theta_cr: 0.7 + epsilon_fmm: 0.001 + comoving_DM_softening: 0.03 + max_physical_DM_softening: 0.01 + comoving_baryon_softening: 0.03 + max_physical_baryon_softening: 0.01 + mesh_side_length: 256 + +# Parameters of the hydro scheme +SPH: + resolution_eta: 1.2348 # "48 Ngb" with the cubic spline kernel + CFL_condition: 0.075 + +MHD: + parabolic_dedner: 2.0 # + hyperbolic_dedner: 0.1 + resistive_eta: 0.00 + hyperbolic_dedner_divv: 0.0 + monopole_subtraction: 1.0 + artificial_diffusion: 1.0 + +# Parameters governing the snapshots +Snapshots: + subdir: snapshots + basename: snap + delta_time: 1.02 + scale_factor_first: 0.05 + compression: 4 + +# Parameters governing the conserved quantities statistics +Statistics: + delta_time: 1.01 + scale_factor_first: 0.016 + +Scheduler: + max_top_level_cells: 16 + +# Parameters related to the initial conditions +InitialConditions: + file_name: SWIFT_ICs_L0150N0256.hdf5 + periodic: 1 + cleanup_h_factors: 0 + cleanup_velocity_factors: 0 + generate_gas_in_ics: 0 + cleanup_smoothing_lengths: 0 + remap_ids: 0 diff --git a/examples/MHDTests/CosmoVolume/makeIC.py b/examples/MHDTests/CosmoVolume/makeIC.py new file mode 100644 index 0000000000000000000000000000000000000000..de3488594fea3529bed5437dc544e06edc6908e0 --- /dev/null +++ b/examples/MHDTests/CosmoVolume/makeIC.py @@ -0,0 +1,56 @@ +############################# +# This file is part of SWIFT +# Copyright +############################# + +import h5py +import numpy as np +import matplotlib.pyplot as plt +import sys +import os + +# Parameters + + +fileOutputName = sys.argv[2] +fileInputName = sys.argv[1] + +###---------------------------### + +infile = h5py.File(fileInputName, "r") + +head = infile["/Header"] +pos_in = infile["/PartType0/Coordinates"][:, :] + +BoxSize = head.attrs["BoxSize"] +afact = head.attrs["Time"] +N_in = head.attrs["NumPart_Total"][0] +print(BoxSize) +print(N_in) +infile.close() + +Bini = 1E-12 +wavelen = BoxSize / 10.0 +wavenum = 2.0 * np.pi / wavelen + +B = np.zeros((N_in,3)) +A = np.zeros((N_in,3)) + +Aini = Bini / wavenum; + +A[:,0] = Aini * (np.sin(pos_in[:,2]*wavenum) + np.cos(pos_in[:,1]*wavenum)) +A[:,1] = Aini * (np.sin(pos_in[:,0]*wavenum) + np.cos(pos_in[:,2]*wavenum)) +A[:,2] = Aini * (np.sin(pos_in[:,1]*wavenum) + np.cos(pos_in[:,0]*wavenum)) +B[:,:] = wavenum * A[:,:] + +os.system("cp "+fileInputName+" "+fileOutputName) +# File +fileOutput = h5py.File(fileOutputName, "a") + + +# Particle group +grp = fileOutput.require_group("/PartType0") +grp.create_dataset("MagneticFluxDensities", data=B, dtype="f") +grp.create_dataset("MagneticVectorPotentials", data=A, dtype="f") + +fileOutput.close() diff --git a/examples/MHDTests/CosmoVolume/plotProjection.py b/examples/MHDTests/CosmoVolume/plotProjection.py new file mode 100755 index 0000000000000000000000000000000000000000..4dc7a8d59e6f810f9dce9d91c8b57c11cb7a2e23 --- /dev/null +++ b/examples/MHDTests/CosmoVolume/plotProjection.py @@ -0,0 +1,178 @@ +#!/usr/bin/env python3 + +# ------------------------------------------------ +# Plots a projection of the DM and baryon mass +# usage: +# $ python3 plotProjection.py <snapnr> +# where <snapnr> is number of snapshot to plot +# ------------------------------------------------ + +import numpy as np +import sys +import os +import matplotlib + +matplotlib.use("Agg") +from matplotlib import pyplot as plt +from matplotlib.colors import LogNorm +from mpl_toolkits.axes_grid1 import make_axes_locatable + +from swiftsimio import load +from swiftsimio.visualisation.projection import project_pixel_grid, project_gas +from swiftsimio.visualisation.smoothing_length_generation import ( + generate_smoothing_lengths, +) + +from unyt import msun, kpc + + +def set_colorbar(ax, im): + """ + Adapt the colorbar a bit for axis object <ax> and + imshow instance <im> + """ + divider = make_axes_locatable(ax) + cax = divider.append_axes("right", size="5%", pad=0.05) + plt.colorbar(im, cax=cax) + return + + +# Grab snapshot + +snapshot_basename = "snapshots/snap" + +try: + snapnr = sys.argv[1] +except IndexError: + print("You need to provide the index of the snapshot to plot.") + print("E.g. python3 plotProjection.py 3") + quit() + +try: + snapnr_int = int(snapnr) +except ValueError: + print("<snapnr> must be an integer.") + print("You provided :'" + snapnr + "'") + +file = snapshot_basename + "_" + str(snapnr_int).zfill(4) + ".hdf5" + +if not os.path.isfile(file): + print("Didn't find snapshot", file) + quit() + + +# Load data +data = load(file) +meta = data.metadata +boxsize = meta.boxsize +extent = [0, boxsize[0].v, 0, boxsize[1].v] + +# Generate smoothing lengths for the dark matter +data.dark_matter.smoothing_length = generate_smoothing_lengths( + data.dark_matter.coordinates, + data.metadata.boxsize, + kernel_gamma=1.8, + neighbours=57, + speedup_fac=2, + dimension=3, +) + +# Project the dark matter mass +dm_mass = project_pixel_grid( + # Note here that we pass in the dark matter dataset not the whole + # data object, to specify what particle type we wish to visualise + data=data.dark_matter, + boxsize=data.metadata.boxsize, + resolution=1024, + project="masses", + parallel=True, + region=None, +) + + +# Project the gas mass +mass_map = project_gas(data, resolution=1024, project="masses", parallel=True) +mass_map.convert_to_units(msun / kpc ** 2) + +######## Magnetic STuff + +B = data.gas.magnetic_flux_densities +normB = np.sqrt(B[:, 0] ** 2 + B[:, 1] ** 2 + B[:, 2] ** 2) +divB = data.gas.magnetic_divergences +h = data.gas.smoothing_lengths + +data.gas.B_Mag = normB +data.gas.DivB_error = np.maximum(h * abs(divB) / normB, 1e-6) + + +divb_map = project_gas(data, resolution=1024, project="DivB_error", parallel=True) +# divb_map.convert_to_units(msun / kpc ** 2) + +bfld_map = project_gas(data, resolution=1024, project="B_Mag", parallel=True) +# bfld_map.convert_to_units(msun / kpc ** 2) + +# Make figure an plot +# fig = plt.figure(figsize=(12, 5), dpi=200) +fig = plt.figure(figsize=(12, 12), dpi=200) + +ax1 = fig.add_subplot(221) +im1 = ax1.imshow( + dm_mass.T, origin="lower", extent=extent, cmap="inferno", norm=LogNorm() +) +ax1.set_title("Dark Matter Mass", usetex=True) +set_colorbar(ax1, im1) + +ax2 = fig.add_subplot(222) +im2 = ax2.imshow( + mass_map.T, + origin="lower", + extent=extent, + cmap="magma", + norm=LogNorm(vmax=5e7, vmin=1e5), +) +ax2.set_title("Baryon Mass", usetex=True) +set_colorbar(ax2, im2) + +ax3 = fig.add_subplot(223) +im3 = ax3.imshow( + divb_map.T, + origin="lower", + extent=extent, + cmap="cividis", + norm=LogNorm(vmax=1e7, vmin=1e5), +) +ax3.set_title("divB", usetex=True) +set_colorbar(ax3, im3) + +ax4 = fig.add_subplot(224) +im4 = ax4.imshow( + bfld_map.T, + origin="lower", + extent=extent, + cmap="magma", + norm=LogNorm(vmax=1e1, vmin=1e-8), +) +ax4.set_title("Magnetic Field", usetex=True) +set_colorbar(ax4, im4) + +# Add xlabels +xunits = data.dark_matter.coordinates.units.latex_representation() +for ax in [ax1, ax2, ax3, ax4]: + ax.set_xlabel("x [" + xunits + "]", usetex=True) + ax.set_ylabel("y [" + xunits + "]", usetex=True) + + +# Add title +title = file.replace("_", r"\_") # exception handle underscore for latex +if meta.cosmology is not None: + title += ", $z$ = {0:.3f}".format(meta.z) +title += ", $t$ = {0:.2e}".format(meta.time.to("Gyr")) +fig.suptitle(title, usetex=True) + +# Save figure +plt.tight_layout() +figname = file[:-5] + ".png" +plt.savefig(figname) +plt.close() + +print("saved", figname) diff --git a/examples/MHDTests/FastRotor/FDI/run.sh b/examples/MHDTests/FastRotor/FDI/run.sh new file mode 100755 index 0000000000000000000000000000000000000000..b692d34e84b22a6de204a20832a02ce4bfca56fe --- /dev/null +++ b/examples/MHDTests/FastRotor/FDI/run.sh @@ -0,0 +1,5 @@ +# Run SWIFT +../../../../sw_FDI --hydro --threads=16 ../FR_schemes.yml 2>&1 > out.log + +# Plot the temperature evolution +python3 ../plot_all.py 0 60 diff --git a/examples/MHDTests/FastRotor/FR_schemes.yml b/examples/MHDTests/FastRotor/FR_schemes.yml new file mode 100644 index 0000000000000000000000000000000000000000..572f3fb6dab86403e8c5a4cd776ddb4e15f1eb9c --- /dev/null +++ b/examples/MHDTests/FastRotor/FR_schemes.yml @@ -0,0 +1,52 @@ +# Define the system of units to use internally. +InternalUnitSystem: + UnitMass_in_cgs: 1 # Grams + UnitLength_in_cgs: 1 # Centimeters + UnitVelocity_in_cgs: 1 # Centimeters per second + UnitCurrent_in_cgs: 1 # Amperes + UnitTemp_in_cgs: 1 # Kelvin + +# Parameters governing the time integration +TimeIntegration: + time_begin: 0. # The starting time of the simulation (in internal units). + time_end: 0.6 # The end time of the simulation (in internal units). + dt_min: 1e-9 # The minimal time-step size of the simulation (in internal units). + dt_max: 1e-3 # The maximal time-step size of the simulation (in internal units). + +# Parameters governing the snapshots +Snapshots: + basename: FastRotor_LR # Common part of the name of output files + time_first: 0. # Time of the first output (in internal units) + delta_time: 0.01 # Time difference between consecutive outputs (in internal units) + compression: 1 + +# Parameters governing the conserved quantities statistics +Statistics: + delta_time: 1e-2 # Time between statistics output + +# Parameters for the hydrodynamics scheme +SPH: + resolution_eta: 1.2348 # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel$ + CFL_condition: 0.075 # Courant-Friedrich-Levy condition for time integration. + viscosity_alpha: 2.0 + +# Parameters for MHD schemes +MHD: + monopole_subtraction: 1.0 + artificial_diffusion: 0.0 + hyperbolic_dedner: 0.1 + hyperbolic_dedner_divv: 0.0 + parabolic_dedner: 2.0 + resistive_eta: 0.01 + +# Physical parameters +PhysicalConstants: + mu_0: 1.0 + +# Parameters related to the initial conditions +InitialConditions: + file_name: ../FastRotor_LR.hdf5 # The file to read + periodic: 1 + +Scheduler: + max_top_level_cells: 256 diff --git a/examples/MHDTests/FastRotor/FastRotor.yml b/examples/MHDTests/FastRotor/FastRotor.yml index 3439619c64435f0a0e6e91bfcce4784242fdc8a3..b57c45a3c728669c9e5711ef07d3dbb5a4dbe840 100644 --- a/examples/MHDTests/FastRotor/FastRotor.yml +++ b/examples/MHDTests/FastRotor/FastRotor.yml @@ -32,9 +32,12 @@ SPH: # Parameters for MHD schemes MHD: - hyperbolic_dedner: 1.0 + monopole_subtraction: 1.0 + artificial_diffusion: 0.0 + hyperbolic_dedner: 0.1 + hyperbolic_dedner_divv: 0.0 parabolic_dedner: 2.0 - resistive_eta: 0.0 + resistive_eta: 0.01 # Physical parameters PhysicalConstants: diff --git a/examples/MHDTests/FastRotor/ODI/run.sh b/examples/MHDTests/FastRotor/ODI/run.sh new file mode 100755 index 0000000000000000000000000000000000000000..508e3a98464bce062f5d9f9957d78cedf2bc082c --- /dev/null +++ b/examples/MHDTests/FastRotor/ODI/run.sh @@ -0,0 +1,5 @@ +# Run SWIFT +../../../../sw_ODI --hydro --threads=16 ../FR_schemes.yml 2>&1 > out.log + +# Plot the temperature evolution +python3 ../plot_all.py 0 60 diff --git a/examples/MHDTests/FastRotor/README.md b/examples/MHDTests/FastRotor/README.md new file mode 100644 index 0000000000000000000000000000000000000000..2fa79cb9fabf8667f2f840556a58191a68f47763 --- /dev/null +++ b/examples/MHDTests/FastRotor/README.md @@ -0,0 +1,21 @@ +Fast Rotor in MHD +Lodrillo, P. & Del Zanna, L., 2000ApJ..530...508 +Stasyszyn et al , 2013MNRAS.428...13S +Following also Seo & Ryu, 2023 +############################## +# Basically the folloing script compiles and run the test in different +# folders for each scheme. usefull for testing +./runSchemes.sh [what Scheme] [FOLDER_TAIL] +[what scheme]: +vep: vector potentials +odi: Oresti's direct induction +fdi: simple direct induction +[FOLDER_TAIL]: +the trailing name for the folders created +############################## + +TODO> +-check for symmetry with the between half of the domain +-make a cut +-upload reference runs or values + diff --git a/examples/MHDTests/FastRotor/VeP/run.sh b/examples/MHDTests/FastRotor/VeP/run.sh new file mode 100755 index 0000000000000000000000000000000000000000..c5f65d08a00a6573f31cdde892ae91c00354b5f8 --- /dev/null +++ b/examples/MHDTests/FastRotor/VeP/run.sh @@ -0,0 +1,5 @@ +# Run SWIFT +../../../../sw_VeP --hydro --threads=16 ../FR_schemes.yml 2>&1 > out.log + +# Plot the temperature evolution +python3 ../plot_all.py 0 60 diff --git a/examples/MHDTests/FastRotor/getGlass.sh b/examples/MHDTests/FastRotor/getGlass.sh old mode 100644 new mode 100755 diff --git a/examples/MHDTests/FastRotor/makeIC.py b/examples/MHDTests/FastRotor/makeIC.py index 7640b8b4bac643356f373eff357b762a1828618b..22f4daafaac0cdc96d0cc8772c6700a7850ba20d 100644 --- a/examples/MHDTests/FastRotor/makeIC.py +++ b/examples/MHDTests/FastRotor/makeIC.py @@ -1,92 +1,218 @@ +############################# +# This file is part of SWIFT +# Copyright +############################# + import h5py import numpy as np import matplotlib.pyplot as plt -###---------------------------### - # Parameters r_in = 0.1 rho_in_0 = 10.0 rho_out_0 = 1.0 P_0 = 1.0 -B_0 = 0.0 # 2.5 / np.sqrt(np.pi) -omega_0 = 20.0 +B_0 = 2.5 / np.sqrt(np.pi) +omega_0 = -20.0 +gamma = 1.4 -fileOutputName = "FastRotor_no_B.hdf5" +fileOutputName = "FastRotor_LR.hdf5" ###---------------------------### -N_side = 8 -N_face = N_side ** 2 -N_unit_cell = N_side ** 3 +glass = h5py.File("glassCube_32.hdf5", "r") -unit_cell = np.zeros((N_unit_cell, 3)) -unit_cell_smoothing_lengths = np.ones(N_unit_cell) +unit_cell = glass["/PartType0/Coordinates"][:, :] +h_unit_cell = glass["/PartType0/SmoothingLength"][:] -for i in range(0, N_unit_cell): - unit_cell[i, 0] = (i % N_side) / N_side - # print(unit_cell[i,0]) - unit_cell[i, 1] = ((i // N_side) % N_side) / N_side - unit_cell[i, 2] = ((i // N_face) % N_side) / N_side +N_unit_cell = len(h_unit_cell) +times = 2 -# plt.scatter(unit_cell[:,0],unit_cell[:,1]) -# plt.savefig("test.png") +ratio = np.cbrt(rho_in_0 / rho_out_0) -# plt.scatter(unit_cell[:,0],unit_cell[:,2]) -# plt.savefig("test2.png") +###---------------------------### -# plt.scatter(unit_cell[:,1],unit_cell[:,2]) -# plt.savefig("test3.png") +cx_out = times +cy_out = times +cz_out = 1 -N_UC_reps = 10 -cx = N_UC_reps -cy = N_UC_reps -cz = 1 +cx_in = int(np.ceil(ratio * cx_out)) +cy_in = int(np.ceil(ratio * cy_out)) +cz_in = int(np.ceil(ratio * cz_out)) lx = 1.0 -ly = 1.0 -lz = 1.0 / float(N_UC_reps) -vol = lx * ly * lz +ly = lx +lz = lx * float(cz_out) / float(cx_out) -pos = np.zeros((int(cx * cy * cz * N_unit_cell), 3)) -h = np.zeros(int(cx * cy * cz * N_unit_cell)) +print(lz) -N_particles = N_unit_cell * cx * cy * cz +lx_c = lx / 2 +ly_c = ly / 2 +lz_c = lz / 2 -for i in range(0, cx): - for j in range(0, cy): - for k in range(0, cz): - ind_low = i * N_unit_cell - ind_upper = ind_low + N_unit_cell +lx_in = lx + (np.ceil(ratio * cx_out) / (ratio * cx_out) - 1.0) * lx +ly_in = ly + (np.ceil(ratio * cy_out) / (ratio * cy_out) - 1.0) * ly +lz_in = lz + (np.ceil(ratio * cz_out) / (ratio * cz_out) - 1.0) * lz - pos[ind_low:ind_upper, 0] = unit_cell[:, 0] + i +print(lz_in) - pos[ind_low:ind_upper, 1] = unit_cell[:, 1] + j +N_out = N_unit_cell * cx_out * cy_out * cz_out +N_in = N_unit_cell * cx_in * cy_in * cz_in - pos[ind_low:ind_upper, 2] = unit_cell[:, 2] + k +pos_out = np.zeros((int(N_out), 3)) +h_out = np.zeros(N_out) - h[ind_low:ind_upper] = unit_cell_smoothing_lengths +pos_in = np.zeros((int(N_in), 3)) +h_in = np.zeros(N_in) -pos[:, 0] = pos[:, 0] * lx / cx -pos[:, 1] = pos[:, 1] * ly / cy -pos[:, 2] = pos[:, 2] * lz / cz +c0 = 0 +c1 = N_unit_cell -fig = plt.figure() -ax = fig.add_subplot(111, projection="3d") -ax.scatter(pos[:, 0], pos[:, 1], pos[:, 2], c="b") +for i in range(0, cx_out): + for j in range(0, cy_out): + for k in range(0, cz_out): + pos_out[c0:c1, 0] = unit_cell[:, 0] + i + pos_out[c0:c1, 1] = unit_cell[:, 1] + j + pos_out[c0:c1, 2] = unit_cell[:, 2] + k + h_out[c0:c1] = h_unit_cell[:] + c0 = c0 + N_unit_cell + c1 = c1 + N_unit_cell -plt.savefig("test.png") +pos_out[:, 0] = pos_out[:, 0] * lx / cx_out +pos_out[:, 1] = pos_out[:, 1] * ly / cy_out +pos_out[:, 2] = pos_out[:, 2] * lz / cz_out +h_out = h_out / cx_out -quit() -###---------------------------### +c0 = 0 +c1 = N_unit_cell + +for i in range(0, cx_in): + for j in range(0, cy_in): + for k in range(0, cz_in): + pos_in[c0:c1, 0] = unit_cell[:, 0] + i + pos_in[c0:c1, 1] = unit_cell[:, 1] + j + pos_in[c0:c1, 2] = unit_cell[:, 2] + k + h_in[c0:c1] = h_unit_cell[:] + c0 = c0 + N_unit_cell + c1 = c1 + N_unit_cell + +pos_in[:, 0] = pos_in[:, 0] * lx_in / cx_in +pos_in[:, 1] = pos_in[:, 1] * ly_in / cy_in +pos_in[:, 2] = pos_in[:, 2] * lz_in / cz_in +h_in = h_in / cx_in + +vol = lx * ly * lz + +pos_out_f = pos_out[ + (pos_out[:, 0] - lx_c) ** 2 + (pos_out[:, 1] - ly_c) ** 2 > r_in ** 2 +] +pos_in_f = pos_in[(pos_in[:, 0] - lx_c) ** 2 + (pos_in[:, 1] - ly_c) ** 2 < r_in ** 2] -# Generate extra arrays +h_out_f = h_out[(pos_out[:, 0] - lx_c) ** 2 + (pos_out[:, 1] - ly_c) ** 2 > r_in ** 2] +h_in_f = h_in[(pos_in[:, 0] - lx_c) ** 2 + (pos_in[:, 1] - ly_c) ** 2 < r_in ** 2] +h_in_f = h_in_f[pos_in_f[:, 2] < lz] +pos_in_f = pos_in_f[pos_in_f[:, 2] < lz] + +pos = np.append(pos_out_f, pos_in_f, axis=0) +h = np.append(h_out_f, h_in_f, axis=0) +N = len(pos) +N_out_f = len(pos_out_f) +N_in_f = len(pos_in_f) + +print(pos_out_f.shape) +print(pos_in_f.shape) +print(pos.shape) + +###---------------------------### + +rot = np.sqrt( + (pos[:, 0] - lx_c) ** 2 + (pos[:, 1] - ly_c) ** 2 +) # + (pos[:,2]-lz_c)**2) +theta = np.arctan2((pos[:, 1] - ly_c), (pos[:, 0] - lx_c)) v = np.zeros((N, 3)) -B = np.zeros((N, 3)) +#B = np.zeros((N, 3)) +#A = np.zeros((N, 3)) ids = np.linspace(1, N, N) -m = np.ones(N) * rho0 * vol / N -u = np.ones(N) * P0 / (rho0 * (gamma - 1.0)) +m = np.ones(N) * rho_out_0 * vol / N_out +u = np.ones(N) +u[:N_out_f] = P_0 / (rho_out_0 * (gamma - 1)) +u[N_out_f:] = P_0 / (rho_in_0 * (gamma - 1)) + +v[N_out_f:, 0] = -rot[N_out_f:] * omega_0 * np.sin(theta[N_out_f:]) +v[N_out_f:, 1] = rot[N_out_f:] * omega_0 * np.cos(theta[N_out_f:]) + +#B[:, 0] = B_0 + + +###---------------------------### +N2=int(2*N) +p=np.zeros((N2, 3)) +p[:N,0]=pos[:,0] +p[N:,0]=pos[:,0] +p[:N,1]=pos[:,1] +p[N:,1]=pos[:,1]+1.0 +p[:N,2]=pos[:,2] +p[N:,2]=pos[:,2] +pos=p +hh =np.zeros(N2) +hh[:N]=h +hh[N:]=h +h=hh +vel=np.zeros((N2,3)) +vel[:N,:]=v[:,:] +vel[N: ,:]=v[:,:] +v=vel +ids = np.linspace(1, N2, N2) +m = np.ones(N2) * rho_out_0 * vol / N_out +uu =np.zeros(N2) +uu[:N]=u +uu[N:]=u +u=uu +B = np.zeros((N2, 3)) +A = np.zeros((N2, 3)) +B[:N, 0] = B_0 +B[N:, 0] = -B_0 +A[:N, 2] = B_0 * pos[:N,1] +A[N:, 2] = B_0 * (2.0-pos[N:,1]) + +# File +fileOutput = h5py.File(fileOutputName, "w") + +# Header +grp = fileOutput.create_group("/Header") +grp.attrs["BoxSize"] = [lx, 2.* ly, lz] ##### +grp.attrs["NumPart_Total"] = [2*N, 0, 0, 0, 0, 0] +grp.attrs["NumPart_Total_HighWord"] = [0, 0, 0, 0, 0, 0] +grp.attrs["NumPart_ThisFile"] = [2*N, 0, 0, 0, 0, 0] +grp.attrs["Time"] = 0.0 +grp.attrs["NumFileOutputsPerSnapshot"] = 1 +grp.attrs["MassTable"] = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0] +grp.attrs["Flag_Entropy_ICs"] = [0, 0, 0, 0, 0, 0] +grp.attrs["Dimension"] = 3 + +# Units +grp = fileOutput.create_group("/Units") +grp.attrs["Unit length in cgs (U_L)"] = 1.0 +grp.attrs["Unit mass in cgs (U_M)"] = 1.0 +grp.attrs["Unit time in cgs (U_t)"] = 1.0 +grp.attrs["Unit current in cgs (U_I)"] = 1.0 +grp.attrs["Unit temperature in cgs (U_T)"] = 1.0 + +# Particle group +grp = fileOutput.create_group("/PartType0") +grp.create_dataset("Coordinates", data=pos, dtype="d") +grp.create_dataset("Velocities", data=v, dtype="f") +grp.create_dataset("Masses", data=m, dtype="f") +grp.create_dataset("SmoothingLength", data=h, dtype="f") +grp.create_dataset("InternalEnergy", data=u, dtype="f") +grp.create_dataset("ParticleIDs", data=ids, dtype="L") +grp.create_dataset("MagneticFluxDensities", data=B, dtype="f") +grp.create_dataset("MagneticVectorPotentials", data = A, dtype = 'f') +# grp.create_dataset("EPalpha", data = epa, dtype = 'f') +# grp.create_dataset("EPbeta" , data = epb, dtype = 'f') + +fileOutput.close() diff --git a/examples/MHDTests/FastRotor/makeICnew.py b/examples/MHDTests/FastRotor/makeICnew.py deleted file mode 100644 index 37c73eb2536cd663c9ebba9195816b4af7d6a87a..0000000000000000000000000000000000000000 --- a/examples/MHDTests/FastRotor/makeICnew.py +++ /dev/null @@ -1,187 +0,0 @@ -############################# -# This file is part of SWIFT -# Copyright -############################# - -import h5py -import numpy as np -import matplotlib.pyplot as plt - -# Parameters - -r_in = 0.1 -rho_in_0 = 10.0 -rho_out_0 = 1.0 -P_0 = 1.0 -B_0 = 2.5 / np.sqrt(np.pi) -omega_0 = -20.0 -gamma = 1.4 - -fileOutputName = "FastRotor_LR.hdf5" - -###---------------------------### - -glass = h5py.File("glassCube_32.hdf5", "r") - -unit_cell = glass["/PartType0/Coordinates"][:, :] -h_unit_cell = glass["/PartType0/SmoothingLength"][:] - -N_unit_cell = len(h_unit_cell) -times = 2 - -ratio = np.cbrt(rho_in_0 / rho_out_0) - -###---------------------------### - -cx_out = times -cy_out = times -cz_out = 1 - -cx_in = int(np.ceil(ratio * cx_out)) -cy_in = int(np.ceil(ratio * cy_out)) -cz_in = int(np.ceil(ratio * cz_out)) - -lx = 1.0 -ly = lx -lz = lx * float(cz_out) / float(cx_out) - -print(lz) - -lx_c = lx / 2 -ly_c = ly / 2 -lz_c = lz / 2 - -lx_in = lx + (np.ceil(ratio * cx_out) / (ratio * cx_out) - 1.0) * lx -ly_in = ly + (np.ceil(ratio * cy_out) / (ratio * cy_out) - 1.0) * ly -lz_in = lz + (np.ceil(ratio * cz_out) / (ratio * cz_out) - 1.0) * lz - -print(lz_in) - -N_out = N_unit_cell * cx_out * cy_out * cz_out -N_in = N_unit_cell * cx_in * cy_in * cz_in - -pos_out = np.zeros((int(N_out), 3)) -h_out = np.zeros(N_out) - -pos_in = np.zeros((int(N_in), 3)) -h_in = np.zeros(N_in) - - -c0 = 0 -c1 = N_unit_cell - -for i in range(0, cx_out): - for j in range(0, cy_out): - for k in range(0, cz_out): - pos_out[c0:c1, 0] = unit_cell[:, 0] + i - pos_out[c0:c1, 1] = unit_cell[:, 1] + j - pos_out[c0:c1, 2] = unit_cell[:, 2] + k - h_out[c0:c1] = h_unit_cell[:] - c0 = c0 + N_unit_cell - c1 = c1 + N_unit_cell - -pos_out[:, 0] = pos_out[:, 0] * lx / cx_out -pos_out[:, 1] = pos_out[:, 1] * ly / cy_out -pos_out[:, 2] = pos_out[:, 2] * lz / cz_out -h_out = h_out / cx_out - - -c0 = 0 -c1 = N_unit_cell - -for i in range(0, cx_in): - for j in range(0, cy_in): - for k in range(0, cz_in): - pos_in[c0:c1, 0] = unit_cell[:, 0] + i - pos_in[c0:c1, 1] = unit_cell[:, 1] + j - pos_in[c0:c1, 2] = unit_cell[:, 2] + k - h_in[c0:c1] = h_unit_cell[:] - c0 = c0 + N_unit_cell - c1 = c1 + N_unit_cell - -pos_in[:, 0] = pos_in[:, 0] * lx_in / cx_in -pos_in[:, 1] = pos_in[:, 1] * ly_in / cy_in -pos_in[:, 2] = pos_in[:, 2] * lz_in / cz_in -h_in = h_in / cx_in - -vol = lx * ly * lz - -pos_out_f = pos_out[ - (pos_out[:, 0] - lx_c) ** 2 + (pos_out[:, 1] - ly_c) ** 2 > r_in ** 2 -] -pos_in_f = pos_in[(pos_in[:, 0] - lx_c) ** 2 + (pos_in[:, 1] - ly_c) ** 2 < r_in ** 2] - -h_out_f = h_out[(pos_out[:, 0] - lx_c) ** 2 + (pos_out[:, 1] - ly_c) ** 2 > r_in ** 2] -h_in_f = h_in[(pos_in[:, 0] - lx_c) ** 2 + (pos_in[:, 1] - ly_c) ** 2 < r_in ** 2] -h_in_f = h_in_f[pos_in_f[:, 2] < lz] - -pos_in_f = pos_in_f[pos_in_f[:, 2] < lz] - -pos = np.append(pos_out_f, pos_in_f, axis=0) -h = np.append(h_out_f, h_in_f, axis=0) -N = len(pos) -N_out_f = len(pos_out_f) -N_in_f = len(pos_in_f) - -print(pos_out_f.shape) -print(pos_in_f.shape) -print(pos.shape) - -###---------------------------### - -rot = np.sqrt( - (pos[:, 0] - lx_c) ** 2 + (pos[:, 1] - ly_c) ** 2 -) # + (pos[:,2]-lz_c)**2) -theta = np.arctan2((pos[:, 1] - ly_c), (pos[:, 0] - lx_c)) -v = np.zeros((N, 3)) -B = np.zeros((N, 3)) -ids = np.linspace(1, N, N) -m = np.ones(N) * rho_out_0 * vol / N_out -u = np.ones(N) -u[:N_out_f] = P_0 / (rho_out_0 * (gamma - 1)) -u[N_out_f:] = P_0 / (rho_in_0 * (gamma - 1)) - -v[N_out_f:, 0] = -rot[N_out_f:] * omega_0 * np.sin(theta[N_out_f:]) -v[N_out_f:, 1] = rot[N_out_f:] * omega_0 * np.cos(theta[N_out_f:]) - -B[:, 0] = B_0 - -###---------------------------### - -# File -fileOutput = h5py.File(fileOutputName, "w") - -# Header -grp = fileOutput.create_group("/Header") -grp.attrs["BoxSize"] = [lx, ly, lz] ##### -grp.attrs["NumPart_Total"] = [N, 0, 0, 0, 0, 0] -grp.attrs["NumPart_Total_HighWord"] = [0, 0, 0, 0, 0, 0] -grp.attrs["NumPart_ThisFile"] = [N, 0, 0, 0, 0, 0] -grp.attrs["Time"] = 0.0 -grp.attrs["NumFileOutputsPerSnapshot"] = 1 -grp.attrs["MassTable"] = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0] -grp.attrs["Flag_Entropy_ICs"] = [0, 0, 0, 0, 0, 0] -grp.attrs["Dimension"] = 3 - -# Units -grp = fileOutput.create_group("/Units") -grp.attrs["Unit length in cgs (U_L)"] = 1.0 -grp.attrs["Unit mass in cgs (U_M)"] = 1.0 -grp.attrs["Unit time in cgs (U_t)"] = 1.0 -grp.attrs["Unit current in cgs (U_I)"] = 1.0 -grp.attrs["Unit temperature in cgs (U_T)"] = 1.0 - -# Particle group -grp = fileOutput.create_group("/PartType0") -grp.create_dataset("Coordinates", data=pos, dtype="d") -grp.create_dataset("Velocities", data=v, dtype="f") -grp.create_dataset("Masses", data=m, dtype="f") -grp.create_dataset("SmoothingLength", data=h, dtype="f") -grp.create_dataset("InternalEnergy", data=u, dtype="f") -grp.create_dataset("ParticleIDs", data=ids, dtype="L") -grp.create_dataset("MagneticFluxDensities", data=B, dtype="f") -# grp.create_dataset("VecPot", data = vp, dtype = 'f') -# grp.create_dataset("EPalpha", data = epa, dtype = 'f') -# grp.create_dataset("EPbeta" , data = epb, dtype = 'f') - -fileOutput.close() diff --git a/examples/MHDTests/FastRotor/makeICnew2.py b/examples/MHDTests/FastRotor/makeICnew2.py deleted file mode 100644 index 26ae4bac1caebdb7b737af3405253e64f5112e95..0000000000000000000000000000000000000000 --- a/examples/MHDTests/FastRotor/makeICnew2.py +++ /dev/null @@ -1,80 +0,0 @@ -############################# -# This file is part of SWIFT -# Copyright -############################# - -import h5py -import numpy as np -import matplotlib.pyplot as plt - -# Parameters - -r_in = 0.1 -rho_in_0 = 10.0 -rho_out_0 = 1.0 -P_0 = 1.0 -B_0 = 2.5 / np.sqrt(np.pi) -omega_0 = 20.0 -gamma = 1.4 - -fileInputName = "FastRotor_temp_0000.hdf5" -fileOutputName = "FastRotor.hdf5" - -###---------------------------### - -particles = h5py.File(fileInputName, "r") - -[lx, ly, lz] = particles["/Header"].attrs.get("BoxSize") -N = particles["/Header"].attrs.get("NumPart_Total")[0] -print(N) -N = int(N) - -rho = particles["/PartType0/Densities"][:] -pos = particles["/PartType0/Coordinates"][:, :] -v = particles["/PartType0/Velocities"][:, :] -B = particles["/PartType0/MagneticFluxDensities"][:, :] -ids = particles["/PartType0/ParticleIDs"][:] -m = particles["/PartType0/Masses"][:] -u = P_0 / (rho * (gamma - 1)) -h = particles["/PartType0/SmoothingLengths"][:] - - -###---------------------------### - -# File -fileOutput = h5py.File(fileOutputName, "w") - -# Header -grp = fileOutput.create_group("/Header") -grp.attrs["BoxSize"] = [lx, ly, lz] ##### -grp.attrs["NumPart_Total"] = [N, 0, 0, 0, 0, 0] -grp.attrs["NumPart_Total_HighWord"] = [0, 0, 0, 0, 0, 0] -grp.attrs["NumPart_ThisFile"] = [N, 0, 0, 0, 0, 0] -grp.attrs["Time"] = 0.0 -grp.attrs["NumFileOutputsPerSnapshot"] = 1 -grp.attrs["MassTable"] = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0] -grp.attrs["Flag_Entropy_ICs"] = [0, 0, 0, 0, 0, 0] -grp.attrs["Dimension"] = 3 - -# Units -grp = fileOutput.create_group("/Units") -grp.attrs["Unit length in cgs (U_L)"] = 1.0 -grp.attrs["Unit mass in cgs (U_M)"] = 1.0 -grp.attrs["Unit time in cgs (U_t)"] = 1.0 -grp.attrs["Unit current in cgs (U_I)"] = 1.0 -grp.attrs["Unit temperature in cgs (U_T)"] = 1.0 - -# Particle group -grp = fileOutput.create_group("/PartType0") -grp.create_dataset("Coordinates", data=pos, dtype="d") -grp.create_dataset("Velocities", data=v, dtype="f") -grp.create_dataset("Masses", data=m, dtype="f") -grp.create_dataset("SmoothingLength", data=h, dtype="f") -grp.create_dataset("InternalEnergy", data=u, dtype="f") -grp.create_dataset("ParticleIDs", data=ids, dtype="L") -grp.create_dataset("MagneticFluxDensities", data=B, dtype="f") -# grp.create_dataset("VecPot", data = vp, dtype = 'f') -# grp.create_dataset("EPalpha", data = epa, dtype = 'f') -# grp.create_dataset("EPbeta" , data = epb, dtype = 'f') - -fileOutput.close() diff --git a/examples/MHDTests/FastRotor/make_gif.py b/examples/MHDTests/FastRotor/make_gif.py deleted file mode 100644 index b592bd754292855aded012cda82a81adc8052250..0000000000000000000000000000000000000000 --- a/examples/MHDTests/FastRotor/make_gif.py +++ /dev/null @@ -1,12 +0,0 @@ -import imageio - -filename_base = "FastRotor_LR_" - -images = [] - -for ii in range(61): - - filename = filename_base + str(ii).zfill(4) + ".png" - images.append(imageio.imread(filename)) - -imageio.mimsave("FastRotor_LR.gif", images) diff --git a/examples/MHDTests/FastRotor/plot.py b/examples/MHDTests/FastRotor/plot.py deleted file mode 100644 index cf0b0b05de4d02afc74cc9dbac7dd1b96a0ab2ec..0000000000000000000000000000000000000000 --- a/examples/MHDTests/FastRotor/plot.py +++ /dev/null @@ -1,164 +0,0 @@ -from swiftsimio import load -from swiftsimio.visualisation.slice import slice_gas -import numpy as np -import sys -import matplotlib.pyplot as plt - -filename = sys.argv[1] -data = load(filename) - -print(data.metadata.gas_properties.field_names) - -# pos = data.gas.coordinates - -# First create a mass-weighted temperature dataset -B = data.gas.magnetic_flux_densities -P_mag = (B[:, 0] ** 2 + B[:, 1] ** 2 + B[:, 2] ** 2) / 2 -print(P_mag) -print(data.gas.masses) -# quit() -data.gas.mass_weighted_magnetic_pressures = data.gas.masses * P_mag - -# Then create a mass-weighted pressure dataset -data.gas.mass_weighted_pressures = data.gas.masses * data.gas.pressures - -# Then create a mass-weighted speed dataset -v = data.gas.velocities -data.gas.mass_weighted_speeds_x = data.gas.masses * v[:, 0] -data.gas.mass_weighted_speeds_y = data.gas.masses * v[:, 1] -data.gas.mass_weighted_speeds = data.gas.masses * np.sqrt( - v[:, 0] ** 2 + v[:, 1] ** 2 + v[:, 2] ** 2 -) - -# Then create a mass-weighted pressure dataset -data.gas.mass_weighted_densities = data.gas.masses * data.gas.densities - -# Map in mass per area -mass_map = slice_gas( - data, - z_slice=0.5 * data.metadata.boxsize[2], - resolution=1024, - project="masses", - parallel=True, -) - -# Map in density per area -rho_pa_map = slice_gas( - data, - z_slice=0.5 * data.metadata.boxsize[2], - resolution=1024, - project="mass_weighted_densities", - parallel=True, -) - -# Map in magnetism squared times mass per area -mass_weighted_magnetic_pressure_map = slice_gas( - data, - z_slice=0.5 * data.metadata.boxsize[2], - resolution=1024, - project="mass_weighted_magnetic_pressures", - parallel=True, -) - -# Map in pressure times mass per area -mass_weighted_pressure_map = slice_gas( - data, - z_slice=0.5 * data.metadata.boxsize[2], - resolution=1024, - project="mass_weighted_pressures", - parallel=True, -) - -# Map in speed squared times mass per area -mass_weighted_speeds_map = slice_gas( - data, - z_slice=0.5 * data.metadata.boxsize[2], - resolution=1024, - project="mass_weighted_speeds", - parallel=True, -) - -# Map in speed squared times mass per area -mass_weighted_speeds_x_map = slice_gas( - data, - z_slice=0.5 * data.metadata.boxsize[2], - resolution=1024, - project="mass_weighted_speeds_x", - parallel=True, -) - -# Map in speed squared times mass per area -mass_weighted_speeds_y_map = slice_gas( - data, - z_slice=0.5 * data.metadata.boxsize[2], - resolution=1024, - project="mass_weighted_speeds_y", - parallel=True, -) - -magnetic_pressure_map = mass_weighted_magnetic_pressure_map / mass_map - -print(magnetic_pressure_map) - -speed_x_map = mass_weighted_speeds_x_map / mass_map -speed_y_map = mass_weighted_speeds_y_map / mass_map -speed_map = mass_weighted_speeds_map / mass_map -rho_map = rho_pa_map / mass_map -pressure_map = mass_weighted_pressure_map / mass_map - -from matplotlib.pyplot import imsave - -# from matplotlib.colors import LogNorm - -# Normalize and save -# imsave(sys.argv[2], np.rot90(magnetic_pressure_map.value), cmap="jet", vmin=0.0177, vmax=2.642) -# imsave(sys.argv[3], speed_map.value, cmap="viridis") - - -levels = 29 - -fig, axs = plt.subplots(2, 2, figsize=(10, 10)) - -im = axs[0, 0].contour( - np.rot90(rho_map.value), levels, colors="k", linewidths=0.5, vmin=0.483, vmax=12.95 -) -# plt.colorbar(im, ax = axs[0,0]) -im = axs[0, 1].contour( - np.rot90(pressure_map.value), - levels, - colors="k", - linewidths=0.5, - vmin=0.0202, - vmax=2.008, -) -# plt.colorbar(im, ax = axs[0,1]) -im = axs[1, 0].contour( - np.rot90(speed_map.value), levels, colors="k", linewidths=0.5, vmin=0.0, vmax=2.6 -) # , vmin=0., vmax=8.18) -# cb = plt.colorbar(im, ax = axs[1,0]) -im = axs[1, 1].contour( - np.rot90(magnetic_pressure_map.value), - levels, - colors="k", - linewidths=0.5, - vmin=0.0177, - vmax=2.642, -) -# plt.colorbar(im, ax = axs[1,1]) - -""" - -im = axs[0,0].contourf(np.rot90(rho_map.value), levels, cmap='viridis') #, vmin=0.483, vmax=12.95) -plt.colorbar(im, ax = axs[0,0]) -im = axs[0,1].contourf(np.rot90(pressure_map.value), levels, cmap='viridis') #, vmin=0.0202, vmax=2.008) -plt.colorbar(im, ax = axs[0,1]) -im = axs[1,0].contourf(np.rot90(speed_map.value), levels, cmap='viridis') #, vmin=0., vmax=2.6) #, vmin=0., vmax=8.18) -cb = plt.colorbar(im, ax = axs[1,0]) -im = axs[1,1].contourf(np.rot90(magnetic_pressure_map.value), levels, cmap='viridis') #, vmin=0.0177, vmax=2.642) -plt.colorbar(im, ax = axs[1,1]) - -""" - -plt.setp(plt.gcf().get_axes(), xticks=[], yticks=[]) - -plt.savefig(sys.argv[2]) diff --git a/examples/MHDTests/FastRotor/plot_all.py b/examples/MHDTests/FastRotor/plot_all.py index 7e94e06a0a8efa2c04482b13d50fbb9ef3368fbb..62466c28968ac376393d8280225e9171990572f1 100644 --- a/examples/MHDTests/FastRotor/plot_all.py +++ b/examples/MHDTests/FastRotor/plot_all.py @@ -2,141 +2,209 @@ from swiftsimio import load from swiftsimio.visualisation.slice import slice_gas import numpy as np import sys +import os +import matplotlib +#matplotlib.use("pdf") import matplotlib.pyplot as plt +from matplotlib.pyplot import imsave +from matplotlib.colors import LogNorm +from matplotlib.colors import Normalize +from mpl_toolkits.axes_grid1 import make_axes_locatable + +def set_colorbar(ax, im): + """ + Adapt the colorbar a bit for axis object <ax> and + imshow instance <im> + """ + divider = make_axes_locatable(ax) + cax = divider.append_axes("right", size="5%", pad=0.01) + plt.colorbar(im, cax=cax) + return filename_base = "FastRotor_LR_" -for ii in range(61): +nini=int(sys.argv[1]) +nfin=int(sys.argv[2]) +#for ii in range(61): +for ii in range(nini,nfin): print(ii) filename = filename_base + str(ii).zfill(4) + ".hdf5" data = load(filename) - + #print(data.metadata.gas_properties.field_names) + boxsize = data.metadata.boxsize + extent = [0, 1.0 , 0, 1.0] + # cut the domian in half + data.metadata.boxsize*=[1.0,0.5,1.0] + + gas_gamma = data.metadata.gas_gamma + print("Gas Gamma:",gas_gamma) + if gas_gamma != 7./5.: + print("WRONG GAS GAMMA") + exit() + + mhdflavour = data.metadata.hydro_scheme["MHD Flavour"] + mhd_scheme = data.metadata.hydro_scheme["MHD Scheme"] + mhdeta = data.metadata.hydro_scheme["Resistive Eta"] + git = data.metadata.code["Git Revision"] + gitBranch = data.metadata.code["Git Branch"] + scheme = data.metadata.hydro_scheme["Scheme"] + kernel = data.metadata.hydro_scheme["Kernel function"] + neighbours = data.metadata.hydro_scheme["Kernel target N_ngb"] + + try: + dedhyp = data.metadata.hydro_scheme["Dedner Hyperbolic Constant"] + dedpar = data.metadata.hydro_scheme["Dedner Parabolic Constant"] + except: + dedhyp = 0.0 + dedpar = 0.0 + + try: + deddivV = data.metadata.hydro_scheme["Dedner Hyperbolic div(v) Constant"] + tensile = data.metadata.hydro_scheme["MHD Tensile Instability Correction Prefactor"] + artdiff = data.metadata.hydro_scheme["Artificial Diffusion Constant"] + except: + deddivV = 0.0 + artdiff = 0.0 + tensile = 1.0 + # First create a mass-weighted temperature dataset + B = data.gas.magnetic_flux_densities + divB = data.gas.magnetic_divergences P_mag = (B[:, 0] ** 2 + B[:, 1] ** 2 + B[:, 2] ** 2) / 2 - + h = data.gas.smoothing_lengths + + normB = np.sqrt(B[:, 0] ** 2 + B[:, 1] ** 2 + B[:, 2] ** 2) + + data.gas.DivB_error = np.maximum(h * abs(divB) / normB, 1e-10) + + # Then create a mass-weighted B error dataset + data.gas.mass_weighted_magnetic_divB_error = data.gas.masses * data.gas.DivB_error + + # Then create a mass-weighted B pressure dataset data.gas.mass_weighted_magnetic_pressures = data.gas.masses * P_mag # Then create a mass-weighted pressure dataset data.gas.mass_weighted_pressures = data.gas.masses * data.gas.pressures + # Then create a mass-weighted Plasma Beta dataset + data.gas.plasma_beta = data.gas.pressures / P_mag + # Then create a mass-weighted speed dataset v = data.gas.velocities - data.gas.mass_weighted_speeds = data.gas.masses * np.sqrt( + data.gas.mass_weighted_speeds = data.gas.masses * ( v[:, 0] ** 2 + v[:, 1] ** 2 + v[:, 2] ** 2 ) - - # Then create a mass-weighted pressure dataset + # Then create a mass-weighted densities dataset data.gas.mass_weighted_densities = data.gas.masses * data.gas.densities # Map in mass per area - mass_map = slice_gas( - data, - z_slice=0.5 * data.metadata.boxsize[2], - resolution=1024, - project="masses", - parallel=True, - ) - + mass_map = slice_gas(data, z_slice=0.5 * data.metadata.boxsize[2], resolution=1024, project="masses", parallel=True ) + # Map in density per area - rho_pa_map = slice_gas( - data, - z_slice=0.5 * data.metadata.boxsize[2], - resolution=1024, - project="mass_weighted_densities", - parallel=True, - ) - + rho_map = slice_gas( + data, z_slice=0.5 * data.metadata.boxsize[2], resolution=1024, + project="mass_weighted_densities", parallel=True) # Map in magnetism squared times mass per area - mass_weighted_magnetic_pressure_map = slice_gas( - data, - z_slice=0.5 * data.metadata.boxsize[2], - resolution=1024, - project="mass_weighted_magnetic_pressures", - parallel=True, - ) + mw_magnetic_pressure_map = slice_gas( + data, z_slice=0.5 * data.metadata.boxsize[2], resolution=1024, + project="mass_weighted_magnetic_pressures", parallel=True) # Map in pressure times mass per area - mass_weighted_pressure_map = slice_gas( - data, - z_slice=0.5 * data.metadata.boxsize[2], - resolution=1024, - project="mass_weighted_pressures", - parallel=True, - ) + mw_pressure_map = slice_gas( + data, z_slice=0.5 * data.metadata.boxsize[2], resolution=1024, + project="mass_weighted_pressures", parallel=True, ) # Map in speed squared times mass per area - mass_weighted_speeds_map = slice_gas( - data, - z_slice=0.5 * data.metadata.boxsize[2], - resolution=1024, - project="mass_weighted_speeds", - parallel=True, + mw_speeds_map = slice_gas( + data, z_slice=0.5 * data.metadata.boxsize[2], resolution=1024, + project="mass_weighted_speeds", parallel=True, ) + + # Map in divB error times mass per area + mw_ErrDivB_map = slice_gas( + data, z_slice=0.5 * data.metadata.boxsize[2], resolution=1024, + project="mass_weighted_magnetic_divB_error", parallel=True, ) + + # Map in Plasma Beta times mass per area + plasma_beta_map = slice_gas( + data, z_slice=0.5 * data.metadata.boxsize[2], resolution=1024, + project="plasma_beta", parallel=True, ) + + rho_map = rho_map / mass_map + magnetic_pressure_map = mw_magnetic_pressure_map / mass_map + speed_map = mw_speeds_map / mass_map + pressure_map = mw_pressure_map / mass_map + ErrDivB_map = mw_ErrDivB_map / mass_map + #plasma_beta_map + + #fig = plt.figure(figsize=(12, 11), dpi=100) + fig = plt.figure(figsize=(12, 8), dpi=100) + + ax1 = fig.add_subplot(231) + im1 = ax1.imshow(rho_map.T, origin="lower", extent=extent, cmap="inferno", norm=LogNorm(vmax=14,vmin=1)) + ax1.set_title("Density") + set_colorbar(ax1, im1) + + ax2 = fig.add_subplot(232) + im2 = ax2.imshow(magnetic_pressure_map.T, origin="lower", extent=extent, cmap="plasma", norm=Normalize(vmax=3,vmin=0.1)) + ax2.set_title("Magnetic Pressure") + set_colorbar(ax2, im2) + + ax3 = fig.add_subplot(233) + im3 = ax3.imshow(speed_map.T, origin="lower", extent=extent, cmap="cividis", norm=Normalize(vmax=1.6,vmin=0.1)) + ax3.set_title("v^2") + set_colorbar(ax3, im3) + + ax4 = fig.add_subplot(234) + im4 = ax4.imshow(pressure_map.T, origin="lower", extent=extent, cmap="viridis", norm=Normalize(vmax=2.1,vmin=0.1)) + ax4.set_title("Internal Pressure") + set_colorbar(ax4, im4) + + #ax5 = fig.add_subplot(235) + #im5 = ax5.imshow(plasma_beta_map.T, origin="lower", extent=extent, cmap="magma", norm=LogNorm()) + #ax5.set_title("plasma beta") + #set_colorbar(ax5, im5) + + ax5 = fig.add_subplot(235) + im5 = ax5.imshow(ErrDivB_map.T, origin="lower", extent=extent, cmap="gray", norm=LogNorm(vmax=1,vmin=0.001)) + ax5.set_title("Err(DivB)") + set_colorbar(ax5, im5) + + for ax in [ax1, ax2, ax3, ax4, ax5]: + ax.set_xlabel("x ") + ax.set_ylabel("y ") + + ax6 = fig.add_subplot(236) + + text_fontsize = 8 + ax6.text( + 0.1, + 0.9, + "Fast Rotor $t=%.2f$" % data.metadata.time, + fontsize=text_fontsize, ) - - magnetic_pressure_map = mass_weighted_magnetic_pressure_map / mass_map - - speed_map = mass_weighted_speeds_map / mass_map - rho_map = rho_pa_map / mass_map - pressure_map = mass_weighted_pressure_map / mass_map - - from matplotlib.pyplot import imsave - - levels = 29 - - fig, axs = plt.subplots(2, 2, figsize=(10, 10)) - - im = axs[0, 0].contour( - np.rot90(rho_map.value), - levels, - colors="k", - linewidths=0.5, - vmin=0.483, - vmax=12.95, - ) - - im = axs[0, 1].contour( - np.rot90(pressure_map.value), - levels, - colors="k", - linewidths=0.5, - vmin=0.0202, - vmax=2.008, + ax6.text(0.1, 0.85, "$SWIFT$ %s" % git.decode("utf-8"), fontsize=text_fontsize) + ax6.text(0.1, 0.8, "$Branch$ %s" % gitBranch.decode("utf-8"), fontsize=text_fontsize) + ax6.text(0.1, 0.75, scheme.decode("utf-8"), fontsize=text_fontsize) + ax6.text(0.1, 0.7, kernel.decode("utf-8"), fontsize=text_fontsize) + ax6.text(0.1, 0.65, "$%.2f$ neighbours" % (neighbours), fontsize=text_fontsize) + ax6.text( + 0.1, + 0.55, + "$Flavour: $ %s" % mhdflavour.decode("utf-8")[0:25], + fontsize=text_fontsize, ) - - im = axs[1, 0].contour( - np.rot90(speed_map.value), - levels, - colors="k", - linewidths=0.5, - vmin=0.0, - vmax=2.6, - ) # , vmin=0., vmax=8.18) - - im = axs[1, 1].contour( - np.rot90(magnetic_pressure_map.value), - levels, - colors="k", - linewidths=0.5, - vmin=0.0177, - vmax=2.642, - ) - - """ - - im = axs[0,0].contourf(np.rot90(rho_map.value), levels, cmap='viridis') #, vmin=0.483, vmax=12.95) - plt.colorbar(im, ax = axs[0,0]) - im = axs[0,1].contourf(np.rot90(pressure_map.value), levels, cmap='viridis') #, vmin=0.0202, vmax=2.008) - plt.colorbar(im, ax = axs[0,1]) - im = axs[1,0].contourf(np.rot90(speed_map.value), levels, cmap='viridis') #, vmin=0., vmax=2.6) #, vmin=0., vmax=8.18) - cb = plt.colorbar(im, ax = axs[1,0]) - im = axs[1,1].contourf(np.rot90(magnetic_pressure_map.value), levels, cmap='viridis') #, vmin=0.0177, vmax=2.642) - plt.colorbar(im, ax = axs[1,1]) - - """ - - plt.setp(plt.gcf().get_axes(), xticks=[], yticks=[]) - - plt.savefig(filename_base + str(ii).zfill(4) + ".png") + ax6.text(0.1, 0.5, "$Resitivity_\\eta:%.4f$ " % (mhdeta), fontsize=text_fontsize) + ax6.text(0.1, 0.45, "$Dedner Parameters: $", fontsize=text_fontsize) + ax6.text(0.1, 0.4, "$[hyp, par, div] [:%.3f,%.3f,%.3f]$ " % (dedhyp,dedpar,deddivV), fontsize=text_fontsize) + ax6.text(0.1, 0.35, "$Tensile Prefactor:%.4f$ " % (tensile), fontsize=text_fontsize) + ax6.text(0.1, 0.3, "$Art. Diffusion:%.4f$ " % (artdiff), fontsize=text_fontsize) + ax6.tick_params(left=False, right=False, labelleft=False, labelbottom=False, bottom=False) + #ax6.set_xlabel("") + #ax6.plot(frameon=False) + + plt.tight_layout() + plt.savefig(filename_base + str(ii).zfill(4) + ".jpg") + plt.close() diff --git a/examples/MHDTests/FastRotor/plot_slice.py b/examples/MHDTests/FastRotor/plot_slice.py deleted file mode 100644 index e42157d503d45a8c4fc4b733bcadcaaf0fab273a..0000000000000000000000000000000000000000 --- a/examples/MHDTests/FastRotor/plot_slice.py +++ /dev/null @@ -1,55 +0,0 @@ -from swiftsimio import load -from swiftsimio.visualisation.slice import slice_gas -import numpy as np -import sys - -filename = sys.argv[1] -data = load(filename) - -# First create a mass-weighted temperature dataset -B = data.gas.magnetic_flux_densities -data.gas.mass_weighted_magnetic_pressures = ( - data.gas.masses * (B[:, 0] ** 2 + B[:, 1] ** 2 + B[:, 2] ** 2) / 2 -) - -# Then create a mass-weighted speed dataset -v = data.gas.velocities -data.gas.mass_weighted_speeds = data.gas.masses * ( - v[:, 0] ** 2 + v[:, 1] ** 2 + v[:, 2] ** 2 -) - -# Map in mass per area -mass_map = slice_gas( - data, - z_slice=0.5 * data.metadata.boxsize[2], - resolution=256, - project="masses", - parallel=True, -) - -# Map in magnetism squared times mass per area -mass_weighted_magnetic_pressure_map = slice_gas( - data, - z_slice=0.5 * data.metadata.boxsize[2], - resolution=256, - project="mass_weighted_magnetic_pressures", - parallel=True, -) -# Map in speed squared times mass per area -mass_weighted_speeds_map = slice_gas( - data, - z_slice=0.5 * data.metadata.boxsize[2], - resolution=256, - project="mass_weighted_speeds", - parallel=True, -) - -magnetic_pressure_map = mass_weighted_magnetic_pressure_map / mass_map -speed_map = mass_weighted_speeds_map / mass_map - -from matplotlib.pyplot import imsave -from matplotlib.colors import LogNorm - -# Normalize and save -imsave(sys.argv[2], magnetic_pressure_map.value, cmap="viridis") -imsave(sys.argv[3], speed_map.value, cmap="viridis") diff --git a/examples/MHDTests/FastRotor/plot_slice_new.py b/examples/MHDTests/FastRotor/plot_slice_new.py deleted file mode 100644 index 76ba090c9a3494a60c880669492aea9c52e97a19..0000000000000000000000000000000000000000 --- a/examples/MHDTests/FastRotor/plot_slice_new.py +++ /dev/null @@ -1,113 +0,0 @@ -from swiftsimio import load -from swiftsimio.visualisation.slice import slice_gas -import numpy as np -import sys -import matplotlib.pyplot as plt - -filename = sys.argv[1] -data = load(filename) - -print(data.metadata.gas_properties.field_names) - -# pos = data.gas.coordinates - -# First create a mass-weighted temperature dataset -B = data.gas.magnetic_flux_densities -P_mag = (B[:, 0] ** 2 + B[:, 1] ** 2 + B[:, 2] ** 2) / 2 -print(P_mag) -print(data.gas.masses) -# quit() -data.gas.mass_weighted_magnetic_pressures = data.gas.masses * P_mag - -# Then create a mass-weighted pressure dataset -data.gas.mass_weighted_pressures = data.gas.masses * data.gas.pressures - -# Then create a mass-weighted speed dataset -v = data.gas.velocities -data.gas.mass_weighted_speeds_x = data.gas.masses * v[:, 0] -data.gas.mass_weighted_speeds_y = data.gas.masses * v[:, 1] -data.gas.mass_weighted_speeds = data.gas.masses * np.sqrt( - v[:, 0] ** 2 + v[:, 1] ** 2 + v[:, 2] ** 2 -) - -# Then create a mass-weighted pressure dataset -data.gas.mass_weighted_densities = data.gas.masses * data.gas.densities - -# Map in mass per area -mass_map = slice_gas(data, resolution=1024, project="masses", parallel=True) - -# Map in density per area -rho_pa_map = slice_gas( - data, resolution=1024, project="mass_weighted_densities", parallel=True -) - -# Map in magnetism squared times mass per area -mass_weighted_magnetic_pressure_map = slice_gas( - data, resolution=1024, project="mass_weighted_magnetic_pressures", parallel=True -) - -# Map in pressure times mass per area -mass_weighted_pressure_map = slice_gas( - data, resolution=1024, project="mass_weighted_pressures", parallel=True -) - -# Map in speed squared times mass per area -mass_weighted_speeds_map = slice_gas( - data, resolution=1024, project="mass_weighted_speeds", parallel=True -) - -# Map in speed squared times mass per area -mass_weighted_speeds_x_map = slice_gas( - data, resolution=1024, project="mass_weighted_speeds_x", parallel=True -) - -# Map in speed squared times mass per area -mass_weighted_speeds_y_map = slice_gas( - data, resolution=1024, project="mass_weighted_speeds_y", parallel=True -) - -magnetic_pressure_map = mass_weighted_magnetic_pressure_map / mass_map - -print(magnetic_pressure_map) - -speed_x_map = mass_weighted_speeds_x_map / mass_map -speed_y_map = mass_weighted_speeds_y_map / mass_map -speed_map = mass_weighted_speeds_map / mass_map -rho_map = rho_pa_map / mass_map -pressure_map = mass_weighted_pressure_map / mass_map - -from matplotlib.pyplot import imsave - -# from matplotlib.colors import LogNorm - -# Normalize and save -# imsave(sys.argv[2], np.rot90(magnetic_pressure_map.value), cmap="jet", vmin=0.0177, vmax=2.642) -# imsave(sys.argv[3], speed_map.value, cmap="viridis") - -levels = 29 - -fig, axs = plt.subplots(2, 2, figsize=(12, 10)) -im = axs[0, 0].contour( - np.rot90(rho_map.value), levels, colors="k", vmin=0.483, vmax=12.95 -) -plt.colorbar(im, ax=axs[0, 0]) -im = axs[0, 1].contour( - np.rot90(pressure_map.value), levels, colors="k", vmin=0.0202, vmax=2.008 -) -plt.colorbar(im, ax=axs[0, 1]) -im = axs[1, 0].contour( - np.rot90(speed_map.value), levels, colors="k", vmin=0.0, vmax=2.6 -) # , vmin=0., vmax=8.18) -cb = plt.colorbar(im, ax=axs[1, 0]) -# vf = axs[1,0].quiver(X, Y, speed_x_map.value, speed_y_map.value, colors='k', scale = 100) - -print(np.rot90(magnetic_pressure_map.value)) - -im = axs[1, 1].contour( - np.rot90(magnetic_pressure_map.value), levels, colors="k", vmin=0.0177, vmax=2.642 -) -plt.colorbar(im, ax=axs[1, 1]) - -# plt.setp(plt.gcf().get_axes(), xticks=[], yticks=[]); - -plt.savefig(sys.argv[2]) diff --git a/examples/MHDTests/FastRotor/run.sh b/examples/MHDTests/FastRotor/run.sh new file mode 100755 index 0000000000000000000000000000000000000000..6a1d6c13d09c7c0c8c6cbc19070a8eb83246cf40 --- /dev/null +++ b/examples/MHDTests/FastRotor/run.sh @@ -0,0 +1,16 @@ +#!/bin/bash + +# Generate the initial conditions if they are not present. +if [ ! -e FastRotor_LR.hdf5 ] +then + echo "Fetching Glass Files..." + ./getGlass.sh + echo "Generating the ICs" + python ./makeIC_VP.py +fi + +# Run SWIFT +../../../swift --hydro --threads=16 ../FastRotor.yml 2>&1 > out.log + +# Plot the temperature evolution +python3 ./plot_all.py 0 60 diff --git a/examples/MHDTests/FastRotor/runSchemes.sh b/examples/MHDTests/FastRotor/runSchemes.sh new file mode 100755 index 0000000000000000000000000000000000000000..d4e99fb498f6bebebd76278b6bb1ca7d3a08cdbc --- /dev/null +++ b/examples/MHDTests/FastRotor/runSchemes.sh @@ -0,0 +1,93 @@ +#!/bin/bash + +# Generate the initial conditions if they are not present. +if [ ! -e FastRotor_LR.hdf5 ] +then + echo "Fetching Glass Files..." + ./getGlass.sh + echo "Generating the ICs" + python ./makeIC.py +fi + +SCHEME_ID=("VeP" "ODI" "FDI") +SCHEME_IDD=("vp" "odi" "fdi") +#################### +case $# in + 1) + echo "Using Dirs as TEST" + DIRS="TEST" + WHAT=$1 + ;; + 2) + DIRS=$2 + WHAT=$1 + echo "Using Dirs as $DIRS" + ;; + *) + echo "Usage $0 [which] [DIRECTORY]" + echo "[what scheme]:" + echo "vep: vector potentials" + echo "odi: Oresti's direct induction" + echo "fdi: simple direct induction" + echo "[FOLDER_TAIL]:" + echo "trailer naming of folders" + echo "" + exit + ;; +esac +##################### +case $WHAT in + vep) + SCHEME_Nr=0 + ;; + odi) + SCHEME_Nr=1 + ;; + fdi) + SCHEME_Nr=2 + ;; + all) + SCHEME_Nr=( 0 1 2 ) + ;; + *) + echo $WHAT" wrong scheme" + exit 2 + ;; +esac + +SCHEME_DIRS=("VeP_$DIRS" "ODI_$DIRS" "FDI_$DIRS") + +for J in ${SCHEME_Nr[@]} +do + echo $J + ID=${SCHEME_ID[$J]} + IDD=${SCHEME_IDD[$J]} + DIR=${SCHEME_DIRS[$J]} + if [ ! -e $DIR ] + then + echo "Folder $DIR does not exist" + mkdir $DIR + fi + cd $DIR + if [ ! -e sw_$ID ] + then + cur_dir=`pwd` + cd ../../../../ + pwd + ./TestAllMHD.sh $IDD "--with-adiabatic-index=7/5" + cd $cur_dir + cp ../../../../sw_$ID . + fi + cat <<-EOF > ./run.sh + #!/bin/bash + # Run SWIFT + ./sw_$ID --hydro --threads=16 ../FR_schemes.yml 2>&1 > out.log + + # Plot the evolution + python3 ../plot_all.py 0 61 2>&1 > plot.log + EOF + chmod u+x ./run.sh + ./run.sh & + cd .. +done + diff --git a/examples/MHDTests/MagneticBlastWave/BW_schemes.yml b/examples/MHDTests/MagneticBlastWave/BW_schemes.yml new file mode 100644 index 0000000000000000000000000000000000000000..0745d7ea17f52db59a64d7e462526877bb34d0d0 --- /dev/null +++ b/examples/MHDTests/MagneticBlastWave/BW_schemes.yml @@ -0,0 +1,52 @@ +# Define the system of units to use internally. +InternalUnitSystem: + UnitMass_in_cgs: 1 # Grams + UnitLength_in_cgs: 1 # Centimeters + UnitVelocity_in_cgs: 1 # Centimeters per second + UnitCurrent_in_cgs: 1 # Amperes + UnitTemp_in_cgs: 1 # Kelvin + +# Parameters governing the time integration +TimeIntegration: + time_begin: 0. # The starting time of the simulation (in internal units). + time_end: 0.4 # The end time of the simulation (in internal units). + dt_min: 1e-9 # The minimal time-step size of the simulation (in internal units). + dt_max: 1e-3 # The maximal time-step size of the simulation (in internal units). + +# Parameters governing the snapshots +Snapshots: + basename: MagneticBlastWave_LR # Common part of the name of output files + time_first: 0. # Time of the first output (in internal units) + delta_time: 0.01 # Time difference between consecutive outputs (in internal units) + compression: 1 + +# Parameters governing the conserved quantities statistics +Statistics: + delta_time: 1e-2 # Time between statistics output + +# Parameters for the hydrodynamics scheme +SPH: + resolution_eta: 1.2348 # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel$ + CFL_condition: 0.075 # Courant-Friedrich-Levy condition for time integration. + viscosity_alpha: 2.0 + +# Parameters for MHD schemes +MHD: + hyperbolic_dedner: 1.0 + parabolic_dedner: 2.0 + hyperbolic_dedner_divv: 0.0 + monopole_subtraction: 1.0 + artificial_diffusion: 0.0 + resistive_eta: 0.01 + +# Physical parameters +PhysicalConstants: + mu_0: 1.0 + +# Parameters related to the initial conditions +InitialConditions: + file_name: ../MagneticBlastWave_LR.hdf5 # The file to read + periodic: 1 + +Scheduler: + max_top_level_cells: 256 diff --git a/examples/MHDTests/MagneticBlastWave/FDI/run.sh b/examples/MHDTests/MagneticBlastWave/FDI/run.sh new file mode 100755 index 0000000000000000000000000000000000000000..c7f6e839f617553ca022e4f0a865046785056161 --- /dev/null +++ b/examples/MHDTests/MagneticBlastWave/FDI/run.sh @@ -0,0 +1,5 @@ +# Run SWIFT +../../../../sw_FDI --hydro --threads=16 ../BW_schemes.yml 2>&1 > out.log + +# Plot the temperature evolution +python3 ../plot_all.py 0 60 > plot.log diff --git a/examples/MHDTests/MagneticBlastWave/ODI/run.sh b/examples/MHDTests/MagneticBlastWave/ODI/run.sh new file mode 100755 index 0000000000000000000000000000000000000000..1e67943dcd986435606636e0b32be7fd0979848b --- /dev/null +++ b/examples/MHDTests/MagneticBlastWave/ODI/run.sh @@ -0,0 +1,5 @@ +# Run SWIFT +../../../../sw_ODI --hydro --threads=16 ../BW_schemes.yml 2>&1 > out.log + +# Plot the temperature evolution +python3 ../plot_all.py 0 60 > plot.log diff --git a/examples/MHDTests/MagneticBlastWave/README.md b/examples/MHDTests/MagneticBlastWave/README.md new file mode 100644 index 0000000000000000000000000000000000000000..e9fecb8cd90602415fd8a7605f8ce10b3c2b6d2c --- /dev/null +++ b/examples/MHDTests/MagneticBlastWave/README.md @@ -0,0 +1,17 @@ +Blast wave folloing Seo & Ryu, 2023. +############################## +# Basically the folloing script compiles and run the test in different +# folders for each scheme. usefull for testing +./runSchemes.sh [what Scheme] [FOLDER_TAIL] +[what scheme]: +vep: vector potentials +odi: Oresti's direct induction +fdi: simple direct induction +[FOLDER_TAIL]: +the trailing name for the folders created +############################## + +TODO> +-check for symmetry with the between half of the domain +-make a cut +-upload reference runs or values diff --git a/examples/MHDTests/MagneticBlastWave/VeP/run.sh b/examples/MHDTests/MagneticBlastWave/VeP/run.sh new file mode 100755 index 0000000000000000000000000000000000000000..3e02d6a037b71312914f449da7f6691ec3428e92 --- /dev/null +++ b/examples/MHDTests/MagneticBlastWave/VeP/run.sh @@ -0,0 +1,5 @@ +# Run SWIFT +../../../../sw_VeP --hydro --threads=16 ../BW_schemes.yml 2>&1 > out.log + +# Plot the temperature evolution +python3 ../plot_all.py 0 60 > plot.log diff --git a/examples/MHDTests/MagneticBlastWave/getGlass.sh b/examples/MHDTests/MagneticBlastWave/getGlass.sh new file mode 100755 index 0000000000000000000000000000000000000000..ffd92e88deae6e91237059adac2a6c2067caee46 --- /dev/null +++ b/examples/MHDTests/MagneticBlastWave/getGlass.sh @@ -0,0 +1,2 @@ +#!/bin/bash +wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/glassCube_32.hdf5 diff --git a/examples/MHDTests/MagneticBlastWave/makeIC.py b/examples/MHDTests/MagneticBlastWave/makeIC.py index 1b24242f94b0746163f57edfb03cc8e01865a602..c9fa05b9c1ce82731e132f0e0f321c971dd22830 100644 --- a/examples/MHDTests/MagneticBlastWave/makeIC.py +++ b/examples/MHDTests/MagneticBlastWave/makeIC.py @@ -9,11 +9,11 @@ import matplotlib.pyplot as plt # Parameters -r_in = 0.1 +r_in = 0.125 rho_0 = 1.0 -P_in_0 = 10.0 -P_out_0 = 0.1 -B_0 = 1.0 / np.sqrt(2.0) +P_in_0 = 100.0 +P_out_0 = 1 +B_0 = 10.0 gamma = 5.0 / 3.0 fileOutputName = "MagneticBlastWave_LR.hdf5" @@ -34,9 +34,9 @@ cx = times cy = times cz = 1 -lx = 2.0 -ly = 2.0 -lz = 2.0 / float(times) +lx = 1.0 +ly = 1.0 +lz = 1.0 / float(times) lx_c = lx / 2 ly_c = ly / 2 @@ -70,28 +70,58 @@ vol = lx * ly * lz ###---------------------------### rot = np.sqrt((pos[:, 0] - lx_c) ** 2 + (pos[:, 1] - ly_c) ** 2) -v = np.zeros((N, 3)) -B = np.zeros((N, 3)) -ids = np.linspace(1, N, N) -m = np.ones(N) * rho_0 * vol / N +#v = np.zeros((N, 3)) +#B = np.zeros((N, 3)) +#ids = np.linspace(1, N, N) +#m = np.ones(N) * rho_0 * vol / N u = np.ones(N) u[rot < r_in] = P_in_0 / (rho_0 * (gamma - 1)) u[rot >= r_in] = P_out_0 / (rho_0 * (gamma - 1)) +#B[:, 0] = B_0 -B[:, 0] = B_0 -B[:, 1] = B_0 - +print("Max Pos: [",max(pos[:,0]),max(pos[:,1]),max(pos[:,2]),"]") +print("Min Pos: [",min(pos[:,0]),min(pos[:,1]),min(pos[:,2]),"]") ###---------------------------### +N2=int(2*N) +p=np.zeros((N2, 3)) +p[:N,0]=pos[:,0] +p[N:,0]=pos[:,0] +p[:N,1]=pos[:,1] +p[N:,1]=pos[:,1]+ly +p[:N,2]=pos[:,2] +p[N:,2]=pos[:,2] +pos=p +hh =np.zeros(N2) +hh[:N]=h +hh[N:]=h +h=hh +v=np.zeros((N2,3)) +ids = np.linspace(1, N2, N2) +m = np.ones(N2) * rho_0 * vol / N +uu =np.zeros(N2) +uu[:N]=u +uu[N:]=u +u=uu +B = np.zeros((N2, 3)) +A = np.zeros((N2, 3)) +B[:N, 0] = B_0 +B[N:, 0] = -B_0 +A[:N, 2] = B_0 * pos[:N,1] +A[N:, 2] = B_0 * (2*ly-pos[N:,1]) + +print("Max Pos: [",max(pos[:,0]),max(pos[:,1]),max(pos[:,2]),"]") +print("Min Pos: [",min(pos[:,0]),min(pos[:,1]),min(pos[:,2]),"]") # File + fileOutput = h5py.File(fileOutputName, "w") # Header grp = fileOutput.create_group("/Header") -grp.attrs["BoxSize"] = [lx, ly, lz] ##### -grp.attrs["NumPart_Total"] = [N, 0, 0, 0, 0, 0] +grp.attrs["BoxSize"] = [lx, 2. * ly, lz] ##### +grp.attrs["NumPart_Total"] = [N2, 0, 0, 0, 0, 0] grp.attrs["NumPart_Total_HighWord"] = [0, 0, 0, 0, 0, 0] -grp.attrs["NumPart_ThisFile"] = [N, 0, 0, 0, 0, 0] +grp.attrs["NumPart_ThisFile"] = [N2, 0, 0, 0, 0, 0] grp.attrs["Time"] = 0.0 grp.attrs["NumFileOutputsPerSnapshot"] = 1 grp.attrs["MassTable"] = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0] @@ -115,8 +145,6 @@ grp.create_dataset("SmoothingLength", data=h, dtype="f") grp.create_dataset("InternalEnergy", data=u, dtype="f") grp.create_dataset("ParticleIDs", data=ids, dtype="L") grp.create_dataset("MagneticFluxDensities", data=B, dtype="f") -# grp.create_dataset("VecPot", data = vp, dtype = 'f') -# grp.create_dataset("EPalpha", data = epa, dtype = 'f') -# grp.create_dataset("EPbeta" , data = epb, dtype = 'f') +grp.create_dataset("MagneticVectorPotentials", data = A, dtype = 'f') fileOutput.close() diff --git a/examples/MHDTests/MagneticBlastWave/plot_all.py b/examples/MHDTests/MagneticBlastWave/plot_all.py index 476e0f94b5dee2f7937db1dfbe2961e119dc97db..3705048ff1fc501f4d7373618cf04f200c3e69fe 100644 --- a/examples/MHDTests/MagneticBlastWave/plot_all.py +++ b/examples/MHDTests/MagneticBlastWave/plot_all.py @@ -1,37 +1,223 @@ from swiftsimio import load +from swiftsimio import mask as sw_mask from swiftsimio.visualisation.projection import project_gas import numpy as np import sys +import matplotlib.pyplot as plt +from matplotlib.pyplot import imsave +from matplotlib.colors import LogNorm +from matplotlib.colors import Normalize +from mpl_toolkits.axes_grid1 import make_axes_locatable -input_filename_base = "MagneticBlastWave_LR_" +def set_colorbar(ax, im): + """ + Adapt the colorbar a bit for axis object <ax> and + imshow instance <im> + """ + divider = make_axes_locatable(ax) + cax = divider.append_axes("right", size="5%", pad=0.01) + plt.colorbar(im, cax=cax) + return -for ii in range(41): +filename_base = "MagneticBlastWave_LR_" + +nini=int(sys.argv[1]) +nfin=int(sys.argv[2]) +#for ii in range(61): +for ii in range(nini,nfin): print(ii) - filename = input_filename_base + str(ii).zfill(4) + ".hdf5" + filename = filename_base + str(ii).zfill(4) + ".hdf5" + #mask=sw_mask(filename) + #bs=mask.metadata.boxsize + #print(bs) + #mask.constrain_spatial([[0.,1.0],[0.,1.0],[0.,0.5]]) + #mask.constrain_spatial([[0.*bs[0],bs[0]],[0.*bs[0], 0.5*bs[1]],[0.*bs[2],bs[2]]]) + #data = load(filename,mask=mask) data = load(filename) + #print(data.metadata.gas_properties.field_names) + #boxsize = data.metadata.boxsize + extent = [0.0, 1.0, 0.0, 1.0] + data.metadata.boxsize*=[1.0,0.5,1.0] + + gas_gamma = data.metadata.gas_gamma + print("Gas Gamma:",gas_gamma) + if gas_gamma != 5./3.: + print("WRONG GAS GAMMA") + exit() + + mhdflavour = data.metadata.hydro_scheme["MHD Flavour"] + mhd_scheme = data.metadata.hydro_scheme["MHD Scheme"] + mhdeta = data.metadata.hydro_scheme["Resistive Eta"] + git = data.metadata.code["Git Revision"] + gitBranch = data.metadata.code["Git Branch"] + scheme = data.metadata.hydro_scheme["Scheme"] + kernel = data.metadata.hydro_scheme["Kernel function"] + neighbours = data.metadata.hydro_scheme["Kernel target N_ngb"] + + try: + dedhyp = data.metadata.hydro_scheme["Dedner Hyperbolic Constant"] + dedpar = data.metadata.hydro_scheme["Dedner Parabolic Constant"] + except: + dedhyp = 0.0 + dedpar = 0.0 + try: + deddivV = data.metadata.hydro_scheme["Dedner Hyperbolic div(v) Constant"] + tensile = data.metadata.hydro_scheme["MHD Tensile Instability Correction Prefactor"] + artdiff = data.metadata.hydro_scheme["Artificial Diffusion Constant"] + except: + deddivV = 0.0 + artdiff = 0.0 + tensile = 1.0 + # First create a mass-weighted temperature dataset - rho = data.gas.densities - data.gas.mass_weighted_densities = data.gas.masses * rho + + B = data.gas.magnetic_flux_densities + divB = data.gas.magnetic_divergences + P_mag = (B[:, 0] ** 2 + B[:, 1] ** 2 + B[:, 2] ** 2) / 2 + h = data.gas.smoothing_lengths + #A = data.gas.magnetic_vector_potentials + + normB = np.sqrt(B[:, 0] ** 2 + B[:, 1] ** 2 + B[:, 2] ** 2) + + DivB_error = np.maximum(h * abs(divB) / normB, 1e-10) + + mmasses=data.gas.masses + + pressure=data.gas.pressures + + # Then create a mass-weighted B error dataset + data.gas.mass_weighted_magnetic_divB_error = mmasses * DivB_error + + # Then create a mass-weighted B pressure dataset + data.gas.mass_weighted_magnetic_pressures = mmasses * P_mag + + # Then create a mass-weighted pressure dataset + data.gas.mass_weighted_pressures = mmasses * data.gas.pressures + + # Then create a mass-weighted Plasma Beta dataset + data.gas.plasma_beta = pressure / P_mag + + # Then create a mass-weighted speed dataset + v = data.gas.velocities + data.gas.mass_weighted_speeds = data.gas.masses * ( + v[:, 0] ** 2 + v[:, 1] ** 2 + v[:, 2] ** 2 + ) + + # Then create a mass-weighted densities dataset + data.gas.mass_weighted_densities = data.gas.masses * data.gas.densities # Map in mass per area mass_map = project_gas(data, resolution=1024, project="masses", parallel=True) - # Map in magnetism squared times mass per area - mass_weighted_density_map = project_gas( + + # Map in density per area + mw_density_map = project_gas( data, resolution=1024, project="mass_weighted_densities", parallel=True ) + # Map in magnetism squared times mass per area + mw_magnetic_pressure_map = project_gas( + data, resolution=1024, project="mass_weighted_magnetic_pressures", parallel=True) - density_map = mass_weighted_density_map / mass_map + # Map in pressure times mass per area + mw_pressure_map = project_gas( + data, resolution=1024, project="mass_weighted_pressures", parallel=True, ) - from matplotlib.pyplot import imsave + # Map in speed squared times mass per area + mw_speeds_map = project_gas( + data, resolution=1024, project="mass_weighted_speeds", parallel=True, ) - # Normalize and save - imsave( - input_filename_base + str(ii).zfill(4) + ".png", - np.rot90(density_map.value), - cmap="jet", - vmin=0.1, - vmax=2.4, + # Map in divB error times mass per area + mw_ErrDivB_map = project_gas( + data, resolution=1024, project="mass_weighted_magnetic_divB_error", parallel=True, ) + + # Map in Plasma Beta times mass per area + plasma_beta_map = project_gas( + data, resolution=1024, project="plasma_beta", parallel=True, ) + + rho_map = mw_density_map / mass_map + magnetic_pressure_map = mw_magnetic_pressure_map / mass_map + speed_map = mw_speeds_map / mass_map + pressure_map = mw_pressure_map / mass_map + ErrDivB_map = np.maximum(mw_ErrDivB_map / mass_map,1E-10) + #plasma_beta_map + + #fig = plt.figure(figsize=(12, 11), dpi=100) + fig = plt.figure(figsize=(12, 8), dpi=100) + + ax1 = fig.add_subplot(231) + im1 = ax1.imshow(rho_map.T, origin="lower", extent=extent, cmap="inferno", norm=Normalize(vmax=3,vmin=0)) + ax1.set_title("Density") + set_colorbar(ax1, im1) + + ax2 = fig.add_subplot(232) + #im2 = ax2.imshow(magnetic_pressure_map.T, origin="lower", extent=extent, cmap="plasma", norm=Normalize(vmax=120,vmin=30)) + im2 = ax2.imshow(magnetic_pressure_map.T, origin="lower", extent=extent, cmap="plasma", norm=Normalize(vmax=65,vmin=25)) + ax2.set_title("Magnetic Pressure") + set_colorbar(ax2, im2) + + ax3 = fig.add_subplot(233) + im3 = ax3.imshow(speed_map.T, origin="lower", extent=extent, cmap="cividis", norm=Normalize(vmax=30,vmin=0.)) + ax3.set_title("Speed") + set_colorbar(ax3, im3) + + ax4 = fig.add_subplot(234) + im4 = ax4.imshow(pressure_map.T, origin="lower", extent=extent, cmap="viridis", norm=Normalize(vmax=50,vmin=0.)) + ax4.set_title("Internal Pressure") + set_colorbar(ax4, im4) + + #ax5 = fig.add_subplot(235) + #im5 = ax5.imshow(plasma_beta_map.T, origin="lower", extent=extent, cmap="magma", norm=LogNorm()) + #ax5.set_title("plasma beta") + #set_colorbar(ax5, im5) + + ax5 = fig.add_subplot(235) + im5 = ax5.imshow(ErrDivB_map.T, origin="lower", extent=extent, cmap="gray", norm=LogNorm(vmax=1,vmin=0.001)) + ax5.set_title("Err(DivB)") + set_colorbar(ax5, im5) + + for ax in [ax1, ax2, ax3, ax4, ax5]: + ax.set_xlabel("x ") + ax.set_ylabel("y ") + + ax6 = fig.add_subplot(236) + + text_fontsize = 8 + ax6.text( + 0.1, + 0.9, + "Blast Wave $t=%.2f$" % data.metadata.time, + fontsize=text_fontsize, ) + ax6.text(0.1, 0.85, "$SWIFT$ %s" % git.decode("utf-8"), fontsize=text_fontsize) + ax6.text(0.1, 0.8, "$Branch$ %s" % gitBranch.decode("utf-8"), fontsize=text_fontsize) + ax6.text(0.1, 0.75, scheme.decode("utf-8"), fontsize=text_fontsize) + ax6.text(0.1, 0.7, kernel.decode("utf-8"), fontsize=text_fontsize) + ax6.text(0.1, 0.65, "$%.2f$ neighbours" % (neighbours), fontsize=text_fontsize) + ax6.text( + 0.1, + 0.55, + "$Flavour: $ %s" % mhdflavour.decode("utf-8")[0:25], + fontsize=text_fontsize, + ) + ax6.text(0.1, 0.5, "$Resitivity_\\eta:%.4f$ " % (mhdeta), fontsize=text_fontsize) + ax6.text(0.1, 0.45, "$Dedner Parameters: $", fontsize=text_fontsize) + ax6.text(0.1, 0.4, "$[hyp, par, div] [:%.3f,%.3f,%.3f]$ " % (dedhyp,dedpar,deddivV), fontsize=text_fontsize) + ax6.text(0.1, 0.35, "$Tensile Prefactor:%.4f$ " % (tensile), fontsize=text_fontsize) + ax6.text(0.1, 0.3, "$Art. Diffusion:%.4f$ " % (artdiff), fontsize=text_fontsize) + ax6.tick_params(left=False, right=False, labelleft=False, labelbottom=False, bottom=False) + + plt.tight_layout() + plt.savefig(filename_base + str(ii).zfill(4) + ".jpg") + plt.close() + #from matplotlib.pyplot import imsave + + # Normalize and save + #imsave( + # input_filename_base + str(ii).zfill(4) + ".png", + # np.rot90(density_map.value), + # cmap="jet", + # vmin=0.1, + # vmax=2.4, + #) diff --git a/examples/MHDTests/MagneticBlastWave/run.sh b/examples/MHDTests/MagneticBlastWave/run.sh new file mode 100755 index 0000000000000000000000000000000000000000..dcea500070739294ddbae29d30d74e1a28034978 --- /dev/null +++ b/examples/MHDTests/MagneticBlastWave/run.sh @@ -0,0 +1,16 @@ +#!/bin/bash + +# Generate the initial conditions if they are not present. +if [ ! -e MagneticBlastWave_LR.hdf5 ] +then + echo "Fetching Glass Files..." + ./getGlass.sh + echo "Generating the ICs" + python ./makeIC.py +fi + +# Run SWIFT +../../../swift --hydro --threads=16 ../FastRotor.yml 2>&1 > out.log + +# Plot the temperature evolution +# python3 ./plot_all.py 0 60 diff --git a/examples/MHDTests/MagneticBlastWave/runSchemes.sh b/examples/MHDTests/MagneticBlastWave/runSchemes.sh new file mode 100755 index 0000000000000000000000000000000000000000..f7fc1a0f1035b4b1005945002e1a76b63e17c859 --- /dev/null +++ b/examples/MHDTests/MagneticBlastWave/runSchemes.sh @@ -0,0 +1,94 @@ +#!/bin/bash + +# Generate the initial conditions if they are not present. +if [ ! -e MagneticBlastWave_LR.hdf5 ] +then + echo "Fetching Glass Files..." + ./getGlass.sh + echo "Generating the ICs" + python ./makeIC.py +fi + +SCHEME_ID=("VeP" "ODI" "FDI") +SCHEME_IDD=("vp" "odi" "fdi") +#################### +case $# in + 1) + echo "Using Dirs as TEST" + DIRS="TEST" + WHAT=$1 + ;; + 2) + DIRS=$2 + WHAT=$1 + echo "Using Dirs as $DIRS" + ;; + *) + echo "Usage $0 [which] [DIRECTORY]" + echo "[what scheme]:" + echo "vep: vector potentials" + echo "odi: Oresti's direct induction" + echo "fdi: simple direct induction" + echo "[FOLDER_TAIL]:" + echo "trailer naming of folders" + echo "" + exit + ;; +esac +##################### +case $WHAT in + vep) + SCHEME_Nr=0 + ;; + odi) + SCHEME_Nr=1 + ;; + fdi) + SCHEME_Nr=2 + ;; + all) + SCHEME_Nr=( 0 1 2 ) + ;; + *) + echo $WHAT" wrong scheme" + exit 2 + ;; +esac + +SCHEME_DIRS=("VeP_$DIRS" "ODI_$DIRS" "FDI_$DIRS") + +for J in ${SCHEME_Nr[@]} +do + echo $J + ID=${SCHEME_ID[$J]} + IDD=${SCHEME_IDD[$J]} + DIR=${SCHEME_DIRS[$J]} + if [ ! -e $DIR ] + then + echo "Folder $DIR does not exist" + mkdir $DIR + fi + cd $DIR + if [ ! -e sw_$ID ] + then + cur_dir=`pwd` + cd ../../../../ + pwd + ./TestAllMHD.sh $IDD + cd $cur_dir + cp ../../../../sw_$ID . + fi + #../../../../sw_VeP --hydro --threads=16 ../BW_schemes.yml 2>&1 > out.log + cat <<-EOF > ./run.sh + #!/bin/bash + # Run SWIFT + ./sw_$ID --hydro --threads=16 ../BW_schemes.yml 2>&1 > out.log + + # Plot the evolution + python3 ../plot_all.py 0 41 2>&1 > plot.log + EOF + chmod u+x ./run.sh + ./run.sh & + cd .. +done + diff --git a/examples/parameter_example.yml b/examples/parameter_example.yml index c94a480fa1658b3b1c541698171dc686379b1d5c..0947c0da43928544de0f2a5f1f973748d8dec09e 100644 --- a/examples/parameter_example.yml +++ b/examples/parameter_example.yml @@ -545,8 +545,9 @@ GrackleCooling: RT_H2_dissociation_rate_cgs: 0 # H2 dissociation rate in cgs [1/s] volumetric_heating_rates_cgs: 0 # Volumetric heating rate in cgs [erg/s/cm3] specific_heating_rates_cgs: 0 # Specific heating rate in cgs [erg/s/g] - - + H2_three_body_rate : 1 # Specific the H2 formation three body rate (0->5,see Grackle documentation) + H2_cie_cooling : 0 # Enable/disable H2 collision-induced emission cooling from Ripamonti & Abel (2004) + cmb_temperature_floor : 1 # Enable/disable an effective CMB temperature floor @@ -788,7 +789,7 @@ SPINJETAGN: include_jets: 1 # Global switch whether to include jet feedback [1] or not [0]. turn_off_radiative_feedback: 0 # Global switch whether to turn off radiative (thermal) feedback [1] or not [0]. This should only be used if 'include_jets' is set to 1, since we want feedback in some form or another. alpha_acc: 0.2 # Viscous alpha of the subgrid accretion disks. Likely to be within the 0.1-0.3 range. The main effect is that it sets the transition accretion rate between the thin and thick disk, as dot(m) = 0.2 * alpha^2. - mdot_crit_ADAF: 0.01 # The transition normalized accretion rate (Eddington ratio) at which the disc goes from thick (low accretion rates) to thin (high accretion rates). The feedback also changes from kinetic jets to thermal isotropic, respectively. + mdot_crit_ADAF: 0.03 # The transition normalized accretion rate (Eddington ratio) at which the disc goes from thick (low accretion rates) to thin (high accretion rates). The feedback also changes from kinetic jets to thermal isotropic, respectively. seed_spin: 0.01 # The (randomly-directed) black hole spin assigned to BHs when they are seeded. Should be strictly between 0 and 1. AGN_jet_velocity_model: Constant # How AGN jet velocities are calculated. If 'Constant', a single value is used. If 'BlackHoleMass', then an empirical relation between halo mass and black hole mass is used to calculate jet velocities. 'HaloMass' is currently not supported. v_jet_km_p_s: 3160. # Jet velocity to use if 'AGN_jet_velocity_model' is 'Constant'. Units are km/s. @@ -807,7 +808,7 @@ SPINJETAGN: fix_jet_efficiency: 0 # Global switch whether to fix jet efficiency to a particular value [1], or use a spin-dependant formula [0]. jet_efficiency: 0.1 # The constant jet efficiency used if 'fix_jet_efficiency' is set to 1. fix_jet_direction: 0 # Global switch whether to fix the jet direction to be along the z-axis, instead of along the spin vector. - accretion_efficiency: 1. # The accretion efficiency (suppression factor of the accretion rate) to use in the thick disc (ADAF), to represent the effects of subgrid winds that take away most of the mass flowing through the accretion disc. + accretion_efficiency: 0.1 # The accretion efficiency (suppression factor of the accretion rate) to use in the thick disc (ADAF), to represent the effects of subgrid winds that take away most of the mass flowing through the accretion disc. fix_radiative_efficiency: 0 # Global switch whether to fix the radiative efficiency to a particular value [1], or use a spin-dependant formula [0]. radiative_efficiency: 0.1 # The constant jet efficiency used if 'fix_radiative_efficiency' is set to 1. Otherwise, this value is used to define the Eddington accretion rate. TD_region: B # How to treat the subgrid accretion disk if it is thin, according to the Shakura & Sunyaev (1973) model. If set to B, region b will be used. If set to C, region c will be used. @@ -815,7 +816,7 @@ SPINJETAGN: turn_off_secondary_feedback: 1 # If set to 1, there will be only radiative (thermal) feedback in the thin disk mode, and only jets in the thick disk mode. jet_h_r_slope: 1. # The slope of the dependence of jet efficiency on aspect ratio of the subgrid accretion disk, H/R. Default value is 1, and another reasonable value is 0 (same jet efficiency for all disks). Reality could be anything in between. This parameter is only used if turn_off_secondary_feedback is set to 0. delta_ADAF: 0.2 # Electron heating parameter, which controls the strength of radiative feedback in thick disks. Should be between 0.1 and 0.5. This parameter is only used if turn_off_secondary_feedback is set to 0. - include_slim_disk: 1 # Global switch whether to include super-Eddington accretion, modeled as the slim disk. If set to 0, disks will be considered thin even at very large accretion rates. + include_slim_disk: 0 # Global switch whether to include super-Eddington accretion, modeled as the slim disk. If set to 0, disks will be considered thin even at very large accretion rates. # Parameters related to the neutrinos -------------------------------------------- diff --git a/src/Makefile.am b/src/Makefile.am index 934ef0508f0ad6541012c7146f648b4b5650ebae..0b8048c51a8aeebd3b86349394c576f1d213dfdc 100644 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -201,7 +201,7 @@ AM_SOURCES += $(GEAR_RT_SOURCES) # Include files for distribution, not installation. nobase_noinst_HEADERS = align.h approx_math.h atomic.h barrier.h cycle.h error.h inline.h kernel_hydro.h kernel_gravity.h -nobase_noinst_HEADERS += gravity_iact.h kernel_long_gravity.h vector.h accumulate.h cache.h exp.h log.h +nobase_noinst_HEADERS += gravity_iact.h kernel_long_gravity.h vector.h accumulate.h cache.h exp.h log.h matrix.h nobase_noinst_HEADERS += runner_doiact_nosort.h runner_doiact_hydro.h runner_doiact_stars.h runner_doiact_black_holes.h runner_doiact_grav.h nobase_noinst_HEADERS += runner_doiact_functions_hydro.h runner_doiact_functions_stars.h runner_doiact_functions_black_holes.h nobase_noinst_HEADERS += runner_doiact_functions_limiter.h runner_doiact_limiter.h units.h intrinsics.h minmax.h @@ -245,6 +245,9 @@ nobase_noinst_HEADERS += hydro/AnarchyPU/hydro_parameters.h nobase_noinst_HEADERS += hydro/SPHENIX/hydro.h hydro/SPHENIX/hydro_iact.h hydro/SPHENIX/hydro_io.h nobase_noinst_HEADERS += hydro/SPHENIX/hydro_debug.h hydro/SPHENIX/hydro_part.h hydro/SPHENIX/hydro_csds.h nobase_noinst_HEADERS += hydro/SPHENIX/hydro_parameters.h +nobase_noinst_HEADERS += hydro/MAGMA/hydro.h hydro/MAGMA/hydro_iact.h hydro/MAGMA/hydro_io.h +nobase_noinst_HEADERS += hydro/MAGMA/hydro_debug.h hydro/MAGMA/hydro_part.h +nobase_noinst_HEADERS += hydro/MAGMA/hydro_parameters.h nobase_noinst_HEADERS += hydro/Gasoline/hydro.h hydro/Gasoline/hydro_iact.h hydro/Gasoline/hydro_io.h nobase_noinst_HEADERS += hydro/Gasoline/hydro_debug.h hydro/Gasoline/hydro_part.h nobase_noinst_HEADERS += hydro/Gasoline/hydro_parameters.h diff --git a/src/black_holes/SPIN_JET/black_holes_properties.h b/src/black_holes/SPIN_JET/black_holes_properties.h index 9a044d297abd22daf8b9390453af56541de0e38d..ff497d45a951a5aef8937c451cac82cd0362f11c 100644 --- a/src/black_holes/SPIN_JET/black_holes_properties.h +++ b/src/black_holes/SPIN_JET/black_holes_properties.h @@ -668,9 +668,10 @@ INLINE static void black_holes_props_init(struct black_holes_props *bp, bp->seed_spin); } - /* Calculate the critical transition accretion rate between the thick and + /* The critical transition accretion rate between the thick and thin disk regimes. */ - bp->mdot_crit_ADAF = 0.2 * bp->alpha_acc_2; + bp->mdot_crit_ADAF = + parser_get_param_float(params, "SPINJETAGN:mdot_crit_ADAF"); /* Calculate the gas-to-total pressure ratio as based on simulations (see Yuan & Narayan 2014) */ diff --git a/src/cooling/grackle/cooling.c b/src/cooling/grackle/cooling.c index 215ec4cd9db8570dc1603a9f5f46ca8bb3afc473..317490ef2f32b7f600f210c8146e6b9feaa1a6db 100644 --- a/src/cooling/grackle/cooling.c +++ b/src/cooling/grackle/cooling.c @@ -179,45 +179,7 @@ void cooling_compute_equilibrium(const struct phys_const* phys_const, const struct cooling_function_data* cooling, const struct part* p, struct xpart* xp) { - /* TODO: this can fail spectacularly and needs to be replaced. */ - - /* get temporary data */ - struct part p_tmp = *p; - struct cooling_function_data cooling_tmp = *cooling; - cooling_tmp.chemistry_data.with_radiative_cooling = 0; - /* need density for computation, therefore quick estimate */ - p_tmp.rho = 0.2387 * p_tmp.mass / pow(p_tmp.h, 3); - - /* compute time step */ - const double alpha = 0.01; - double dt = fabs(cooling_time(phys_const, us, hydro_properties, cosmo, - &cooling_tmp, &p_tmp, xp)); - cooling_new_energy(phys_const, us, cosmo, hydro_properties, &cooling_tmp, - &p_tmp, xp, dt, dt); - dt = alpha * fabs(cooling_time(phys_const, us, hydro_properties, cosmo, - &cooling_tmp, &p_tmp, xp)); - - /* init simple variables */ - int step = 0; - const int max_step = cooling_tmp.max_step; - const float conv_limit = cooling_tmp.convergence_limit; - struct xpart old; - - do { - /* update variables */ - step += 1; - old = *xp; - - /* update chemistry */ - cooling_new_energy(phys_const, us, cosmo, hydro_properties, &cooling_tmp, - &p_tmp, xp, dt, dt); - } while (step < max_step && !cooling_converged(xp, &old, conv_limit)); - - if (step == max_step) - error( - "A particle element fraction failed to converge." - "You can change 'GrackleCooling:MaxSteps' or " - "'GrackleCooling:ConvergenceLimit' to avoid this problem"); + return; } /** @@ -249,15 +211,13 @@ void cooling_first_init_part(const struct phys_const* phys_const, * a better determination will be done in cooling_post_init_part */ /* primordial chemistry >= 1 */ + /* assume neutral gas */ xp->cooling_data.HI_frac = grackle_data->HydrogenFractionByMass; xp->cooling_data.HII_frac = zero; - xp->cooling_data.HeI_frac = zero; + xp->cooling_data.HeI_frac = 1. - grackle_data->HydrogenFractionByMass; xp->cooling_data.HeII_frac = zero; - xp->cooling_data.HeIII_frac = 1. - grackle_data->HydrogenFractionByMass; - xp->cooling_data.e_frac = xp->cooling_data.HII_frac + - 0.25 * xp->cooling_data.HeII_frac + - 0.5 * xp->cooling_data.HeIII_frac; - + xp->cooling_data.HeIII_frac = zero; + xp->cooling_data.e_frac = zero; #endif // MODE >= 1 #if COOLING_GRACKLE_MODE >= 2 @@ -301,9 +261,9 @@ void cooling_post_init_part(const struct phys_const* phys_const, // message("rho = %g energy = %g",rho,energy); #if COOLING_GRACKLE_MODE > 0 - /* TODO: this can fail spectacularly and needs to be replaced. */ - // cooling_compute_equilibrium(phys_const, us, hydro_props, cosmo, cooling, - // p,xp); + /* The function below currently does nothing. Will have to be updated. */ + cooling_compute_equilibrium(phys_const, us, hydro_props, cosmo, cooling, p, + xp); #endif } @@ -366,6 +326,12 @@ void cooling_print_backend(const struct cooling_function_data* cooling) { cooling->chemistry_data.with_radiative_cooling); message("grackle_chemistry_data.primordial_chemistry = %d", cooling->chemistry_data.primordial_chemistry); + message("grackle_chemistry_data.three_body_rate = %d", + cooling->chemistry_data.three_body_rate); + message("grackle_chemistry_data.cmb_temperature_floor = %d", + cooling->chemistry_data.cmb_temperature_floor); + message("grackle_chemistry_data.cie_cooling = %d", + cooling->chemistry_data.cie_cooling); message("grackle_chemistry_data.dust_chemistry = %d", cooling->chemistry_data.dust_chemistry); message("grackle_chemistry_data.metal_cooling = %d", @@ -1149,9 +1115,19 @@ void cooling_init_grackle(struct cooling_function_data* cooling) { chemistry->primordial_chemistry = cooling->primordial_chemistry; chemistry->metal_cooling = cooling->with_metal_cooling; chemistry->UVbackground = cooling->with_uv_background; + chemistry->three_body_rate = cooling->H2_three_body_rate; + chemistry->cmb_temperature_floor = cooling->cmb_temperature_floor; + chemistry->cie_cooling = cooling->H2_cie_cooling; chemistry->grackle_data_file = cooling->cloudy_table; /* radiative transfer */ +#if COOLING_GRACKLE_MODE == 0 + if (cooling->use_radiative_transfer) + error( + "The parameter use_radiative_transfer cannot be set to 1 in Grackle " + "mode 0 !"); +#endif + chemistry->use_radiative_transfer = cooling->use_radiative_transfer; if (cooling->volumetric_heating_rates > 0) @@ -1173,7 +1149,10 @@ void cooling_init_grackle(struct cooling_function_data* cooling) { "heating rates, not both"); /* self shielding */ - chemistry->self_shielding_method = cooling->self_shielding_method; + if (cooling->self_shielding_method <= 0) + chemistry->self_shielding_method = 0; + else + chemistry->self_shielding_method = cooling->self_shielding_method; if (local_initialize_chemistry_data(&cooling->chemistry_data, &cooling->chemistry_rates, diff --git a/src/cooling/grackle/cooling_io.h b/src/cooling/grackle/cooling_io.h index fb15b730f9c75260df0a56fca0d2e19c60e748df..1ecdad94589e6e3c97946f23fcd6b94a39329e14 100644 --- a/src/cooling/grackle/cooling_io.h +++ b/src/cooling/grackle/cooling_io.h @@ -98,17 +98,15 @@ __attribute__((always_inline)) INLINE static int cooling_write_particles( #endif #if COOLING_GRACKLE_MODE >= 2 - list += num; - - list[0] = + list[6] = io_make_output_field("HM", FLOAT, 1, UNIT_CONV_NO_UNITS, 0.f, xparts, cooling_data.HM_frac, "H- mass fraction"); - list[1] = + list[7] = io_make_output_field("H2I", FLOAT, 1, UNIT_CONV_NO_UNITS, 0.f, xparts, cooling_data.H2I_frac, "H2I mass fraction"); - list[2] = + list[8] = io_make_output_field("H2II", FLOAT, 1, UNIT_CONV_NO_UNITS, 0.f, xparts, cooling_data.H2II_frac, "H2II mass fraction"); @@ -116,20 +114,17 @@ __attribute__((always_inline)) INLINE static int cooling_write_particles( #endif #if COOLING_GRACKLE_MODE >= 3 - list += num; - - list[0] = + list[9] = io_make_output_field("DI", FLOAT, 1, UNIT_CONV_NO_UNITS, 0.f, xparts, cooling_data.DI_frac, "DI mass fraction"); - list[1] = + list[10] = io_make_output_field("DII", FLOAT, 1, UNIT_CONV_NO_UNITS, 0.f, xparts, cooling_data.DII_frac, "DII mass fraction"); - list[2] = + list[11] = io_make_output_field("HDI", FLOAT, 1, UNIT_CONV_NO_UNITS, 0.f, xparts, cooling_data.HDI_frac, "HDI mass fraction"); - num += 3; #endif @@ -161,6 +156,15 @@ __attribute__((always_inline)) INLINE static void cooling_read_parameters( error("Cannot run primordial chemistry %i when compiled with %i", cooling->primordial_chemistry, COOLING_GRACKLE_MODE); + cooling->H2_three_body_rate = parser_get_opt_param_int( + parameter_file, "GrackleCooling:H2_three_body_rate", 0); + + cooling->H2_cie_cooling = parser_get_opt_param_int( + parameter_file, "GrackleCooling:H2_cie_cooling", 0); + + cooling->cmb_temperature_floor = parser_get_opt_param_int( + parameter_file, "GrackleCooling:cmb_temperature_floor", 1); + cooling->with_uv_background = parser_get_param_int(parameter_file, "GrackleCooling:with_UV_background"); diff --git a/src/cooling/grackle/cooling_properties.h b/src/cooling/grackle/cooling_properties.h index 93b602cffdb4b6130bb773ffa66b7c0c0a6dd394..f08083b15b3a14f6b49dfc6a6df15dd6084518d7 100644 --- a/src/cooling/grackle/cooling_properties.h +++ b/src/cooling/grackle/cooling_properties.h @@ -44,6 +44,16 @@ struct cooling_function_data { /*! Chemistry network */ int primordial_chemistry; + /*! Set the three-body reaction rate (see grackle documentation) */ + int H2_three_body_rate; + + /*! Enable/disable H2 collision-induced emission cooling from Ripamonti & Abel + * (2004) */ + int H2_cie_cooling; + + /*! Enable/disable CMB temperature floor */ + int cmb_temperature_floor; + /*! Redshift to use for the UV backgroud (-1 to use cosmological one) */ double redshift; diff --git a/src/debug.c b/src/debug.c index 6282f8d11840548933425a24f163405082a57598..b6a89b825ab38f064c140f6cca6b632407d4f8d8 100644 --- a/src/debug.c +++ b/src/debug.c @@ -76,6 +76,8 @@ #include "./hydro/Planetary/hydro_debug.h" #elif defined(SPHENIX_SPH) #include "./hydro/SPHENIX/hydro_debug.h" +#elif defined(MAGMA_SPH) +#include "./hydro/MAGMA/hydro_debug.h" #elif defined(GASOLINE_SPH) #include "./hydro/Gasoline/hydro_debug.h" #elif defined(ANARCHY_PU_SPH) diff --git a/src/forcing/roberts_flow_acceleration/forcing.h b/src/forcing/roberts_flow_acceleration/forcing.h index 9a9d33d3c88e8d855081474eca40b9429ce8eb30..5310dfb59a532df74dbff2f5d153550c3688ef11 100644 --- a/src/forcing/roberts_flow_acceleration/forcing.h +++ b/src/forcing/roberts_flow_acceleration/forcing.h @@ -65,10 +65,10 @@ __attribute__((always_inline)) INLINE static void forcing_terms_apply( const double time, const struct forcing_terms* terms, const struct space* s, const struct phys_const* phys_const, struct part* p, struct xpart* xp) { - const float v_sig = hydro_get_signal_velocity(p); + const float c_s = hydro_get_comoving_soundspeed(p); const double L = s->dim[0]; const float u0 = terms->u0; - const float nu = terms->nu * p->viscosity.alpha * v_sig * p->h; + const float nu = terms->nu * p->viscosity.alpha * c_s * p->h; const float Vz_factor = terms->Vz_factor; const double k0 = 2. * M_PI / L; const double kf = M_SQRT2 * k0; @@ -90,13 +90,6 @@ __attribute__((always_inline)) INLINE static void forcing_terms_apply( p->a_hydro[0] += f[0]; p->a_hydro[1] += f[1]; p->a_hydro[2] += f[2]; - - if (time == 0.f) { - /* Force the velocity */ - xp->v_full[0] = v_Rob[0]; - xp->v_full[1] = v_Rob[1]; - xp->v_full[2] = v_Rob[2]; - } } /** diff --git a/src/hydro.c b/src/hydro.c index ae4cfd9742a2426ceaf46fd19984d10a9f3dd74c..e9c6649a8004faf50b53d266f03ead91553a832b 100644 --- a/src/hydro.c +++ b/src/hydro.c @@ -373,7 +373,7 @@ void hydro_exact_density_check(struct space *s, const struct engine *e, h_max_limited); wrong_rho++; } -#if defined(SPHENIX_SPH) +#if defined(SPHENIX_SPH) || defined(MAGMA_SPH) if (check_force && !found_inhibited && (fabsf(pi->n_gradient / pi->n_gradient_exact - 1.f) > rel_tol || fabsf(pi->n_gradient_exact / pi->n_gradient - 1.f) > rel_tol)) { diff --git a/src/hydro.h b/src/hydro.h index 1f96d1f3031033c85d5915ecb3ca1d0959738d70..644c3ef6cfb6bebe4a2286ccff199969def91eb2 100644 --- a/src/hydro.h +++ b/src/hydro.h @@ -75,6 +75,10 @@ #include "./hydro/SPHENIX/hydro.h" #include "./hydro/SPHENIX/hydro_iact.h" #define SPH_IMPLEMENTATION "SPHENIX (Borrow+ 2020)" +#elif defined(MAGMA_SPH) +#include "./hydro/MAGMA/hydro.h" +#include "./hydro/MAGMA/hydro_iact.h" +#define SPH_IMPLEMENTATION "MAGMA-2 (Rosswog+ 2020)" #elif defined(GASOLINE_SPH) #include "./hydro/Gasoline/hydro.h" #include "./hydro/Gasoline/hydro_iact.h" diff --git a/src/hydro/MAGMA/hydro.h b/src/hydro/MAGMA/hydro.h new file mode 100644 index 0000000000000000000000000000000000000000..7e228fc5f8ead21fcacff3fe90f776db8e17fe99 --- /dev/null +++ b/src/hydro/MAGMA/hydro.h @@ -0,0 +1,1043 @@ +/******************************************************************************* + * This file is part of SWIFT. + * Copyright (c) 2023 Matthieu Schaller (schaller@strw.leidenuniv.nl) + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU Lesser General Public License as published + * by the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. + * + ******************************************************************************/ +#ifndef SWIFT_MAGMA_HYDRO_H +#define SWIFT_MAGMA_HYDRO_H + +/** + * @file MAGMA/hydro.h + * @brief MAGMA2 implementation of SPH (Non-neighbour loop + * equations) + * + * The thermal variable is the internal energy (u). + */ + +#include "adiabatic_index.h" +#include "approx_math.h" +#include "cosmology.h" +#include "dimension.h" +#include "entropy_floor.h" +#include "equation_of_state.h" +#include "hydro_parameters.h" +#include "hydro_properties.h" +#include "hydro_space.h" +#include "kernel_hydro.h" +#include "minmax.h" +#include "pressure_floor.h" + +/* System includes */ +#include <string.h> + +/** + * @brief Returns the comoving internal energy of a particle at the last + * time the particle was kicked. + * + * @param p The particle of interest + * @param xp The extended data of the particle of interest. + */ +__attribute__((always_inline)) INLINE static float +hydro_get_comoving_internal_energy(const struct part *restrict p, + const struct xpart *restrict xp) { + + return xp->u_full; +} + +/** + * @brief Returns the physical internal energy of a particle at the last + * time the particle was kicked. + * + * @param p The particle of interest. + * @param xp The extended data of the particle of interest. + * @param cosmo The cosmological model. + */ +__attribute__((always_inline)) INLINE static float +hydro_get_physical_internal_energy(const struct part *restrict p, + const struct xpart *restrict xp, + const struct cosmology *cosmo) { + + return xp->u_full * cosmo->a_factor_internal_energy; +} + +/** + * @brief Returns the comoving internal energy of a particle drifted to the + * current time. + * + * @param p The particle of interest + */ +__attribute__((always_inline)) INLINE static float +hydro_get_drifted_comoving_internal_energy(const struct part *restrict p) { + + return p->u; +} + +/** + * @brief Returns the physical internal energy of a particle drifted to the + * current time. + * + * @param p The particle of interest. + * @param cosmo The cosmological model. + */ +__attribute__((always_inline)) INLINE static float +hydro_get_drifted_physical_internal_energy(const struct part *restrict p, + const struct cosmology *cosmo) { + + return p->u * cosmo->a_factor_internal_energy; +} + +/** + * @brief Returns the comoving pressure of a particle + * + * Computes the pressure based on the particle's properties. + * + * @param p The particle of interest + */ +__attribute__((always_inline)) INLINE static float hydro_get_comoving_pressure( + const struct part *restrict p) { + + return gas_pressure_from_internal_energy(p->rho, p->u); +} + +/** + * @brief Returns the physical pressure of a particle + * + * Computes the pressure based on the particle's properties and + * convert it to physical coordinates. + * + * @param p The particle of interest + * @param cosmo The cosmological model. + */ +__attribute__((always_inline)) INLINE static float hydro_get_physical_pressure( + const struct part *restrict p, const struct cosmology *cosmo) { + + return cosmo->a_factor_pressure * + gas_pressure_from_internal_energy(p->rho, p->u); +} + +/** + * @brief Returns the comoving entropy of a particle at the last + * time the particle was kicked. + * + * @param p The particle of interest. + * @param xp The extended data of the particle of interest. + */ +__attribute__((always_inline)) INLINE static float hydro_get_comoving_entropy( + const struct part *restrict p, const struct xpart *restrict xp) { + + return gas_entropy_from_internal_energy(p->rho, xp->u_full); +} + +/** + * @brief Returns the physical entropy of a particle at the last + * time the particle was kicked. + * + * @param p The particle of interest. + * @param xp The extended data of the particle of interest. + * @param cosmo The cosmological model. + */ +__attribute__((always_inline)) INLINE static float hydro_get_physical_entropy( + const struct part *restrict p, const struct xpart *restrict xp, + const struct cosmology *cosmo) { + + /* Note: no cosmological conversion required here with our choice of + * coordinates. */ + return gas_entropy_from_internal_energy(p->rho, xp->u_full); +} + +/** + * @brief Returns the comoving entropy of a particle drifted to the + * current time. + * + * @param p The particle of interest. + */ +__attribute__((always_inline)) INLINE static float +hydro_get_drifted_comoving_entropy(const struct part *restrict p) { + + return gas_entropy_from_internal_energy(p->rho, p->u); +} + +/** + * @brief Returns the physical entropy of a particle drifted to the + * current time. + * + * @param p The particle of interest. + * @param cosmo The cosmological model. + */ +__attribute__((always_inline)) INLINE static float +hydro_get_drifted_physical_entropy(const struct part *restrict p, + const struct cosmology *cosmo) { + + /* Note: no cosmological conversion required here with our choice of + * coordinates. */ + return gas_entropy_from_internal_energy(p->rho, p->u); +} + +/** + * @brief Returns the comoving sound speed of a particle + * + * @param p The particle of interest + */ +__attribute__((always_inline)) INLINE static float +hydro_get_comoving_soundspeed(const struct part *restrict p) { + + return p->force.soundspeed; +} + +/** + * @brief Returns the physical sound speed of a particle + * + * @param p The particle of interest + * @param cosmo The cosmological model. + */ +__attribute__((always_inline)) INLINE static float +hydro_get_physical_soundspeed(const struct part *restrict p, + const struct cosmology *cosmo) { + + return cosmo->a_factor_sound_speed * p->force.soundspeed; +} + +/** + * @brief Returns the comoving density of a particle + * + * @param p The particle of interest + */ +__attribute__((always_inline)) INLINE static float hydro_get_comoving_density( + const struct part *restrict p) { + + return p->rho; +} + +/** + * @brief Returns the comoving density of a particle. + * + * @param p The particle of interest + * @param cosmo The cosmological model. + */ +__attribute__((always_inline)) INLINE static float hydro_get_physical_density( + const struct part *restrict p, const struct cosmology *cosmo) { + + return cosmo->a3_inv * p->rho; +} + +/** + * @brief Returns the mass of a particle + * + * @param p The particle of interest + */ +__attribute__((always_inline)) INLINE static float hydro_get_mass( + const struct part *restrict p) { + + return p->mass; +} + +/** + * @brief Sets the mass of a particle + * + * @param p The particle of interest + * @param m The mass to set. + */ +__attribute__((always_inline)) INLINE static void hydro_set_mass( + struct part *restrict p, float m) { + + p->mass = m; +} + +/** + * @brief Returns the time derivative of co-moving internal energy of a particle + * + * We assume a constant density. + * + * @param p The particle of interest + */ +__attribute__((always_inline)) INLINE static float +hydro_get_comoving_internal_energy_dt(const struct part *restrict p) { + + return p->u_dt; +} + +/** + * @brief Returns the time derivative of internal energy of a particle + * + * We assume a constant density. + * + * @param p The particle of interest + * @param cosmo Cosmology data structure + */ +__attribute__((always_inline)) INLINE static float +hydro_get_physical_internal_energy_dt(const struct part *restrict p, + const struct cosmology *cosmo) { + + return p->u_dt * cosmo->a_factor_internal_energy; +} + +/** + * @brief Sets the time derivative of the co-moving internal energy of a + * particle + * + * We assume a constant density for the conversion to entropy. + * + * @param p The particle of interest. + * @param du_dt The new time derivative of the internal energy. + */ +__attribute__((always_inline)) INLINE static void +hydro_set_comoving_internal_energy_dt(struct part *restrict p, float du_dt) { + + p->u_dt = du_dt; +} + +/** + * @brief Returns the time derivative of internal energy of a particle + * + * We assume a constant density. + * + * @param p The particle of interest. + * @param cosmo Cosmology data structure + * @param du_dt The new time derivative of the internal energy. + */ +__attribute__((always_inline)) INLINE static void +hydro_set_physical_internal_energy_dt(struct part *restrict p, + const struct cosmology *cosmo, + float du_dt) { + + p->u_dt = du_dt / cosmo->a_factor_internal_energy; +} + +/** + * @brief Sets the physical entropy of a particle + * + * @param p The particle of interest. + * @param xp The extended particle data. + * @param cosmo Cosmology data structure + * @param entropy The physical entropy + */ +__attribute__((always_inline)) INLINE static void hydro_set_physical_entropy( + struct part *p, struct xpart *xp, const struct cosmology *cosmo, + const float entropy) { + + /* Note there is no conversion from physical to comoving entropy */ + const float comoving_entropy = entropy; + xp->u_full = gas_internal_energy_from_entropy(p->rho, comoving_entropy); +} + +/** + * @brief Sets the physical internal energy of a particle + * + * @param p The particle of interest. + * @param xp The extended particle data. + * @param cosmo Cosmology data structure + * @param u The physical internal energy + */ +__attribute__((always_inline)) INLINE static void +hydro_set_physical_internal_energy(struct part *p, struct xpart *xp, + const struct cosmology *cosmo, + const float u) { + + xp->u_full = u / cosmo->a_factor_internal_energy; +} + +/** + * @brief Sets the drifted physical internal energy of a particle + * + * @param p The particle of interest. + * @param cosmo Cosmology data structure + * @param u The physical internal energy + */ +__attribute__((always_inline)) INLINE static void +hydro_set_drifted_physical_internal_energy( + struct part *p, const struct cosmology *cosmo, + const struct pressure_floor_props *pressure_floor, const float u) { + + p->u = u / cosmo->a_factor_internal_energy; + + /* Now recompute the extra quantities */ + + /* Compute the sound speed */ + const float pressure = gas_pressure_from_internal_energy(p->rho, p->u); + const float soundspeed = gas_soundspeed_from_pressure(p->rho, pressure); + + /* Update variables. */ + p->force.pressure = pressure; + p->force.soundspeed = soundspeed; + + p->force.v_sig = max(p->force.v_sig, 2.f * soundspeed); +} + +/** + * @brief Correct the signal velocity of the particle partaking in + * supernova (kinetic) feedback based on the velocity kick the particle receives + * + * @param p The particle of interest. + * @param cosmo Cosmology data structure + * @param dv_phys The velocity kick received by the particle expressed in + * physical units (note that dv_phys must be positive or equal to zero) + */ +__attribute__((always_inline)) INLINE static void +hydro_set_v_sig_based_on_velocity_kick(struct part *p, + const struct cosmology *cosmo, + const float dv_phys) { + + /* Compute the velocity kick in comoving coordinates */ + const float dv = dv_phys / cosmo->a_factor_sound_speed; + + /* Sound speed */ + const float soundspeed = hydro_get_comoving_soundspeed(p); + + /* Update the signal velocity */ + p->force.v_sig = + max(2.f * soundspeed, p->force.v_sig + const_viscosity_beta * dv); +} + +/** + * @brief Update the value of the viscosity alpha for the scheme. + * + * @param p the particle of interest + * @param alpha the new value for the viscosity coefficient. + */ +__attribute__((always_inline)) INLINE static void hydro_set_viscosity_alpha( + struct part *restrict p, float alpha) { + /* This scheme has fixed alpha */ +} + +/** + * @brief Update the value of the diffusive coefficients to the + * feedback reset value for the scheme. + * + * @param p the particle of interest + */ +__attribute__((always_inline)) INLINE static void +hydro_diffusive_feedback_reset(struct part *restrict p) { + /* This scheme has fixed alpha */ +} + +/** + * @brief Computes the hydro time-step of a given particle + * + * This function returns the time-step of a particle given its hydro-dynamical + * state. A typical time-step calculation would be the use of the CFL condition. + * + * @param p Pointer to the particle data + * @param xp Pointer to the extended particle data + * @param hydro_properties The SPH parameters + * @param cosmo The cosmological model. + */ +__attribute__((always_inline)) INLINE static float hydro_compute_timestep( + const struct part *restrict p, const struct xpart *restrict xp, + const struct hydro_props *restrict hydro_properties, + const struct cosmology *restrict cosmo) { + + const float CFL_condition = hydro_properties->CFL_condition; + + /* CFL condition */ + const float dt_cfl = 2.f * kernel_gamma * CFL_condition * cosmo->a * p->h / + (cosmo->a_factor_sound_speed * p->force.v_sig); + + return dt_cfl; +} + +/** + * @brief Compute the signal velocity between two gas particles + * + * This is eq. (103) of Price D., JCoPh, 2012, Vol. 231, Issue 3. + * + * @param dx Comoving vector separating both particles (pi - pj). + * @brief pi The first #part. + * @brief pj The second #part. + * @brief mu_ij The velocity on the axis linking the particles, or zero if the + * particles are moving away from each other, + * @brief beta The non-linear viscosity constant. + */ +__attribute__((always_inline)) INLINE static float hydro_signal_velocity( + const float dx[3], const struct part *restrict pi, + const struct part *restrict pj, const float mu_ij, const float beta) { + + const float ci = pi->force.soundspeed; + const float cj = pj->force.soundspeed; + + return ci + cj - beta * mu_ij; +} + +/** + * @brief returns the signal velocity + * + * @brief p the particle + */ +__attribute__((always_inline)) INLINE static float hydro_get_signal_velocity( + const struct part *restrict p) { + + return p->force.v_sig; +} +/** + * @brief returns the div_v + * + * @brief p the particle + */ +__attribute__((always_inline)) INLINE static float hydro_get_div_v( + const struct part *restrict p) { + + return 0.; + ; +} + +/** + * @brief Does some extra hydro operations once the actual physical time step + * for the particle is known. + * + * @param p The particle to act upon. + * @param dt Physical time step of the particle during the next step. + */ +__attribute__((always_inline)) INLINE static void hydro_timestep_extra( + struct part *p, float dt) {} + +/** + * @brief Prepares a particle for the density calculation. + * + * Zeroes all the relevant arrays in preparation for the sums taking place in + * the various density loop over neighbours. Typically, all fields of the + * density sub-structure of a particle get zeroed in here. + * + * @param p The particle to act upon + * @param hs #hydro_space containing hydro specific space information. + */ +__attribute__((always_inline)) INLINE static void hydro_init_part( + struct part *restrict p, const struct hydro_space *hs) { + + p->density.wcount = 0.f; + p->density.wcount_dh = 0.f; + p->rho = 0.f; + p->density.rho_dh = 0.f; + p->force.f = 0.f; +} + +/** + * @brief Finishes the density calculation. + * + * Multiplies the density and number of neighbours by the appropiate constants + * and add the self-contribution term. + * + * Also adds/multiplies the cosmological terms if need be. + * + * @param p The particle to act upon + * @param cosmo The cosmological model. + */ +__attribute__((always_inline)) INLINE static void hydro_end_density( + struct part *restrict p, const struct cosmology *cosmo) { + + /* Some smoothing length multiples. */ + const float h = p->h; + const float h_inv = 1.0f / h; /* 1/h */ + const float h_inv_dim = pow_dimension(h_inv); /* 1/h^d */ + const float h_inv_dim_plus_one = h_inv_dim * h_inv; /* 1/h^(d+1) */ + + /* Final operation on the density (add self-contribution). */ + p->rho += p->mass * kernel_root; + p->density.rho_dh -= hydro_dimension * p->mass * kernel_root; + p->density.wcount += kernel_root; + p->density.wcount_dh -= hydro_dimension * kernel_root; + + /* Finish the calculation by inserting the missing h-factors */ + p->rho *= h_inv_dim; + p->density.rho_dh *= h_inv_dim_plus_one; + p->density.wcount *= h_inv_dim; + p->density.wcount_dh *= h_inv_dim_plus_one; +} + +/** + * @brief Prepare a particle for the gradient calculation. + * + * This function is called after the density loop and before the gradient loop. + * Nothing to do in this scheme as there are no terms needing a preparation + * after the density loop. + * + * @param p The particle to act upon. + * @param xp The extended particle data to act upon. + * @param cosmo The cosmological model. + * @param hydro_props Hydrodynamic properties. + */ +__attribute__((always_inline)) INLINE static void hydro_prepare_gradient( + struct part *restrict p, struct xpart *restrict xp, + const struct cosmology *cosmo, const struct hydro_props *hydro_props, + const struct pressure_floor_props *pressure_floor) {} + +/** + * @brief Resets the variables that are required for a gradient calculation. + * + * This function is called after hydro_prepare_gradient. + * + * @param p The particle to act upon. + * @param xp The extended particle data to act upon. + * @param cosmo The cosmological model. + */ +__attribute__((always_inline)) INLINE static void hydro_reset_gradient( + struct part *restrict p) { + + zero_sym_matrix(&p->gradient.c_matrix_inv); + for (int i = 0; i < 3; ++i) p->gradient.gradient_vx[i] = 0.f; + for (int i = 0; i < 3; ++i) p->gradient.gradient_vy[i] = 0.f; + for (int i = 0; i < 3; ++i) p->gradient.gradient_vz[i] = 0.f; + for (int i = 0; i < 3; ++i) p->gradient.gradient_u[i] = 0.f; +} + +/** + * @brief Finishes the gradient calculation. + * + * Multiplies the C-matrix by the appropiate constants. + * + * Also adds/multiplies the cosmological terms if need be. + * Nothing to do in this scheme as the gradient loop is not used. + * + * @param p The particle to act upon. + */ +__attribute__((always_inline)) INLINE static void hydro_end_gradient( + struct part *p) { + + /* Some smoothing length multiples. */ + const float h = p->h; + const float h_inv = 1.0f / h; /* 1/h */ + const float h_inv_dim = pow_dimension(h_inv); /* 1/h^d */ + + /* Finish the construction of the inverse of the c-matrix by + * multiplying in the factors of h coming from W */ + for (int i = 0; i < 6; ++i) { + p->gradient.c_matrix_inv.elements[i] *= h_inv_dim; + } + + /* Finish the construction of the inverse of the velocity gradient + * multiplying in the factors of h coming from W */ + for (int i = 0; i < 3; ++i) p->gradient.gradient_vx[i] *= h_inv_dim; + for (int i = 0; i < 3; ++i) p->gradient.gradient_vy[i] *= h_inv_dim; + for (int i = 0; i < 3; ++i) p->gradient.gradient_vz[i] *= h_inv_dim; +} + +/** + * @brief Sets all particle fields to sensible values when the #part has 0 ngbs. + * + * In the desperate case where a particle has no neighbours (likely because + * of the h_max ceiling), set the particle fields to something sensible to avoid + * NaNs in the next calculations. + * + * @param p The particle to act upon + * @param xp The extended particle data to act upon + * @param cosmo The cosmological model. + */ +__attribute__((always_inline)) INLINE static void hydro_part_has_no_neighbours( + struct part *restrict p, struct xpart *restrict xp, + const struct cosmology *cosmo) { + + /* Some smoothing length multiples. */ + const float h = p->h; + const float h_inv = 1.0f / h; /* 1/h */ + const float h_inv_dim = pow_dimension(h_inv); /* 1/h^d */ + + warning( + "Gas particle with ID %lld treated as having no neighbours (h: %g, " + "wcount: %g).", + p->id, h, p->density.wcount); + + /* Re-set problematic values */ + p->rho = p->mass * kernel_root * h_inv_dim; + p->density.wcount = kernel_root * h_inv_dim; + p->density.rho_dh = 0.f; + p->density.wcount_dh = 0.f; +} + +/** + * @brief Prepare a particle for the force calculation. + * + * This function is called in the ghost task to convert some quantities coming + * from the density loop over neighbours into quantities ready to be used in the + * force loop over neighbours. Quantities are typically read from the density + * sub-structure and written to the force sub-structure. + * Examples of calculations done here include the calculation of viscosity term + * constants, thermal conduction terms, hydro conversions, etc. + * + * @param p The particle to act upon + * @param xp The extended particle data to act upon + * @param cosmo The current cosmological model. + * @param hydro_props Hydrodynamic properties. + * @param dt_alpha The time-step used to evolve non-cosmological quantities such + * as the artificial viscosity. + * @param dt_therm The time-step used to evolve hydrodynamical quantities. + */ +__attribute__((always_inline)) INLINE static void hydro_prepare_force( + struct part *restrict p, struct xpart *restrict xp, + const struct cosmology *cosmo, const struct hydro_props *hydro_props, + const struct pressure_floor_props *pressure_floor, const float dt_alpha, + const float dt_therm) { + + /* Compute the pressure */ + const float pressure = gas_pressure_from_internal_energy(p->rho, p->u); + + /* Compute the sound speed */ + const float soundspeed = gas_soundspeed_from_pressure(p->rho, pressure); + + /* Invert the c-matrix */ + float c_matrix_temp[3][3]; + get_matrix_from_sym_matrix(c_matrix_temp, &p->gradient.c_matrix_inv); + int res = invert_dimension_by_dimension_matrix(c_matrix_temp); + if (res) { + sym_matrix_print(&p->gradient.c_matrix_inv); + error("Error inverting matrix"); + } + + /* TODO: Write a routine to invert symmetric matrices */ + + /* Finish computation of velocity gradient (eq. 18) */ + const float gradient_vx[3] = {p->gradient.gradient_vx[0], + p->gradient.gradient_vx[1], + p->gradient.gradient_vx[2]}; + const float gradient_vy[3] = {p->gradient.gradient_vy[0], + p->gradient.gradient_vy[1], + p->gradient.gradient_vy[2]}; + const float gradient_vz[3] = {p->gradient.gradient_vz[0], + p->gradient.gradient_vz[1], + p->gradient.gradient_vz[2]}; + + p->force.gradient_vx[0] = c_matrix_temp[0][0] * gradient_vx[0] + + c_matrix_temp[0][1] * gradient_vx[1] + + c_matrix_temp[0][2] * gradient_vx[2]; + p->force.gradient_vx[1] = c_matrix_temp[1][0] * gradient_vx[0] + + c_matrix_temp[1][1] * gradient_vx[1] + + c_matrix_temp[1][2] * gradient_vx[2]; + p->force.gradient_vx[2] = c_matrix_temp[2][0] * gradient_vx[0] + + c_matrix_temp[2][1] * gradient_vx[1] + + c_matrix_temp[2][2] * gradient_vx[2]; + + p->force.gradient_vy[0] = c_matrix_temp[0][0] * gradient_vy[0] + + c_matrix_temp[0][1] * gradient_vy[1] + + c_matrix_temp[0][2] * gradient_vy[2]; + p->force.gradient_vy[1] = c_matrix_temp[1][0] * gradient_vy[0] + + c_matrix_temp[1][1] * gradient_vy[1] + + c_matrix_temp[1][2] * gradient_vy[2]; + p->force.gradient_vy[2] = c_matrix_temp[2][0] * gradient_vy[0] + + c_matrix_temp[2][1] * gradient_vy[1] + + c_matrix_temp[2][2] * gradient_vy[2]; + + p->force.gradient_vz[0] = c_matrix_temp[0][0] * gradient_vz[0] + + c_matrix_temp[0][1] * gradient_vz[1] + + c_matrix_temp[0][2] * gradient_vz[2]; + p->force.gradient_vz[1] = c_matrix_temp[1][0] * gradient_vz[0] + + c_matrix_temp[1][1] * gradient_vz[1] + + c_matrix_temp[1][2] * gradient_vz[2]; + p->force.gradient_vz[2] = c_matrix_temp[2][0] * gradient_vz[0] + + c_matrix_temp[2][1] * gradient_vz[1] + + c_matrix_temp[2][2] * gradient_vz[2]; + + /* Finish computation of u gradient (same as eq. 18) */ + const float gradient_u[3] = {p->gradient.gradient_u[0], + p->gradient.gradient_u[1], + p->gradient.gradient_u[2]}; + + p->force.gradient_u[0] = c_matrix_temp[0][0] * gradient_u[0] + + c_matrix_temp[0][1] * gradient_u[1] + + c_matrix_temp[0][2] * gradient_u[2]; + p->force.gradient_u[1] = c_matrix_temp[1][0] * gradient_u[0] + + c_matrix_temp[1][1] * gradient_u[1] + + c_matrix_temp[1][2] * gradient_u[2]; + p->force.gradient_u[2] = c_matrix_temp[2][0] * gradient_u[0] + + c_matrix_temp[2][1] * gradient_u[1] + + c_matrix_temp[2][2] * gradient_u[2]; + + /* Update other variables. */ + p->force.pressure = pressure; + p->force.soundspeed = soundspeed; + get_sym_matrix_from_matrix(&p->force.c_matrix, c_matrix_temp); +} + +/** + * @brief Reset acceleration fields of a particle + * + * Resets all hydro acceleration and time derivative fields in preparation + * for the sums taking place in the various force tasks. + * + * @param p The particle to act upon + */ +__attribute__((always_inline)) INLINE static void hydro_reset_acceleration( + struct part *restrict p) { + + /* Reset the acceleration. */ + p->a_hydro[0] = 0.0f; + p->a_hydro[1] = 0.0f; + p->a_hydro[2] = 0.0f; + + /* Reset the time derivatives. */ + p->u_dt = 0.0f; + p->force.h_dt = 0.0f; + p->force.v_sig = 2.f * p->force.soundspeed; +} + +/** + * @brief Sets the values to be predicted in the drifts to their values at a + * kick time + * + * @param p The particle. + * @param xp The extended data of this particle. + * @param cosmo The cosmological model + */ +__attribute__((always_inline)) INLINE static void hydro_reset_predicted_values( + struct part *restrict p, const struct xpart *restrict xp, + const struct cosmology *cosmo, + const struct pressure_floor_props *pressure_floor) { + + /* Re-set the predicted velocities */ + p->v[0] = xp->v_full[0]; + p->v[1] = xp->v_full[1]; + p->v[2] = xp->v_full[2]; + + /* Re-set the entropy */ + p->u = xp->u_full; + + /* Re-compute the pressure */ + const float pressure = gas_pressure_from_internal_energy(p->rho, p->u); + + /* Compute the new sound speed */ + const float soundspeed = gas_soundspeed_from_pressure(p->rho, pressure); + + /* Update variables */ + p->force.pressure = pressure; + p->force.soundspeed = soundspeed; + + p->force.v_sig = max(p->force.v_sig, 2.f * soundspeed); +} + +/** + * @brief Predict additional particle fields forward in time when drifting + * + * Additional hydrodynamic quantites are drifted forward in time here. These + * include thermal quantities (thermal energy or total energy or entropy, ...). + * + * Note the different time-step sizes used for the different quantities as they + * include cosmological factors. + * + * @param p The particle. + * @param xp The extended data of the particle. + * @param dt_drift The drift time-step for positions. + * @param dt_therm The drift time-step for thermal quantities. + * @param dt_kick_grav The time-step for gravity quantities. + * @param cosmo The cosmological model. + * @param hydro_props The properties of the hydro scheme. + * @param floor_props The properties of the entropy floor. + */ +__attribute__((always_inline)) INLINE static void hydro_predict_extra( + struct part *restrict p, const struct xpart *restrict xp, float dt_drift, + float dt_therm, float dt_kick_grav, const struct cosmology *cosmo, + const struct hydro_props *hydro_props, + const struct entropy_floor_properties *floor_props, + const struct pressure_floor_props *pressure_floor) { + + /* Predict the internal energy */ + p->u += p->u_dt * dt_therm; + + const float h_inv = 1.f / p->h; + + /* Predict smoothing length */ + const float w1 = p->force.h_dt * h_inv * dt_drift; + if (fabsf(w1) < 0.2f) { + p->h *= approx_expf(w1); /* 4th order expansion of exp(w) */ + } else { + p->h *= expf(w1); + } + + /* Predict density */ + const float w2 = -hydro_dimension * w1; + if (fabsf(w2) < 0.2f) { + p->rho *= approx_expf(w2); /* 4th order expansion of exp(w) */ + } else { + p->rho *= expf(w2); + } + + /* Check against entropy floor */ + const float floor_A = entropy_floor(p, cosmo, floor_props); + const float floor_u = gas_internal_energy_from_entropy(p->rho, floor_A); + + /* Check against absolute minimum */ + const float min_u = + hydro_props->minimal_internal_energy / cosmo->a_factor_internal_energy; + + p->u = max(p->u, floor_u); + p->u = max(p->u, min_u); + + /* Compute the new pressure */ + const float pressure = gas_pressure_from_internal_energy(p->rho, p->u); + + /* Compute the new sound speed */ + const float soundspeed = gas_soundspeed_from_pressure(p->rho, pressure); + + p->force.pressure = pressure; + p->force.soundspeed = soundspeed; + + p->force.v_sig = max(p->force.v_sig, 2.f * soundspeed); +} + +/** + * @brief Finishes the force calculation. + * + * Multiplies the force and accelerations by the appropiate constants + * and add the self-contribution term. In most cases, there is little + * to do here. + * + * Cosmological terms are also added/multiplied here. + * + * @param p The particle to act upon + * @param cosmo The current cosmological model. + */ +__attribute__((always_inline)) INLINE static void hydro_end_force( + struct part *restrict p, const struct cosmology *cosmo) { + + p->force.h_dt *= p->h * hydro_dimension_inv; +} + +/** + * @brief Kick the additional variables + * + * Additional hydrodynamic quantites are kicked forward in time here. These + * include thermal quantities (thermal energy or total energy or entropy, ...). + * + * @param p The particle to act upon. + * @param xp The particle extended data to act upon. + * @param dt_therm The time-step for this kick (for thermodynamic quantities). + * @param dt_grav The time-step for this kick (for gravity quantities). + * @param dt_grav_mesh The time-step for this kick (mesh gravity). + * @param dt_hydro The time-step for this kick (for hydro quantities). + * @param dt_kick_corr The time-step for this kick (for gravity corrections). + * @param cosmo The cosmological model. + * @param hydro_props The constants used in the scheme. + * @param floor_props The properties of the entropy floor. + */ +__attribute__((always_inline)) INLINE static void hydro_kick_extra( + struct part *restrict p, struct xpart *restrict xp, float dt_therm, + float dt_grav, float dt_grav_mesh, float dt_hydro, float dt_kick_corr, + const struct cosmology *cosmo, const struct hydro_props *hydro_props, + const struct entropy_floor_properties *floor_props) { + + /* Integrate the internal energy forward in time */ + const float delta_u = p->u_dt * dt_therm; + + /* Do not decrease the energy by more than a factor of 2*/ + xp->u_full = max(xp->u_full + delta_u, 0.5f * xp->u_full); + + /* Check against entropy floor */ + const float floor_A = entropy_floor(p, cosmo, floor_props); + const float floor_u = gas_internal_energy_from_entropy(p->rho, floor_A); + + /* Check against absolute minimum */ + const float min_u = + hydro_props->minimal_internal_energy / cosmo->a_factor_internal_energy; + + /* Take highest of both limits */ + const float energy_min = max(min_u, floor_u); + + if (xp->u_full < energy_min) { + xp->u_full = energy_min; + p->u_dt = 0.f; + } +} + +/** + * @brief Converts hydro quantity of a particle at the start of a run + * + * This function is called once at the end of the engine_init_particle() + * routine (at the start of a calculation) after the densities of + * particles have been computed. + * This can be used to convert internal energy into entropy for instance. + * + * @param p The particle to act upon + * @param xp The extended particle to act upon + * @param cosmo The cosmological model. + * @param hydro_props The constants used in the scheme. + */ +__attribute__((always_inline)) INLINE static void hydro_convert_quantities( + struct part *restrict p, struct xpart *restrict xp, + const struct cosmology *cosmo, const struct hydro_props *hydro_props, + const struct pressure_floor_props *pressure_floor) { + + /* Convert the physcial internal energy to the comoving one. */ + /* u' = a^(3(g-1)) u */ + const float factor = 1.f / cosmo->a_factor_internal_energy; + p->u *= factor; + xp->u_full = p->u; + + /* Apply the minimal energy limit */ + const float min_comoving_energy = + hydro_props->minimal_internal_energy / cosmo->a_factor_internal_energy; + if (xp->u_full < min_comoving_energy) { + xp->u_full = min_comoving_energy; + p->u = min_comoving_energy; + p->u_dt = 0.f; + } + + /* Compute the pressure */ + const float pressure = gas_pressure_from_internal_energy(p->rho, p->u); + + /* Compute the sound speed */ + const float soundspeed = gas_soundspeed_from_internal_energy(p->rho, p->u); + + p->force.pressure = pressure; + p->force.soundspeed = soundspeed; +} + +/** + * @brief Initialises the particles for the first time + * + * This function is called only once just after the ICs have been + * read in to do some conversions or assignments between the particle + * and extended particle fields. + * + * @param p The particle to act upon + * @param xp The extended particle data to act upon + */ +__attribute__((always_inline)) INLINE static void hydro_first_init_part( + struct part *restrict p, struct xpart *restrict xp) { + + p->time_bin = 0; + xp->v_full[0] = p->v[0]; + xp->v_full[1] = p->v[1]; + xp->v_full[2] = p->v[2]; + xp->u_full = p->u; + + hydro_reset_acceleration(p); + hydro_init_part(p, NULL); +} + +/** + * @brief Overwrite the initial internal energy of a particle. + * + * Note that in the cases where the thermodynamic variable is not + * internal energy but gets converted later, we must overwrite that + * field. The conversion to the actual variable happens later after + * the initial fake time-step. + * + * @param p The #part to write to. + * @param u_init The new initial internal energy. + */ +__attribute__((always_inline)) INLINE static void +hydro_set_init_internal_energy(struct part *p, float u_init) { + + p->u = u_init; +} + +/** + * @brief Operations performed when a particle gets removed from the + * simulation volume. + * + * @param p The particle. + * @param xp The extended particle data. + * @param time The simulation time. + */ +__attribute__((always_inline)) INLINE static void hydro_remove_part( + const struct part *p, const struct xpart *xp, const double time) {} + +#endif /* SWIFT_MAGMA_HYDRO_H */ diff --git a/src/hydro/MAGMA/hydro_debug.h b/src/hydro/MAGMA/hydro_debug.h new file mode 100644 index 0000000000000000000000000000000000000000..d96c0b471954f9f988a62761ded5d678e2bde712 --- /dev/null +++ b/src/hydro/MAGMA/hydro_debug.h @@ -0,0 +1,58 @@ +/******************************************************************************* + * This file is part of SWIFT. + * Copyright (c) 2016 Matthieu Schaller (schaller@strw.leidenuniv.nl) + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU Lesser General Public License as published + * by the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. + * + ******************************************************************************/ +#ifndef SWIFT_MAGMA_HYDRO_DEBUG_H +#define SWIFT_MAGMA_HYDRO_DEBUG_H + +/** + * @file Minimal/hydro_debug.h + * @brief Minimal conservative implementation of SPH (Debugging routines) + * + * The thermal variable is the internal energy (u). Simple constant + * viscosity term with the Balsara (1995) switch. No thermal conduction + * term is implemented. + * + * This corresponds to equations (43), (44), (45), (101), (103) and (104) with + * \f$\beta=3\f$ and \f$\alpha_u=0\f$ of + * Price, D., Journal of Computational Physics, 2012, Volume 231, Issue 3, + * pp. 759-794. + */ + +__attribute__((always_inline)) INLINE static void hydro_debug_particle( + const struct part* p, const struct xpart* xp) { + warning("[PID%lld] part:", p->id); + warning( + "[PID%lld] " + "x=[%.6g, %.6g, %.6g], v=[%.3g, %.3g, %.3g], " + "a=[%.3g, %.3g, %.3g], " + "m=%.3g, u=%.3g, du/dt=%.3g, P=%.3g, c_s=%.3g, " + "v_sig=%.3g, h=%.3g, dh/dt=%.3g, wcount=%.3g, rho=%.3g, " + "dh_drho=%.3g, time_bin=%d wakeup=%d", + p->id, p->x[0], p->x[1], p->x[2], p->v[0], p->v[1], p->v[2], + p->a_hydro[0], p->a_hydro[1], p->a_hydro[2], p->mass, p->u, p->u_dt, + hydro_get_comoving_pressure(p), p->force.soundspeed, p->force.v_sig, p->h, + p->force.h_dt, p->density.wcount, p->rho, p->density.rho_dh, p->time_bin, + p->limiter_data.wakeup); + if (xp != NULL) { + warning("[PID%lld] xpart:", p->id); + warning("[PID%lld] v_full=[%.3g, %.3g, %.3g]", p->id, xp->v_full[0], + xp->v_full[1], xp->v_full[2]); + } +} + +#endif /* SWIFT_MAGMA_HYDRO_DEBUG_H */ diff --git a/src/hydro/MAGMA/hydro_iact.h b/src/hydro/MAGMA/hydro_iact.h new file mode 100644 index 0000000000000000000000000000000000000000..d9da3766cc2fb993712f77c4d79458b7d6637ef0 --- /dev/null +++ b/src/hydro/MAGMA/hydro_iact.h @@ -0,0 +1,511 @@ +/******************************************************************************* + * This file is part of SWIFT. + * Copyright (c) 2016 Matthieu Schaller (schaller@strw.leidenuniv.nl) + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU Lesser General Public License as published + * by the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. + * + ******************************************************************************/ +#ifndef SWIFT_MAGMA_HYDRO_IACT_H +#define SWIFT_MAGMA_HYDRO_IACT_H + +/** + * @file Minimal/hydro_iact.h + * @brief Minimal conservative implementation of SPH (Neighbour loop equations) + * + * The thermal variable is the internal energy (u). Simple constant + * viscosity term with the Balsara (1995) switch. No thermal conduction + * term is implemented. + * + * This corresponds to equations (43), (44), (45), (101), (103) and (104) with + * \f$\beta=3\f$ and \f$\alpha_u=0\f$ of Price, D., Journal of Computational + * Physics, 2012, Volume 231, Issue 3, pp. 759-794. + */ + +#include "adiabatic_index.h" +#include "hydro_parameters.h" +#include "minmax.h" +#include "signal_velocity.h" + +/** + * @brief Density interaction between two particles (non-symmetric). + * + * @param r2 Comoving square distance between the two particles. + * @param dx Comoving vector separating both particles (pi - pj). + * @param hi Comoving smoothing-length of particle i. + * @param hj Comoving smoothing-length of particle j. + * @param pi First particle. + * @param pj Second particle (not updated). + * @param a Current scale factor. + * @param H Current Hubble parameter. + */ +__attribute__((always_inline)) INLINE static void runner_iact_nonsym_density( + const float r2, const float dx[3], const float hi, const float hj, + struct part *restrict pi, const struct part *restrict pj, const float mu_0, + const float a, const float H) { + + float wi, wi_dx; + +#ifdef SWIFT_DEBUG_CHECKS + if (pi->time_bin >= time_bin_inhibited) + error("Inhibited pi in interaction function!"); + if (pj->time_bin >= time_bin_inhibited) + error("Inhibited pj in interaction function!"); +#endif + + /* Get the masses. */ + const float mj = pj->mass; + + /* Get r and 1/r. */ + const float r = sqrtf(r2); + + const float h_inv = 1.f / hi; + const float ui = r * h_inv; + kernel_deval(ui, &wi, &wi_dx); + + pi->rho += mj * wi; + pi->density.rho_dh -= mj * (hydro_dimension * wi + ui * wi_dx); + pi->density.wcount += wi; + pi->density.wcount_dh -= (hydro_dimension * wi + ui * wi_dx); +} + +/** + * @brief Density interaction between two particles. + * + * @param r2 Comoving square distance between the two particles. + * @param dx Comoving vector separating both particles (pi - pj). + * @param hi Comoving smoothing-length of particle i. + * @param hj Comoving smoothing-length of particle j. + * @param pi First particle. + * @param pj Second particle. + * @param a Current scale factor. + * @param H Current Hubble parameter. + */ +__attribute__((always_inline)) INLINE static void runner_iact_density( + const float r2, const float dx[3], const float hi, const float hj, + struct part *restrict pi, struct part *restrict pj, const float mu_0, + const float a, const float H) { + + runner_iact_nonsym_density(r2, dx, hi, hj, pi, pj, mu_0, a, H); + const float dx_rev[3] = {-dx[0], -dx[1], -dx[2]}; + runner_iact_nonsym_density(r2, dx_rev, hj, hi, pj, pi, mu_0, a, H); +} + +/** + * @brief Calculate the gradient interaction between particle i and particle j: + * non-symmetric version + * + * Nothing to do here in this scheme. + * + * @param r2 Comoving squared distance between particle i and particle j. + * @param dx Comoving distance vector between the particles (dx = pi->x - + * pj->x). + * @param hi Comoving smoothing-length of particle i. + * @param hj Comoving smoothing-length of particle j. + * @param pi Particle i. + * @param pj Particle j. + * @param a Current scale factor. + * @param H Current Hubble parameter. + */ +__attribute__((always_inline)) INLINE static void runner_iact_nonsym_gradient( + const float r2, const float dx[3], const float hi, const float hj, + struct part *restrict pi, struct part *restrict pj, const float mu_0, + const float a, const float H) { + + /* Get r. */ + const float r = sqrtf(r2); + + /* Compute the kernel function */ + const float h_inv = 1.f / hi; + const float ui = r * h_inv; + float w; + kernel_eval(ui, &w); + + /* Get the mass and density. */ + const float mj = pj->mass; + const float rhoj = pj->rho; + + /* Velocity difference */ + const float vij[3] = {pj->v[0] - pi->v[0], pj->v[1] - pi->v[1], + pj->v[2] - pi->v[2]}; + + /* Internal energy difference */ + const float uij = pj->u - pi->u; + + const float common_term = w * mj / rhoj; + + /* The inverse of the C-matrix. eq. 6 + * It's symmetric so recall we only store the 6 useful terms. */ + pi->gradient.c_matrix_inv.xx += common_term * dx[0] * dx[0]; + pi->gradient.c_matrix_inv.yy += common_term * dx[1] * dx[1]; + pi->gradient.c_matrix_inv.zz += common_term * dx[2] * dx[2]; + pi->gradient.c_matrix_inv.xy += common_term * dx[0] * dx[1]; + pi->gradient.c_matrix_inv.xz += common_term * dx[0] * dx[2]; + pi->gradient.c_matrix_inv.yz += common_term * dx[1] * dx[2]; + + /* Gradient of v (recall dx is pi - pj), eq. 18 */ + pi->gradient.gradient_vx[0] -= common_term * vij[0] * dx[0]; + pi->gradient.gradient_vx[1] -= common_term * vij[0] * dx[1]; + pi->gradient.gradient_vx[2] -= common_term * vij[0] * dx[2]; + + pi->gradient.gradient_vy[0] -= common_term * vij[1] * dx[0]; + pi->gradient.gradient_vy[1] -= common_term * vij[1] * dx[1]; + pi->gradient.gradient_vy[2] -= common_term * vij[1] * dx[2]; + + pi->gradient.gradient_vz[0] -= common_term * vij[2] * dx[0]; + pi->gradient.gradient_vz[1] -= common_term * vij[2] * dx[1]; + pi->gradient.gradient_vz[2] -= common_term * vij[2] * dx[2]; + + pi->gradient.gradient_u[0] -= common_term * uij * dx[0]; + pi->gradient.gradient_u[1] -= common_term * uij * dx[1]; + pi->gradient.gradient_u[2] -= common_term * uij * dx[2]; +} + +/** + * @brief Calculate the gradient interaction between particle i and particle j + * + * Nothing to do here in this scheme. + * + * @param r2 Comoving squared distance between particle i and particle j. + * @param dx Comoving distance vector between the particles (dx = pi->x - + * pj->x). + * @param hi Comoving smoothing-length of particle i. + * @param hj Comoving smoothing-length of particle j. + * @param pi Particle i. + * @param pj Particle j. + * @param a Current scale factor. + * @param H Current Hubble parameter. + */ +__attribute__((always_inline)) INLINE static void runner_iact_gradient( + const float r2, const float dx[3], const float hi, const float hj, + struct part *restrict pi, struct part *restrict pj, const float mu_0, + const float a, const float H) { + + runner_iact_nonsym_gradient(r2, dx, hi, hj, pi, pj, mu_0, a, H); + const float dx_rev[3] = {-dx[0], -dx[1], -dx[2]}; + runner_iact_nonsym_gradient(r2, dx_rev, hj, hi, pj, pi, mu_0, a, H); +} + +/** + * @brief Force interaction between two particles (non-symmetric). + * + * @param r2 Comoving square distance between the two particles. + * @param dx Comoving vector separating both particles (pi - pj). + * @param hi Comoving smoothing-length of particle i. + * @param hj Comoving smoothing-length of particle j. + * @param pi First particle. + * @param pj Second particle (not updated). + * @param a Current scale factor. + * @param H Current Hubble parameter. + */ +__attribute__((always_inline)) INLINE static void runner_iact_nonsym_force( + const float r2, const float dx[3], const float hi, const float hj, + struct part *restrict pi, const struct part *restrict pj, const float mu_0, + const float a, const float H) { + +#ifdef SWIFT_DEBUG_CHECKS + if (pi->time_bin >= time_bin_inhibited) + error("Inhibited pi in interaction function!"); + if (pj->time_bin >= time_bin_inhibited) + error("Inhibited pj in interaction function!"); +#endif + + /* Cosmological factors entering the EoMs */ + const float fac_mu = pow_three_gamma_minus_five_over_two(a); + const float a2_Hubble = a * a * H; + + /* Get r and 1/r. */ + const float r = sqrtf(r2); + const float r_inv = r ? 1.0f / r : 0.0f; + + /* Recover some data */ + // const float mi = pi->mass; + const float mj = pj->mass; + const float rhoi = pi->rho; + const float rhoj = pj->rho; + const float pressurei = pi->force.pressure; + const float pressurej = pj->force.pressure; + const float ci = pi->force.soundspeed; + const float cj = pj->force.soundspeed; + + /* Get the kernel for hi. */ + const float hi_inv = 1.0f / hi; + const float hid_inv = pow_dimension(hi_inv); /* 1/h^d */ + const float ui = r * hi_inv; + float wi, wi_dx; + kernel_deval(ui, &wi, &wi_dx); + + /* Get the kernel for hj. */ + const float hj_inv = 1.0f / hj; + const float hjd_inv = pow_dimension(hj_inv); /* 1/h^d */ + const float uj = r * hj_inv; + float wj, wj_dx; + kernel_deval(uj, &wj, &wj_dx); + + /* Velocity difference */ + const float v_ij[3] = {pi->v[0] - pj->v[0], /* x */ + pi->v[1] - pj->v[1], /* y */ + pi->v[2] - pj->v[2]}; /* z */ + + /* Compute dv dot r. */ + const float dvdr = v_ij[0] * dx[0] + v_ij[1] * dx[1] + v_ij[2] * dx[2]; + + /* Add Hubble flow */ + const float dvdr_Hubble = dvdr + a2_Hubble * r2; + + /* Are the particles moving towards each others ? */ + const float omega_ij = min(dvdr_Hubble, 0.f); + const float mu_ij = fac_mu * r_inv * omega_ij; /* This is 0 or negative */ + + /* Compute signal velocity */ + const float v_sig = + signal_velocity(dx, pi, pj, mu_ij, const_viscosity_beta, a, mu_0); + + /* De-dimentionalised distances (eq. 16, recall dx = xi - xj)*/ + const float eta_i[3] = {dx[0] / hi, dx[1] / hi, dx[2] / hi}; + const float eta_j[3] = {-dx[0] / hj, -dx[1] / hj, -dx[2] / hj}; + + /* Norms of the eta vectors (eq. 16) */ + const float eta_square_i = + eta_i[0] * eta_i[0] + eta_i[1] * eta_i[1] + eta_i[2] * eta_i[2]; + const float eta_square_j = + eta_j[0] * eta_j[0] + eta_j[1] * eta_j[1] + eta_j[2] * eta_j[2]; + + /* Reconstructed velocities at the mid-point (before reconstruction) */ + float v_rec_i[3] = {pi->v[0], pi->v[1], pi->v[2]}; + float v_rec_j[3] = {pj->v[0], pj->v[1], pj->v[2]}; + +#ifndef USE_ZEROTH_ORDER_VELOCITIES + + /* Vectors from the particles to the mid-point */ + const float delta_i[3] = {0.5f * (pj->x[0] - pi->x[0]), + 0.5f * (pj->x[1] - pi->x[1]), + 0.5f * (pj->x[2] - pi->x[2])}; + const float delta_j[3] = {-delta_i[0], -delta_i[1], -delta_i[2]}; + + /* Terms entering the limiter (eq. 23) */ + const float eta_ij = sqrtf(fminf(eta_square_i, eta_square_j)); + const float eta_crit = 0.5f; + + /* Van Leer limiter fraction (eq. 22) */ + const float A_ij_num = pi->gradient.gradient_vx[0] * dx[0] * dx[0] + + pi->gradient.gradient_vx[1] * dx[0] * dx[1] + + pi->gradient.gradient_vx[2] * dx[0] * dx[2] + + pi->gradient.gradient_vy[0] * dx[1] * dx[0] + + pi->gradient.gradient_vy[1] * dx[1] * dx[1] + + pi->gradient.gradient_vy[2] * dx[1] * dx[2] + + pi->gradient.gradient_vz[0] * dx[2] * dx[0] + + pi->gradient.gradient_vz[1] * dx[2] * dx[1] + + pi->gradient.gradient_vz[2] * dx[2] * dx[2]; + + const float A_ij_den = pj->gradient.gradient_vx[0] * dx[0] * dx[0] + + pj->gradient.gradient_vx[1] * dx[0] * dx[1] + + pj->gradient.gradient_vx[2] * dx[0] * dx[2] + + pj->gradient.gradient_vy[0] * dx[1] * dx[0] + + pj->gradient.gradient_vy[1] * dx[1] * dx[1] + + pj->gradient.gradient_vy[2] * dx[1] * dx[2] + + pj->gradient.gradient_vz[0] * dx[2] * dx[0] + + pj->gradient.gradient_vz[1] * dx[2] * dx[1] + + pj->gradient.gradient_vz[2] * dx[2] * dx[2]; + + const float A_ij = A_ij_den != 0.f ? A_ij_num / A_ij_den : 0.f; + + /* Slope limiter exponential term (eq. 21, right term) */ + const float exp_term = + eta_ij < eta_crit + ? expf(-25.f * (eta_ij - eta_crit) * (eta_ij - eta_crit)) + : 1.f; + + /* Van Leer limiter (eq. 21) */ + const float fraction = + (A_ij != -1.f) ? 4.f * A_ij / ((1.f + A_ij) * (1.f + A_ij)) : 1.f; + const float Phi_ij = fmaxf(0.f, fminf(1.f, fraction)) * exp_term; + + /* Mid-point reconstruction, first order (eq. 17) */ + v_rec_i[0] += Phi_ij * pi->gradient.gradient_vx[0] * delta_i[0]; + v_rec_i[0] += Phi_ij * pi->gradient.gradient_vx[1] * delta_i[1]; + v_rec_i[0] += Phi_ij * pi->gradient.gradient_vx[2] * delta_i[2]; + v_rec_i[1] += Phi_ij * pi->gradient.gradient_vy[0] * delta_i[0]; + v_rec_i[1] += Phi_ij * pi->gradient.gradient_vy[1] * delta_i[1]; + v_rec_i[1] += Phi_ij * pi->gradient.gradient_vy[2] * delta_i[2]; + v_rec_i[2] += Phi_ij * pi->gradient.gradient_vz[0] * delta_i[0]; + v_rec_i[2] += Phi_ij * pi->gradient.gradient_vz[1] * delta_i[1]; + v_rec_i[2] += Phi_ij * pi->gradient.gradient_vz[2] * delta_i[2]; + + v_rec_j[0] += Phi_ij * pj->gradient.gradient_vx[0] * delta_j[0]; + v_rec_j[0] += Phi_ij * pj->gradient.gradient_vx[1] * delta_j[1]; + v_rec_j[0] += Phi_ij * pj->gradient.gradient_vx[2] * delta_j[2]; + v_rec_j[1] += Phi_ij * pj->gradient.gradient_vy[0] * delta_j[0]; + v_rec_j[1] += Phi_ij * pj->gradient.gradient_vy[1] * delta_j[1]; + v_rec_j[1] += Phi_ij * pj->gradient.gradient_vy[2] * delta_j[2]; + v_rec_j[2] += Phi_ij * pj->gradient.gradient_vz[0] * delta_j[0]; + v_rec_j[2] += Phi_ij * pj->gradient.gradient_vz[1] * delta_j[1]; + v_rec_j[2] += Phi_ij * pj->gradient.gradient_vz[2] * delta_j[2]; + +#endif + + /* Difference in velocity at the mid-point */ + const float v_rec_ij[3] = {v_rec_i[0] - v_rec_j[0], v_rec_i[1] - v_rec_j[1], + v_rec_i[2] - v_rec_j[2]}; + + /* Normalised relative velocity (eq. 15) */ + const float vel_rel_i = + eta_i[0] * v_rec_ij[0] + eta_i[1] * v_rec_ij[1] + eta_i[2] * v_rec_ij[2]; + const float vel_rel_j = eta_j[0] * -v_rec_ij[0] + eta_j[1] * -v_rec_ij[1] + + eta_j[2] * -v_rec_ij[2]; + + /* Terms entering the viscosity (eq. 15) */ + const float eps_squared = const_viscosity_epsilon * const_viscosity_epsilon; + const float mu_i = fminf(0.f, vel_rel_i / (eta_square_i + eps_squared)); + const float mu_j = fminf(0.f, vel_rel_j / (eta_square_j + eps_squared)); + + /* Eq. 14 */ + const float Qi = + 1.f * rhoi * + (-const_viscosity_alpha * ci * mu_i + const_viscosity_beta * mu_i * mu_i); + const float Qj = + 1.f * rhoj * + (-const_viscosity_alpha * cj * mu_j + const_viscosity_beta * mu_j * mu_j); + + /* Construct the gradient functions (eq. 4 and 5) */ + float G_i[3], G_j[3]; + sym_matrix_multiply_by_vector(G_i, &pi->force.c_matrix, dx); + sym_matrix_multiply_by_vector(G_j, &pj->force.c_matrix, dx); + + /* Note we multiply by -1 as dx is (pi - pj) and not (pj - pi) */ + G_i[0] *= -wi * hid_inv; + G_i[1] *= -wi * hid_inv; + G_i[2] *= -wi * hid_inv; + G_j[0] *= -wj * hjd_inv; + G_j[1] *= -wj * hjd_inv; + G_j[2] *= -wj * hjd_inv; + +#ifdef USE_STANDARD_KERNEL_GRADIENTS + const float wi_dr = hid_inv * hi_inv * wi_dx; + const float wj_dr = hjd_inv * hj_inv * wj_dx; + G_i[0] = wi_dr * r_inv * dx[0]; + G_i[1] = wi_dr * r_inv * dx[1]; + G_i[2] = wi_dr * r_inv * dx[2]; + G_j[0] = wj_dr * r_inv * dx[0]; + G_j[1] = wj_dr * r_inv * dx[1]; + G_j[2] = wj_dr * r_inv * dx[2]; +#endif + +#ifdef I_LOVE_GASOLINE + + const float G_ij[3] = {0.5f * (G_i[0] + G_j[0]), 0.5f * (G_i[1] + G_j[1]), + 0.5f * (G_i[2] + G_j[2])}; + + const float acc_term = (pressurei + pressurej + Qi + Qj) / (rhoi * rhoj); + const float du_term = (pressurei + Qi) / (rhoi * rhoj); + + /* Raw fluid acceleration (eq. 10) */ + pi->a_hydro[0] -= mj * acc_term * G_ij[0]; + pi->a_hydro[1] -= mj * acc_term * G_ij[1]; + pi->a_hydro[2] -= mj * acc_term * G_ij[2]; + + /* Equivalent of div v */ + const float v_ij_dot_G_ij = + v_ij[0] * G_ij[0] + v_ij[1] * G_ij[1] + v_ij[2] * G_ij[2]; + + /* Raw change in internal energy (eq. 11) */ + pi->u_dt += du_term * mj * v_ij_dot_G_ij; + +#else + + /* Compute pressure terms */ + const float P_over_rho2_i = (pressurei + Qi) / (rhoi * rhoi); + const float P_over_rho2_j = (pressurej + Qj) / (rhoj * rhoj); + + /* Raw fluid acceleration (eq. 2) */ + pi->a_hydro[0] -= mj * (P_over_rho2_i * G_i[0] + P_over_rho2_j * G_j[0]); + pi->a_hydro[1] -= mj * (P_over_rho2_i * G_i[1] + P_over_rho2_j * G_j[1]); + pi->a_hydro[2] -= mj * (P_over_rho2_i * G_i[2] + P_over_rho2_j * G_j[2]); + + /* Equivalent of div v */ + const float v_ij_dot_G_i = + v_ij[0] * G_i[0] + v_ij[1] * G_i[1] + v_ij[2] * G_i[2]; + + /* Raw change in internal energy (eq. 3) */ + pi->u_dt += P_over_rho2_i * mj * v_ij_dot_G_i; + +#endif + + /* Reconstructed internal energies at the mid-point (before reconstruction) */ + float u_rec_i = pi->u; + float u_rec_j = pj->u; + +#ifndef USE_ZEROTH_ORDER_VELOCITIES + + /* Mid-point reconstruction, first order (eq. 17) */ + u_rec_i += pi->gradient.gradient_u[0] * delta_i[0]; + u_rec_i += pi->gradient.gradient_u[1] * delta_i[1]; + u_rec_i += pi->gradient.gradient_u[2] * delta_i[2]; + + u_rec_j += pj->gradient.gradient_u[0] * delta_j[0]; + u_rec_j += pj->gradient.gradient_u[1] * delta_j[1]; + u_rec_j += pj->gradient.gradient_u[2] * delta_j[2]; + +#endif + + /* Difference in internal energy */ + const float delta_u = u_rec_i - u_rec_j; + + /* Norm of the G vectors */ + const float sum_G[3] = {G_i[0] + G_j[0], G_i[1] + G_j[1], G_i[2] + G_j[2]}; + const float norm_sum_G = + sqrtf(sum_G[0] * sum_G[0] + sum_G[1] * sum_G[1] + sum_G[2] * sum_G[2]); + + /* Diffusion signal velocity (eq. 26) */ +#ifdef GRAVITY_DIFF_VELOCITY + const float v_sig_u = + sqrtf(v_ij[0] * v_ij[0] + v_ij[1] * v_ij[1] + v_ij[2] * v_ij[2]); +#else + const float v_sig_u = + sqrtf(2. * fabsf(pressurei - pressurej) / (rhoi + rhoj)); +#endif + + /* Diffusion term (eq. 24) */ + pi->u_dt += -const_diffusion_alpha * mj * delta_u * v_sig_u * norm_sum_G / + (rhoi + rhoj); + + /* Get the time derivative for h. */ + pi->force.h_dt -= mj * dvdr * r_inv / rhoj * wi_dx * hi_inv * hid_inv; + + /* Update the signal velocity. */ + pi->force.v_sig = max(pi->force.v_sig, v_sig); +} + +/** + * @brief Force interaction between two particles. + * + * @param r2 Comoving square distance between the two particles. + * @param dx Comoving vector separating both particles (pi - pj). + * @param hi Comoving smoothing-length of particle i. + * @param hj Comoving smoothing-length of particle j. + * @param pi First particle. + * @param pj Second particle. + * @param a Current scale factor. + * @param H Current Hubble parameter. + */ +__attribute__((always_inline)) INLINE static void runner_iact_force( + const float r2, const float dx[3], const float hi, const float hj, + struct part *restrict pi, struct part *restrict pj, const float mu_0, + const float a, const float H) { + + runner_iact_nonsym_force(r2, dx, hi, hj, pi, pj, mu_0, a, H); + const float dx_rev[3] = {-dx[0], -dx[1], -dx[2]}; + runner_iact_nonsym_force(r2, dx_rev, hj, hi, pj, pi, mu_0, a, H); +} + +#endif /* SWIFT_MAGMA_HYDRO_IACT_H */ diff --git a/src/hydro/MAGMA/hydro_io.h b/src/hydro/MAGMA/hydro_io.h new file mode 100644 index 0000000000000000000000000000000000000000..fff5eb78ae00e0905f6704c7cc737139d81760c9 --- /dev/null +++ b/src/hydro/MAGMA/hydro_io.h @@ -0,0 +1,250 @@ +/******************************************************************************* + * This file is part of SWIFT. + * Copyright (c) 2016 Matthieu Schaller (schaller@strw.leidenuniv.nl) + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU Lesser General Public License as published + * by the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. + * + ******************************************************************************/ +#ifndef SWIFT_MAGMA_HYDRO_IO_H +#define SWIFT_MAGMA_HYDRO_IO_H + +/** + * @file Minimal/hydro_io.h + * @brief Minimal conservative implementation of SPH (i/o routines) + * + * The thermal variable is the internal energy (u). Simple constant + * viscosity term with the Balsara (1995) switch. No thermal conduction + * term is implemented. + * + * This corresponds to equations (43), (44), (45), (101), (103) and (104) with + * \f$\beta=3\f$ and \f$\alpha_u=0\f$ of + * Price, D., Journal of Computational Physics, 2012, Volume 231, Issue 3, + * pp. 759-794. + */ + +#include "adiabatic_index.h" +#include "hydro.h" +#include "hydro_parameters.h" +#include "io_properties.h" +#include "kernel_hydro.h" + +/** + * @brief Specifies which particle fields to read from a dataset + * + * @param parts The particle array. + * @param list The list of i/o properties to read. + * @param num_fields The number of i/o fields to read. + */ +INLINE static void hydro_read_particles(struct part* parts, + struct io_props* list, + int* num_fields) { + + *num_fields = 8; + + /* List what we want to read */ + list[0] = io_make_input_field("Coordinates", DOUBLE, 3, COMPULSORY, + UNIT_CONV_LENGTH, parts, x); + list[1] = io_make_input_field("Velocities", FLOAT, 3, COMPULSORY, + UNIT_CONV_SPEED, parts, v); + list[2] = io_make_input_field("Masses", FLOAT, 1, COMPULSORY, UNIT_CONV_MASS, + parts, mass); + list[3] = io_make_input_field("SmoothingLength", FLOAT, 1, COMPULSORY, + UNIT_CONV_LENGTH, parts, h); + list[4] = io_make_input_field("InternalEnergy", FLOAT, 1, COMPULSORY, + UNIT_CONV_ENERGY_PER_UNIT_MASS, parts, u); + list[5] = io_make_input_field("ParticleIDs", ULONGLONG, 1, COMPULSORY, + UNIT_CONV_NO_UNITS, parts, id); + list[6] = io_make_input_field("Accelerations", FLOAT, 3, OPTIONAL, + UNIT_CONV_ACCELERATION, parts, a_hydro); + list[7] = io_make_input_field("Density", FLOAT, 1, OPTIONAL, + UNIT_CONV_DENSITY, parts, rho); +} + +INLINE static void convert_S(const struct engine* e, const struct part* p, + const struct xpart* xp, float* ret) { + + ret[0] = hydro_get_comoving_entropy(p, xp); +} + +INLINE static void convert_P(const struct engine* e, const struct part* p, + const struct xpart* xp, float* ret) { + + ret[0] = hydro_get_comoving_pressure(p); +} + +INLINE static void convert_part_pos(const struct engine* e, + const struct part* p, + const struct xpart* xp, double* ret) { + + const struct space* s = e->s; + if (s->periodic) { + ret[0] = box_wrap(p->x[0], 0.0, s->dim[0]); + ret[1] = box_wrap(p->x[1], 0.0, s->dim[1]); + ret[2] = box_wrap(p->x[2], 0.0, s->dim[2]); + } else { + ret[0] = p->x[0]; + ret[1] = p->x[1]; + ret[2] = p->x[2]; + } + if (e->snapshot_use_delta_from_edge) { + ret[0] = min(ret[0], s->dim[0] - e->snapshot_delta_from_edge); + ret[1] = min(ret[1], s->dim[1] - e->snapshot_delta_from_edge); + ret[2] = min(ret[2], s->dim[2] - e->snapshot_delta_from_edge); + } +} + +INLINE static void convert_part_vel(const struct engine* e, + const struct part* p, + const struct xpart* xp, float* ret) { + + const int with_cosmology = (e->policy & engine_policy_cosmology); + const struct cosmology* cosmo = e->cosmology; + const integertime_t ti_current = e->ti_current; + const double time_base = e->time_base; + const float dt_kick_grav_mesh = e->dt_kick_grav_mesh_for_io; + + const integertime_t ti_beg = get_integer_time_begin(ti_current, p->time_bin); + const integertime_t ti_end = get_integer_time_end(ti_current, p->time_bin); + + /* Get time-step since the last kick */ + float dt_kick_grav, dt_kick_hydro; + if (with_cosmology) { + dt_kick_grav = cosmology_get_grav_kick_factor(cosmo, ti_beg, ti_current); + dt_kick_grav -= + cosmology_get_grav_kick_factor(cosmo, ti_beg, (ti_beg + ti_end) / 2); + dt_kick_hydro = cosmology_get_hydro_kick_factor(cosmo, ti_beg, ti_current); + dt_kick_hydro -= + cosmology_get_hydro_kick_factor(cosmo, ti_beg, (ti_beg + ti_end) / 2); + } else { + dt_kick_grav = (ti_current - ((ti_beg + ti_end) / 2)) * time_base; + dt_kick_hydro = (ti_current - ((ti_beg + ti_end) / 2)) * time_base; + } + + /* Extrapolate the velocites to the current time (hydro term)*/ + ret[0] = xp->v_full[0] + p->a_hydro[0] * dt_kick_hydro; + ret[1] = xp->v_full[1] + p->a_hydro[1] * dt_kick_hydro; + ret[2] = xp->v_full[2] + p->a_hydro[2] * dt_kick_hydro; + + /* Add the gravity term */ + if (p->gpart != NULL) { + ret[0] += p->gpart->a_grav[0] * dt_kick_grav; + ret[1] += p->gpart->a_grav[1] * dt_kick_grav; + ret[2] += p->gpart->a_grav[2] * dt_kick_grav; + } + + /* And the mesh gravity term */ + if (p->gpart != NULL) { + ret[0] += p->gpart->a_grav_mesh[0] * dt_kick_grav_mesh; + ret[1] += p->gpart->a_grav_mesh[1] * dt_kick_grav_mesh; + ret[2] += p->gpart->a_grav_mesh[2] * dt_kick_grav_mesh; + } + + /* Conversion from internal units to peculiar velocities */ + ret[0] *= cosmo->a_inv; + ret[1] *= cosmo->a_inv; + ret[2] *= cosmo->a_inv; +} + +INLINE static void convert_part_potential(const struct engine* e, + const struct part* p, + const struct xpart* xp, float* ret) { + + if (p->gpart != NULL) + ret[0] = gravity_get_comoving_potential(p->gpart); + else + ret[0] = 0.f; +} + +/** + * @brief Specifies which particle fields to write to a dataset + * + * @param parts The particle array. + * @param xparts The extended particle array. + * @param list The list of i/o properties to write. + * @param num_fields The number of i/o fields to write. + */ +INLINE static void hydro_write_particles(const struct part* parts, + const struct xpart* xparts, + struct io_props* list, + int* num_fields) { + + *num_fields = 10; + + /* List what we want to write */ + list[0] = io_make_output_field_convert_part( + "Coordinates", DOUBLE, 3, UNIT_CONV_LENGTH, 1.f, parts, xparts, + convert_part_pos, "Co-moving positions of the particles"); + + list[1] = io_make_output_field_convert_part( + "Velocities", FLOAT, 3, UNIT_CONV_SPEED, 0.f, parts, xparts, + convert_part_vel, + "Peculiar velocities of the stars. This is (a * dx/dt) where x is the " + "co-moving positions of the particles"); + + list[2] = io_make_output_field("Masses", FLOAT, 1, UNIT_CONV_MASS, 0.f, parts, + mass, "Masses of the particles"); + + list[3] = io_make_output_field( + "SmoothingLengths", FLOAT, 1, UNIT_CONV_LENGTH, 1.f, parts, h, + "Co-moving smoothing lengths (FWHM of the kernel) of the particles"); + + list[4] = io_make_output_field( + "InternalEnergies", FLOAT, 1, UNIT_CONV_ENERGY_PER_UNIT_MASS, + -3.f * hydro_gamma_minus_one, parts, u, + "Co-moving thermal energies per unit mass of the particles"); + + list[5] = + io_make_output_field("ParticleIDs", ULONGLONG, 1, UNIT_CONV_NO_UNITS, 0.f, + parts, id, "Unique IDs of the particles"); + + list[6] = io_make_output_field("Densities", FLOAT, 1, UNIT_CONV_DENSITY, -3.f, + parts, rho, + "Co-moving mass densities of the particles"); + + list[7] = io_make_output_field_convert_part( + "Entropies", FLOAT, 1, UNIT_CONV_ENTROPY_PER_UNIT_MASS, 0.f, parts, + xparts, convert_S, "Co-moving entropies per unit mass of the particles"); + + list[8] = io_make_output_field_convert_part( + "Pressures", FLOAT, 1, UNIT_CONV_PRESSURE, -3.f * hydro_gamma, parts, + xparts, convert_P, "Co-moving pressures of the particles"); + + list[9] = io_make_output_field_convert_part( + "Potentials", FLOAT, 1, UNIT_CONV_POTENTIAL, -1.f, parts, xparts, + convert_part_potential, + "Co-moving gravitational potential at position of the particles"); +} + +/** + * @brief Writes the current model of SPH to the file + * @param h_grpsph The HDF5 group in which to write + */ +INLINE static void hydro_write_flavour(hid_t h_grpsph) { + + /* Viscosity and thermal conduction */ + /* Nothing in this minimal model... */ + io_write_attribute_s(h_grpsph, "Thermal Conductivity Model", "No treatment"); + io_write_attribute_s( + h_grpsph, "Viscosity Model", + "as in Springel (2005), i.e. Monaghan (1992) with Balsara (1995) switch"); +} + +/** + * @brief Are we writing entropy in the internal energy field ? + * + * @return 1 if entropy is in 'internal energy', 0 otherwise. + */ +INLINE static int writeEntropyFlag(void) { return 0; } + +#endif /* SWIFT_MAGMA_HYDRO_IO_H */ diff --git a/src/hydro/MAGMA/hydro_parameters.h b/src/hydro/MAGMA/hydro_parameters.h new file mode 100644 index 0000000000000000000000000000000000000000..afeb90dd792903ee5c37474e1e5561e6b9f50c22 --- /dev/null +++ b/src/hydro/MAGMA/hydro_parameters.h @@ -0,0 +1,194 @@ +/******************************************************************************* + * This file is part of SWIFT. + * Copyright (c) 2019 Josh Borrow (joshua.borrow@durham.ac.uk) + * + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU Lesser General Public License as published + * by the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. + * + ******************************************************************************/ + +#ifndef SWIFT_MAGMA_HYDRO_PARAMETERS_H +#define SWIFT_MAGMA_HYDRO_PARAMETERS_H + +/* Configuration file */ +#include <config.h> + +/* Global headers */ +#if defined(HAVE_HDF5) +#include <hdf5.h> +#endif + +/* Local headers */ +#include "common_io.h" +#include "error.h" +#include "inline.h" + +/** + * @file Minimal/hydro_parameters.h + * @brief Minimal conservative implementation of SPH . (default parameters) + * + * This file defines a number of things that are used in + * hydro_properties.c as defaults for run-time parameters + * as well as a number of compile-time parameters. + */ + +/* Viscosity parameters -- FIXED -- MUST BE DEFINED AT COMPILE-TIME */ + +/* Cosmology default beta=3.0. + * Alpha can be set in the parameter file. + * Beta is defined as in e.g. Price (2010) Eqn (103) */ +#define const_viscosity_alpha 1.0f +#define const_viscosity_beta 2.0f +#define const_viscosity_epsilon 0.1f + +#define const_diffusion_alpha 0.05f + +/* The viscosity that the particles are reset to after being hit by a + * feedback event. This should be set to the same value as the + * hydro_props_default_viscosity_alpha in fixed schemes, and likely + * to hydro_props_default_viscosity_alpha_max in variable schemes. */ +#define hydro_props_default_viscosity_alpha_feedback_reset 0.8f + +/* Viscosity paramaters -- Defaults; can be changed at run-time */ + +/* The "initial" hydro viscosity, or the fixed value for non-variable + * schemes. This usually takes the value 0.8. */ +#define hydro_props_default_viscosity_alpha 0.8f + +/* Structs that store the relevant variables */ + +/*! Artificial viscosity parameters */ +struct viscosity_global_data { + /*! For the fixed, simple case. Also used to set the initial AV + coefficient for variable schemes. */ + float alpha; +}; + +/*! Thermal diffusion parameters */ +struct diffusion_global_data {}; + +/* Functions for reading from parameter file */ + +/* Forward declartions */ +struct swift_params; +struct phys_const; +struct unit_system; + +/* Viscosity */ + +/** + * @brief Initialises the viscosity parameters in the struct from + * the parameter file, or sets them to defaults. + * + * @param params: the pointer to the swift_params file + * @param us: pointer to the internal unit system + * @param phys_const: pointer to the physical constants system + * @param viscosity: pointer to the viscosity_global_data struct to be filled. + **/ +static INLINE void viscosity_init(struct swift_params* params, + const struct unit_system* us, + const struct phys_const* phys_const, + struct viscosity_global_data* viscosity) { + + /* Read the artificial viscosity parameters from the file, if they exist, + * otherwise set them to the defaults defined above. */ + + viscosity->alpha = parser_get_opt_param_float( + params, "SPH:viscosity_alpha", hydro_props_default_viscosity_alpha); +} + +/** + * @brief Initialises a viscosity struct to sensible numbers for mocking + * purposes. + * + * @param viscosity: pointer to the viscosity_global_data struct to be filled. + **/ +static INLINE void viscosity_init_no_hydro( + struct viscosity_global_data* viscosity) { + viscosity->alpha = hydro_props_default_viscosity_alpha; +} + +/** + * @brief Prints out the viscosity parameters at the start of a run. + * + * @param viscosity: pointer to the viscosity_global_data struct found in + * hydro_properties + **/ +static INLINE void viscosity_print( + const struct viscosity_global_data* viscosity) { + message("Artificial viscosity parameters set to alpha: %.3f", + viscosity->alpha); +} + +#if defined(HAVE_HDF5) +/** + * @brief Prints the viscosity information to the snapshot when writing. + * + * @param h_grpsph: the SPH group in the ICs to write attributes to. + * @param viscosity: pointer to the viscosity_global_data struct. + **/ +static INLINE void viscosity_print_snapshot( + hid_t h_grpsph, const struct viscosity_global_data* viscosity) { + + io_write_attribute_f(h_grpsph, "Alpha viscosity", viscosity->alpha); + io_write_attribute_f(h_grpsph, "Beta viscosity", const_viscosity_beta); +} +#endif + +/* Diffusion */ + +/** + * @brief Initialises the diffusion parameters in the struct from + * the parameter file, or sets them to defaults. + * + * @param params: the pointer to the swift_params file + * @param us: pointer to the internal unit system + * @param phys_const: pointer to the physical constants system + * @param diffusion: pointer to the diffusion struct to be filled. + **/ +static INLINE void diffusion_init(struct swift_params* params, + const struct unit_system* us, + const struct phys_const* phys_const, + struct diffusion_global_data* diffusion) {} + +/** + * @brief Initialises a diffusion struct to sensible numbers for mocking + * purposes. + * + * @param diffusion: pointer to the diffusion_global_data struct to be filled. + **/ +static INLINE void diffusion_init_no_hydro( + struct diffusion_global_data* diffusion) {} + +/** + * @brief Prints out the diffusion parameters at the start of a run. + * + * @param diffusion: pointer to the diffusion_global_data struct found in + * hydro_properties + **/ +static INLINE void diffusion_print( + const struct diffusion_global_data* diffusion) {} + +#ifdef HAVE_HDF5 +/** + * @brief Prints the diffusion information to the snapshot when writing. + * + * @param h_grpsph: the SPH group in the ICs to write attributes to. + * @param diffusion: pointer to the diffusion_global_data struct. + **/ +static INLINE void diffusion_print_snapshot( + hid_t h_grpsph, const struct diffusion_global_data* diffusion) {} +#endif + +#endif /* SWIFT_MAGMA_HYDRO_PARAMETERS_H */ diff --git a/src/hydro/MAGMA/hydro_part.h b/src/hydro/MAGMA/hydro_part.h new file mode 100644 index 0000000000000000000000000000000000000000..933cefa53487f709b0e3a81f2603badadae84ba4 --- /dev/null +++ b/src/hydro/MAGMA/hydro_part.h @@ -0,0 +1,261 @@ +/******************************************************************************* + * This file is part of SWIFT. + * Copyright (c) 2016 Matthieu Schaller (schaller@strw.leidenuniv.nl) + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU Lesser General Public License as published + * by the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. + * + ******************************************************************************/ +#ifndef SWIFT_MAGMA_HYDRO_PART_H +#define SWIFT_MAGMA_HYDRO_PART_H + +/** + * @file Minimal/hydro_part.h + * @brief Minimal conservative implementation of SPH (Particle definition) + * + * The thermal variable is the internal energy (u). Simple constant + * viscosity term with the Balsara (1995) switch. No thermal conduction + * term is implemented. + * + * This corresponds to equations (43), (44), (45), (101), (103) and (104) with + * \f$\beta=3\f$ and \f$\alpha_u=0\f$ of Price, D., Journal of Computational + * Physics, 2012, Volume 231, Issue 3, pp. 759-794. + */ + +#include "black_holes_struct.h" +#include "chemistry_struct.h" +#include "cooling_struct.h" +#include "feedback_struct.h" +#include "matrix.h" +#include "mhd_struct.h" +#include "particle_splitting_struct.h" +#include "rt_struct.h" +#include "sink_struct.h" +#include "star_formation_struct.h" +#include "timestep_limiter_struct.h" +#include "tracers_struct.h" + +/** + * @brief Particle fields not needed during the SPH loops over neighbours. + * + * This structure contains the particle fields that are not used in the + * density or force loops. Quantities should be used in the kick, drift and + * potentially ghost tasks only. + */ +struct xpart { + + /*! Offset between current position and position at last tree rebuild. */ + float x_diff[3]; + + /*! Offset between the current position and position at the last sort. */ + float x_diff_sort[3]; + + /*! Velocity at the last full step. */ + float v_full[3]; + + /*! Gravitational acceleration at the end of the last step */ + float a_grav[3]; + + /*! Internal energy at the last full step. */ + float u_full; + + /*! Additional data used to record particle splits */ + struct particle_splitting_data split_data; + + /*! Additional data used to record cooling information */ + struct cooling_xpart_data cooling_data; + + /* Additional data used by the tracers */ + struct tracers_xpart_data tracers_data; + + /* Additional data used by the tracers */ + struct star_formation_xpart_data sf_data; + + /*! Additional data used by the feedback */ + struct feedback_xpart_data feedback_data; + + /*! Additional data used by the MHD scheme */ + struct mhd_xpart_data mhd_data; + +} SWIFT_STRUCT_ALIGN; + +/** + * @brief Particle fields for the SPH particles + * + * The density and force substructures are used to contain variables only used + * within the density and force loops over neighbours. All more permanent + * variables should be declared in the main part of the part structure, + */ +struct part { + + /*! Particle unique ID. */ + long long id; + + /*! Pointer to corresponding gravity part. */ + struct gpart* gpart; + + /*! Particle position. */ + double x[3]; + + /*! Particle predicted velocity. */ + float v[3]; + + /*! Particle acceleration. */ + float a_hydro[3]; + + /*! Particle mass. */ + float mass; + + /*! Particle smoothing length. */ + float h; + + /*! Particle internal energy. */ + float u; + + /*! Time derivative of the internal energy. */ + float u_dt; + + /*! Particle density. */ + float rho; + + /** + * @brief Structure for the variables only used in the density loop over + * neighbours. + * + * Quantities in this sub-structure should only be accessed in the density + * loop over neighbours and the ghost task. + */ + struct { + + /*! Neighbour number count. */ + float wcount; + + /*! Derivative of the neighbour number with respect to h. */ + float wcount_dh; + + /*! Derivative of density with respect to h */ + float rho_dh; + + } density; + + /** + * @brief Structure for the variables only used in the gradient loop over + * neighbours. + * + * Quantities should only be accessed in the gradient loop over neighbours + * and the extra ghost task. + */ + struct { + + /*! The inverse of 'correction matrix' (e.q. 6) - It's symmetric */ + struct sym_matrix c_matrix_inv; + + /*! Gradient of the x-component of the velocity */ + float gradient_vx[3]; + + /*! Gradient of the y-component of the velocity */ + float gradient_vy[3]; + + /*! Gradient of the z-component of the velocity */ + float gradient_vz[3]; + + /*! Gradient of the internal energy */ + float gradient_u[3]; + + } gradient; + + /** + * @brief Structure for the variables only used in the force loop over + * neighbours. + * + * Quantities in this sub-structure should only be accessed in the force + * loop over neighbours, the extra ghost, drift and kick tasks. + */ + struct { + + /*! Particle pressure. */ + float pressure; + + /*! Particle soundspeed. */ + float soundspeed; + + /*! Particle signal velocity */ + float v_sig; + + /*! Time derivative of smoothing length */ + float h_dt; + + /*! The 'correction matrix' (e.q. 6) - It's symmetric */ + struct sym_matrix c_matrix; + + /*! Gradient of the x-component of the velocity */ + float gradient_vx[3]; + + /*! Gradient of the y-component of the velocity */ + float gradient_vy[3]; + + /*! Gradient of the z-component of the velocity */ + float gradient_vz[3]; + + /*! Gradient of the internal energy */ + float gradient_u[3]; + + float f; + + float balsara; + + } force; + + /*! Additional data used by the MHD scheme */ + struct mhd_part_data mhd_data; + + /*! Chemistry information */ + struct chemistry_part_data chemistry_data; + + /*! Cooling information */ + struct cooling_part_data cooling_data; + + /*! Additional data used by the feedback */ + struct feedback_part_data feedback_data; + + /*! Black holes information (e.g. swallowing ID) */ + struct black_holes_part_data black_holes_data; + + /*! Sink information (e.g. swallowing ID) */ + struct sink_part_data sink_data; + + /*! Additional Radiative Transfer Data */ + struct rt_part_data rt_data; + + /*! RT sub-cycling time stepping data */ + struct rt_timestepping_data rt_time_data; + + /*! Time-step length */ + timebin_t time_bin; + + /*! Time-step limiter information */ + struct timestep_limiter_data limiter_data; + +#ifdef SWIFT_DEBUG_CHECKS + + /* Time of the last drift */ + integertime_t ti_drift; + + /* Time of the last kick */ + integertime_t ti_kick; + +#endif + +} SWIFT_STRUCT_ALIGN; + +#endif /* SWIFT_MAGMA_HYDRO_PART_H */ diff --git a/src/hydro_io.h b/src/hydro_io.h index c0f82c2e7b4f52adc870b59f41d12c0e14e84a46..60d2af7963793cfbec090eceecca2e9afc1942f5 100644 --- a/src/hydro_io.h +++ b/src/hydro_io.h @@ -45,6 +45,8 @@ #include "./hydro/Planetary/hydro_io.h" #elif defined(SPHENIX_SPH) #include "./hydro/SPHENIX/hydro_io.h" +#elif defined(MAGMA_SPH) +#include "./hydro/MAGMA/hydro_io.h" #elif defined(GASOLINE_SPH) #include "./hydro/Gasoline/hydro_io.h" #elif defined(ANARCHY_PU_SPH) diff --git a/src/hydro_parameters.h b/src/hydro_parameters.h index 10cdb1b4cdc7f548bcbe2f9744302e63057691b3..0ab0dc63e0b8c316e21a3a1325e42cd738be7858 100644 --- a/src/hydro_parameters.h +++ b/src/hydro_parameters.h @@ -54,6 +54,8 @@ #include "./hydro/Planetary/hydro_parameters.h" #elif defined(SPHENIX_SPH) #include "./hydro/SPHENIX/hydro_parameters.h" +#elif defined(MAGMA_SPH) +#include "./hydro/MAGMA/hydro_parameters.h" #elif defined(GASOLINE_SPH) #include "./hydro/Gasoline/hydro_parameters.h" #elif defined(ANARCHY_PU_SPH) diff --git a/src/lightcone/lightcone_particle_io.c b/src/lightcone/lightcone_particle_io.c index 9d3584412a56ca60e077fb9f3d8d59e24f7276a7..7b1ab838c0155d27338c22b71682fdaf4b06df03 100644 --- a/src/lightcone/lightcone_particle_io.c +++ b/src/lightcone/lightcone_particle_io.c @@ -642,34 +642,32 @@ void append_dataset(const struct unit_system *snapshot_units, hid_t mem_type_id, hsize_t chunk_size, int lossy_compression, enum lossy_compression_schemes compression_scheme, - int gzip_level, const int rank, const hsize_t *dims, + int gzip_level, const int rank, const hsize_t dims[2], const hsize_t num_written, const void *data) { - const int max_rank = 2; - if (rank > max_rank) - error("HDF5 dataset has too may dimensions. Increase max_rank."); + if (rank > 2) error("HDF5 dataset has too may dimensions."); if (rank < 1) error("HDF5 dataset must be at least one dimensional"); /* If we have zero elements to append, there's nothing to do */ if (dims[0] == 0) return; /* Determine size of the dataset after we append our data */ - hsize_t full_dims[max_rank]; + hsize_t full_dims[rank]; for (int i = 0; i < rank; i += 1) full_dims[i] = dims[i]; full_dims[0] += num_written; /* Determine maximum size in each dimension */ - hsize_t max_dims[max_rank]; + hsize_t max_dims[rank]; for (int i = 1; i < rank; i += 1) max_dims[i] = full_dims[i]; max_dims[0] = H5S_UNLIMITED; /* Determine chunk size in each dimension */ - hsize_t chunk_dims[max_rank]; + hsize_t chunk_dims[rank]; for (int i = 1; i < rank; i += 1) chunk_dims[i] = full_dims[i]; chunk_dims[0] = (hsize_t)chunk_size; /* Find offset to region to write in each dimension */ - hsize_t offset[max_rank]; + hsize_t offset[rank]; for (int i = 1; i < rank; i += 1) offset[i] = 0; offset[0] = num_written; diff --git a/src/matrix.h b/src/matrix.h new file mode 100644 index 0000000000000000000000000000000000000000..c0d217fe84ca140dc743bec4749e58a8065fe602 --- /dev/null +++ b/src/matrix.h @@ -0,0 +1,90 @@ +/******************************************************************************* + * This file is part of SWIFT. + * Copyright (c) 2023 Matthieu Schaller (schaller@strw.leidenuniv.nl) + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU Lesser General Public License as published + * by the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. + * + ******************************************************************************/ +#ifndef SWIFT_MATRIX_H +#define SWIFT_MATRIX_H + +/* Local includes */ +#include "error.h" + +#if !defined(HYDRO_DIMENSION_3D) +#error "Not yet defined!" +#endif + +struct sym_matrix { + + union { + struct { + float elements[6]; + }; + struct { + float xx; + float yy; + float zz; + float xy; + float xz; + float yz; + }; + }; +}; + +__attribute__((always_inline)) INLINE static void zero_sym_matrix( + struct sym_matrix *M) { + for (int i = 0; i < 6; ++i) M->elements[i] = 0.f; +} + +__attribute__((always_inline)) INLINE static void get_matrix_from_sym_matrix( + float out[3][3], const struct sym_matrix *in) { + + out[0][0] = in->xx; + out[0][1] = in->xy; + out[0][2] = in->xz; + out[1][0] = in->xy; + out[1][1] = in->yy; + out[1][2] = in->yz; + out[2][0] = in->xz; + out[2][1] = in->yz; + out[2][2] = in->zz; +} + +__attribute__((always_inline)) INLINE static void get_sym_matrix_from_matrix( + struct sym_matrix *out, const float in[3][3]) { + out->xx = in[0][0]; + out->yy = in[1][1]; + out->zz = in[2][2]; + out->xy = in[0][1]; + out->xz = in[0][2]; + out->yz = in[1][2]; +} + +__attribute__((always_inline)) INLINE static void sym_matrix_multiply_by_vector( + float out[3], const struct sym_matrix *M, const float v[3]) { + + out[0] = M->xx * v[0] + M->xy * v[1] + M->xz * v[2]; + out[1] = M->xy * v[0] + M->yy * v[1] + M->yz * v[2]; + out[2] = M->xz * v[0] + M->yz * v[1] + M->zz * v[2]; +} + +__attribute__((always_inline)) INLINE static void sym_matrix_print( + const struct sym_matrix *M) { + message("|%.4f %.4f %.4f|", M->xx, M->xy, M->xz); + message("|%.4f %.4f %.4f|", M->xy, M->yy, M->yz); + message("|%.4f %.4f %.4f|", M->xz, M->yz, M->zz); +} + +#endif /* SWIFT_MATRIX_H */ diff --git a/src/mhd/DInduction/mhd.h b/src/mhd/DInduction/mhd.h index c3d2f4272ccf0a1f754ef981f1fad944ad67a0f9..9cf201c6d58dcabef04a7154814c71a7f28e9804 100644 --- a/src/mhd/DInduction/mhd.h +++ b/src/mhd/DInduction/mhd.h @@ -115,17 +115,19 @@ __attribute__((always_inline)) INLINE static float mhd_compute_timestep( const struct hydro_props *hydro_properties, const struct cosmology *cosmo, const float mu_0) { + const float afac_divB = pow(cosmo->a, -mhd_comoving_factor - 0.5f); + const float afac_resistive = cosmo->a * cosmo->a; /* Dt from 1/DivOperator(Alfven speed) */ float dt_divB = p->mhd_data.divB != 0.0f - ? cosmo->a * hydro_properties->CFL_condition * + ? afac_divB * hydro_properties->CFL_condition * sqrtf(p->rho / (p->mhd_data.divB * p->mhd_data.divB) * mu_0) : FLT_MAX; const float resistive_eta = p->mhd_data.resistive_eta; const float dt_eta = resistive_eta != 0.0f - ? cosmo->a * hydro_properties->CFL_condition * p->h * - p->h / resistive_eta * 0.5 + ? afac_resistive * hydro_properties->CFL_condition * + p->h * p->h / resistive_eta * 0.5 : FLT_MAX; return fminf(dt_eta, dt_divB); @@ -542,6 +544,14 @@ __attribute__((always_inline)) INLINE static void mhd_convert_quantities( const struct hydro_props *hydro_props) { /* Set Restitivity Eta */ p->mhd_data.resistive_eta = hydro_props->mhd.mhd_eta; + + p->mhd_data.BPred[0] *= pow(cosmo->a, -mhd_comoving_factor); + p->mhd_data.BPred[1] *= pow(cosmo->a, -mhd_comoving_factor); + p->mhd_data.BPred[2] *= pow(cosmo->a, -mhd_comoving_factor); + + xp->mhd_data.Bfld_full[0] = p->mhd_data.BPred[0]; + xp->mhd_data.Bfld_full[1] = p->mhd_data.BPred[1]; + xp->mhd_data.Bfld_full[2] = p->mhd_data.BPred[2]; } /** @@ -558,22 +568,6 @@ __attribute__((always_inline)) INLINE static void mhd_first_init_part( struct part *restrict p, struct xpart *restrict xp, const struct mhd_global_data *mhd_data, const double Lsize) { - const float define_Bfield_in_ics = mhd_data->define_Bfield_in_ics; - const float Nvort = 10; - const float Bini = define_Bfield_in_ics; - if (define_Bfield_in_ics) { - p->mhd_data.BPred[0] = Bini * (sin(2 * M_PI * p->x[2] / Lsize * Nvort) + - cos(2 * M_PI * p->x[1] / Lsize * Nvort)); - p->mhd_data.BPred[1] = Bini * (sin(2 * M_PI * p->x[0] / Lsize * Nvort) + - cos(2 * M_PI * p->x[2] / Lsize * Nvort)); - p->mhd_data.BPred[2] = Bini * (sin(2 * M_PI * p->x[1] / Lsize * Nvort) + - cos(2 * M_PI * p->x[0] / Lsize * Nvort)); - } - - xp->mhd_data.Bfld_full[0] = p->mhd_data.BPred[0]; - xp->mhd_data.Bfld_full[1] = p->mhd_data.BPred[1]; - xp->mhd_data.Bfld_full[2] = p->mhd_data.BPred[2]; - p->mhd_data.phi = 0.f; mhd_reset_acceleration(p); diff --git a/src/mhd/DInduction/mhd_parameters.h b/src/mhd/DInduction/mhd_parameters.h index 9a48b731255cd9de482706c06ac461c3afc67422..e8969ad4912df2c5a0f8b50f8c3f7d5749323a35 100644 --- a/src/mhd/DInduction/mhd_parameters.h +++ b/src/mhd/DInduction/mhd_parameters.h @@ -74,7 +74,6 @@ struct mhd_global_data { float hyp_dedner; float par_dedner; float mhd_eta; - float define_Bfield_in_ics; }; /* Functions for reading from parameter file */ @@ -106,14 +105,6 @@ static INLINE void mhd_init(struct swift_params* params, mhd->hyp_dedner = parser_get_param_float(params, "MHD:hyperbolic_dedner"); mhd->par_dedner = parser_get_param_float(params, "MHD:parabolic_dedner"); mhd->mhd_eta = parser_get_param_float(params, "MHD:resistive_eta"); - mhd->define_Bfield_in_ics = - parser_get_opt_param_float(params, "MHD:define_B_in_ics", 0.f); - // calculate the comoving seed field - if (mhd->define_Bfield_in_ics != 0.f) { - float a_beg = parser_get_param_float(params, "Cosmology:a_begin"); - mhd->define_Bfield_in_ics = - mhd->define_Bfield_in_ics * pow(a_beg, -mhd_comoving_factor); - } } /** @@ -138,9 +129,6 @@ static INLINE void mhd_print(const struct mhd_global_data* mhd) { message("Dedner Hyperbolic/Parabolic: %.3f, %.3f ", mhd->hyp_dedner, mhd->par_dedner); message("MHD global dissipation Eta: %.3f", mhd->mhd_eta); - if (mhd->define_Bfield_in_ics) - message("Starting with a Initial co-moving Bfield: %4.3e Gauss", - mhd->define_Bfield_in_ics); } #if defined(HAVE_HDF5) @@ -158,8 +146,6 @@ static INLINE void mhd_print_snapshot(hid_t h_grpsph, io_write_attribute_f(h_grpsph, "Dedner Parabolic Constant", mhd_data->par_dedner); io_write_attribute_f(h_grpsph, "Resistive Eta", mhd_data->mhd_eta); - io_write_attribute_f(h_grpsph, "Generate comoving BField in ICs", - mhd_data->define_Bfield_in_ics); io_write_attribute_f(h_grpsph, "Comoving exponent", mhd_comoving_factor); } #endif diff --git a/src/mhd/DirectInduction/mhd.h b/src/mhd/DirectInduction/mhd.h index 1414c75a83a47be55ec6e2e516b7f1f9762a7d14..a69beaa5344338a9e289170ca3c0958a88107228 100644 --- a/src/mhd/DirectInduction/mhd.h +++ b/src/mhd/DirectInduction/mhd.h @@ -134,12 +134,14 @@ __attribute__((always_inline)) INLINE static float mhd_compute_timestep( const float dt_B_derivatives = dt_B_factor != 0.f - ? hydro_properties->CFL_condition * + ? hydro_properties->CFL_condition * cosmo->a / + cosmo->a_factor_sound_speed * sqrtf(p->rho * mu_0 / (dt_B_factor * dt_B_factor)) : FLT_MAX; const float dt_eta = p->mhd_data.resistive_eta != 0.f - ? hydro_properties->CFL_condition * p->h * p->h / + ? hydro_properties->CFL_condition * cosmo->a * + cosmo->a * p->h * p->h / p->mhd_data.resistive_eta : FLT_MAX; @@ -230,10 +232,17 @@ __attribute__((always_inline)) INLINE static float mhd_signal_velocity( const struct part *restrict pj, const float mu_ij, const float beta, const float a, const float mu_0) { + const float v_sigi = + mhd_get_magnetosonic_speed(pi, a, mu_0); + const float v_sigj = + mhd_get_magnetosonic_speed(pj, a, mu_0); + + /* const float v_sigi = mhd_get_fast_magnetosonic_wave_phase_velocity(dx, pi, a, mu_0); const float v_sigj = mhd_get_fast_magnetosonic_wave_phase_velocity(dx, pj, a, mu_0); + */ const float v_sig = v_sigi + v_sigj - beta * mu_ij; @@ -386,8 +395,8 @@ __attribute__((always_inline)) INLINE static void mhd_prepare_force( * @param a The current value of the cosmological scale factor */ __attribute__((always_inline)) INLINE static float mhd_get_psi_over_ch_dt( - struct part *p, const float a, const struct hydro_props *hydro_props, - const float mu_0) { + struct part *p, const float a, const float a_factor_sound_speed, + const float H, const struct hydro_props *hydro_props, const float mu_0) { /* Retrieve inverse of smoothing length. */ const float h = p->h; @@ -402,11 +411,17 @@ __attribute__((always_inline)) INLINE static float mhd_get_psi_over_ch_dt( const float par = hydro_props->mhd.par_dedner; const float div_B = p->mhd_data.divB; - const float div_v = hydro_get_div_v(p); + const float div_v = a * a * hydro_get_div_v(p) - 3.0f * a * a * H; const float psi_over_ch = p->mhd_data.psi_over_ch; - return -hyp * ch * div_B - hyp_divv * psi_over_ch * div_v - - par * psi_over_ch * ch * h_inv; + const float hyperbolic_term = + -hyp * a * a * a_factor_sound_speed * a_factor_sound_speed * ch * div_B; + const float hyperbolic_divv_term = -hyp_divv * psi_over_ch * div_v; + const float parabolic_term = + -par * a * a_factor_sound_speed * psi_over_ch * ch * h_inv; + const float Hubble_term = a * a * H * psi_over_ch; + + return hyperbolic_term + hyperbolic_divv_term + parabolic_term + Hubble_term; } /** @@ -469,10 +484,12 @@ __attribute__((always_inline)) INLINE static void mhd_predict_extra( p->mhd_data.B_over_rho[1] += p->mhd_data.B_over_rho_dt[1] * dt_therm; p->mhd_data.B_over_rho[2] += p->mhd_data.B_over_rho_dt[2] * dt_therm; - p->mhd_data.psi_over_ch_dt = - mhd_get_psi_over_ch_dt(p, cosmo->a, hydro_props, mu_0); + p->mhd_data.psi_over_ch_dt = mhd_get_psi_over_ch_dt( + p, cosmo->a, cosmo->a_factor_sound_speed, cosmo->H, hydro_props, mu_0); p->mhd_data.psi_over_ch += - mhd_get_psi_over_ch_dt(p, cosmo->a, hydro_props, mu_0) * dt_therm; + mhd_get_psi_over_ch_dt(p, cosmo->a, cosmo->a_factor_sound_speed, cosmo->H, + hydro_props, mu_0) * + dt_therm; } /** @@ -489,7 +506,19 @@ __attribute__((always_inline)) INLINE static void mhd_predict_extra( */ __attribute__((always_inline)) INLINE static void mhd_end_force( struct part *p, const struct cosmology *cosmo, - const struct hydro_props *hydro_props, const float mu_0) {} + const struct hydro_props *hydro_props, const float mu_0) { + + /* Hubble expansion contribution to induction equation */ + + const float Hubble_induction_pref = + cosmo->a * cosmo->a * cosmo->H * (1.5f * hydro_gamma - 2.f); + p->mhd_data.B_over_rho_dt[0] += + Hubble_induction_pref * p->mhd_data.B_over_rho[0]; + p->mhd_data.B_over_rho_dt[1] += + Hubble_induction_pref * p->mhd_data.B_over_rho[1]; + p->mhd_data.B_over_rho_dt[2] += + Hubble_induction_pref * p->mhd_data.B_over_rho[2]; +} /** * @brief Kick the additional variables @@ -553,6 +582,11 @@ __attribute__((always_inline)) INLINE static void mhd_convert_quantities( p->mhd_data.B_over_rho[1] /= p->rho; p->mhd_data.B_over_rho[2] /= p->rho; + /* Convert to co-moving B/rho */ + p->mhd_data.B_over_rho[0] *= pow(cosmo->a, 1.5f * hydro_gamma); + p->mhd_data.B_over_rho[1] *= pow(cosmo->a, 1.5f * hydro_gamma); + p->mhd_data.B_over_rho[2] *= pow(cosmo->a, 1.5f * hydro_gamma); + xp->mhd_data.B_over_rho_full[0] = p->mhd_data.B_over_rho[0]; xp->mhd_data.B_over_rho_full[1] = p->mhd_data.B_over_rho[1]; xp->mhd_data.B_over_rho_full[2] = p->mhd_data.B_over_rho[2]; diff --git a/src/mhd/DirectInduction/mhd_iact.h b/src/mhd/DirectInduction/mhd_iact.h index 9530cab1ceb9cb93a64fc791b31b34e3c029bc3e..ec101145d904cc940b904796fc082bd4521b0588 100644 --- a/src/mhd/DirectInduction/mhd_iact.h +++ b/src/mhd/DirectInduction/mhd_iact.h @@ -606,7 +606,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_mhd_force( const float corr_ratio_i = fabsf(normBi * psi_over_ch_i_inv); const float corr_ratio_j = fabsf(normBj * psi_over_ch_j_inv); - + const float Qi = corr_ratio_i < 2 ? 0.5 * corr_ratio_i : 1.0f; const float Qj = corr_ratio_j < 2 ? 0.5 * corr_ratio_j : 1.0f; @@ -897,7 +897,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_mhd_force( psi_over_ch_i != 0.f ? 1.f / psi_over_ch_i : 0.; const float corr_ratio_i = fabsf(normBi * psi_over_ch_i_inv); - + const float Qi = corr_ratio_i < 2 ? 0.5 * corr_ratio_i : 1.0f; pi->mhd_data.B_over_rho_dt[0] -= mj * Qi * grad_psi_i * dx[0]; diff --git a/src/mhd/DirectInduction/mhd_io.h b/src/mhd/DirectInduction/mhd_io.h index 367d25d7b17b2805a71d824761231119ebf4b37a..30ccb7d30e1c919678e52d6a20cb443110bf4f9c 100644 --- a/src/mhd/DirectInduction/mhd_io.h +++ b/src/mhd/DirectInduction/mhd_io.h @@ -19,6 +19,8 @@ #ifndef SWIFT_DIRECT_INDUCTION_MHD_IO_H #define SWIFT_DIRECT_INDUCTION_MHD_IO_H +#include "adiabatic_index.h" + /** * @brief Specifies which particle fields to read from a dataset * @@ -57,16 +59,19 @@ INLINE static int mhd_write_particles(const struct part* parts, struct io_props* list) { list[0] = io_make_output_field_convert_part( - "MagneticFluxDensities", FLOAT, 3, UNIT_CONV_MAGNETIC_FIELD, 1.f, parts, - xparts, convert_B, "Magnetic flux densities of the particles"); + "MagneticFluxDensities", FLOAT, 3, UNIT_CONV_MAGNETIC_FIELD, + -1.5f * hydro_gamma, parts, xparts, convert_B, + "Magnetic flux densities of the particles"); list[1] = io_make_output_field( - "MagneticDivergences", FLOAT, 1, UNIT_CONV_MAGNETIC_DIVERGENCE, 1.f, - parts, mhd_data.divB, "co-moving DivB of the particle"); + "MagneticDivergences", FLOAT, 1, UNIT_CONV_MAGNETIC_DIVERGENCE, + -1.5f * hydro_gamma - 1.f, parts, mhd_data.divB, + "co-moving DivB of the particle"); list[2] = io_make_output_field( - "DednerScalars", FLOAT, 1, UNIT_CONV_ELECTRIC_CHARGE_FIELD_STRENGTH, 1.f, - parts, mhd_data.psi_over_ch, "Dedner scalar associated to the particle"); + "DednerScalars", FLOAT, 1, UNIT_CONV_ELECTRIC_CHARGE_FIELD_STRENGTH, + -1.5f * hydro_gamma - 1.f, parts, mhd_data.psi_over_ch, + "Dedner scalar associated to the particle"); list[3] = io_make_output_field( "DednerScalarsdt", FLOAT, 1, @@ -80,8 +85,9 @@ INLINE static int mhd_write_particles(const struct part* parts, "Time derivative of Magnetic flux densities of the particles"); list[5] = io_make_output_field( - "MagneticFluxCurl", FLOAT, 3, UNIT_CONV_MAGNETIC_CURL, 1.f, parts, - mhd_data.curl_B, "The curl of Magnetic flux densities of the particles"); + "MagneticFluxCurl", FLOAT, 3, UNIT_CONV_MAGNETIC_CURL, + -1.5f * hydro_gamma - 1.f, parts, mhd_data.curl_B, + "The curl of Magnetic flux densities of the particles"); list[6] = io_make_output_field( "AlphaAR", FLOAT, 1, UNIT_CONV_NO_UNITS, 1.f, parts, mhd_data.alpha_AR, @@ -99,7 +105,7 @@ INLINE static void mhd_write_flavour(hid_t h_grpsph) { io_write_attribute_s( h_grpsph, "MHD Flavour", "Orestis - Direct Induction, divB Subtraction, " - "Artificial Resistivity & Dedner cleaning. Price et al. (2018)."); + "Artificial Resistivity & Dedner Cleaning. Price et al. (2018)."); } #endif /* SWIFT_DIRECT_INDUCTION_MHD_IO_H */ diff --git a/src/mhd/DirectInduction/mhd_parameters.h b/src/mhd/DirectInduction/mhd_parameters.h index a0d2ba4e9ca710a53640c8e9de17c720d2194268..e10aace0eb9041cfa5d7e6fc323a88a47bfc6b61 100644 --- a/src/mhd/DirectInduction/mhd_parameters.h +++ b/src/mhd/DirectInduction/mhd_parameters.h @@ -84,7 +84,6 @@ struct mhd_global_data { float hyp_dedner_divv; float par_dedner; float mhd_eta; - float define_Bfield_in_ics; }; /* Functions for reading from parameter file */ @@ -122,8 +121,6 @@ static INLINE void mhd_init(struct swift_params* params, parser_get_param_float(params, "MHD:artificial_diffusion"); mhd->hyp_dedner_divv = parser_get_param_float(params, "MHD:hyperbolic_dedner_divv"); - mhd->define_Bfield_in_ics = - parser_get_opt_param_float(params, "MHD:define_B_in_ics", 0.f); } /** @@ -147,15 +144,10 @@ static INLINE void mhd_print(const struct mhd_global_data* mhd) { message("MHD tensile instability correction prefactor: %.3f ", mhd->monopole_subtraction); - message("Artificial resistivity: %.3f ", mhd->art_diffusion); + message("Artificial diffusion: %.3f ", mhd->art_diffusion); message("Dedner Hyperbolic/Hyperbolic div(v)/Parabolic: %.3f, %.3f, %.3f ", mhd->hyp_dedner, mhd->hyp_dedner_divv, mhd->par_dedner); - message("MHD global dissipation Eta: %.3f", mhd->mhd_eta); - if (mhd->define_Bfield_in_ics) - message( - "NOT IMPLEMENTED YET!Starting with a Initial co-moving Bfield: %4.3e " - "Gauss", - mhd->define_Bfield_in_ics); + message("MHD global Resistive Eta: %.3f", mhd->mhd_eta); } #if defined(HAVE_HDF5) @@ -178,8 +170,6 @@ static INLINE void mhd_print_snapshot(hid_t h_grpsph, io_write_attribute_f(h_grpsph, "Dedner Parabolic Constant", mhd_data->par_dedner); io_write_attribute_f(h_grpsph, "Resistive Eta", mhd_data->mhd_eta); - // io_write_attribute_f(h_grpsph, "Generate comoving BField in ICs", - // mhd_data->define_Bfield_in_ics); // io_write_attribute_f(h_grpsph, "Comoving exponent", mhd_comoving_factor); } #endif diff --git a/src/mhd/None/mhd.h b/src/mhd/None/mhd.h index 4147f48d7b23ab19e62bf38de87c1b720afa74a6..bb84990e64001b88a2db06331db1070e331f1ead 100644 --- a/src/mhd/None/mhd.h +++ b/src/mhd/None/mhd.h @@ -27,10 +27,9 @@ */ __attribute__((always_inline)) INLINE static float mhd_get_magnetic_energy( const struct part *p, const struct xpart *xp, const float mu_0) { - - error("Calling mhd_get_magnetic_energy when compiling without MHD!"); return 0.f; } + /** * @brief Returns the magnetic field squared contained in the particle. * @@ -40,7 +39,6 @@ __attribute__((always_inline)) INLINE static float mhd_get_magnetic_energy( __attribute__((always_inline)) INLINE static float mhd_get_Bms( const struct part *p, const struct xpart *xp) { - error("Calling mhd_get_Bms when compiling without MHD!"); return 0.f; } @@ -66,7 +64,6 @@ __attribute__((always_inline)) INLINE static float mhd_get_magnetic_divergence( __attribute__((always_inline)) INLINE static float mhd_get_magnetic_helicity( const struct part *p, const struct xpart *xp) { - error("Calling mhd_get_magnetic_helicity when compiling without MHD!"); return 0.f; } @@ -79,7 +76,6 @@ __attribute__((always_inline)) INLINE static float mhd_get_magnetic_helicity( __attribute__((always_inline)) INLINE static float mhd_get_cross_helicity( const struct part *p, const struct xpart *xp) { - error("Calling mhd_get_cross_helicity when compiling without MHD!"); return 0.f; } @@ -94,8 +90,7 @@ __attribute__((always_inline)) INLINE static float mhd_get_cross_helicity( __attribute__((always_inline)) INLINE static float mhd_get_divB_error( const struct part *p, const struct xpart *xp) { - error("Calling mhd_get_fast_divB_error when compiling without MHD!"); - return -1.f; + return 0.f; } /** diff --git a/src/mhd/VPotential/mhd.h b/src/mhd/VPotential/mhd.h index fe390a7b3723a83291bb6d70a4725a9b8e769ead..814c26ffb227e7cb0aaf6f26faf1395954c61863 100644 --- a/src/mhd/VPotential/mhd.h +++ b/src/mhd/VPotential/mhd.h @@ -121,17 +121,19 @@ __attribute__((always_inline)) INLINE static float mhd_compute_timestep( const struct hydro_props *hydro_properties, const struct cosmology *cosmo, const float mu_0) { + const float afac_divB = pow(cosmo->a, -mhd_comoving_factor - 0.5f); + const float afac_resistive = cosmo->a * cosmo->a; /* Dt from 1/DivOperator(Alfven speed) */ float dt_divB = p->mhd_data.divB != 0.0f - ? cosmo->a * hydro_properties->CFL_condition * + ? afac_divB * hydro_properties->CFL_condition * sqrtf(p->rho / (p->mhd_data.divB * p->mhd_data.divB) * mu_0) : FLT_MAX; const float resistive_eta = p->mhd_data.resistive_eta; const float dt_eta = resistive_eta != 0.0f - ? cosmo->a * hydro_properties->CFL_condition * p->h * - p->h / resistive_eta * 0.5 + ? afac_resistive * hydro_properties->CFL_condition * + p->h * p->h / resistive_eta * 0.5 : FLT_MAX; return fminf(dt_eta, dt_divB); @@ -270,12 +272,12 @@ __attribute__((always_inline)) INLINE static float hydro_get_dGau_dt( const struct cosmology *c) { const float v_sig = hydro_get_signal_velocity(p); - const float afac1 = pow(c->a, 2.f * mhd_comoving_factor); - const float afac2 = pow(c->a, 1.f + mhd_comoving_factor); + const float afac1 = pow(c->a, 2.f * mhd_comoving_factor - 1.f); + const float afac2 = pow(c->a, mhd_comoving_factor); return (-p->mhd_data.divA * v_sig * v_sig * 0.1 * afac1 - 2.0f * v_sig * Gauge / p->h * afac2 - - (2.f + mhd_comoving_factor) * c->a * c->a * c->H * Gauge) / + (1.f + mhd_comoving_factor) * c->a * c->a * c->H * Gauge) / 2.f; } @@ -571,6 +573,14 @@ __attribute__((always_inline)) INLINE static void mhd_convert_quantities( const struct hydro_props *hydro_props) { /* Set Restitivity Eta */ p->mhd_data.resistive_eta = hydro_props->mhd.mhd_eta; + + p->mhd_data.APred[0] *= pow(cosmo->a, -mhd_comoving_factor + 1.f); + p->mhd_data.APred[1] *= pow(cosmo->a, -mhd_comoving_factor + 1.f); + p->mhd_data.APred[2] *= pow(cosmo->a, -mhd_comoving_factor + 1.f); + + xp->mhd_data.APot_full[0] = p->mhd_data.APred[0]; + xp->mhd_data.APot_full[1] = p->mhd_data.APred[1]; + xp->mhd_data.APot_full[2] = p->mhd_data.APred[2]; } /** @@ -586,39 +596,7 @@ __attribute__((always_inline)) INLINE static void mhd_convert_quantities( __attribute__((always_inline)) INLINE static void mhd_first_init_part( struct part *restrict p, struct xpart *restrict xp, const struct mhd_global_data *mhd_data, const double Lsize) { - - // const float Lsize = s->dims[0]; - const float define_Bfield_in_ics = mhd_data->define_Bfield_in_ics; - const float define_Afield_in_ics = mhd_data->define_Afield_in_ics; - const float Nvort = 10; - const float Bini = define_Bfield_in_ics; - const float Aini = define_Afield_in_ics / (2 * M_PI * Nvort) * Lsize; - if (define_Bfield_in_ics) { - p->mhd_data.BPred[0] = Bini * (sin(2 * M_PI * p->x[2] / Lsize * Nvort) + - cos(2 * M_PI * p->x[1] / Lsize * Nvort)); - p->mhd_data.BPred[1] = Bini * (sin(2 * M_PI * p->x[0] / Lsize * Nvort) + - cos(2 * M_PI * p->x[2] / Lsize * Nvort)); - p->mhd_data.BPred[2] = Bini * (sin(2 * M_PI * p->x[1] / Lsize * Nvort) + - cos(2 * M_PI * p->x[0] / Lsize * Nvort)); - p->mhd_data.APred[0] = Aini * (sin(2 * M_PI * p->x[2] / Lsize * Nvort) + - cos(2 * M_PI * p->x[1] / Lsize * Nvort)); - p->mhd_data.APred[1] = Aini * (sin(2 * M_PI * p->x[0] / Lsize * Nvort) + - cos(2 * M_PI * p->x[2] / Lsize * Nvort)); - p->mhd_data.APred[2] = Aini * (sin(2 * M_PI * p->x[1] / Lsize * Nvort) + - cos(2 * M_PI * p->x[0] / Lsize * Nvort)); - } - - // p->mhd_data.BPred[0] = 0.0; - // p->mhd_data.BPred[1] = 0.0; - // p->mhd_data.BPred[2] = 0.0; - - // p->mhd_data.APred[0] = 0.0; - // p->mhd_data.APred[1] = 0.0; - // p->mhd_data.APred[2] = 0.0; - xp->mhd_data.APot_full[0] = p->mhd_data.APred[0]; - xp->mhd_data.APot_full[1] = p->mhd_data.APred[1]; - xp->mhd_data.APot_full[2] = p->mhd_data.APred[2]; - xp->mhd_data.Gau_full = 0.0f * p->mhd_data.Gau; + xp->mhd_data.Gau_full = 0.0f; p->mhd_data.divB = 0.0f; mhd_reset_acceleration(p); diff --git a/src/mhd/VPotential/mhd_parameters.h b/src/mhd/VPotential/mhd_parameters.h index f98bffc4eeb6d7b96c469b2b038d93d438d9bd29..c3d4c1ce5c1da2de0394a2615f3820412c1e1626 100644 --- a/src/mhd/VPotential/mhd_parameters.h +++ b/src/mhd/VPotential/mhd_parameters.h @@ -64,8 +64,6 @@ struct mhd_global_data { /*! For the fixed, simple case of direct induction. */ float mhd_eta; - float define_Bfield_in_ics; - float define_Afield_in_ics; }; /* Functions for reading from parameter file */ @@ -95,18 +93,6 @@ static INLINE void mhd_init(struct swift_params* params, * otherwise set them to the defaults defined above. */ mhd->mhd_eta = parser_get_param_float(params, "MHD:resistive_eta"); - mhd->define_Bfield_in_ics = - parser_get_opt_param_float(params, "MHD:define_B_in_ics", 0.f); - // calculate the comoving seed field - if (mhd->define_Bfield_in_ics != 0.f) { - float a_beg = parser_get_param_float(params, "Cosmology:a_begin"); - mhd->define_Afield_in_ics = - mhd->define_Bfield_in_ics * pow(a_beg, -mhd_comoving_factor); - mhd->define_Bfield_in_ics = - mhd->define_Bfield_in_ics * pow(a_beg, -mhd_comoving_factor); - } else - mhd->define_Afield_in_ics = 0.f; - // mhd->define_Afield_in_ics = mhd->define_Bfield_in_ics; } /** @@ -129,9 +115,6 @@ static INLINE void mhd_init(struct swift_params* params, static INLINE void mhd_print(const struct mhd_global_data* mhd) { message("MHD global dissipation Eta: %.3f", mhd->mhd_eta); - if (mhd->define_Bfield_in_ics) - message("Starting with a Initial co-moving Bfield: %4.3e Gauss", - mhd->define_Bfield_in_ics); } #if defined(HAVE_HDF5) @@ -145,8 +128,6 @@ static INLINE void mhd_print_snapshot(hid_t h_grpsph, const struct mhd_global_data* mhd_data) { io_write_attribute_f(h_grpsph, "Resistive Eta", mhd_data->mhd_eta); - io_write_attribute_f(h_grpsph, "Generate comoving BField in ICs", - mhd_data->define_Bfield_in_ics); io_write_attribute_f(h_grpsph, "Comoving exponent", mhd_comoving_factor); } #endif diff --git a/src/part.h b/src/part.h index dce15ece9df269b51f5ceebefa32e3f9b52abc0b..de8ab7ea39a99d9d72c4b1539fb0cb4062420107 100644 --- a/src/part.h +++ b/src/part.h @@ -84,6 +84,10 @@ struct threadpool; #include "./hydro/SPHENIX/hydro_part.h" #define hydro_need_extra_init_loop 0 #define EXTRA_HYDRO_LOOP +#elif defined(MAGMA_SPH) +#include "./hydro/MAGMA/hydro_part.h" +#define hydro_need_extra_init_loop 0 +#define EXTRA_HYDRO_LOOP #elif defined(GASOLINE_SPH) #include "./hydro/Gasoline/hydro_part.h" #define hydro_need_extra_init_loop 0 diff --git a/src/rt/GEAR/rt_thermochemistry.h b/src/rt/GEAR/rt_thermochemistry.h index 273a240e3d7ed0ec89f6e6928224ff5d25451409..40aafd2b91f7e7f7c4a9a628d93c3b354bbcc352 100644 --- a/src/rt/GEAR/rt_thermochemistry.h +++ b/src/rt/GEAR/rt_thermochemistry.h @@ -100,7 +100,7 @@ __attribute__((always_inline)) INLINE static void rt_tchem_first_init_part( * @param phys_const The physical constants in internal units. * @param us The internal system of units. * @param dt The time-step of this particle. - * @depth recursion depth + * @param depth recursion depth */ INLINE static void rt_do_thermochemistry( struct part* restrict p, struct xpart* restrict xp, diff --git a/tests/test125cells.c b/tests/test125cells.c index 9d2c3638668d5714478fa3a36d95ed24793a235c..00b6b2ac6efb7cfbbfd3b99dce6e7bc4c647940d 100644 --- a/tests/test125cells.c +++ b/tests/test125cells.c @@ -111,9 +111,10 @@ void set_energy_state(struct part *part, enum pressure_field press, float size, part->entropy = pressure / pow_gamma(density); #elif defined(PHANTOM_SPH) part->u = pressure / (hydro_gamma_minus_one * density); -#elif defined(MINIMAL_SPH) || defined(HOPKINS_PU_SPH) || \ - defined(HOPKINS_PU_SPH_MONAGHAN) || defined(ANARCHY_PU_SPH) || \ - defined(SPHENIX_SPH) || defined(PHANTOM_SPH) || defined(GASOLINE_SPH) +#elif defined(MINIMAL_SPH) || defined(HOPKINS_PU_SPH) || \ + defined(HOPKINS_PU_SPH_MONAGHAN) || defined(ANARCHY_PU_SPH) || \ + defined(SPHENIX_SPH) || defined(PHANTOM_SPH) || defined(GASOLINE_SPH) || \ + defined(MAGMA_SPH) part->u = pressure / (hydro_gamma_minus_one * density); #elif defined(PLANETARY_SPH) part->u = pressure / (hydro_gamma_minus_one * density); @@ -390,10 +391,11 @@ void dump_particle_fields(char *fileName, struct cell *main_cell, main_cell->hydro.parts[pid].v[0], main_cell->hydro.parts[pid].v[1], main_cell->hydro.parts[pid].v[2], main_cell->hydro.parts[pid].h, hydro_get_comoving_density(&main_cell->hydro.parts[pid]), -#if defined(MINIMAL_SPH) || defined(PLANETARY_SPH) || \ - defined(GIZMO_MFV_SPH) || defined(GIZMO_MFM_SPH) || \ - defined(SHADOWFAX_SPH) || defined(HOPKINS_PU_SPH) || \ - defined(HOPKINS_PU_SPH_MONAGHAN) || defined(GASOLINE_SPH) +#if defined(MINIMAL_SPH) || defined(PLANETARY_SPH) || \ + defined(GIZMO_MFV_SPH) || defined(GIZMO_MFM_SPH) || \ + defined(SHADOWFAX_SPH) || defined(HOPKINS_PU_SPH) || \ + defined(HOPKINS_PU_SPH_MONAGHAN) || defined(GASOLINE_SPH) || \ + defined(MAGMA_SPH) 0.f, #elif defined(ANARCHY_PU_SPH) || defined(SPHENIX_SPH) || defined(PHANTOM_SPH) main_cell->hydro.parts[pid].viscosity.div_v, diff --git a/tests/testRandomCone.c b/tests/testRandomCone.c index 08c953c00acf2adb6997be010b44ed045229b5ab..a7e4568de213770cbd6ce4a7afd4c41cb6244100 100644 --- a/tests/testRandomCone.c +++ b/tests/testRandomCone.c @@ -18,6 +18,9 @@ ******************************************************************************/ #include <config.h> +/* System includes. */ +#include <fenv.h> + /* Local headers. */ #include "swift.h" @@ -47,10 +50,9 @@ const int N_cube = 5; * the uniformity of the distribution. * @param tolerance The tolerance of each bin relative to the expected value. */ -void test_cone(int64_t id_bh, const integertime_t ti_current, - const enum random_number_type type, double opening_angle, - float unit_vector[3], const int N_test, const int N_bins, - const double tolerance) { +float test_cone(int64_t id_bh, const integertime_t ti_current, + const enum random_number_type type, double opening_angle, + float unit_vector[3]) { /* Compute cosine that corresponds to the maximum opening angle */ const double cos_theta_max = cos(opening_angle); @@ -58,56 +60,41 @@ void test_cone(int64_t id_bh, const integertime_t ti_current, /* Initialize an array that will hold a random vector every step */ float rand_vector[3]; - /* Initialize an array that will hold the binned number of drawn cosines, - i.e. this is the probability density function that we wish to test. */ - double binned_cosines[N_bins]; - for (int j = 0; j < N_bins; ++j) { - binned_cosines[j] = 0.; + /* Generate a random unit vector within a cone around unit_vector */ + random_direction_in_cone(id_bh, ti_current, type, opening_angle, unit_vector, + rand_vector); + + /* Check that this vector is actually within the cone we want */ + const double cos_rand_unit = rand_vector[0] * unit_vector[0] + + rand_vector[1] * unit_vector[1] + + rand_vector[2] * unit_vector[2]; + if (cos_rand_unit < 0.99999 * cos_theta_max) { + printf("Cos_opening_angle is: %f, Random cos is: %f\n", cos_theta_max, + cos_rand_unit); + error("Generated random unit vector is outside cone."); } - for (int k = 0; k < N_test; ++k) { + return cos_rand_unit; +} - /* Generate a random unit vector within a cone around unit_vector */ - random_direction_in_cone(id_bh, ti_current, type, opening_angle, - unit_vector, rand_vector); - - /* Check that this vector is actually within the cone we want */ - const double cos_rand_unit = rand_vector[0] * unit_vector[0] + - rand_vector[1] * unit_vector[1] + - rand_vector[2] * unit_vector[2]; - if (cos_rand_unit < 0.99999 * cos_theta_max) { - printf("Cos_opening_angle is: %f, Random cos is: %f\n", cos_theta_max, - cos_rand_unit); - error("Generated random unit vector is outside cone."); - } +int main(int argc, char *argv[]) { - /* Add the unit vector to the probability density function array. The solid - * angle subtended by some angle theta grows as (1-cos(theta)). Furthermore, - * we are limited to the spherical cap defined by the angles [0, theta_max]. - * Therefore the variable which we expect to be uniformly distributed is (1 - * - cos(theta)) / (1 - cos(theta_max)). */ - double uniform_variable = (1. - cos_rand_unit) / (1 - cos_theta_max); - for (int j = 0; j < N_bins; ++j) { - if ((uniform_variable > (double)j / (double)N_bins) && - (uniform_variable < (double)(j + 1) / (double)N_bins)) { - binned_cosines[j] = binned_cosines[j] + 1. / (double)N_test; - } - } - } + /* Initialize CPU frequency, this also starts time. */ + unsigned long long cpufreq = 0; + clocks_set_cpufreq(cpufreq); - /* Check whether the binned quantity is really uniformly distributed. If it - * is, the density (value) of each bin should be 1/N_bin. */ - for (int j = 0; j < N_bins; ++j) { - if ((binned_cosines[j] < (1. - tolerance) / (double)N_bins) || - (binned_cosines[j] > (1. + tolerance) / (double)N_bins)) { - error( - "Generated distribution of random unit vectors within a cone exceeds " - "the limit imposed by the tolerance."); - } - } -} +/* Choke on FPEs */ +#ifdef HAVE_FE_ENABLE_EXCEPT + feenableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW); +#endif -int main(int argc, char *argv[]) { + /* Get some randomness going */ + const int seed = time(NULL); + message("Seed = %d", seed); + srand(seed); + + /* Log the swift random seed */ + message("SWIFT random seed = %d", SWIFT_RANDOM_SEED_XOR); /* Test the random-vector-in-cone function, for different values of opening * angle from 0 to pi/2 (in radians). For each of these opening angles we draw @@ -132,15 +119,16 @@ int main(int argc, char *argv[]) { const double cos_unit = random_unit_interval(id_bh, ti_current, random_number_BH_kick); const double sin_unit = sqrtf(max(0., (1. - cos_unit) * (1. + cos_unit))); - const double phi_unit = random_unit_interval(id_bh * id_bh, ti_current, - random_number_BH_kick); + const double phi_unit = + (2. * M_PI) * random_unit_interval(id_bh * id_bh, ti_current, + random_number_BH_kick); unit_vector[0] = sin_unit * cos(phi_unit); unit_vector[1] = sin_unit * sin(phi_unit); unit_vector[2] = cos_unit; /* Do the test. */ test_cone(id_bh, ti_current, random_number_BH_kick, opening_angle, - unit_vector, 100000, 10, 0.1); + unit_vector); } } @@ -148,26 +136,72 @@ int main(int argc, char *argv[]) { * bins, but for just one opening angle and one randomly generated cone */ const double opening_angle_0 = 0.2; - /* Generate an id for the bh and a time. We do this for every opening - * angle and every cone. */ - const long long id_bh_0 = rand() * (1LL << 31) + rand(); - const integertime_t ti_current_0 = rand() * (1LL << 31) + rand(); + /* Compute cosine that corresponds to the maximum opening angle */ + const double cos_theta_max = cos(opening_angle_0); /* Generate a random unit vector that defines a cone, along with the * opening angle. */ + const long long id_bh_0 = rand() * (1LL << 31) + rand(); + const integertime_t ti_current_0 = rand() * (1LL << 31) + rand(); + float unit_vector_0[3]; const double cos_unit = random_unit_interval(id_bh_0, ti_current_0, random_number_BH_kick); const double sin_unit = sqrtf(max(0., (1. - cos_unit) * (1. + cos_unit))); - const double phi_unit = random_unit_interval(id_bh_0 * id_bh_0, ti_current_0, - random_number_BH_kick); + const double phi_unit = + (2. * M_PI) * random_unit_interval(id_bh_0 * id_bh_0, ti_current_0, + random_number_BH_kick); unit_vector_0[0] = sin_unit * cos(phi_unit); unit_vector_0[1] = sin_unit * sin(phi_unit); unit_vector_0[2] = cos_unit; - /* Do the test. */ - test_cone(id_bh_0, ti_current_0, random_number_BH_kick, opening_angle_0, - unit_vector_0, 100000000, 100, 0.01); + /* Some parameters to test the uniformity of drawn vectors */ + int N_test = 10000000; + int N_bins = 100; + float tolerance = 0.05; + + /* Initialize an array that will hold the binned number of drawn cosines, + i.e. this is the probability density function that we wish to test. */ + double binned_cosines[N_bins]; + for (int j = 0; j < N_bins; ++j) { + binned_cosines[j] = 0.; + } + + /* Draw N_test vectors and bin them to test uniformity */ + for (int k = 0; k < N_test; ++k) { + + const long long id_bh = rand() * (1LL << 31) + rand(); + const integertime_t ti_current = rand() * (1LL << 31) + rand(); + + /* Do the test, with a newly generated BH id and time */ + const float cos_rand_unit = + test_cone(id_bh, ti_current, random_number_BH_kick, opening_angle_0, + unit_vector_0); + + /* Add the unit vector to the probability density function array. The solid + * angle subtended by some angle theta grows as (1-cos(theta)). Furthermore, + * we are limited to the spherical cap defined by the angles [0, theta_max]. + * Therefore the variable which we expect to be uniformly distributed is (1 + * - cos(theta)) / (1 - cos(theta_max)). */ + double uniform_variable = (1. - cos_rand_unit) / (1 - cos_theta_max); + for (int j = 0; j < N_bins; ++j) { + if ((uniform_variable > (double)j / (double)N_bins) && + (uniform_variable < (double)(j + 1) / (double)N_bins)) { + binned_cosines[j] = binned_cosines[j] + 1. / (double)N_test; + } + } + } + + /* Check whether the binned quantity is really uniformly distributed. If it + * is, the density (value) of each bin should be 1/N_bin. */ + for (int j = 0; j < N_bins; ++j) { + if ((binned_cosines[j] < (1. - tolerance) / (double)N_bins) || + (binned_cosines[j] > (1. + tolerance) / (double)N_bins)) { + error( + "Generated distribution of random unit vectors within a cone exceeds " + "the limit imposed by the tolerance."); + } + } /* We now repeat the same process, but we do not generate random unit vectors * to define the cones. Instead, we sample unit vectors along the grid @@ -204,7 +238,7 @@ int main(int argc, char *argv[]) { /* Do the test. */ test_cone(id_bh, ti_current, random_number_BH_kick, opening_angle, - unit_vector, 100000, 10, 0.1); + unit_vector); } } }