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

Merged master to acc_parameter_overhaul

Browse Source
This commit is contained in:
jpekkila
2019-07-03 17:37:37 +03:00
parent d9b07a59b5 e8a5579b50
commit acc53e1c2b
13 changed files with 527 additions and 145 deletions
Download Patch File Download Diff File Expand all files Collapse all files

17
CMakeLists.txt
View File

@@ -128,19 +128,14 @@ set(CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS_RELEASE}\
set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG}\ set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG}\
-O0 -g") -O0 -g")
set (CXX_FLAGS_WARNING "-Wall -Wextra -Werror -Wno-error=unused-parameter\ if ("${CMAKE_CXX_COMPILER_ID}" STREQUAL "GNU")
-Wno-error=unused-function -Wno-error=unknown-pragmas") set (CXX_FLAGS_WARNING "-Wall -Wextra -Werror")# -Wdouble-promotion -Wfloat-conversion")
elseif ("${CMAKE_CXX_COMPILER_ID}" STREQUAL "Intel")
# Also warn about implicit conversions if the compiler supports it set (CXX_FLAGS_WARNING "-Wall -Wextra -Werror-all -Wsign-conversion")
if (${CMAKE_CXX_COMPILER_ID} STREQUAL "GNU") else()
set (CXX_FLAGS_WARNING "${CXX_FLAGS_WARNING} -Wdouble-promotion -Wfloat-conversion") message(WARNING "Using an unknown compiler. Compilation warning flags were not set.")
endif() endif()
# Other flags. -D_FORCE_INLINES is a workaround to some CUDA/C++ "feature"
# which botches the compilation ("memcpy was not declared in this scope")
# (Not required with cc >= 3.0)
#set(CXX_FLAGS_ETC "-D_FORCE_INLINES")
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS}\ set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS}\
${CXX_FLAGS_WARNING}\ ${CXX_FLAGS_WARNING}\
${CXX_FLAGS_ETC}\ ${CXX_FLAGS_ETC}\

100
acc/mhd_solver/stencil_process.sps
View File

@@ -12,6 +12,9 @@ uniform Scalar dsx;
uniform Scalar dsy; uniform Scalar dsy;
uniform Scalar dsz; uniform Scalar dsz;
uniform Scalar lnT0;
uniform Scalar lnrho0;
uniform int nx_min; uniform int nx_min;
uniform int ny_min; uniform int ny_min;
uniform int nz_min; uniform int nz_min;
@@ -58,9 +61,10 @@ continuity(in Vector uu, in Scalar lnrho) {
Vector Vector
momentum(in Vector uu, in Scalar lnrho, in Scalar ss, in Vector aa) { momentum(in Vector uu, in Scalar lnrho, in Scalar ss, in Vector aa) {
const Matrix S = stress_tensor(uu); const Matrix S = stress_tensor(uu);
const Scalar cs2 = cs2_sound * exp(gamma * value(ss) / cp_sound + (gamma - 1) * (value(lnrho) - LNRHO0)); const Scalar cs2 = cs2_sound * exp(gamma * value(ss) / cp_sound + (gamma - 1) * (value(lnrho) - lnrho0));
const Vector j = (Scalar(1.) / mu0) * (gradient_of_divergence(aa) - laplace_vec(aa)); // Current density const Vector j = (Scalar(1.) / mu0) * (gradient_of_divergence(aa) - laplace_vec(aa)); // Current density
const Vector B = curl(aa); const Vector B = curl(aa);
//TODO: DOES INTHERMAL VERSTION INCLUDE THE MAGNETIC FIELD?
const Scalar inv_rho = Scalar(1.) / exp(value(lnrho)); const Scalar inv_rho = Scalar(1.) / exp(value(lnrho));
// Regex replace CPU constants with get\(AC_([a-zA-Z_0-9]*)\) // Regex replace CPU constants with get\(AC_([a-zA-Z_0-9]*)\)
@@ -100,20 +104,22 @@ momentum(in Vector uu, in Scalar lnrho, in Scalar tt) {
#else #else
Vector Vector
momentum(in Vector uu, in Scalar lnrho) { momentum(in Vector uu, in Scalar lnrho) {
// !!!!!!!!!!!!!!!!%JP: NOTE TODO IMPORTANT!!!!!!!!!!!!!!!!!!!!!!!! Vector mom;
// NOT CHECKED FOR CORRECTNESS: USE AT YOUR OWN RISK
const Matrix S = stress_tensor(uu); const Matrix S = stress_tensor(uu);
const Scalar cs2 = cs2_sound * exp((gamma - 1) * (value(lnrho) - LNRHO0));
// Regex replace CPU constants with get\(AC_([a-zA-Z_0-9]*)\) // Isothermal: we have constant speed of sound
// \1
const Vector mom = - mul(gradients(uu), value(uu)) mom = -mul(gradients(uu), value(uu)) -
- cs2 * gradient(lnrho) cs2_sound * gradient(lnrho) +
+ nu_visc * ( nu_visc *
laplace_vec(uu) (laplace_vec(uu) + Scalar(1. / 3.) * gradient_of_divergence(uu) +
+ Scalar(1. / 3.) * gradient_of_divergence(uu) Scalar(2.) * mul(S, gradient(lnrho))) + zeta * gradient_of_divergence(uu);
+ Scalar(2.) * mul(S, gradient(lnrho))
) #if LGRAVITY
+ zeta * gradient_of_divergence(uu); mom = mom - (Vector){0, 0, -10.0};
#endif
return mom; return mom;
} }
#endif #endif
@@ -139,8 +145,8 @@ induction(in Vector uu, in Vector aa) {
#if LENTROPY #if LENTROPY
Scalar Scalar
lnT( in Scalar ss, in Scalar lnrho) { lnT( in Scalar ss, in Scalar lnrho) {
const Scalar lnT = LNT0 + gamma * value(ss) / cp_sound + const Scalar lnT = lnT0 + gamma * value(ss) / cp_sound +
(gamma - Scalar(1.)) * (value(lnrho) - LNRHO0); (gamma - Scalar(1.)) * (value(lnrho) - lnrho0);
return lnT; return lnT;
} }
@@ -193,8 +199,12 @@ heat_transfer(in Vector uu, in Scalar lnrho, in Scalar tt)
} }
#endif #endif
#if LFORCING #if LFORCING
Vector Vector
simple_vortex_forcing(Vector a, Vector b, Scalar magnitude) simple_vortex_forcing(Vector a, Vector b, Scalar magnitude)
{ {
@@ -207,23 +217,67 @@ simple_outward_flow_forcing(Vector a, Vector b, Scalar magnitude)
return magnitude * (1 / length(b - a)) * normalized(b - a); // Outward flow return magnitude * (1 / length(b - a)) * normalized(b - a); // Outward flow
} }
// The Pencil Code forcing_hel_noshear(), manual Eq. 222, inspired forcing function with adjustable helicity
Vector Vector
forcing(int3 globalVertexIdx) helical_forcing(Scalar magnitude, Vector k_force, Vector xx, Vector ff_re, Vector ff_im, Scalar phi)
{
xx.x = xx.x*(2.0*M_PI/(dsx*(DCONST_INT(AC_ny_max) - DCONST_INT(AC_ny_min))));
xx.y = xx.y*(2.0*M_PI/(dsy*(DCONST_INT(AC_ny_max) - DCONST_INT(AC_ny_min))));
xx.z = xx.z*(2.0*M_PI/(dsz*(DCONST_INT(AC_ny_max) - DCONST_INT(AC_ny_min))));
Scalar cos_phi = cos(phi);
Scalar sin_phi = sin(phi);
Scalar cos_k_dox_x = cos(dot(k_force, xx));
Scalar sin_k_dox_x = sin(dot(k_force, xx));
// Phase affect only the x-component
//Scalar real_comp = cos_k_dox_x;
//Scalar imag_comp = sin_k_dox_x;
Scalar real_comp_phase = cos_k_dox_x*cos_phi - sin_k_dox_x*sin_phi;
Scalar imag_comp_phase = cos_k_dox_x*sin_phi + sin_k_dox_x*cos_phi;
Vector force = (Vector){ ff_re.x*real_comp_phase - ff_im.x*imag_comp_phase,
ff_re.y*real_comp_phase - ff_im.y*imag_comp_phase,
ff_re.z*real_comp_phase - ff_im.z*imag_comp_phase};
return force;
}
Vector
forcing(int3 globalVertexIdx, Scalar dt)
{ {
Vector a = Scalar(.5) * (Vector){globalGrid.n.x * dsx, Vector a = Scalar(.5) * (Vector){globalGrid.n.x * dsx,
globalGrid.n.y * dsy, globalGrid.n.y * dsy,
globalGrid.n.z * dsz}; // source (origin) globalGrid.n.z * dsz}; // source (origin)
Vector b = (Vector){(globalVertexIdx.x - nx_min) * dsx, Vector xx = (Vector){(globalVertexIdx.x - nx_min) * dsx,
(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)
const Scalar cs2 = cs2_sound;
const Scalar cs = sqrt(cs2);
//Placeholders until determined properly
Scalar magnitude = DCONST_REAL(AC_forcing_magnitude);
Scalar phase = DCONST_REAL(AC_forcing_phase);
Vector k_force = (Vector){ DCONST_REAL(AC_k_forcex), DCONST_REAL(AC_k_forcey), DCONST_REAL(AC_k_forcez)};
Vector ff_re = (Vector){DCONST_REAL(AC_ff_hel_rex), DCONST_REAL(AC_ff_hel_rey), DCONST_REAL(AC_ff_hel_rez)};
Vector ff_im = (Vector){DCONST_REAL(AC_ff_hel_imx), DCONST_REAL(AC_ff_hel_imy), DCONST_REAL(AC_ff_hel_imz)};
Scalar magnitude = 0.05;
//Determine that forcing funtion type at this point. //Determine that forcing funtion type at this point.
Vector c = simple_vortex_forcing(a, b, magnitude); //Vector force = simple_vortex_forcing(a, xx, magnitude);
//Vector c = simple_outward_flow_forcing(a, b, magnitude); //Vector force = simple_outward_flow_forcing(a, xx, magnitude);
Vector force = helical_forcing(magnitude, k_force, xx, ff_re,ff_im, phase);
if (is_valid(c)) { return c; } //Scaling N = magnitude*cs*sqrt(k*cs/dt) * dt
const Scalar NN = cs*sqrt(DCONST_REAL(AC_kaver)*cs);
//MV: Like in the Pencil Code. I don't understandf the logic here.
force.x = sqrt(dt)*NN*force.x;
force.y = sqrt(dt)*NN*force.y;
force.z = sqrt(dt)*NN*force.z;
if (is_valid(force)) { return force; }
else { return (Vector){0, 0, 0}; } else { return (Vector){0, 0, 0}; }
} }
#endif // LFORCING #endif // LFORCING
@@ -272,7 +326,7 @@ solve(Scalar dt) {
#if LFORCING #if LFORCING
if (step_number == 2) { if (step_number == 2) {
out_uu = out_uu + dt * forcing(globalVertexIdx); out_uu = out_uu + forcing(globalVertexIdx, dt);
} }
#endif #endif
} }

14
analysis/python/jupyter/notebook_example.ipynb
View File

@@ -24,7 +24,8 @@
"metadata": {}, "metadata": {},
"outputs": [], "outputs": [],
"source": [ "source": [
"meshdir = \"/scratch/data/mvaisala/forcingtest/\"" "meshdir = \"/scratch/data/mvaisala/forcingtest/\"\n",
"#meshdir = \"/scratch/data/mvaisala/normaltest/\""
] ]
}, },
{ {
@@ -34,7 +35,7 @@
"outputs": [], "outputs": [],
"source": [ "source": [
"#imesh = 30000\n", "#imesh = 30000\n",
"imesh = 30000\n", "imesh = 100\n",
"mesh = ad.read.Mesh(imesh, fdir=meshdir)" "mesh = ad.read.Mesh(imesh, fdir=meshdir)"
] ]
}, },
@@ -101,15 +102,8 @@
"metadata": {}, "metadata": {},
"outputs": [], "outputs": [],
"source": [ "source": [
"vis.lineplot.plot_ts(ts, lnrho=1, uutot=1, uux=1, uuy=1, ss=1)" "vis.lineplot.plot_ts(ts, lnrho=1, uutot=1, uux=1, uuy=1, uuz=1, ss=1)"
] ]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
} }
], ],
"metadata": { "metadata": {

4
config/astaroth.conf
View File

@@ -38,6 +38,10 @@ AC_chi = 0.0001
// Forcing // Forcing
AC_relhel = 0.0 AC_relhel = 0.0
AC_forcing_magnitude = 1e-5
AC_kmin = 0.8
AC_kmax = 1.2
// Entropy // Entropy
AC_cp_sound = 1.0 AC_cp_sound = 1.0

27
include/astaroth.h
View File

@@ -136,7 +136,6 @@ extern "C" {
FUNC(AC_cs_sound), \ FUNC(AC_cs_sound), \
FUNC(AC_eta), \ FUNC(AC_eta), \
FUNC(AC_mu0), \ FUNC(AC_mu0), \
FUNC(AC_relhel), \
FUNC(AC_cp_sound), \ FUNC(AC_cp_sound), \
FUNC(AC_gamma), \ FUNC(AC_gamma), \
FUNC(AC_cv_sound), \ FUNC(AC_cv_sound), \
@@ -152,6 +151,23 @@ extern "C" {
FUNC(AC_angl_uu), \ FUNC(AC_angl_uu), \
FUNC(AC_lnrho_edge),\ FUNC(AC_lnrho_edge),\
FUNC(AC_lnrho_out),\ FUNC(AC_lnrho_out),\
/* Forcing parameters. User configured. */\
FUNC(AC_forcing_magnitude),\
FUNC(AC_relhel), \
FUNC(AC_kmin), \
FUNC(AC_kmax), \
/* Forcing parameters. Set by the generator. */\
FUNC(AC_forcing_phase),\
FUNC(AC_k_forcex),\
FUNC(AC_k_forcey),\
FUNC(AC_k_forcez),\
FUNC(AC_kaver),\
FUNC(AC_ff_hel_rex),\
FUNC(AC_ff_hel_rey),\
FUNC(AC_ff_hel_rez),\
FUNC(AC_ff_hel_imx),\
FUNC(AC_ff_hel_imy),\
FUNC(AC_ff_hel_imz),\
/* Additional helper params */\ /* Additional helper params */\
/* (deduced from other params do not set these directly!) */\ /* (deduced from other params do not set these directly!) */\
FUNC(AC_G_CONST),\ FUNC(AC_G_CONST),\
@@ -386,6 +402,15 @@ AcResult acQuit(void);
unless otherwise stated. */ unless otherwise stated. */
AcResult acSynchronize(void); AcResult acSynchronize(void);
/** */
AcResult acLoadDeviceConstant(const AcRealParam param, const AcReal value);
/** Tool for loading forcing vector information into the device memory
*/
AcResult acForcingVec(const AcReal forcing_magnitude, const AcReal3 k_force,
const AcReal3 ff_hel_re, const AcReal3 ff_hel_im, const AcReal forcing_phase,
const AcReal kaver);
/* End extern "C" */ /* End extern "C" */
#ifdef __cplusplus #ifdef __cplusplus
} }

41
src/core/astaroth.cu
View File

@@ -505,3 +505,44 @@ acSynchronize(void)
return AC_SUCCESS; return AC_SUCCESS;
} }
AcResult
acLoadDeviceConstant(const AcRealParam param, const AcReal value)
{
for (int i = 0; i < num_devices; ++i) {
loadDeviceConstant(devices[i], param, value);
}
return AC_SUCCESS;
}
// Tool for loading forcing vector information into the device memory
// %JP: Added a generic function for loading device constants (acLoadDeviceConstant).
// This acForcingVec should go outside the core library since it references user-defined
// parameters such as AC_forcing_magnitude which may not be defined in all projects.
// host_forcing.cc is probably a good place for this.
AcResult
acForcingVec(const AcReal forcing_magnitude, const AcReal3 k_force, const AcReal3 ff_hel_re,
const AcReal3 ff_hel_im, const AcReal forcing_phase, const AcReal kaver)
{
for (int i = 0; i < num_devices; ++i) {
loadDeviceConstant(devices[i], AC_forcing_magnitude, forcing_magnitude);
loadDeviceConstant(devices[i], AC_forcing_phase, forcing_phase);
loadDeviceConstant(devices[i], AC_k_forcex, k_force.x);
loadDeviceConstant(devices[i], AC_k_forcey, k_force.y);
loadDeviceConstant(devices[i], AC_k_forcez, k_force.z);
loadDeviceConstant(devices[i], AC_ff_hel_rex, ff_hel_re.x);
loadDeviceConstant(devices[i], AC_ff_hel_rey, ff_hel_re.y);
loadDeviceConstant(devices[i], AC_ff_hel_rez, ff_hel_re.z);
loadDeviceConstant(devices[i], AC_ff_hel_imx, ff_hel_im.x);
loadDeviceConstant(devices[i], AC_ff_hel_imy, ff_hel_im.y);
loadDeviceConstant(devices[i], AC_ff_hel_imz, ff_hel_im.z);
loadDeviceConstant(devices[i], AC_kaver, kaver);
}
return AC_SUCCESS;
}

53
src/core/kernels/kernels.cuh
View File

@@ -543,39 +543,6 @@ normalized(const AcReal3& vec)
return inv_len * vec; return inv_len * vec;
} }
// Sinusoidal forcing
// https://arxiv.org/pdf/1704.04676.pdf
// NOTE: This method of forcing is depracated. However, it will remain in here
// until a corresponding scheme exists in the new code.
__constant__ AcReal3 forcing_vec;
__constant__ AcReal forcing_phi;
static __device__ __forceinline__ AcReal3
DEPRECATED_forcing(const int i, const int j, const int k)
{
#define DOMAIN_SIZE_X (DCONST_INT(AC_nx) * DCONST_REAL(AC_dsx))
#define DOMAIN_SIZE_Y (DCONST_INT(AC_ny) * DCONST_REAL(AC_dsy))
#define DOMAIN_SIZE_Z (DCONST_INT(AC_nz) * DCONST_REAL(AC_dsz))
const AcReal3 k_vec = (AcReal3){
(i - DCONST_INT(AC_nx_min)) * DCONST_REAL(AC_dsx) - AcReal(.5) * DOMAIN_SIZE_X,
(j - DCONST_INT(AC_ny_min)) * DCONST_REAL(AC_dsy) - AcReal(.5) * DOMAIN_SIZE_Y,
(k - DCONST_INT(AC_nz_min)) * DCONST_REAL(AC_dsz) - AcReal(.5) * DOMAIN_SIZE_Z};
AcReal inv_len = reciprocal_len(k_vec);
if (isnan(inv_len) || isinf(inv_len))
inv_len = 0;
if (inv_len > 2) // hack to make it cool
inv_len = 2;
const AcReal k_dot_x = dot(k_vec, forcing_vec);
const AcReal waves = cos(k_dot_x) * cos(forcing_phi) - sin(k_dot_x) * sin(forcing_phi);
return inv_len * inv_len * waves * forcing_vec;
}
// Note: LNT0 and LNRHO0 must be set very carefully: if the magnitude is different that other values
// in the mesh, then we will inherently lose precision
#define LNT0 (AcReal(0.0))
#define LNRHO0 (AcReal(0.0))
#define H_CONST (AcReal(0.0)) #define H_CONST (AcReal(0.0))
#define C_CONST (AcReal(0.0)) #define C_CONST (AcReal(0.0))
@@ -739,31 +706,11 @@ read_out(const int idx, AcReal* __restrict__ field[], const int3 handle)
* ============================================================================= * =============================================================================
*/ */
static AcReal
randf(void)
{
return AcReal(rand()) / AcReal(RAND_MAX);
}
AcResult AcResult
rk3_step_async(const cudaStream_t stream, const int& step_number, const int3& start, rk3_step_async(const cudaStream_t stream, const int& step_number, const int3& start,
const int3& end, const AcReal dt, VertexBufferArray* buffer) const int3& end, const AcReal dt, VertexBufferArray* buffer)
{ {
const dim3 tpb(32, 1, 4); const dim3 tpb(32, 1, 4);
/////////////////// Forcing
#if LFORCING
const AcReal ff_scale = AcReal(.2);
static AcReal3 ff = ff_scale * (AcReal3){1, 0, 0};
const AcReal radians = randf() * AcReal(2 * M_PI) / 360 / 8;
const AcMatrix rotz = create_rotz(radians);
ff = mul(rotz, ff);
cudaMemcpyToSymbolAsync(forcing_vec, &ff, sizeof(ff), 0, cudaMemcpyHostToDevice, stream);
const AcReal ff_phi = AcReal(M_PI); // AcReal(2 * M_PI) * randf();
cudaMemcpyToSymbolAsync(forcing_phi, &ff_phi, sizeof(ff_phi), 0, cudaMemcpyHostToDevice,
stream);
#endif // LFORCING
//////////////////////////
const int nx = end.x - start.x; const int nx = end.x - start.x;
const int ny = end.y - start.y; const int ny = end.y - start.y;

10
src/standalone/autotest.cc
View File

@@ -1,5 +1,5 @@
/* /*
Copyright (C) 2014-2018, Johannes Pekkilae, Miikka Vaeisalae. Copyright (C) 2014-2019, Johannes Pekkilae, Miikka Vaeisalae.
This file is part of Astaroth. This file is part of Astaroth.
@@ -32,6 +32,7 @@
#include "core/math_utils.h" #include "core/math_utils.h"
#include "model/host_memory.h" #include "model/host_memory.h"
#include "model/host_timestep.h" #include "model/host_timestep.h"
#include "model/host_forcing.h"
#include "model/model_boundconds.h" #include "model/model_boundconds.h"
#include "model/model_reduce.h" #include "model/model_reduce.h"
#include "model/model_rk3.h" #include "model/model_rk3.h"
@@ -699,6 +700,13 @@ run_autotest(void)
VTXBUF_UUX, VTXBUF_UUY, VTXBUF_UUZ); VTXBUF_UUX, VTXBUF_UUY, VTXBUF_UUZ);
const AcReal dt = host_timestep(umax, config); const AcReal dt = host_timestep(umax, config);
#if LFORCING
//CURRENTLY AUTOTEST NOT SUPPORTED WITH FORCING!!!
#endif
// Host integration step // Host integration step
model_rk3(dt, model_mesh); model_rk3(dt, model_mesh);
boundconds(config, model_mesh); boundconds(config, model_mesh);

177
src/standalone/model/host_forcing.cc Normal file
View File

@@ -0,0 +1,177 @@
/*
Copyright (C) 2014-2019, Johannes Pekkilae, Miikka Vaeisalae.
This file is part of Astaroth.
Astaroth is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
Astaroth 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 General Public License
along with Astaroth. If not, see <http://www.gnu.org/licenses/>.
*/
/**
* @file
* \brief Brief info.
*
* Detailed info.
*
*/
#include "host_forcing.h"
#include "core/math_utils.h"
// The is a wrapper for genering random numbers with a chosen system.
AcReal
get_random_number_01()
{
// TODO: Implement better randon number generator http://www.cplusplus.com/reference/random/
return AcReal(rand()) / AcReal(RAND_MAX);
}
AcReal3
cross(const AcReal3& a, const AcReal3& b)
{
AcReal3 c;
c.x = a.y * b.z - a.z * b.y;
c.y = a.z * b.x - a.x * b.z;
c.z = a.x * b.y - a.y * b.x;
return c;
}
AcReal
dot(const AcReal3& a, const AcReal3& b)
{
return a.x * b.x + a.y * b.y + a.z * b.z;
}
AcReal3
vec_norm(const AcReal3& a)
{
AcReal3 c;
AcReal norm = dot(a, a);
c.x = a.x / sqrt(norm);
c.y = a.y / sqrt(norm);
c.z = a.z / sqrt(norm);
return c;
}
AcReal3
vec_multi_scal(const AcReal scal, const AcReal3& a)
{
AcReal3 c;
c.x = a.x * scal;
c.y = a.y * scal;
c.z = a.z * scal;
return c;
}
// Generate forcing wave vector k_force
AcReal3
helical_forcing_k_generator(const AcReal kmax, const AcReal kmin)
{
AcReal phi, theta, kk; // Spherical direction coordinates
AcReal3 k_force; // forcing wave vector
AcReal delta_k = kmax - kmin;
// Generate vector in spherical coordinates
phi = AcReal(2.0) * AcReal(M_PI) * get_random_number_01();
theta = AcReal(M_PI) * get_random_number_01();
kk = delta_k * get_random_number_01() + kmin;
// Cast into Cartesian form
k_force = (AcReal3){AcReal(kk * sin(theta) * cos(phi)), //
AcReal(kk * sin(theta) * sin(phi)), //
AcReal(kk * cos(theta))};
// printf("k_force.x %f, k_force.y %f, k_force.z %f \n", k_force.x, k_force.y, k_force.z);
// Round the numbers. In that way k(x/y/z) will get complete waves.
k_force.x = round(k_force.x);
k_force.y = round(k_force.y);
k_force.z = round(k_force.z);
// printf("After rounding --> k_force.x %f, k_force.y %f, k_force.z %f \n", k_force.x,
// k_force.y, k_force.z);
return k_force;
}
// Generate the unit perpendicular unit vector e required for helical forcing
// Addapted from Pencil code forcing.f90 hel_vec() subroutine.
void
helical_forcing_e_generator(AcReal3* e_force, const AcReal3 k_force)
{
AcReal3 k_cross_e = cross(k_force, *e_force);
k_cross_e = vec_norm(k_cross_e);
AcReal3 k_cross_k_cross_e = cross(k_force, k_cross_e);
k_cross_k_cross_e = vec_norm(k_cross_k_cross_e);
AcReal phi = AcReal(2.0) * AcReal(M_PI) * get_random_number_01();
AcReal3 ee_tmp1 = vec_multi_scal(cos(phi), k_cross_e);
AcReal3 ee_tmp2 = vec_multi_scal(sin(phi), k_cross_k_cross_e);
*e_force = (AcReal3){ee_tmp1.x + ee_tmp2.x, ee_tmp1.y + ee_tmp2.y, ee_tmp1.z + ee_tmp2.z};
}
// PC Manual Eq. 223
void
helical_forcing_special_vector(AcReal3* ff_hel_re, AcReal3* ff_hel_im, const AcReal3 k_force,
const AcReal3 e_force, const AcReal relhel)
{
// k dot e
AcReal3 kdote;
kdote.x = k_force.x * e_force.x;
kdote.y = k_force.y * e_force.y;
kdote.z = k_force.z * e_force.z;
// k cross e
AcReal3 k_cross_e;
k_cross_e.x = k_force.y * e_force.z - k_force.z * e_force.y;
k_cross_e.y = k_force.z * e_force.x - k_force.x * e_force.z;
k_cross_e.z = k_force.x * e_force.y - k_force.y * e_force.x;
// k cross k cross e
AcReal3 k_cross_k_cross_e;
k_cross_k_cross_e.x = k_force.y * k_cross_e.z - k_force.z * k_cross_e.y;
k_cross_k_cross_e.y = k_force.z * k_cross_e.x - k_force.x * k_cross_e.z;
k_cross_k_cross_e.z = k_force.x * k_cross_e.y - k_force.y * k_cross_e.x;
// abs(k)
AcReal kabs = sqrt(k_force.x * k_force.x + k_force.y * k_force.y + k_force.z * k_force.z);
AcReal denominator = sqrt(AcReal(1.0) + relhel * relhel) * kabs *
sqrt(kabs * kabs -
(kdote.x * kdote.x + kdote.y * kdote.y + kdote.z * kdote.z));
// MV: I suspect there is a typo in the Pencil Code manual!
//*ff_hel_re = (AcReal3){-relhel*kabs*k_cross_e.x/denominator,
// -relhel*kabs*k_cross_e.y/denominator,
// -relhel*kabs*k_cross_e.z/denominator};
//*ff_hel_im = (AcReal3){k_cross_k_cross_e.x/denominator,
// k_cross_k_cross_e.y/denominator,
// k_cross_k_cross_e.z/denominator};
// See PC forcing.f90 forcing_hel_both()
*ff_hel_re = (AcReal3){kabs * k_cross_e.x / denominator, kabs * k_cross_e.y,
kabs * k_cross_e.z};
*ff_hel_im = (AcReal3){relhel * k_cross_k_cross_e.x / denominator, relhel * k_cross_k_cross_e.y,
relhel * k_cross_k_cross_e.z};
}

45
src/standalone/model/host_forcing.h Normal file
View File

@@ -0,0 +1,45 @@
/*
Copyright (C) 2014-2019, Johannes Pekkilae, Miikka Vaeisalae.
This file is part of Astaroth.
Astaroth is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
Astaroth 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 General Public License
along with Astaroth. If not, see <http://www.gnu.org/licenses/>.
*/
/**
* @file
* \brief Brief info.
*
* Detailed info.
*
*/
#pragma once
#include "astaroth.h"
AcReal get_random_number_01();
AcReal3 cross(const AcReal3& a, const AcReal3& b);
AcReal dot(const AcReal3& a, const AcReal3& b);
AcReal3 vec_norm(const AcReal3& a);
AcReal3 vec_multi_scal(const AcReal scal, const AcReal3& a);
AcReal3 helical_forcing_k_generator(const AcReal kmax, const AcReal kmin);
void helical_forcing_e_generator(AcReal3* e_force, const AcReal3 k_force);
void helical_forcing_special_vector(AcReal3* ff_hel_re, AcReal3* ff_hel_im, const AcReal3 k_force, const AcReal3 e_force, const AcReal relhel);

3
src/standalone/model/host_timestep.cc
View File

@@ -50,7 +50,8 @@ host_timestep(const AcReal& umax, const AcMeshInfo& mesh_info)
// New, closer to the actual Courant timestep // New, closer to the actual Courant timestep
// See Pencil Code user manual p. 38 (timestep section) // See Pencil Code user manual p. 38 (timestep section)
const long double uu_dt = cdt * dsmin / (fabsl(umax) + sqrtl(cs2_sound + 0.0l)); const long double uu_dt = cdt * dsmin / (fabsl(umax) + sqrtl(cs2_sound + 0.0l));
const long double visc_dt = cdtv * dsmin * dsmin / max(max(nu_visc, eta), max(gamma, chi)) + 1; // TODO NOTE: comment the +1 out to get scientifically accurate results const long double visc_dt = cdtv * dsmin * dsmin / max(max(nu_visc, eta), max(gamma, chi));// + 1; // TODO NOTE: comment the +1 out to get scientifically accurate results
//MV: White the +1? It was messing up my computations!
const long double dt = min(uu_dt, visc_dt); const long double dt = min(uu_dt, visc_dt);
return AcReal(timescale) * AcReal(dt); return AcReal(timescale) * AcReal(dt);

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

@@ -1,5 +1,5 @@
/* /*
Copyright (C) 2014-2018, Johannes Pekkilae, Miikka Vaeisalae. Copyright (C) 2014-2019, Johannes Pekkilae, Miikka Vaeisalae.
This file is part of Astaroth. This file is part of Astaroth.
@@ -553,11 +553,6 @@ normalized(const ModelVector& vec)
return inv_len * vec; return inv_len * vec;
} }
// Note: LNT0 and LNRHO0 must be set very carefully: if the magnitude is different that other values
// in the mesh, then we will inherently lose precision
#define LNT0 (ModelScalar(0.0))
#define LNRHO0 (ModelScalar(0.0))
#define H_CONST (ModelScalar(0.0)) #define H_CONST (ModelScalar(0.0))
#define C_CONST (ModelScalar(0.0)) #define C_CONST (ModelScalar(0.0))
@@ -571,8 +566,9 @@ momentum(const ModelVectorData& uu, const ModelScalarData& lnrho
{ {
#if LENTROPY #if LENTROPY
const ModelMatrix S = stress_tensor(uu); const ModelMatrix S = stress_tensor(uu);
const ModelScalar cs2 = get(AC_cs2_sound) * expl(get(AC_gamma) * value(ss) / get(AC_cp_sound) + const ModelScalar cs2 = get(AC_cs2_sound) *
(get(AC_gamma) - 1) * (value(lnrho) - LNRHO0)); expl(get(AC_gamma) * value(ss) / get(AC_cp_sound) +
(get(AC_gamma) - 1) * (value(lnrho) - get(AC_lnrho0)));
const ModelVector j = (ModelScalar(1.) / get(AC_mu0)) * const ModelVector j = (ModelScalar(1.) / get(AC_mu0)) *
(gradient_of_divergence(aa) - laplace_vec(aa)); // Current density (gradient_of_divergence(aa) - laplace_vec(aa)); // Current density
const ModelVector B = curl(aa); const ModelVector B = curl(aa);
@@ -591,7 +587,8 @@ momentum(const ModelVectorData& uu, const ModelScalarData& lnrho
// !!!!!!!!!!!!!!!!%JP: NOTE TODO IMPORTANT!!!!!!!!!!!!!!!!!!!!!!!! // !!!!!!!!!!!!!!!!%JP: NOTE TODO IMPORTANT!!!!!!!!!!!!!!!!!!!!!!!!
// NOT CHECKED FOR CORRECTNESS: USE AT YOUR OWN RISK // NOT CHECKED FOR CORRECTNESS: USE AT YOUR OWN RISK
const ModelMatrix S = stress_tensor(uu); const ModelMatrix S = stress_tensor(uu);
const ModelScalar cs2 = get(AC_cs2_sound) * expl((get(AC_gamma) - 1) * (value(lnrho) - LNRHO0)); const ModelScalar cs2 = get(AC_cs2_sound) *
expl((get(AC_gamma) - 1) * (value(lnrho) - get(AC_lnrho0)));
const ModelVector mom = -mul(gradients(uu), value(uu)) - cs2 * gradient(lnrho) + const ModelVector mom = -mul(gradients(uu), value(uu)) - cs2 * gradient(lnrho) +
get(AC_nu_visc) * (laplace_vec(uu) + get(AC_nu_visc) * (laplace_vec(uu) +
@@ -623,8 +620,8 @@ induction(const ModelVectorData& uu, const ModelVectorData& aa)
static inline ModelScalar static inline ModelScalar
lnT(const ModelScalarData& ss, const ModelScalarData& lnrho) lnT(const ModelScalarData& ss, const ModelScalarData& lnrho)
{ {
const ModelScalar lnT = LNT0 + get(AC_gamma) * value(ss) / get(AC_cp_sound) + const ModelScalar lnT = get(AC_lnT0) + get(AC_gamma) * value(ss) / get(AC_cp_sound) +
(get(AC_gamma) - ModelScalar(1.)) * (value(lnrho) - LNRHO0); (get(AC_gamma) - ModelScalar(1.)) * (value(lnrho) - get(AC_lnrho0));
return lnT; return lnT;
} }
@@ -678,40 +675,95 @@ entropy(const ModelScalarData& ss, const ModelVectorData& uu, const ModelScalarD
*/ */
} }
static bool static inline bool
is_valid(const ModelScalar a) is_valid(const ModelScalar a)
{ {
return !isnan(a) && !isinf(a); return !isnan(a) && !isinf(a);
} }
static bool static inline bool
is_valid(const ModelVector& a) is_valid(const ModelVector& a)
{ {
return is_valid(a.x) && is_valid(a.y) && is_valid(a.z); return is_valid(a.x) && is_valid(a.y) && is_valid(a.z);
} }
#if 0
//FORCING NOT SUPPORTED FOR AUTOTEST
static inline ModelVector static inline ModelVector
forcing(int3 globalVertexIdx) simple_vortex_forcing(ModelVector a, ModelVector b, ModelScalar magnitude)
{
return magnitude * cross(normalized(b - a), (ModelVector){0, 0, 1}); // Vortex
}
static inline ModelVector
simple_outward_flow_forcing(ModelVector a, ModelVector b, ModelScalar magnitude)
{
return magnitude * (1 / length(b - a)) * normalized(b - a); // Outward flow
}
// 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 ff_im, ModelScalar phi)
{
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 = cos(phi);
ModelScalar sin_phi = sin(phi);
ModelScalar cos_k_dox_x = cos(dot(k_force, xx));
ModelScalar sin_k_dox_x = sin(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;
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,
ff_re.z*real_comp_phase - ff_im.z*imag_comp_phase};
return force;
}
static inline ModelVector
forcing(int3 globalVertexIdx, ModelScalar dt)
{ {
// source (origin)
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_ny) * get(AC_dsy),
get(AC_nz) * get(AC_dsz)}; get(AC_nz) * get(AC_dsz)}; // source (origin)
// sink (current index) ModelVector xx = (ModelVector){(globalVertexIdx.x - get(AC_nx_min)) * get(AC_dsx),
ModelVector b = (ModelVector){(globalVertexIdx.x - get(AC_nx_min)) * get(AC_dsx), (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))}; // sink (current index)
(globalVertexIdx.z - get(AC_nz_min)) * get(AC_dsz)}; const ModelScalar cs2 = get(AC_cs2_sound);
const ModelScalar cs = sqrt(cs2);
ModelScalar magnitude = 0.05; //Placeholders until determined properly
// Vector c = magnitude * (1 / length(b - a)) * normalized(b - a); // Outward flow ModelScalar magnitude = get(AC_forcing_magnitude);
ModelVector c = magnitude * cross(normalized(b - a), (ModelVector){0, 0, 1}); // Vortex ModelScalar phase = get(AC_forcing_phase);
if (is_valid(c)) { ModelVector k_force = (ModelVector){ get(AC_k_forcex), get(AC_k_forcey), get(AC_k_forcez)};
return c; 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)};
else {
return (ModelVector){0, 0, 0};
} //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);
ModelVector force = helical_forcing(magnitude, k_force, xx, ff_re,ff_im, phase);
//Scaling N = magnitude*cs*sqrt(k*cs/dt) * dt
const ModelScalar NN = cs*sqrt(get(AC_kaver)*cs);
//MV: Like in the Pencil Code. I don't understandf the logic here.
force.x = sqrt(dt)*NN*force.x;
force.y = sqrt(dt)*NN*force.y;
force.z = sqrt(dt)*NN*force.z;
if (is_valid(force)) { return force; }
else { return (ModelVector){0, 0, 0}; }
} }
#endif
static void static void
solve_alpha_step(const int step_number, const ModelScalar dt, const int i, const int j, const int k, solve_alpha_step(const int step_number, const ModelScalar dt, const int i, const int j, const int k,
@@ -774,12 +826,13 @@ solve_beta_step(const int step_number, const ModelScalar dt, const int i, const
for (int w = 0; w < NUM_VTXBUF_HANDLES; ++w) for (int w = 0; w < NUM_VTXBUF_HANDLES; ++w)
out->vertex_buffer[w][idx] += beta[step_number] * in.vertex_buffer[w][idx]; out->vertex_buffer[w][idx] += beta[step_number] * in.vertex_buffer[w][idx];
(void)dt; // Suppress unused variable warning if forcing not used
#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}, dt);
out->vertex_buffer[VTXBUF_UUX][idx] += force.x * dt; out->vertex_buffer[VTXBUF_UUX][idx] += force.x;
out->vertex_buffer[VTXBUF_UUY][idx] += force.y * dt; out->vertex_buffer[VTXBUF_UUY][idx] += force.y;
out->vertex_buffer[VTXBUF_UUZ][idx] += force.z * dt; out->vertex_buffer[VTXBUF_UUZ][idx] += force.z;
} }
#endif #endif
} }

38
src/standalone/simulation.cc
View File

@@ -29,6 +29,7 @@
#include "config_loader.h" #include "config_loader.h"
#include "core/errchk.h" #include "core/errchk.h"
#include "core/math_utils.h" #include "core/math_utils.h"
#include "model/host_forcing.h"
#include "model/host_memory.h" #include "model/host_memory.h"
#include "model/host_timestep.h" #include "model/host_timestep.h"
#include "model/model_reduce.h" #include "model/model_reduce.h"
@@ -236,10 +237,47 @@ run_simulation(void)
AcReal bin_save_t = mesh_info.real_params[AC_bin_save_t]; AcReal bin_save_t = mesh_info.real_params[AC_bin_save_t];
AcReal bin_crit_t = bin_save_t; AcReal bin_crit_t = bin_save_t;
// Forcing properties
// AcReal relhel = mesh_info.real_params[AC_relhel]; // JP: gcc warning: unused
// AcReal magnitude = mesh_info.real_params[AC_forcing_magnitude]; // JP: gcc warning: unused
// AcReal kmin = mesh_info.real_params[AC_kmin]; // JP: gcc warning: unused
// AcReal kmax = mesh_info.real_params[AC_kmax]; // JP: gcc warning: unused
// AcReal kaver = (kmax - kmin) / AcReal(2.0); // JP: gcc warning: unused
/* initialize random seed: */
srand(312256655);
/* Step the simulation */ /* Step the simulation */
for (int i = 1; i < max_steps; ++i) { for (int i = 1; i < max_steps; ++i) {
const AcReal umax = acReduceVec(RTYPE_MAX, VTXBUF_UUX, VTXBUF_UUY, VTXBUF_UUZ); const AcReal umax = acReduceVec(RTYPE_MAX, VTXBUF_UUX, VTXBUF_UUY, VTXBUF_UUZ);
const AcReal dt = host_timestep(umax, mesh_info); const AcReal dt = host_timestep(umax, mesh_info);
#if LFORCING
// Generate a forcing vector before canculating an integration step.
// Generate forcing wave vector k_force
AcReal3 k_force;
k_force = helical_forcing_k_generator(kmax, kmin);
// Randomize the phase
AcReal phase = AcReal(2.0) * AcReal(M_PI) * get_random_number_01();
// Generate e for k. Needed for the sake of isotrophy.
AcReal3 e_force;
if ((k_force.y == 0.0) && (k_force.z == 0.0)) {
e_force = (AcReal3){0.0, 1.0, 0.0};
}
else {
e_force = (AcReal3){1.0, 0.0, 0.0};
}
helical_forcing_e_generator(&e_force, k_force);
AcReal3 ff_hel_re, ff_hel_im;
helical_forcing_special_vector(&ff_hel_re, &ff_hel_im, k_force, e_force, relhel);
acForcingVec(magnitude, k_force, ff_hel_re, ff_hel_im, phase, kaver);
#endif
acIntegrate(dt); acIntegrate(dt);
t_step += dt; t_step += dt;