Skip to content
Snippets Groups Projects
Commit 5fd77838 authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Port fixes to the forcing terms from MHD_canvas

parent 43c583a7
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
src/feedback/EAGLE_thermal/feedback.h
+ 2
− 1
View file @ 5fd77838
......@@ -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
......
src/forcing/roberts_flow/forcing.h
+ 94
− 7
View file @ 5fd77838
......@@ -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 */
src/forcing/roberts_flow_acceleration/forcing.h
+ 94
− 20
View file @ 5fd77838
......@@ -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 */
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment