cwpearson/astaroth
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Made the simple forcing to scale with dt s.t. it does not explode so easily

Browse Source
This commit is contained in:
jpekkila
2019-06-19 16:34:23 +03:00
parent f3cbc4984c
commit a7515fbbd7
2 changed files with 10 additions and 10 deletions
Download Patch File Download Diff File Expand all files Collapse all files

4
acc/mhd_solver/stencil_process.sps
View File

@@ -196,7 +196,7 @@ forcing(int3 globalVertexIdx)
(globalVertexIdx.y - ny_min) * dsy, (globalVertexIdx.y - ny_min) * dsy,
(globalVertexIdx.z - nz_min) * dsz}; // sink (current index) (globalVertexIdx.z - nz_min) * dsz}; // sink (current index)
Scalar magnitude = .1; Scalar magnitude = 0.05;
// Vector c = magnitude * (1 / length(b - a)) * normalized(b - a); // Outward flow // Vector c = magnitude * (1 / length(b - a)) * normalized(b - a); // Outward flow
Vector c = magnitude * cross(normalized(b - a), (Vector){0, 0, 1}); // Vortex Vector c = magnitude * cross(normalized(b - a), (Vector){0, 0, 1}); // Vortex
if (is_valid(c)) { return c; } if (is_valid(c)) { return c; }
@@ -248,7 +248,7 @@ solve(Scalar dt) {
#if LFORCING #if LFORCING
if (step_number == 2) { if (step_number == 2) {
out_uu = out_uu + forcing(globalVertexIdx); out_uu = out_uu + dt * forcing(globalVertexIdx);
} }
#endif #endif
} }

16
src/standalone/model/model_rk3.cc
View File

@@ -720,7 +720,7 @@ forcing(int3 globalVertexIdx)
(globalVertexIdx.y - get(AC_ny_min)) * get(AC_dsy), (globalVertexIdx.y - get(AC_ny_min)) * get(AC_dsy),
(globalVertexIdx.z - get(AC_nz_min)) * get(AC_dsz)}; (globalVertexIdx.z - get(AC_nz_min)) * get(AC_dsz)};
ModelScalar magnitude = .1; ModelScalar magnitude = 0.05;
// Vector c = magnitude * (1 / length(b - a)) * normalized(b - a); // Outward flow // Vector c = magnitude * (1 / length(b - a)) * normalized(b - a); // Outward flow
ModelVector c = magnitude * cross(normalized(b - a), (ModelVector){0, 0, 1}); // Vortex ModelVector c = magnitude * cross(normalized(b - a), (ModelVector){0, 0, 1}); // Vortex
if (is_valid(c)) { if (is_valid(c)) {
@@ -780,8 +780,8 @@ solve_alpha_step(const int step_number, const ModelScalar dt, const int i, const
} }
static void static void
solve_beta_step(const int step_number, const int i, const int j, const int k, const ModelMesh& in, solve_beta_step(const int step_number, const ModelScalar dt, const int i, const int j, const int k,
ModelMesh* out) const ModelMesh& in, ModelMesh* out)
{ {
const int idx = AC_VTXBUF_IDX(i, j, k, in.info); const int idx = AC_VTXBUF_IDX(i, j, k, in.info);
@@ -795,9 +795,9 @@ solve_beta_step(const int step_number, const int i, const int j, const int k, co
#if LFORCING #if LFORCING
if (step_number == 2) { if (step_number == 2) {
ModelVector force = forcing((int3){i, j, k}); ModelVector force = forcing((int3){i, j, k});
out->vertex_buffer[VTXBUF_UUX][idx] += force.x; out->vertex_buffer[VTXBUF_UUX][idx] += force.x * dt;
out->vertex_buffer[VTXBUF_UUY][idx] += force.y; out->vertex_buffer[VTXBUF_UUY][idx] += force.y * dt;
out->vertex_buffer[VTXBUF_UUZ][idx] += force.z; out->vertex_buffer[VTXBUF_UUZ][idx] += force.z * dt;
} }
#endif #endif
} }
@@ -822,7 +822,7 @@ model_rk3_step(const int step_number, const ModelScalar dt, ModelMesh* mesh)
for (int k = get(AC_nz_min); k < get(AC_nz_max); ++k) { for (int k = get(AC_nz_min); k < get(AC_nz_max); ++k) {
for (int j = get(AC_ny_min); j < get(AC_ny_max); ++j) { for (int j = get(AC_ny_min); j < get(AC_ny_max); ++j) {
for (int i = get(AC_nx_min); i < get(AC_nx_max); ++i) { for (int i = get(AC_nx_min); i < get(AC_nx_max); ++i) {
solve_beta_step(step_number, i, j, k, *tmp, mesh); solve_beta_step(step_number, dt, i, j, k, *tmp, mesh);
} }
} }
} }
@@ -852,7 +852,7 @@ model_rk3(const ModelScalar dt, ModelMesh* mesh)
for (int k = get(AC_nz_min); k < get(AC_nz_max); ++k) { for (int k = get(AC_nz_min); k < get(AC_nz_max); ++k) {
for (int j = get(AC_ny_min); j < get(AC_ny_max); ++j) { for (int j = get(AC_ny_min); j < get(AC_ny_max); ++j) {
for (int i = get(AC_nx_min); i < get(AC_nx_max); ++i) { for (int i = get(AC_nx_min); i < get(AC_nx_max); ++i) {
solve_beta_step(step_number, i, j, k, *tmp, mesh); solve_beta_step(step_number, dt, i, j, k, *tmp, mesh);
} }
} }
} }