MAGMA2 SPH

.gitignore

@@ -72,6 +72,7 @@ examples/*/*/restart
 examples/Cooling/CoolingRates/cooling_rates
 examples/Cooling/CoolingRates/cooling_element_*.dat
 examples/Cooling/CoolingRates/cooling_output.dat
+examples/Cooling/CoolingSedovBlast_3D/run/*
 examples/SubgridTests/StellarEvolution/StellarEvolutionSolution*
 examples/SubgridTests/CosmologicalStellarEvolution/StellarEvolutionSolution*
 examples/SmallCosmoVolume/SmallCosmoVolume_DM/power_spectra

configure.ac

@@ -331,6 +331,18 @@ if test "$enable_debugging_checks" = "yes"; then
    AC_DEFINE([SWIFT_DEBUG_CHECKS],1,[Enable expensive debugging])
 fi
 
+# Check if expensive debugging is on.
+AC_ARG_ENABLE([magma2-debugging-checks],
+   [AS_HELP_STRING([--enable-magma2-debugging-checks],
+     [Activate expensive MAGMA2 debugging checks @<:@yes/no@:>@]
+   )],
+   [enable_magma2_debugging_checks="$enableval"],
+   [enable_magma2_debugging_checks="no"]
+)
+if test "$enable_magma2_debugging_checks" = "yes"; then
+   AC_DEFINE([MAGMA2_DEBUG_CHECKS],1,[Enable expensive MAGMA2 debugging])
+fi
+
 # Check if cell graph is on.
 AC_ARG_ENABLE([cell-graph],
    [AS_HELP_STRING([--enable-cell-graph],
@@ -2223,7 +2235,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, shadowswift, planetary, remix, sphenix, gasoline, anarchy-pu default: sphenix@:>@]
+      [Hydro dynamics to use @<:@gadget2, minimal, magma2, pressure-entropy, pressure-energy, pressure-energy-monaghan, phantom, gizmo-mfv, gizmo-mfm, shadowswift, planetary, remix, sphenix, gasoline, anarchy-pu default: sphenix@:>@]
    )],
    [with_hydro="$withval"],
    [with_hydro="sphenix"]
@@ -2279,6 +2291,9 @@ case "$with_hydro" in
    gasoline)
       AC_DEFINE([GASOLINE_SPH], [1], [Gasoline SPH])
    ;;
+   magma2)
+      AC_DEFINE([MAGMA2_SPH], [1], [MAGMA2 SPH])
+   ;;
    anarchy-du)
       AC_DEFINE([SPHENIX_SPH], [1], [SPHENIX SPH])
    ;;

examples/Cooling/CoolingSedovBlast_3D/makeIC.py

@@ -38,8 +38,8 @@
 print(rho0, "cm s^-3")
 
 uL = 1.0e3 * pc
-uM = 1.99e30
-uv = 1.0e10
+uM = 1.99e43
+uv = 1.0e5
 ut = uL / uv
 uU = uv ** 2
 print("ut:", ut / 3.154e7, "yr")

examples/Cooling/CoolingSedovBlast_3D/sedov.yml

@@ -1,6 +1,6 @@
 # Define the system of units to use internally. 
 InternalUnitSystem:
-  UnitMass_in_cgs:     1.99e33   # g
+  UnitMass_in_cgs:     1.99e43   # g
   UnitLength_in_cgs:   3.086e21   # cm
   UnitVelocity_in_cgs: 1.e5   # Centimeters per second
   UnitCurrent_in_cgs:  1   # Amperes

examples/Cosmology/ZeldovichPancake_3D/makeIC.py

@@ -43,7 +43,7 @@
 
 # Some useful variables in h-full units
 H_0 = 1.0 / Mpc_in_m * 10 ** 5  # h s^-1
-rho_0 = 3.0 * H_0 ** 2 / (8 * math.pi * G_in_SI)  # h^2 kg m^-3
+rho_0 = 3.0 * H_0 ** 2 / (8 * pi * G_in_SI)  # h^2 kg m^-3
 lambda_i = 64.0 / H_0 * 10 ** 5  # h^-1 m (= 64 h^-1 Mpc)
 x_min = -0.5 * lambda_i
 x_max = 0.5 * lambda_i

examples/HydroTests/KelvinHelmholtzGrowthRate_3D/run.sh

@@ -14,7 +14,7 @@ then
 fi
 
 # Run SWIFT
-../../../swift --hydro --threads=1 kelvinHelmholtzGrowthRate.yml 2>&1 | tee output.log
+../../../swift --hydro --threads=16 kelvinHelmholtzGrowthRate.yml 2>&1 | tee output.log
 
 # Plot the solution
 python3 plotSolution.py 100

examples/HydroTests/Rayleigh-Taylor_2D/makeMovie.py

@@ -71,7 +71,7 @@ def project(data, m_res, property, ylim):
     x = x[mask]
     y = y[mask] - np.float64(ylim[0])
 
-    h = data.gas.smoothing_length[mask]
+    h = data.gas.smoothing_lengths[mask]
 
     if property == "density":
         property = "masses"

src/debug.c

@@ -80,6 +80,8 @@
 #include "./hydro/SPHENIX/hydro_debug.h"
 #elif defined(GASOLINE_SPH)
 #include "./hydro/Gasoline/hydro_debug.h"
+#elif defined(MAGMA2_SPH)
+#include "./hydro/MAGMA2/hydro_debug.h"
 #elif defined(ANARCHY_PU_SPH)
 #include "./hydro/AnarchyPU/hydro_debug.h"
 #else

src/forcing/roberts_flow_acceleration/forcing.h

@@ -87,8 +87,14 @@ __attribute__((always_inline)) INLINE static void forcing_terms_apply(
 
   /* Effective viscosity from artificial viscosity, as in eq. 100 from
    * arXiv:1012.1885 */
+#ifndef MAGMA2_SPH
   const float nu = terms->nu * p->viscosity.alpha * c_s * p->h /
                    (2.f * (hydro_dimension + 2.f));
+#else
+  /* TODO: MAGMA2 does not track alpha per particle */
+  const float nu = terms->nu * const_viscosity_alpha * c_s * p->h /
+                   (2.f * (hydro_dimension + 1.f));
+#endif
 
   const float Vz_factor = terms->Vz_factor;
   const double k0 = (2. * M_PI / L) * terms->kv;

src/hydro.h

@@ -82,6 +82,10 @@
 #include "./hydro/Gasoline/hydro.h"
 #include "./hydro/Gasoline/hydro_iact.h"
 #define SPH_IMPLEMENTATION "Gasoline-2 (Wadsley+ 2017)"
+#elif defined(MAGMA2_SPH)
+#include "./hydro/MAGMA2/hydro.h"
+#include "./hydro/MAGMA2/hydro_iact.h"
+#define SPH_IMPLEMENTATION "MAGMA2 (Rosswog 2020)"
 #elif defined(ANARCHY_PU_SPH)
 #include "./hydro/AnarchyPU/hydro.h"
 #include "./hydro/AnarchyPU/hydro_iact.h"