diff --git a/src/feedback/EAGLE_thermal/feedback.h b/src/feedback/EAGLE_thermal/feedback.h
index fa546289fca1e3921f8fb1e87996fb78eac23dd8..27bb0c78f13d9d256537f0bc0e081d17736f8e08 100644
--- a/src/feedback/EAGLE_thermal/feedback.h
+++ b/src/feedback/EAGLE_thermal/feedback.h
@@ -217,7 +217,8 @@ __attribute__((always_inline)) INLINE static void feedback_prepare_feedback(
const float h_inv_dim = pow_dimension(h_inv); /* 1/h^d */
sp->feedback_data.to_collect.ngb_rho *= h_inv_dim;
- const float rho_inv = 1.f / sp->feedback_data.to_collect.ngb_rho;
+ const float rho = sp->feedback_data.to_collect.ngb_rho;
+ const float rho_inv = rho != 0.f ? 1.f / rho : 0.f;
sp->feedback_data.to_collect.ngb_Z *= h_inv_dim * rho_inv;
/* Compute amount of enrichment and feedback that needs to be done in this
diff --git a/src/forcing/roberts_flow/forcing.h b/src/forcing/roberts_flow/forcing.h
index 5c8f38e147fae55396fe41aca9cd302a1b21c263..73dc4c161fea3416a7b9519b50338ca3729adfbc 100644
--- a/src/forcing/roberts_flow/forcing.h
+++ b/src/forcing/roberts_flow/forcing.h
@@ -33,6 +33,15 @@
#include "space.h"
#include "units.h"
+/* Type of flow */
+enum flow {
+ Brandenburg_flow,
+ Roberts_flow_1,
+ Roberts_flow_2,
+ Roberts_flow_3,
+ Roberts_flow_4
+};
+
/**
* @brief Forcing Term Properties
*/
@@ -43,6 +52,12 @@ struct forcing_terms {
/*! Velocity scaling along the z direction */
float Vz_factor;
+
+ /*! Kind of RobertsFlow */
+ enum flow Flow_kind;
+
+ /*! Wavenumber of the flow*/
+ float kv;
};
/**
@@ -63,19 +78,82 @@ __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) {
+ enum flow Flow_kind = terms->Flow_kind;
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 k0 = (2. * M_PI / L) * terms->kv;
const double kf = M_SQRT2 * k0;
+ double v_Rob[3];
+ double Psi;
+
+ /* Switching between different kinds of flows */
+ /* Note: Roberts worked in yz plane, we work in xy plane just for convenience
+ * and also because A.Brandenburg has flows in xy plane, so our formulas
+ * differ from Robets article by several rotations. Theese rotations are
+ * equivalent to yzx -> xyz permutation*/
+
+ switch (Flow_kind) {
+
+ case Brandenburg_flow:
+
+ /* Eq. 8 of Tilgner & Brandenburg, 2008, MNRAS, 391, 1477 */
+ // Psi = (u0 / k0) * cos(k0 * p->x[0]) * cos(k0 * p->x[1]);
+
+ /* Eq. 7 of Tilgner & Brandenburg, 2008, MNRAS, 391, 1477 */
+ // v_Rob[0] = u0 * cos(k0 * p->x[0]) * sin(k0 * p->x[1]);
+ // v_Rob[1] = -u0 * sin(k0 * p->x[0]) * cos(k0 * p->x[1]);
+ // v_Rob[2] = kf * Psi;
+
+ // Velocity used to compare with A.B. runs (from overleaf)
+ Psi = (u0 / k0) * sin(k0 * p->x[0]) * sin(k0 * p->x[1]);
+ v_Rob[0] = u0 * sin(k0 * p->x[0]) * cos(k0 * p->x[1]);
+ v_Rob[1] = -u0 * cos(k0 * p->x[0]) * sin(k0 * p->x[1]);
+ v_Rob[2] = kf * Psi;
+
+ break;
+
+ case Roberts_flow_1:
+ /* Eq. 5.1 of Roberts, Feb. 3, 1972, Vol. 271, No. 1216 (Feb. 3, 1972),
+ * pp. 411-454.*/
+ v_Rob[0] = u0 * sin(k0 * p->x[0]);
+ v_Rob[1] = u0 * sin(k0 * p->x[1]);
+ v_Rob[2] = u0 * (cos(k0 * p->x[0]) - cos(k0 * p->x[1]));
+ break;
+
+ case Roberts_flow_2:
+
+ /* Eq. 6.1 of Roberts, Feb. 3, 1972, Vol. 271, No. 1216 (Feb. 3, 1972),
+ pp. 411-454.*/
+ v_Rob[0] = u0 * sin(k0 * p->x[0]);
+ v_Rob[1] = u0 * sin(k0 * p->x[1]);
+ v_Rob[2] = u0 * (cos(k0 * p->x[0]) + cos(k0 * p->x[1]));
+ break;
+
+ case Roberts_flow_3:
+
+ /* Eq. 6.2 of Roberts, Feb. 3, 1972, Vol. 271, No. 1216 (Feb. 3, 1972),
+ pp. 411-454.*/
+ v_Rob[0] = u0 * sin(k0 * p->x[0]);
+ v_Rob[1] = u0 * sin(k0 * p->x[1]);
+ v_Rob[2] = u0 * 2 * cos(k0 * p->x[0]) * cos(k0 * p->x[1]);
+ break;
+
+ case Roberts_flow_4:
+
+ /* Eq. 6.3 of Roberts, Feb. 3, 1972, Vol. 271, No. 1216 (Feb. 3, 1972),
+ pp. 411-454.*/
+ v_Rob[0] = u0 * sin(k0 * p->x[0]);
+ v_Rob[1] = u0 * sin(k0 * p->x[1]);
+ v_Rob[2] = u0 * sin(k0 * (p->x[0] + p->x[1]));
+ break;
- /* Eq. 8 of Tilgner & Brandenburg, 2008, MNRAS, 391, 1477 */
- const double Psi = (u0 / k0) * cos(k0 * p->x[0]) * cos(k0 * p->x[1]);
+ default:
- /* 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};
+ v_Rob[0] = 0.f;
+ v_Rob[1] = 0.f;
+ v_Rob[2] = 0.f;
+ }
/* Force the velocity and possibly scale the z-direction */
xp->v_full[0] = v_Rob[0];
@@ -111,6 +189,7 @@ 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);
+ message("Forcing 'Roberts flow' Kind: %i .", terms->Flow_kind);
}
/**
@@ -131,6 +210,14 @@ static INLINE void forcing_terms_init(struct swift_params* parameter_file,
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);
+ terms->Flow_kind =
+ parser_get_param_int(parameter_file, "RobertsFlowForcing:Flow_kind");
+ terms->kv = parser_get_param_double(parameter_file, "RobertsFlowForcing:kv");
+
+ if (terms->Flow_kind > 4 || terms->Flow_kind < 0)
+ error(
+ "Error: Flow_kind variable can take integer values from [0,4] "
+ "interval.");
}
#endif /* SWIFT_FORCING_ROBERTS_FLOW_H */
diff --git a/src/forcing/roberts_flow_acceleration/forcing.h b/src/forcing/roberts_flow_acceleration/forcing.h
index 9a9d33d3c88e8d855081474eca40b9429ce8eb30..d739e7e3a9ed3538dd84c1fc59f933fec80d6353 100644
--- a/src/forcing/roberts_flow_acceleration/forcing.h
+++ b/src/forcing/roberts_flow_acceleration/forcing.h
@@ -34,6 +34,15 @@
#include "space.h"
#include "units.h"
+/* Type of flow */
+enum flow {
+ Brandenburg_flow,
+ Roberts_flow_1,
+ Roberts_flow_2,
+ Roberts_flow_3,
+ Roberts_flow_4
+};
+
/**
* @brief Forcing Term Properties
*/
@@ -47,6 +56,12 @@ struct forcing_terms {
/*! Velocity scaling along the z direction */
float Vz_factor;
+
+ /*! Kind of RobertsFlow */
+ enum flow Flow_kind;
+
+ /*! Wavenumber of the flow*/
+ float kv;
};
/**
@@ -65,22 +80,79 @@ __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);
+ enum flow Flow_kind = terms->Flow_kind;
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;
+ const float c_s = hydro_get_comoving_soundspeed(p);
- /* Eq. 8 of Tilgner & Brandenburg, 2008, MNRAS, 391, 1477 */
- const double Psi = (u0 / k0) * cos(k0 * p->x[0]) * cos(k0 * p->x[1]);
+ /* Effective viscosity from artificial viscosity, as in eq. 100 from
+ * arXiv:1012.1885 */
+ const float nu = terms->nu * p->viscosity.alpha * c_s * p->h /
+ (2.f * (hydro_dimension + 2.f));
- /* 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};
+ const float Vz_factor = terms->Vz_factor;
+ const double k0 = (2. * M_PI / L) * terms->kv;
+ const double kf = M_SQRT2 * k0;
+ double v_Rob[3];
+ double Psi;
+
+ switch (Flow_kind) {
+
+ case Brandenburg_flow:
+ /* Eq. 8 of Tilgner & Brandenburg, 2008, MNRAS, 391, 1477 */
+ // Psi = (u0 / k0) * cos(k0 * p->x[0]) * cos(k0 * p->x[1]);
+
+ /* Eq. 7 of Tilgner & Brandenburg, 2008, MNRAS, 391, 1477 */
+ // v_Rob[0] = u0 * cos(k0 * p->x[0]) * sin(k0 * p->x[1]);
+ // v_Rob[1] = -u0 * sin(k0 * p->x[0]) * cos(k0 * p->x[1]);
+ // v_Rob[2] = kf * Psi;
+
+ // Velocity used to compare with A.B. runs (from overleaf)
+ Psi = (u0 / k0) * sin(k0 * p->x[0]) * sin(k0 * p->x[1]);
+ v_Rob[0] = u0 * sin(k0 * p->x[0]) * cos(k0 * p->x[1]);
+ v_Rob[1] = -u0 * cos(k0 * p->x[0]) * sin(k0 * p->x[1]);
+ v_Rob[2] = kf * Psi;
+
+ break;
+
+ case Roberts_flow_1:
+ /* Eq. 5.1 of Roberts, Feb. 3, 1972, Vol. 271, No. 1216 (Feb. 3, 1972),
+ * pp. 411-454.*/
+ v_Rob[0] = u0 * sin(k0 * p->x[0]);
+ v_Rob[1] = u0 * sin(k0 * p->x[1]);
+ v_Rob[2] = u0 * (cos(k0 * p->x[0]) - cos(k0 * p->x[1]));
+ break;
+
+ case Roberts_flow_2:
+ /* Eq. 6.1 of Roberts, Feb. 3, 1972, Vol. 271, No. 1216 (Feb. 3, 1972),
+ * pp. 411-454.*/
+ v_Rob[0] = u0 * sin(k0 * p->x[0]);
+ v_Rob[1] = u0 * sin(k0 * p->x[1]);
+ v_Rob[2] = u0 * (cos(k0 * p->x[0]) + cos(k0 * p->x[1]));
+ break;
+
+ case Roberts_flow_3:
+ /* Eq. 6.2 of Roberts, Feb. 3, 1972, Vol. 271, No. 1216 (Feb. 3, 1972),
+ * pp. 411-454.*/
+ v_Rob[0] = u0 * sin(k0 * p->x[0]);
+ v_Rob[1] = u0 * sin(k0 * p->x[1]);
+ v_Rob[2] = u0 * 2 * cos(k0 * p->x[0]) * cos(k0 * p->x[1]);
+ break;
+
+ case Roberts_flow_4:
+ /* Eq. 6.3 of Roberts, Feb. 3, 1972, Vol. 271, No. 1216 (Feb. 3, 1972),
+ * pp. 411-454.*/
+ v_Rob[0] = u0 * sin(k0 * p->x[0]);
+ v_Rob[1] = u0 * sin(k0 * p->x[1]);
+ v_Rob[2] = u0 * sin(k0 * (p->x[0] + p->x[1]));
+ break;
+
+ default:
+
+ v_Rob[0] = 0.f;
+ v_Rob[1] = 0.f;
+ v_Rob[2] = 0.f;
+ }
/* 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],
@@ -89,14 +161,7 @@ __attribute__((always_inline)) INLINE static void forcing_terms_apply(
/* 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];
- }
+ p->a_hydro[2] += f[2] * Vz_factor;
}
/**
@@ -129,6 +194,7 @@ static INLINE void forcing_terms_print(const struct forcing_terms* terms) {
"Forcing terms is 'Roberts flow using accelerations'. U0: %.5f. nu: "
"%.5f. Vz factor: %.5f.",
terms->u0, terms->nu, terms->Vz_factor);
+ message("Forcing 'Roberts flow' Kind: %i .", terms->Flow_kind);
}
/**
@@ -152,6 +218,14 @@ static INLINE void forcing_terms_init(struct swift_params* parameter_file,
"RobertsFlowAccelerationForcing:nu");
terms->Vz_factor = parser_get_opt_param_float(
parameter_file, "RobertsFlowAccelerationForcing:Vz_factor", 1.f);
+ terms->Flow_kind =
+ parser_get_param_int(parameter_file, "RobertsFlowForcing:Flow_kind");
+ terms->kv = parser_get_param_double(parameter_file, "RobertsFlowForcing:kv");
+
+ if (terms->Flow_kind > 4 || terms->Flow_kind < 0)
+ error(
+ "Error: Flow_kind variable can take integer values from [0,4] "
+ "interval.");
}
#endif /* SWIFT_FORCING_ROBERTS_FLOW_ACCELERATION_H */