From 156600bff2381f723c029a1707a043dfed56c43d Mon Sep 17 00:00:00 2001
From: Matthieu Schaller <schaller@strw.leidenuniv.nl>
Date: Thu, 5 Oct 2023 14:34:48 +0000
Subject: [PATCH] Add forcing terms to the hydro equations
---
configure.ac | 25 +++
examples/parameter_example.yml | 13 ++
src/Makefile.am | 3 +
src/engine.c | 10 ++
src/engine.h | 5 +
src/forcing.c | 51 ++++++
src/forcing.h | 45 +++++
src/forcing/none/forcing.h | 109 ++++++++++++
src/forcing/roberts_flow/forcing.h | 136 +++++++++++++++
.../roberts_flow_acceleration/forcing.h | 157 ++++++++++++++++++
src/runner_others.c | 7 +
src/swift.h | 1 +
src/timestep.h | 9 +-
swift.c | 12 +-
swift_fof.c | 1 +
15 files changed, 579 insertions(+), 5 deletions(-)
create mode 100644 src/forcing.c
create mode 100644 src/forcing.h
create mode 100644 src/forcing/none/forcing.h
create mode 100644 src/forcing/roberts_flow/forcing.h
create mode 100644 src/forcing/roberts_flow_acceleration/forcing.h
diff --git a/configure.ac b/configure.ac
index 136d607645..0eac6e9da2 100644
--- a/configure.ac
+++ b/configure.ac
@@ -2621,6 +2621,30 @@ case "$with_sink" in
;;
esac
+# Forcing terms
+AC_ARG_WITH([forcing],
+ [AS_HELP_STRING([--with-forcing=<term>],
+ [Hydrodynamics forcing terms @<:@none, roberts-flow, roberts-flow-acceleration default: none@:>@]
+ )],
+ [with_forcing="$withval"],
+ [with_forcing="none"]
+)
+case "$with_forcing" in
+ none)
+ AC_DEFINE([FORCING_NONE], [1], [No external forcing terms])
+ ;;
+ roberts-flow)
+ AC_DEFINE([FORCING_ROBERTS_FLOW], [1], [Roberts' flow external forcing terms])
+ ;;
+ roberts-flow-acceleration)
+ AC_DEFINE([FORCING_ROBERTS_FLOW_ACCELERATION], [1], [Roberts' flow external forcing terms entering the equations as an acceleration term])
+ ;;
+ *)
+ AC_MSG_ERROR([Unknown external forcing term: $with_forcing])
+ ;;
+esac
+
+
# External potential
AC_ARG_WITH([ext-potential],
[AS_HELP_STRING([--with-ext-potential=<pot>],
@@ -3041,6 +3065,7 @@ AC_MSG_RESULT([
No gravity below ID : $no_gravity_below_id
Make gravity glass : $gravity_glass_making
External potential : $with_potential
+ Forcing terms : $with_forcing
Pressure floor : $with_pressure_floor
Entropy floor : $with_entropy_floor
diff --git a/examples/parameter_example.yml b/examples/parameter_example.yml
index 0e5cac61fc..5e1cdcd036 100644
--- a/examples/parameter_example.yml
+++ b/examples/parameter_example.yml
@@ -438,6 +438,19 @@ MWPotential2014Potential:
r_c_kpc: 1.9 # Cut-off radius of the bulge (in kpc)
potential_factors: [0.4367419745056084, 1.002641971008805, 0.022264787598364262] # Coefficients that adjust the strength of the halo (1st component), the disk (2nd component) and the bulge (3rd component)
+# Parameters related to forcing terms ----------------------------------------------
+
+# Roberts' flow (Tilgner & Brandenburg, 2008, MNRAS, 391, 1477) where the velocity is imposed
+RobertsFlowForcing:
+ u0: 1.5 # Forcing velocity (internal units)
+ Vz_factor: 1.0 # (optional) Scaling of the velocity along the z axis
+
+# Roberts' flow (Tilgner & Brandenburg, 2008, MNRAS, 391, 1477) where the forcing is an acceleration
+RobertsFlowAccelerationForcing:
+ u0: 1.5 # Forcing velocity (internal units)
+ nu: 1.0 # Viscosity (internal units) to convert velocities to accelerations.
+ Vz_factor: 1.0 # (optional) Scaling of the velocity along the z axis
+
# Parameters related to entropy floors ----------------------------------------------
EAGLEEntropyFloor:
diff --git a/src/Makefile.am b/src/Makefile.am
index b0ad938f55..ac7ea7da7f 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -74,6 +74,7 @@ include_HEADERS += lightcone/lightcone.h lightcone/lightcone_particle_io.h light
include_HEADERS += lightcone/lightcone_crossing.h lightcone/lightcone_array.h lightcone/lightcone_map.h
include_HEADERS += lightcone/lightcone_map_types.h lightcone/projected_kernel.h lightcone/lightcone_shell.h
include_HEADERS += lightcone/healpix_util.h lightcone/pixel_index.h
+include_HEADERS += forcing.h
include_HEADERS += power_spectrum.h
include_HEADERS += ghost_stats.h
@@ -186,6 +187,7 @@ AM_SOURCES += lightcone/lightcone.c lightcone/lightcone_particle_io.c lightcone/
AM_SOURCES += lightcone/healpix_util.c lightcone/lightcone_array.c lightcone/lightcone_map.c
AM_SOURCES += lightcone/lightcone_map_types.c lightcone/projected_kernel.c lightcone/lightcone_shell.c
AM_SOURCES += power_spectrum.c
+AM_SOURCES += forcing.c
AM_SOURCES += ghost_stats.c
AM_SOURCES += $(EAGLE_EXTRA_IO_SOURCES)
AM_SOURCES += $(QLA_COOLING_SOURCES) $(QLA_EAGLE_COOLING_SOURCES)
@@ -363,6 +365,7 @@ nobase_noinst_HEADERS += stars/EAGLE/stars.h stars/EAGLE/stars_iact.h stars/EAGL
nobase_noinst_HEADERS += stars/EAGLE/stars_debug.h stars/EAGLE/stars_part.h
nobase_noinst_HEADERS += stars/GEAR/stars.h stars/GEAR/stars_iact.h stars/GEAR/stars_io.h
nobase_noinst_HEADERS += stars/GEAR/stars_debug.h stars/GEAR/stars_csds.h stars/GEAR/stars_part.h
+nobase_noinst_HEADERS += forcing/none/forcing.h forcing/roberts_flow/forcing.h forcing/roberts_flow_acceleration/forcing.h
nobase_noinst_HEADERS += potential/none/potential.h potential/point_mass/potential.h
nobase_noinst_HEADERS += potential/isothermal/potential.h potential/disc_patch/potential.h
nobase_noinst_HEADERS += potential/sine_wave/potential.h potential/constant/potential.h
diff --git a/src/engine.c b/src/engine.c
index 0873b83322..b545fd66d1 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -70,6 +70,7 @@
#include "extra_io.h"
#include "feedback.h"
#include "fof.h"
+#include "forcing.h"
#include "gravity.h"
#include "gravity_cache.h"
#include "hydro.h"
@@ -3065,6 +3066,7 @@ void engine_init(
struct pressure_floor_props *pressure_floor, struct rt_props *rt,
struct pm_mesh *mesh, struct power_spectrum_data *pow_data,
const struct external_potential *potential,
+ const struct forcing_terms *forcing_terms,
struct cooling_function_data *cooling_func,
const struct star_formation *starform,
const struct chemistry_global_data *chemistry,
@@ -3205,6 +3207,7 @@ void engine_init(
e->mesh = mesh;
e->power_data = pow_data;
e->external_potential = potential;
+ e->forcing_terms = forcing_terms;
e->cooling_func = cooling_func;
e->star_formation = starform;
e->feedback_props = feedback;
@@ -3644,6 +3647,7 @@ void engine_clean(struct engine *e, const int fof, const int restart) {
free((void *)e->parameter_file);
free((void *)e->output_options);
free((void *)e->external_potential);
+ free((void *)e->forcing_terms);
free((void *)e->black_holes_properties);
free((void *)e->pressure_floor_props);
free((void *)e->rt_props);
@@ -3718,6 +3722,7 @@ void engine_struct_dump(struct engine *e, FILE *stream) {
pm_mesh_struct_dump(e->mesh, stream);
power_spectrum_struct_dump(e->power_data, stream);
potential_struct_dump(e->external_potential, stream);
+ forcing_terms_struct_dump(e->forcing_terms, stream);
cooling_struct_dump(e->cooling_func, stream);
starformation_struct_dump(e->star_formation, stream);
feedback_struct_dump(e->feedback_props, stream);
@@ -3833,6 +3838,11 @@ void engine_struct_restore(struct engine *e, FILE *stream) {
potential_struct_restore(external_potential, stream);
e->external_potential = external_potential;
+ struct forcing_terms *forcing_terms =
+ (struct forcing_terms *)malloc(sizeof(struct forcing_terms));
+ forcing_terms_struct_restore(forcing_terms, stream);
+ e->forcing_terms = forcing_terms;
+
struct cooling_function_data *cooling_func =
(struct cooling_function_data *)malloc(
sizeof(struct cooling_function_data));
diff --git a/src/engine.h b/src/engine.h
index 11394dca3e..2760fc4c23 100644
--- a/src/engine.h
+++ b/src/engine.h
@@ -54,6 +54,7 @@
struct black_holes_properties;
struct extra_io_properties;
struct external_potential;
+struct forcing_terms;
/**
* @brief The different policies the #engine can follow.
@@ -536,6 +537,9 @@ struct engine {
/* Properties of external gravitational potential */
const struct external_potential *external_potential;
+ /* Properties of the hydrodynamics forcing terms */
+ const struct forcing_terms *forcing_terms;
+
/* Properties of the cooling scheme */
struct cooling_function_data *cooling_func;
@@ -713,6 +717,7 @@ void engine_init(
struct pressure_floor_props *pressure_floor, struct rt_props *rt,
struct pm_mesh *mesh, struct power_spectrum_data *pow_data,
const struct external_potential *potential,
+ const struct forcing_terms *forcing_terms,
struct cooling_function_data *cooling_func,
const struct star_formation *starform,
const struct chemistry_global_data *chemistry,
diff --git a/src/forcing.c b/src/forcing.c
new file mode 100644
index 0000000000..2a1cfbea97
--- /dev/null
+++ b/src/forcing.c
@@ -0,0 +1,51 @@
+/*******************************************************************************
+ * 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/>.
+ *
+ ******************************************************************************/
+
+/* Config parameters. */
+#include <config.h>
+
+/* This object's header. */
+#include "forcing.h"
+#include "restart.h"
+
+/**
+ * @brief Write an forcing_terms struct to the given FILE as a stream of
+ * bytes.
+ *
+ * @param terms the struct
+ * @param stream the file stream
+ */
+void forcing_terms_struct_dump(const struct forcing_terms* terms,
+ FILE* stream) {
+ restart_write_blocks((void*)terms, sizeof(struct forcing_terms), 1, stream,
+ "forcingterms", "forcing terms");
+}
+
+/**
+ * @brief Restore a forcing_terms struct from the given FILE as a stream of
+ * bytes.
+ *
+ * @param terms the struct
+ * @param stream the file stream
+ */
+void forcing_terms_struct_restore(const struct forcing_terms* terms,
+ FILE* stream) {
+ restart_read_blocks((void*)terms, sizeof(struct forcing_terms), 1, stream,
+ NULL, "forcing terms");
+}
diff --git a/src/forcing.h b/src/forcing.h
new file mode 100644
index 0000000000..8b1d748658
--- /dev/null
+++ b/src/forcing.h
@@ -0,0 +1,45 @@
+/*******************************************************************************
+ * 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_FORCING_H
+#define SWIFT_FORCING_H
+
+/**
+ * @file src/potential.h
+ * @brief Branches between the different external gravitational potentials.
+ */
+
+/* Config parameters. */
+#include <config.h>
+
+/* Import the right external potential definition */
+#if defined(FORCING_NONE)
+#include "./forcing/none/forcing.h"
+#elif defined(FORCING_ROBERTS_FLOW)
+#include "./forcing/roberts_flow/forcing.h"
+#elif defined(FORCING_ROBERTS_FLOW_ACCELERATION)
+#include "./forcing/roberts_flow_acceleration/forcing.h"
+#else
+#error "Invalid choice of forcing terms"
+#endif
+
+void forcing_terms_struct_dump(const struct forcing_terms* terms, FILE* stream);
+void forcing_terms_struct_restore(const struct forcing_terms* terms,
+ FILE* stream);
+
+#endif /* SWIFT_FORCING_H */
diff --git a/src/forcing/none/forcing.h b/src/forcing/none/forcing.h
new file mode 100644
index 0000000000..52b811ae90
--- /dev/null
+++ b/src/forcing/none/forcing.h
@@ -0,0 +1,109 @@
+/*******************************************************************************
+ * 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_FORCING_NONE_H
+#define SWIFT_FORCING_NONE_H
+
+/* Config parameters. */
+#include <config.h>
+
+/* Standard includes. */
+#include <float.h>
+
+/* Local includes. */
+#include "error.h"
+#include "parser.h"
+#include "part.h"
+#include "physical_constants.h"
+#include "space.h"
+#include "units.h"
+
+/**
+ * @brief Forcing Term Properties
+ */
+struct forcing_terms {};
+
+/**
+ * @brief Computes the forcing terms.
+ *
+ * We do nothing in this 'none' scheme.
+ *
+ * @param time The current time.
+ * @param terms The properties of the forcing terms.
+ * @param s The #space we act on.
+ * @param phys_const The physical constants in internal units.
+ * @param p Pointer to the particle data.
+ * @param xp Pointer to the extended particle data.
+ */
+__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) {
+ /* Nothing to do here */
+}
+
+/**
+ * @brief Computes the time-step condition due to the forcing terms.
+ *
+ * Nothing to do here. --> Return FLT_MAX.
+ *
+ * @param time The current time.
+ * @param terms The properties of the forcing terms.
+ * @param phys_const The physical constants in internal units.
+ * @param p Pointer to the particle data.
+ * @param xp Pointer to the extended particle data.
+ */
+__attribute__((always_inline)) INLINE static float forcing_terms_timestep(
+ double time, const struct forcing_terms* terms,
+ const struct phys_const* phys_const, const struct part* p,
+ const struct xpart* xp) {
+
+ /* No time-step size limit */
+ return FLT_MAX;
+}
+
+/**
+ * @brief Prints the properties of the forcing terms to stdout.
+ *
+ * @param terms The #forcing_terms properties of the run.
+ */
+static INLINE void forcing_terms_print(const struct forcing_terms* terms) {
+
+ message("Forcing terms is 'No forcing terms'.");
+}
+
+/**
+ * @brief Initialises the forcing term properties
+ *
+ * Nothing to do here.
+ *
+ * @param parameter_file The parsed parameter file
+ * @param phys_const Physical constants in internal units
+ * @param us The current internal system of units
+ * @param s The #space object.
+ * @param terms The forcing term properties to initialize
+ */
+static INLINE void forcing_terms_init(struct swift_params* parameter_file,
+ const struct phys_const* phys_const,
+ const struct unit_system* us,
+ const struct space* s,
+ struct forcing_terms* terms) {
+
+ /* Nothing to do here */
+}
+
+#endif /* SWIFT_FORCING_NONE_H */
diff --git a/src/forcing/roberts_flow/forcing.h b/src/forcing/roberts_flow/forcing.h
new file mode 100644
index 0000000000..5c8f38e147
--- /dev/null
+++ b/src/forcing/roberts_flow/forcing.h
@@ -0,0 +1,136 @@
+/*******************************************************************************
+ * 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_FORCING_ROBERTS_FLOW_H
+#define SWIFT_FORCING_ROBERTS_FLOW_H
+
+/* Config parameters. */
+#include <config.h>
+
+/* Standard includes. */
+#include <float.h>
+
+/* Local includes. */
+#include "error.h"
+#include "parser.h"
+#include "part.h"
+#include "physical_constants.h"
+#include "space.h"
+#include "units.h"
+
+/**
+ * @brief Forcing Term Properties
+ */
+struct forcing_terms {
+
+ /*! Reference velocity (internal units) */
+ float u0;
+
+ /*! Velocity scaling along the z direction */
+ float Vz_factor;
+};
+
+/**
+ * @brief Computes the forcing terms.
+ *
+ * Based on Tilgner & Brandenburg, 2008, MNRAS, 391, 1477.
+ * This version differs from the paper by imposing the velocity directly rather
+ * than by giving the particles an acceleration.
+ *
+ * @param time The current time.
+ * @param terms The properties of the forcing terms.
+ * @param s The #space we act on.
+ * @param phys_const The physical constants in internal units.
+ * @param p Pointer to the particle data.
+ * @param xp Pointer to the extended particle data.
+ */
+__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 double L = s->dim[0];
+ const float u0 = terms->u0;
+ const float Vz_factor = terms->Vz_factor;
+ const double k0 = 2. * M_PI / L;
+ const double kf = M_SQRT2 * k0;
+
+ /* Eq. 8 of Tilgner & Brandenburg, 2008, MNRAS, 391, 1477 */
+ const double Psi = (u0 / k0) * cos(k0 * p->x[0]) * cos(k0 * p->x[1]);
+
+ /* Eq. 7 of Tilgner & Brandenburg, 2008, MNRAS, 391, 1477 */
+ const double v_Rob[3] = {u0 * cos(k0 * p->x[0]) * sin(k0 * p->x[1]),
+ -u0 * sin(k0 * p->x[0]) * cos(k0 * p->x[1]),
+ kf * Psi};
+
+ /* Force the velocity and possibly scale the z-direction */
+ xp->v_full[0] = v_Rob[0];
+ xp->v_full[1] = v_Rob[1];
+ xp->v_full[2] = v_Rob[2] * Vz_factor;
+}
+
+/**
+ * @brief Computes the time-step condition due to the forcing terms.
+ *
+ * Nothing to do here. --> Return FLT_MAX.
+ *
+ * @param time The current time.
+ * @param terms The properties of the forcing terms.
+ * @param phys_const The physical constants in internal units.
+ * @param p Pointer to the particle data.
+ * @param xp Pointer to the extended particle data.
+ */
+__attribute__((always_inline)) INLINE static float forcing_terms_timestep(
+ double time, const struct forcing_terms* terms,
+ const struct phys_const* phys_const, const struct part* p,
+ const struct xpart* xp) {
+
+ return FLT_MAX;
+}
+
+/**
+ * @brief Prints the properties of the forcing terms to stdout.
+ *
+ * @param terms The #forcing_terms properties of the run.
+ */
+static INLINE void forcing_terms_print(const struct forcing_terms* terms) {
+
+ message("Forcing terms is 'Roberts flow'. U0: %.5f / Vz factor: %.5f.",
+ terms->u0, terms->Vz_factor);
+}
+
+/**
+ * @brief Initialises the forcing term properties
+ *
+ * @param parameter_file The parsed parameter file
+ * @param phys_const Physical constants in internal units
+ * @param us The current internal system of units
+ * @param s The #space object.
+ * @param terms The forcing term properties to initialize
+ */
+static INLINE void forcing_terms_init(struct swift_params* parameter_file,
+ const struct phys_const* phys_const,
+ const struct unit_system* us,
+ const struct space* s,
+ struct forcing_terms* terms) {
+
+ terms->u0 = parser_get_param_double(parameter_file, "RobertsFlowForcing:u0");
+ terms->Vz_factor = parser_get_opt_param_float(
+ parameter_file, "RobertsFlowForcing:Vz_factor", 1.f);
+}
+
+#endif /* SWIFT_FORCING_ROBERTS_FLOW_H */
diff --git a/src/forcing/roberts_flow_acceleration/forcing.h b/src/forcing/roberts_flow_acceleration/forcing.h
new file mode 100644
index 0000000000..9a9d33d3c8
--- /dev/null
+++ b/src/forcing/roberts_flow_acceleration/forcing.h
@@ -0,0 +1,157 @@
+/*******************************************************************************
+ * 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_FORCING_ROBERTS_FLOW_ACCELERATION_H
+#define SWIFT_FORCING_ROBERTS_FLOW_ACCELERATION_H
+
+/* Config parameters. */
+#include <config.h>
+
+/* Standard includes. */
+#include <float.h>
+
+/* Local includes. */
+#include "error.h"
+#include "hydro.h"
+#include "parser.h"
+#include "part.h"
+#include "physical_constants.h"
+#include "space.h"
+#include "units.h"
+
+/**
+ * @brief Forcing Term Properties
+ */
+struct forcing_terms {
+
+ /*! Reference velocity (internal units) */
+ float u0;
+
+ /*! 'viscosity' to convert from velocity to acceleration (internal units) */
+ float nu;
+
+ /*! Velocity scaling along the z direction */
+ float Vz_factor;
+};
+
+/**
+ * @brief Computes the forcing terms.
+ *
+ * Based on Tilgner & Brandenburg, 2008, MNRAS, 391, 1477
+ *
+ * @param time The current time.
+ * @param terms The properties of the forcing terms.
+ * @param s The #space we act on.
+ * @param phys_const The physical constants in internal units.
+ * @param p Pointer to the particle data.
+ * @param xp Pointer to the extended particle data.
+ */
+__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 double L = s->dim[0];
+ const float u0 = terms->u0;
+ const float nu = terms->nu * p->viscosity.alpha * v_sig * p->h;
+ const float Vz_factor = terms->Vz_factor;
+ const double k0 = 2. * M_PI / L;
+ const double kf = M_SQRT2 * k0;
+
+ /* Eq. 8 of Tilgner & Brandenburg, 2008, MNRAS, 391, 1477 */
+ const double Psi = (u0 / k0) * cos(k0 * p->x[0]) * cos(k0 * p->x[1]);
+
+ /* Eq. 7 of Tilgner & Brandenburg, 2008, MNRAS, 391, 1477
+ * (with optional scaling along z) */
+ const double v_Rob[3] = {u0 * cos(k0 * p->x[0]) * sin(k0 * p->x[1]),
+ -u0 * sin(k0 * p->x[0]) * cos(k0 * p->x[1]),
+ kf * Psi * Vz_factor};
+
+ /* Eq. 6 of Tilgner & Brandenburg, 2008, MNRAS, 391, 1477 */
+ const double f[3] = {nu * kf * kf * v_Rob[0], nu * kf * kf * v_Rob[1],
+ nu * kf * kf * v_Rob[2]};
+
+ /* Update the accelerations */
+ 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];
+ }
+}
+
+/**
+ * @brief Computes the time-step condition due to the forcing terms.
+ *
+ * Nothing to do here. --> Return FLT_MAX.
+ *
+ * @param time The current time.
+ * @param terms The properties of the forcing terms.
+ * @param phys_const The physical constants in internal units.
+ * @param p Pointer to the particle data.
+ * @param xp Pointer to the extended particle data.
+ */
+__attribute__((always_inline)) INLINE static float forcing_terms_timestep(
+ double time, const struct forcing_terms* terms,
+ const struct phys_const* phys_const, const struct part* p,
+ const struct xpart* xp) {
+
+ return FLT_MAX;
+}
+
+/**
+ * @brief Prints the properties of the forcing terms to stdout.
+ *
+ * @param terms The #forcing_terms properties of the run.
+ */
+static INLINE void forcing_terms_print(const struct forcing_terms* terms) {
+
+ message(
+ "Forcing terms is 'Roberts flow using accelerations'. U0: %.5f. nu: "
+ "%.5f. Vz factor: %.5f.",
+ terms->u0, terms->nu, terms->Vz_factor);
+}
+
+/**
+ * @brief Initialises the forcing term properties
+ *
+ * @param parameter_file The parsed parameter file
+ * @param phys_const Physical constants in internal units
+ * @param us The current internal system of units
+ * @param s The #space object.
+ * @param terms The forcing term properties to initialize
+ */
+static INLINE void forcing_terms_init(struct swift_params* parameter_file,
+ const struct phys_const* phys_const,
+ const struct unit_system* us,
+ const struct space* s,
+ struct forcing_terms* terms) {
+
+ terms->u0 = parser_get_param_double(parameter_file,
+ "RobertsFlowAccelerationForcing:u0");
+ terms->nu = parser_get_param_double(parameter_file,
+ "RobertsFlowAccelerationForcing:nu");
+ terms->Vz_factor = parser_get_opt_param_float(
+ parameter_file, "RobertsFlowAccelerationForcing:Vz_factor", 1.f);
+}
+
+#endif /* SWIFT_FORCING_ROBERTS_FLOW_ACCELERATION_H */
diff --git a/src/runner_others.c b/src/runner_others.c
index 2112a7a2b4..c70bc0ccb5 100644
--- a/src/runner_others.c
+++ b/src/runner_others.c
@@ -48,6 +48,7 @@
#include "error.h"
#include "feedback.h"
#include "fof.h"
+#include "forcing.h"
#include "gravity.h"
#include "hydro.h"
#include "potential.h"
@@ -634,12 +635,14 @@ void runner_do_end_hydro_force(struct runner *r, struct cell *c, int timer) {
const struct cosmology *cosmo = e->cosmology;
const int count = c->hydro.count;
struct part *restrict parts = c->hydro.parts;
+ struct xpart *restrict xparts = c->hydro.xparts;
/* Loop over the gas particles in this cell. */
for (int k = 0; k < count; k++) {
/* Get a handle on the part. */
struct part *restrict p = &parts[k];
+ struct xpart *restrict xp = &xparts[k];
double dt = 0;
if (part_is_active(p, e)) {
@@ -661,6 +664,10 @@ void runner_do_end_hydro_force(struct runner *r, struct cell *c, int timer) {
timestep_limiter_end_force(p);
chemistry_end_force(p, cosmo, with_cosmology, e->time, dt);
+ /* Apply the forcing terms (if any) */
+ forcing_terms_apply(e->time, e->forcing_terms, e->s,
+ e->physical_constants, p, xp);
+
#ifdef SWIFT_BOUNDARY_PARTICLES
/* Get the ID of the part */
diff --git a/src/swift.h b/src/swift.h
index 9cee769a6e..d0dc2e890f 100644
--- a/src/swift.h
+++ b/src/swift.h
@@ -46,6 +46,7 @@
#include "feedback.h"
#include "feedback_properties.h"
#include "fof.h"
+#include "forcing.h"
#include "gravity.h"
#include "gravity_derivatives.h"
#include "gravity_properties.h"
diff --git a/src/timestep.h b/src/timestep.h
index 851befc3b5..a77a2ebdc4 100644
--- a/src/timestep.h
+++ b/src/timestep.h
@@ -25,6 +25,7 @@
/* Local headers. */
#include "cooling.h"
#include "debug.h"
+#include "forcing.h"
#include "potential.h"
#include "rt.h"
#include "timeline.h"
@@ -173,14 +174,18 @@ __attribute__((always_inline)) INLINE static integertime_t get_part_timestep(
new_dt_grav = min(new_dt_self_grav, new_dt_ext_grav);
}
+ /* Compute the next timestep (forcing terms condition) */
+ const float new_dt_forcing = forcing_terms_timestep(
+ e->time, e->forcing_terms, e->physical_constants, p, xp);
+
/* Compute the next timestep (chemistry condition, e.g. diffusion) */
const float new_dt_chemistry =
chemistry_timestep(e->physical_constants, e->cosmology, e->internal_units,
e->hydro_properties, e->chemistry, p);
/* Take the minimum of all */
- float new_dt = min5(new_dt_hydro, new_dt_cooling, new_dt_grav, new_dt_mhd,
- new_dt_chemistry);
+ float new_dt = min3(new_dt_hydro, new_dt_cooling, new_dt_grav);
+ new_dt = min4(new_dt, new_dt_mhd, new_dt_chemistry, new_dt_forcing);
/* Limit change in smoothing length */
const float dt_h_change =
diff --git a/swift.c b/swift.c
index f877b0b16f..4507ac055c 100644
--- a/swift.c
+++ b/swift.c
@@ -90,6 +90,7 @@ int main(int argc, char *argv[]) {
struct cooling_function_data cooling_func;
struct cosmology cosmo;
struct external_potential potential;
+ struct forcing_terms forcing_terms;
struct extra_io_properties extra_io_props;
struct star_formation starform;
struct pm_mesh mesh;
@@ -1416,6 +1417,11 @@ int main(int argc, char *argv[]) {
potential_init(params, &prog_const, &us, &s, &potential);
if (myrank == 0) potential_print(&potential);
+ /* Initialise the forcing terms */
+ bzero(&forcing_terms, sizeof(struct forcing_terms));
+ forcing_terms_init(params, &prog_const, &us, &s, &forcing_terms);
+ if (myrank == 0) forcing_terms_print(&forcing_terms);
+
/* Initialise the long-range gravity mesh */
if (with_self_gravity && periodic) {
#ifdef HAVE_FFTW
@@ -1534,9 +1540,9 @@ int main(int argc, char *argv[]) {
&gravity_properties, &stars_properties, &black_holes_properties,
&sink_properties, &neutrino_properties, &neutrino_response,
&feedback_properties, &pressure_floor_props, &rt_properties,
- &mesh, &pow_data, &potential, &cooling_func, &starform,
- &chemistry, &extra_io_props, &fof_properties, &los_properties,
- &lightcone_array_properties, &ics_metadata);
+ &mesh, &pow_data, &potential, &forcing_terms, &cooling_func,
+ &starform, &chemistry, &extra_io_props, &fof_properties,
+ &los_properties, &lightcone_array_properties, &ics_metadata);
engine_config(/*restart=*/0, /*fof=*/0, &e, params, nr_nodes, myrank,
nr_threads, nr_pool_threads, with_aff, talking, restart_dir,
restart_file, &reparttype);
diff --git a/swift_fof.c b/swift_fof.c
index 996015bae3..af7aefd784 100644
--- a/swift_fof.c
+++ b/swift_fof.c
@@ -660,6 +660,7 @@ int main(int argc, char *argv[]) {
/*neutrino_response=*/NULL, /*feedback_properties=*/NULL,
/*pressure_floor_properties=*/NULL,
/*rt_properties=*/NULL, &mesh, /*pow_data=*/NULL, /*potential=*/NULL,
+ /*forcing_terms=*/NULL,
/*cooling_func=*/NULL, /*starform=*/NULL, /*chemistry=*/NULL,
/*extra_io_props=*/NULL, &fof_properties, /*los_properties=*/NULL,
/*lightcone_properties=*/NULL, &ics_metadata);
--
GitLab