Add ANARCHY-DU and refactor `hydro_properties.c`

configure.ac

@@ -1371,7 +1371,7 @@ esac
 # 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, default, gizmo-mfv, gizmo-mfm, shadowfax, planetary, anarchy-pu default: gadget2@:>@]
+      [Hydro dynamics to use @<:@gadget2, minimal, pressure-entropy, pressure-energy, pressure-energy-monaghan, default, gizmo-mfv, gizmo-mfm, shadowfax, planetary, anarchy-du, anarchy-pu default: gadget2@:>@]
    )],
    [with_hydro="$withval"],
    [with_hydro="gadget2"]
@@ -1408,6 +1408,9 @@ case "$with_hydro" in
    planetary)
       AC_DEFINE([PLANETARY_SPH], [1], [Planetary SPH])
    ;;
+   anarchy-du)
+      AC_DEFINE([ANARCHY_DU_SPH], [1], [ANARCHY (DU) SPH])
+   ;;
    anarchy-pu)
       AC_DEFINE([ANARCHY_PU_SPH], [1], [ANARCHY (PU) SPH])
    ;;

examples/HydroTests/SodShock_1D/plotSolution.py

 ###############################################################################
  # This file is part of SWIFT.
  # Copyright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ #               2019 Josh Borrow (josh.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
@@ -270,7 +271,7 @@ ylim(0.01, 1.1)
 
 # Internal energy profile -------------------------
 subplot(234)
-scatter(x, u, marker='.', c=alpha_diff, s=4.0)
+plot(x, u, '.', color='r', ms=4.0)
 plot(x_s, u_s, '--', color='k', alpha=0.8, lw=1.2)
 xlabel("${\\rm{Position}}~x$", labelpad=0)
 ylabel("${\\rm{Internal~Energy}}~u$", labelpad=0)

examples/HydroTests/SodShock_1D/sodShock.yml

@@ -27,10 +27,6 @@ Statistics:
 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.
-  viscosity_alpha_min:   0.01
-  viscosity_alpha:       0.01
-  viscosity_alpha_max:   2.0
-  viscosity_length:      0.02
 
 # Parameters related to the initial conditions
 InitialConditions:

src/Makefile.am

@@ -94,15 +94,31 @@ nobase_noinst_HEADERS = align.h approx_math.h atomic.h barrier.h cycle.h error.h
 		 gravity/Potential/gravity_debug.h gravity/Potential/gravity_part.h  \
 		 equation_of_state.h \
 		 equation_of_state/ideal_gas/equation_of_state.h equation_of_state/isothermal/equation_of_state.h \
-	 	 hydro.h hydro_io.h \
+	 	 hydro.h hydro_io.h hydro_parameters.h \
 		 hydro/Minimal/hydro.h hydro/Minimal/hydro_iact.h hydro/Minimal/hydro_io.h \
                  hydro/Minimal/hydro_debug.h hydro/Minimal/hydro_part.h \
+				 hydro/Minimal/hydro_parameters.h \
 		 hydro/Default/hydro.h hydro/Default/hydro_iact.h hydro/Default/hydro_io.h \
                  hydro/Default/hydro_debug.h hydro/Default/hydro_part.h \
+				 hydro/Default/hydro_parameters.h \
 		 hydro/Gadget2/hydro.h hydro/Gadget2/hydro_iact.h hydro/Gadget2/hydro_io.h \
                  hydro/Gadget2/hydro_debug.h hydro/Gadget2/hydro_part.h \
+				 hydro/Gadget2/hydro_parameters.h \
 		 hydro/PressureEntropy/hydro.h hydro/PressureEntropy/hydro_iact.h hydro/PressureEntropy/hydro_io.h \
                  hydro/PressureEntropy/hydro_debug.h hydro/PressureEntropy/hydro_part.h \
+				 hydro/PressureEntropy/hydro_parameters.h \
+		 hydro/PressureEnergy/hydro.h hydro/PressureEnergy/hydro_iact.h hydro/PressureEnergy/hydro_io.h \
+                 hydro/PressureEnergy/hydro_debug.h hydro/PressureEnergy/hydro_part.h \
+				 hydro/PressureEnergy/hydro_parameters.h \
+		 hydro/PressureEnergyMorrisMonaghanAV/hydro.h hydro/PressureEnergyMorrisMonaghanAV/hydro_iact.h hydro/PressureEnergyMorrisMonaghanAV/hydro_io.h \
+                 hydro/PressureEnergyMorrisMonaghanAV/hydro_debug.h hydro/PressureEnergyMorrisMonaghanAV/hydro_part.h \
+				 hydro/PressureEnergyMorrisMonaghanAV/hydro_parameters.h \
+		 hydro/AnarchyPU/hydro.h hydro/PressureEnergy/hydro_iact.h hydro/PressureEnergy/hydro_io.h \
+                 hydro/AnarchyPU/hydro_debug.h hydro/PressureEnergy/hydro_part.h \
+				 hydro/AnarchyPU/hydro_parameters.h \
+		 hydro/AnarchyDU/hydro.h hydro/PressureEnergy/hydro_iact.h hydro/PressureEnergy/hydro_io.h \
+                 hydro/AnarchyDU/hydro_debug.h hydro/PressureEnergy/hydro_part.h \
+				 hydro/AnarchyDU/hydro_parameters.h \
 		 hydro/GizmoMFV/hydro.h hydro/GizmoMFV/hydro_iact.h \
                  hydro/GizmoMFV/hydro_io.h hydro/GizmoMFV/hydro_debug.h \
                  hydro/GizmoMFV/hydro_part.h \
@@ -114,6 +130,7 @@ nobase_noinst_HEADERS = align.h approx_math.h atomic.h barrier.h cycle.h error.h
                  hydro/GizmoMFV/hydro_slope_limiters.h \
                  hydro/GizmoMFV/hydro_unphysical.h \
                  hydro/GizmoMFV/hydro_velocities.h \
+				 hydro/GizmoMFV/hydro_parameters.h \
 		 hydro/GizmoMFM/hydro.h hydro/GizmoMFM/hydro_iact.h \
                  hydro/GizmoMFM/hydro_io.h hydro/GizmoMFM/hydro_debug.h \
                  hydro/GizmoMFM/hydro_part.h \
@@ -124,6 +141,7 @@ nobase_noinst_HEADERS = align.h approx_math.h atomic.h barrier.h cycle.h error.h
                  hydro/GizmoMFM/hydro_slope_limiters_face.h \
                  hydro/GizmoMFM/hydro_slope_limiters.h \
                  hydro/GizmoMFM/hydro_unphysical.h \
+				 hydro/GizmoMFM/hydro_parameters.h \
                  hydro/Shadowswift/hydro_debug.h \
                  hydro/Shadowswift/hydro_gradients.h hydro/Shadowswift/hydro.h \
                  hydro/Shadowswift/hydro_iact.h \
@@ -140,6 +158,7 @@ nobase_noinst_HEADERS = align.h approx_math.h atomic.h barrier.h cycle.h error.h
                  hydro/Shadowswift/voronoi3d_cell.h \
                  hydro/Shadowswift/voronoi_algorithm.h \
                  hydro/Shadowswift/voronoi_cell.h \
+                 hydro/Shadowswift/hydro_parameters.h \
 	         riemann/riemann_hllc.h riemann/riemann_trrs.h \
 		 riemann/riemann_exact.h riemann/riemann_vacuum.h \
                  riemann/riemann_checks.h \

src/const.h

@@ -20,16 +20,6 @@
 #ifndef SWIFT_CONST_H
 #define SWIFT_CONST_H
 
-/* SPH Viscosity constants. */
-/* Cosmology default beta=3.0. Planetary default beta=4.0
- * Alpha can be set in the parameter file.
- * Beta is defined as in e.g. Price (2010) Eqn (103) */
-#define const_viscosity_beta 3.0f
-
-/* SPH Thermal conductivity constants. */
-#define const_conductivity_alpha \
-  1.f /* Value taken from (Price,2008), not used in legacy gadget mode */
-
 /* Time integration constants. */
 #define const_max_u_change 0.1f
 

src/debug.c

@@ -60,6 +60,8 @@
 #include "./hydro/Shadowswift/hydro_debug.h"
 #elif defined(PLANETARY_SPH)
 #include "./hydro/Planetary/hydro_debug.h"
+#elif defined(ANARCHY_DU_SPH)
+#include "./hydro/AnarchyDU/hydro_debug.h"
 #elif defined(ANARCHY_PU_SPH)
 #include "./hydro/AnarchyPU/hydro_debug.h"
 #else

src/hydro.h

@@ -72,6 +72,10 @@
 #include "./hydro/Planetary/hydro.h"
 #include "./hydro/Planetary/hydro_iact.h"
 #define SPH_IMPLEMENTATION "Minimal version of SPH with multiple materials"
+#elif defined(ANARCHY_DU_SPH)
+#include "./hydro/AnarchyDU/hydro.h"
+#include "./hydro/AnarchyDU/hydro_iact.h"
+#define SPH_IMPLEMENTATION "ANARCHY (Density-Energy) SPH (Borrow+ in prep)"
 #elif defined(ANARCHY_PU_SPH)
 #include "./hydro/AnarchyPU/hydro.h"
 #include "./hydro/AnarchyPU/hydro_iact.h"

src/hydro/AnarchyDU/hydro_debug.h

+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Coypright (c) 2019 Josh Borrow (joshua.borrow@durham.ac.uk) &
+ *                    Matthieu Schaller (matthieu.schaller@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_ANARCHY_DU_HYDRO_DEBUG_H
+#define SWIFT_ANARCHY_DU_HYDRO_DEBUG_H
+
+/**
+ * @file AnarchyDU/hydro_debug.h
+ * @brief Density-Energy conservative implementation of SPH,
+ *        with added ANARCHY physics (Cullen & Denhen 2011 AV,
+ *        Price 2008 thermal diffusion (Debugging routines)
+ */
+
+__attribute__((always_inline)) INLINE static void hydro_debug_particle(
+    const struct part* p, const struct xpart* xp) {
+  printf(
+      "x=[%.3e,%.3e,%.3e], "
+      "v=[%.3e,%.3e,%.3e],v_full=[%.3e,%.3e,%.3e] \n a=[%.3e,%.3e,%.3e], "
+      "u=%.3e, du/dt=%.3e v_sig=%.3e, P=%.3e\n"
+      "h=%.3e, dh/dt=%.3e wcount=%d, m=%.3e, dh_drho=%.3e, rho=%.3e, \n"
+      "alpha=%.3e \n"
+      "time_bin=%d\n",
+      p->x[0], p->x[1], p->x[2], p->v[0], p->v[1], p->v[2], xp->v_full[0],
+      xp->v_full[1], xp->v_full[2], p->a_hydro[0], p->a_hydro[1], p->a_hydro[2],
+      p->u, p->u_dt, p->viscosity.v_sig, hydro_get_comoving_pressure(p), p->h,
+      p->force.h_dt, (int)p->density.wcount, p->mass, p->density.rho_dh, p->rho,
+      p->viscosity.alpha, p->time_bin);
+}
+
+#endif /* SWIFT_ANARCHY_DU_HYDRO_DEBUG_H */

src/hydro/AnarchyDU/hydro_io.h

+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Coypright (c) 2019 Josh Borrow (joshua.borrow@durham.ac.uk) &
+ *                    Matthieu Schaller (matthieu.schaller@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 DURPOSE.  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_ANARCHY_DU_HYDRO_IO_H
+#define SWIFT_ANARCHY_DU_HYDRO_IO_H
+
+/**
+ * @file AnarchyDU/hydro_io.h
+ * @brief Density-Energy conservative implementation of SPH,
+ *        with added ANARCHY physics (Cullen & Denhen 2011 AV,
+ *        Price 2008 thermal diffusion (i/o routines)
+ */
+
+#include "adiabatic_index.h"
+#include "hydro.h"
+#include "io_properties.h"
+#include "kernel_hydro.h"
+
+#include "./hydro_parameters.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) {
+
+  if (e->s->periodic) {
+    ret[0] = box_wrap(p->x[0], 0.0, e->s->dim[0]);
+    ret[1] = box_wrap(p->x[1], 0.0, e->s->dim[1]);
+    ret[2] = box_wrap(p->x[2], 0.0, e->s->dim[2]);
+  } else {
+    ret[0] = p->x[0];
+    ret[1] = p->x[1];
+    ret[2] = p->x[2];
+  }
+}
+
+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 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_get_drifted_velocities(p, xp, dt_kick_hydro, dt_kick_grav, ret);
+
+  /* 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;
+}
+
+INLINE static void convert_viscosity(const struct engine* e,
+                                     const struct part* p,
+                                     const struct xpart* xp, float* ret) {
+  ret[0] = p->viscosity.alpha;
+}
+
+INLINE static void convert_diffusion(const struct engine* e,
+                                     const struct part* p,
+                                     const struct xpart* xp, float* ret) {
+  ret[0] = p->diffusion.alpha;
+}
+
+/**
+ * @brief Specifies which particle fields to write to a dataset
+ *
+ * @param parts The 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 = 12;
+
+  /* List what we want to write */
+  list[0] = io_make_output_field_convert_part("Coordinates", DOUBLE, 3,
+                                              UNIT_CONV_LENGTH, parts, xparts,
+                                              convert_part_pos);
+  list[1] = io_make_output_field_convert_part(
+      "Velocities", FLOAT, 3, UNIT_CONV_SPEED, parts, xparts, convert_part_vel);
+  list[2] =
+      io_make_output_field("Masses", FLOAT, 1, UNIT_CONV_MASS, parts, mass);
+  list[3] = io_make_output_field("SmoothingLength", FLOAT, 1, UNIT_CONV_LENGTH,
+                                 parts, h);
+  list[4] = io_make_output_field("InternalEnergy", FLOAT, 1,
+                                 UNIT_CONV_ENERGY_PER_UNIT_MASS, parts, u);
+  list[5] = io_make_output_field("ParticleIDs", ULONGLONG, 1,
+                                 UNIT_CONV_NO_UNITS, parts, id);
+  list[6] =
+      io_make_output_field("Density", FLOAT, 1, UNIT_CONV_DENSITY, parts, rho);
+  list[7] = io_make_output_field_convert_part(
+      "Pressure", FLOAT, 1, UNIT_CONV_PRESSURE, parts, xparts, convert_P);
+  list[8] = io_make_output_field_convert_part("Entropy", FLOAT, 1,
+                                              UNIT_CONV_ENTROPY_PER_UNIT_MASS,
+                                              parts, xparts, convert_S);
+  list[9] = io_make_output_field_convert_part("Potential", FLOAT, 1,
+                                              UNIT_CONV_POTENTIAL, parts,
+                                              xparts, convert_part_potential);
+  list[10] = io_make_output_field_convert_part("Viscosity", FLOAT, 1,
+                                               UNIT_CONV_NO_UNITS, parts,
+                                               xparts, convert_viscosity);
+  list[11] = io_make_output_field_convert_part("Diffusion", FLOAT, 1,
+                                               UNIT_CONV_NO_UNITS, parts,
+                                               xparts, convert_diffusion);
+}
+
+/**
+ * @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",
+                       "Simple treatment as in Price (2008)");
+  io_write_attribute_s(h_grpsph, "Viscosity Model",
+                       "Simplified version of Cullen & Denhen (2011)");
+
+  /* Time integration properties */
+  io_write_attribute_f(h_grpsph, "Maximal Delta u change over dt",
+                       const_max_u_change);
+}
+
+/**
+ * @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_ANARCHY_DU_HYDRO_IO_H */

src/hydro/AnarchyDU/hydro_parameters.h

+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Coypright (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_ANARCHY_DU_HYDRO_PARAMETERS_H
+#define SWIFT_ANARCHY_DU_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 AnarchyDU/hydro_parameters.h
+ * @brief Density-Energy conservative implementation of SPH,
+ *        with added ANARCHY physics (Cullen & Denhen 2011 AV,
+ *        Price 2008 thermal diffusion (default compile-time
+ *        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_beta 3.0f
+
+/*! 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 2.0f
+
+/* 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.1f
+
+/*! Minimal value for the viscosity alpha in variable schemes. */
+#define hydro_props_default_viscosity_alpha_min 0.0f
+
+/*! Maximal value for the viscosity alpha in variable schemes. */
+#define hydro_props_default_viscosity_alpha_max 2.0f
+
+/*! Decay length for the viscosity scheme. This is scheme dependent. In
+ * non-variable schemes this must be defined but is not used. This also
+ * sets the decay length for the diffusion. */
+#define hydro_props_default_viscosity_length 0.25f
+
+/* Diffusion parameters -- Defaults; can be changed at run-time */
+
+/*! The "initial" diffusion, or the fixed value for non-variable
+ * schemes. This usually takes the value 0.0. */
+#define hydro_props_default_diffusion_alpha 0.0f
+
+/*! Beta coefficient for the diffusion. This controls how fast the
+ * diffusion coefficient peaks, and how high it can get. Chosen to be
+ * very small in schemes where little diffusion is needed, 0.2-1.0 in
+ * schemes (e.g. density-energy) where diffusion is needed to solve
+ * the contact discontinuity problem. */
+#define hydro_props_default_diffusion_beta 0.25f
+
+/*! Maximal value for the diffusion alpha in variable schemes. */
+#define hydro_props_default_diffusion_alpha_max 1.0f
+
+/*! Minimal value for the diffusion alpha in variable schemes. */
+#define hydro_props_default_diffusion_alpha_min 0.0f
+
+/* 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;
+
+  /*! Artificial viscosity (max) for the variable case (e.g. M&M) */
+  float alpha_max;
+
+  /*! Artificial viscosity (min) for the variable case (e.g. M&M) */
+  float alpha_min;
+
+  /*! The decay length of the artificial viscosity (used in M&M, etc.) */
+  float length;
+};
+
+/*! Thermal diffusion parameters */
+struct diffusion_global_data {
+
+  /*! Initialisation value, or the case for constant thermal diffusion coeffs
+   */
+  float alpha;
+
+  /*! Tuning parameter for speed of ramp up/down */
+  float beta;
+
+  /*! Maximal value for alpha_diff */
+  float alpha_max;
+
+  /*! Minimal value for alpha_diff */
+  float alpha_min;
+};
+
+/* 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 unit_system: pointer to the 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);
+
+  viscosity->alpha_max =
+      parser_get_opt_param_float(params, "SPH:viscosity_alpha_max",
+                                 hydro_props_default_viscosity_alpha_max);
+
+  viscosity->alpha_min =
+      parser_get_opt_param_float(params, "SPH:viscosity_alpha_min",
+                                 hydro_props_default_viscosity_alpha_min);
+
+  viscosity->length = parser_get_opt_param_float(
+      params, "SPH:viscosity_length", hydro_props_default_viscosity_length);
+}
+
+/**
+ * @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;
+  viscosity->alpha_max = hydro_props_default_viscosity_alpha_max;
+  viscosity->alpha_min = hydro_props_default_viscosity_alpha_min;
+  viscosity->length = hydro_props_default_viscosity_length;
+}
+
+/**
+ * @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) {