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

Added forcing from stencil_process.sps to autotests. 3 Tests fail.

Browse Source
This commit is contained in:
jpekkila
2019-08-06 19:15:28 +03:00
parent 0b7f43da91
commit d7e26e8f21
7 changed files with 70 additions and 26 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
acc/mhd_solver/stencil_defines.h
View File

@@ -30,7 +30,7 @@
#define LMAGNETIC (1)
#define LENTROPY (1)
#define LTEMPERATURE (0)
#define LFORCING (0)
#define LFORCING (1)
#define LUPWD (0)
#define AC_THERMAL_CONDUCTIVITY (AcReal(0.001)) // TODO: make an actual config parameter

2
config/astaroth.conf
View File

@@ -40,7 +40,7 @@ AC_chi = 0.0001
AC_relhel = 0.0
AC_forcing_magnitude = 1e-5
AC_kmin = 0.8
AC_kmax = 1.2
AC_kmax = 1.2
// Entropy

6
src/standalone/autotest.cc
View File

@@ -422,6 +422,12 @@ check_rk3(const AcMeshInfo& mesh_info)
// const AcReal dt = host_timestep(umax, mesh_info);
const AcReal dt = AcReal(1e-2); // Use a small constant timestep to avoid instabilities
#if LFORCING
const ForcingParams forcing_params = generateForcingParams(model_mesh->info);
loadForcingParamsToHost(forcing_params, model_mesh);
loadForcingParamsToDevice(forcing_params);
#endif
acIntegrate(dt);
model_rk3(dt, model_mesh);

24
src/standalone/model/host_forcing.cc
View File

@@ -255,6 +255,30 @@ loadForcingParamsToHost(const ForcingParams& forcing_params, AcMesh* mesh)
mesh->info.real_params[AC_kaver] = forcing_params.kaver;
}
void
loadForcingParamsToHost(const ForcingParams& forcing_params, ModelMesh* mesh)
{
// %JP: Left some regex magic here in case we need to modify the ForcingParams struct
// acLoadDeviceConstant\(([A-Za-z_]*), ([a-z_.]*)\);
// mesh->info.real_params[$1] = $2;
mesh->info.real_params[AC_forcing_magnitude] = forcing_params.magnitude;
mesh->info.real_params[AC_forcing_phase] = forcing_params.phase;
mesh->info.real_params[AC_k_forcex] = forcing_params.k_force.x;
mesh->info.real_params[AC_k_forcey] = forcing_params.k_force.y;
mesh->info.real_params[AC_k_forcez] = forcing_params.k_force.z;
mesh->info.real_params[AC_ff_hel_rex] = forcing_params.ff_hel_re.x;
mesh->info.real_params[AC_ff_hel_rey] = forcing_params.ff_hel_re.y;
mesh->info.real_params[AC_ff_hel_rez] = forcing_params.ff_hel_re.z;
mesh->info.real_params[AC_ff_hel_imx] = forcing_params.ff_hel_im.x;
mesh->info.real_params[AC_ff_hel_imy] = forcing_params.ff_hel_im.y;
mesh->info.real_params[AC_ff_hel_imz] = forcing_params.ff_hel_im.z;
mesh->info.real_params[AC_kaver] = forcing_params.kaver;
}
ForcingParams
generateForcingParams(const AcMeshInfo& mesh_info)
{

3
src/standalone/model/host_forcing.h
View File

@@ -28,6 +28,8 @@
#pragma once
#include "astaroth.h"
#include "modelmesh.h"
AcReal get_random_number_01();
AcReal3 cross(const AcReal3& a, const AcReal3& b);
@@ -64,5 +66,6 @@ typedef struct {
void loadForcingParamsToDevice(const ForcingParams& forcing_params);
void loadForcingParamsToHost(const ForcingParams& forcing_params, AcMesh* mesh);
void loadForcingParamsToHost(const ForcingParams& forcing_params, ModelMesh* mesh);
ForcingParams generateForcingParams(const AcMeshInfo& mesh_info);

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

@@ -660,15 +660,13 @@ is_valid(const ModelVector& a)
}
#if LFORCING
// FORCING NOT SUPPORTED FOR AUTOTEST
static inline ModelVector
ModelVector
simple_vortex_forcing(ModelVector a, ModelVector b, ModelScalar magnitude)
{
return magnitude * cross(normalized(b - a), (ModelVector){0, 0, 1}); // Vortex
}
static inline ModelVector
ModelVector
simple_outward_flow_forcing(ModelVector a, ModelVector b, ModelScalar magnitude)
{
return magnitude * (1 / length(b - a)) * normalized(b - a); // Outward flow
@@ -676,22 +674,24 @@ simple_outward_flow_forcing(ModelVector a, ModelVector b, ModelScalar magnitude)
// The Pencil Code forcing_hel_noshear(), manual Eq. 222, inspired forcing function with adjustable
// helicity
static inline ModelVector
helical_forcing(ModelScalar /* magnitude */, ModelVector k_force, ModelVector xx, ModelVector ff_re,
ModelVector
helical_forcing(ModelScalar magnitude, ModelVector k_force, ModelVector xx, ModelVector ff_re,
ModelVector ff_im, ModelScalar phi)
{
(void)magnitude; // WARNING: unused
xx.x = xx.x * (2.0l * M_PI / (get(AC_dsx) * (get(AC_ny_max) - get(AC_ny_min))));
xx.y = xx.y * (2.0l * M_PI / (get(AC_dsy) * (get(AC_ny_max) - get(AC_ny_min))));
xx.z = xx.z * (2.0l * M_PI / (get(AC_dsz) * (get(AC_ny_max) - get(AC_ny_min))));
xx.x = xx.x * (2.0 * M_PI / (get(AC_dsx) * (get(AC_ny_max) - get(AC_ny_min))));
xx.y = xx.y * (2.0 * M_PI / (get(AC_dsy) * (get(AC_ny_max) - get(AC_ny_min))));
xx.z = xx.z * (2.0 * M_PI / (get(AC_dsz) * (get(AC_ny_max) - get(AC_ny_min))));
ModelScalar cos_phi = cosl(phi);
ModelScalar sin_phi = sinl(phi);
ModelScalar cos_k_dox_x = cosl(dot(k_force, xx));
ModelScalar sin_k_dox_x = sinl(dot(k_force, xx));
ModelScalar cosl_phi = cosl(phi);
ModelScalar sinl_phi = sinl(phi);
ModelScalar cosl_k_dox_x = cosl(dot(k_force, xx));
ModelScalar sinl_k_dox_x = sinl(dot(k_force, xx));
// Phase affect only the x-component
ModelScalar real_comp_phase = cos_k_dox_x * cos_phi - sin_k_dox_x * sin_phi;
ModelScalar imag_comp_phase = cos_k_dox_x * sin_phi + sin_k_dox_x * cos_phi;
// ModelScalar real_comp = cosl_k_dox_x;
// ModelScalar imag_comp = sinl_k_dox_x;
ModelScalar real_comp_phase = cosl_k_dox_x * cosl_phi - sinl_k_dox_x * sinl_phi;
ModelScalar imag_comp_phase = cosl_k_dox_x * sinl_phi + sinl_k_dox_x * cosl_phi;
ModelVector force = (ModelVector){ff_re.x * real_comp_phase - ff_im.x * imag_comp_phase,
ff_re.y * real_comp_phase - ff_im.y * imag_comp_phase,
@@ -700,18 +700,17 @@ helical_forcing(ModelScalar /* magnitude */, ModelVector k_force, ModelVector xx
return force;
}
static inline ModelVector
ModelVector
forcing(int3 globalVertexIdx, ModelScalar dt)
{
/*
ModelVector a = ModelScalar(.5) * (ModelVector){get(AC_nx) * get(AC_dsx),
ModelVector a = ModelScalar(.5) * (ModelVector){get(AC_nx) * get(AC_dsx),
get(AC_ny) * get(AC_dsy),
get(AC_nz) * get(AC_dsz)}; // source (origin)
*/
ModelVector xx = (ModelVector){
(globalVertexIdx.x - get(AC_nx_min)) * get(AC_dsx),
(globalVertexIdx.y - get(AC_ny_min) * get(AC_dsy)),
(globalVertexIdx.z - get(AC_nz_min) * get(AC_dsz))}; // sink (current index)
(void)a; // WARNING: not used
ModelVector xx = (ModelVector){(globalVertexIdx.x - get(AC_nx_min)) * get(AC_dsx),
(globalVertexIdx.y - get(AC_ny_min)) * get(AC_dsy),
(globalVertexIdx.z - get(AC_nz_min)) *
get(AC_dsz)}; // sink (current index)
const ModelScalar cs2 = get(AC_cs2_sound);
const ModelScalar cs = sqrtl(cs2);
@@ -722,6 +721,11 @@ forcing(int3 globalVertexIdx, ModelScalar dt)
ModelVector ff_re = (ModelVector){get(AC_ff_hel_rex), get(AC_ff_hel_rey), get(AC_ff_hel_rez)};
ModelVector ff_im = (ModelVector){get(AC_ff_hel_imx), get(AC_ff_hel_imy), get(AC_ff_hel_imz)};
(void)phase; // WARNING: unused with simple forcing. Should be defined in helical_forcing
(void)k_force; // WARNING: unused with simple forcing. Should be defined in helical_forcing
(void)ff_re; // WARNING: unused with simple forcing. Should be defined in helical_forcing
(void)ff_im; // WARNING: unused with simple forcing. Should be defined in helical_forcing
// Determine that forcing funtion type at this point.
// ModelVector force = simple_vortex_forcing(a, xx, magnitude);
// ModelVector force = simple_outward_flow_forcing(a, xx, magnitude);

7
src/standalone/renderer.cc
View File

@@ -34,6 +34,7 @@
#include "config_loader.h"
#include "core/errchk.h"
#include "core/math_utils.h"
#include "model/host_forcing.h"
#include "model/host_memory.h"
#include "model/host_timestep.h"
#include "model/model_reduce.h"
@@ -383,6 +384,12 @@ run_renderer(void)
#if 1
const AcReal umax = acReduceVec(RTYPE_MAX, VTXBUF_UUX, VTXBUF_UUY, VTXBUF_UUZ);
const AcReal dt = host_timestep(umax, mesh_info);
#if LFORCING
const ForcingParams forcing_params = generateForcingParams(mesh->info);
loadForcingParamsToDevice(forcing_params);
#endif
acIntegrate(dt);
#else
ModelMesh* model_mesh = modelmesh_create(mesh->info);