diff --git a/acc/compile.sh b/acc/compile.sh
index 7f7c143..0fe6b4c 100755
--- a/acc/compile.sh
+++ b/acc/compile.sh
@@ -8,17 +8,21 @@ FILENAME="${FULL_NAME%.*}"
EXTENSION="${FULL_NAME##*.}"
if [ "${EXTENSION}" = "sas" ]; then
- echo "Generating stencil assembly stage ${FILENAME}.sas -> stencil_assembly.cuh"
COMPILE_FLAGS="-sas" # Generate stencil assembly stage
CUH_FILENAME="stencil_assembly.cuh"
+ echo "Generating stencil assembly stage ${FILENAME}.sas -> ${CUH_FILENAME}"
elif [ "${EXTENSION}" = "sps" ]; then
- echo "Generating stencil processing stage: ${FILENAME}.sps -> stencil_process.cuh"
COMPILE_FLAGS="-sps" # Generate stencil processing stage
CUH_FILENAME="stencil_process.cuh"
+ echo "Generating stencil processing stage: ${FILENAME}.sps -> ${CUH_FILENAME}"
+elif [ "${EXTENSION}" = "sdh" ]; then
+ COMPILE_FLAGS="-sdh" # Generate stencil definition header
+ CUH_FILENAME="stencil_defines.h"
+ echo "Generating stencil definition header: ${FILENAME}.sdh -> ${CUH_FILENAME}"
else
echo "Error: unknown extension" ${EXTENSION} "of file" ${FULL_NAME}
- echo "Extension should be either .sas or .sps"
+ echo "Extension should be either .sas, .sps or .sdh"
exit
fi
-${ACC_DIR}/preprocess.sh $2 $1 | ${ACC_DIR}/build/acc ${COMPILE_FLAGS} > ${CUH_FILENAME}
+${ACC_DIR}/preprocess.sh $1 | ${ACC_DIR}/build/acc ${COMPILE_FLAGS} > ${CUH_FILENAME}
diff --git a/acc/mhd_solver/.gitignore b/acc/mhd_solver/.gitignore
new file mode 100644
index 0000000..bc4b7d8
--- /dev/null
+++ b/acc/mhd_solver/.gitignore
@@ -0,0 +1,5 @@
+build
+testbin
+
+# Except this file
+!.gitignore
diff --git a/acc/mhd_solver/stencil_assembly.sas b/acc/mhd_solver/stencil_assembly.sas
index 5c4f14a..85886e6 100644
--- a/acc/mhd_solver/stencil_assembly.sas
+++ b/acc/mhd_solver/stencil_assembly.sas
@@ -1,3 +1,5 @@
+#include "stencil_definition.sdh"
+
Preprocessed Scalar
value(in ScalarField vertex)
{
@@ -7,9 +9,7 @@ value(in ScalarField vertex)
Preprocessed Vector
gradient(in ScalarField vertex)
{
- return (Vector){derx(vertexIdx, vertex),
- dery(vertexIdx, vertex),
- derz(vertexIdx, vertex)};
+ return (Vector){derx(vertexIdx, vertex), dery(vertexIdx, vertex), derz(vertexIdx, vertex)};
}
#if LUPWD
@@ -17,46 +17,46 @@ gradient(in ScalarField vertex)
Preprocessed Scalar
der6x_upwd(in ScalarField vertex)
{
- Scalar inv_ds = DCONST_REAL(AC_inv_dsx);
+ Scalar inv_ds = AC_inv_dsx;
- return (Scalar){ Scalar(1.0/60.0)*inv_ds* (
- - Scalar(20.0)* vertex[vertexIdx.x, vertexIdx.y, vertexIdx.z]
- + Scalar(15.0)*(vertex[vertexIdx.x+1, vertexIdx.y, vertexIdx.z]
- + vertex[vertexIdx.x-1, vertexIdx.y, vertexIdx.z])
- - Scalar( 6.0)*(vertex[vertexIdx.x+2, vertexIdx.y, vertexIdx.z]
- + vertex[vertexIdx.x-2, vertexIdx.y, vertexIdx.z])
- + vertex[vertexIdx.x+3, vertexIdx.y, vertexIdx.z]
- + vertex[vertexIdx.x-3, vertexIdx.y, vertexIdx.z])};
+ return (Scalar){Scalar(1.0 / 60.0) * inv_ds *
+ (-Scalar(20.0) * vertex[vertexIdx.x, vertexIdx.y, vertexIdx.z] +
+ Scalar(15.0) * (vertex[vertexIdx.x + 1, vertexIdx.y, vertexIdx.z] +
+ vertex[vertexIdx.x - 1, vertexIdx.y, vertexIdx.z]) -
+ Scalar(6.0) * (vertex[vertexIdx.x + 2, vertexIdx.y, vertexIdx.z] +
+ vertex[vertexIdx.x - 2, vertexIdx.y, vertexIdx.z]) +
+ vertex[vertexIdx.x + 3, vertexIdx.y, vertexIdx.z] +
+ vertex[vertexIdx.x - 3, vertexIdx.y, vertexIdx.z])};
}
Preprocessed Scalar
der6y_upwd(in ScalarField vertex)
{
- Scalar inv_ds = DCONST_REAL(AC_inv_dsy);
+ Scalar inv_ds = AC_inv_dsy;
- return (Scalar){ Scalar(1.0/60.0)*inv_ds* (
- -Scalar( 20.0)* vertex[vertexIdx.x, vertexIdx.y, vertexIdx.z]
- +Scalar( 15.0)*(vertex[vertexIdx.x, vertexIdx.y+1, vertexIdx.z]
- + vertex[vertexIdx.x, vertexIdx.y-1, vertexIdx.z])
- -Scalar( 6.0)*(vertex[vertexIdx.x, vertexIdx.y+2, vertexIdx.z]
- + vertex[vertexIdx.x, vertexIdx.y-2, vertexIdx.z])
- + vertex[vertexIdx.x, vertexIdx.y+3, vertexIdx.z]
- + vertex[vertexIdx.x, vertexIdx.y-3, vertexIdx.z])};
+ return (Scalar){Scalar(1.0 / 60.0) * inv_ds *
+ (-Scalar(20.0) * vertex[vertexIdx.x, vertexIdx.y, vertexIdx.z] +
+ Scalar(15.0) * (vertex[vertexIdx.x, vertexIdx.y + 1, vertexIdx.z] +
+ vertex[vertexIdx.x, vertexIdx.y - 1, vertexIdx.z]) -
+ Scalar(6.0) * (vertex[vertexIdx.x, vertexIdx.y + 2, vertexIdx.z] +
+ vertex[vertexIdx.x, vertexIdx.y - 2, vertexIdx.z]) +
+ vertex[vertexIdx.x, vertexIdx.y + 3, vertexIdx.z] +
+ vertex[vertexIdx.x, vertexIdx.y - 3, vertexIdx.z])};
}
Preprocessed Scalar
der6z_upwd(in ScalarField vertex)
{
- Scalar inv_ds = DCONST_REAL(AC_inv_dsz);
+ Scalar inv_ds = AC_inv_dsz;
- return (Scalar){ Scalar(1.0/60.0)*inv_ds* (
- -Scalar( 20.0)* vertex[vertexIdx.x, vertexIdx.y, vertexIdx.z]
- +Scalar( 15.0)*(vertex[vertexIdx.x, vertexIdx.y, vertexIdx.z+1]
- + vertex[vertexIdx.x, vertexIdx.y, vertexIdx.z-1])
- -Scalar( 6.0)*(vertex[vertexIdx.x, vertexIdx.y, vertexIdx.z+2]
- + vertex[vertexIdx.x, vertexIdx.y, vertexIdx.z-2])
- + vertex[vertexIdx.x, vertexIdx.y, vertexIdx.z+3]
- + vertex[vertexIdx.x, vertexIdx.y, vertexIdx.z-3])};
+ return (Scalar){Scalar(1.0 / 60.0) * inv_ds *
+ (-Scalar(20.0) * vertex[vertexIdx.x, vertexIdx.y, vertexIdx.z] +
+ Scalar(15.0) * (vertex[vertexIdx.x, vertexIdx.y, vertexIdx.z + 1] +
+ vertex[vertexIdx.x, vertexIdx.y, vertexIdx.z - 1]) -
+ Scalar(6.0) * (vertex[vertexIdx.x, vertexIdx.y, vertexIdx.z + 2] +
+ vertex[vertexIdx.x, vertexIdx.y, vertexIdx.z - 2]) +
+ vertex[vertexIdx.x, vertexIdx.y, vertexIdx.z + 3] +
+ vertex[vertexIdx.x, vertexIdx.y, vertexIdx.z - 3])};
}
#endif
@@ -66,9 +66,10 @@ hessian(in ScalarField vertex)
{
Matrix hessian;
- hessian.row[0] = (Vector){derxx(vertexIdx, vertex), derxy(vertexIdx, vertex), derxz(vertexIdx, vertex)};
- hessian.row[1] = (Vector){hessian.row[0].y, deryy(vertexIdx, vertex), deryz(vertexIdx, vertex)};
- hessian.row[2] = (Vector){hessian.row[0].z, hessian.row[1].z, derzz(vertexIdx, vertex)};
+ hessian.row[0] = (Vector){derxx(vertexIdx, vertex), derxy(vertexIdx, vertex),
+ derxz(vertexIdx, vertex)};
+ hessian.row[1] = (Vector){hessian.row[0].y, deryy(vertexIdx, vertex), deryz(vertexIdx, vertex)};
+ hessian.row[2] = (Vector){hessian.row[0].z, hessian.row[1].z, derzz(vertexIdx, vertex)};
return hessian;
}
diff --git a/acc/mhd_solver/stencil_defines.h b/acc/mhd_solver/stencil_defines.h
deleted file mode 100644
index b4bc622..0000000
--- a/acc/mhd_solver/stencil_defines.h
+++ /dev/null
@@ -1,163 +0,0 @@
-/*
- 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 .
-*/
-#pragma once
-
-/*
- * =============================================================================
- * Logical switches
- * =============================================================================
- */
-#define STENCIL_ORDER (6)
-#define NGHOST (STENCIL_ORDER / 2)
-#define LDENSITY (1)
-#define LHYDRO (1)
-#define LMAGNETIC (1)
-#define LENTROPY (1)
-#define LTEMPERATURE (0)
-#define LFORCING (1)
-#define LUPWD (1)
-
-#define AC_THERMAL_CONDUCTIVITY (AcReal(0.001)) // TODO: make an actual config parameter
-
-/*
- * =============================================================================
- * User-defined parameters
- * =============================================================================
- */
-// clang-format off
-#define AC_FOR_USER_INT_PARAM_TYPES(FUNC)\
- /* Other */\
- FUNC(AC_max_steps), \
- FUNC(AC_save_steps), \
- FUNC(AC_bin_steps), \
- FUNC(AC_bc_type),
-
-#define AC_FOR_USER_INT3_PARAM_TYPES(FUNC)
-
-#define AC_FOR_USER_REAL_PARAM_TYPES(FUNC)\
- /* cparams */\
- FUNC(AC_dsx), \
- FUNC(AC_dsy), \
- FUNC(AC_dsz), \
- FUNC(AC_dsmin), \
- /* physical grid*/\
- FUNC(AC_xlen), \
- FUNC(AC_ylen), \
- FUNC(AC_zlen), \
- FUNC(AC_xorig), \
- FUNC(AC_yorig), \
- FUNC(AC_zorig), \
- /*Physical units*/\
- FUNC(AC_unit_density),\
- FUNC(AC_unit_velocity),\
- FUNC(AC_unit_length),\
- /* properties of gravitating star*/\
- FUNC(AC_star_pos_x),\
- FUNC(AC_star_pos_y),\
- FUNC(AC_star_pos_z),\
- FUNC(AC_M_star),\
- /* Run params */\
- FUNC(AC_cdt), \
- FUNC(AC_cdtv), \
- FUNC(AC_cdts), \
- FUNC(AC_nu_visc), \
- FUNC(AC_cs_sound), \
- FUNC(AC_eta), \
- FUNC(AC_mu0), \
- FUNC(AC_cp_sound), \
- FUNC(AC_gamma), \
- FUNC(AC_cv_sound), \
- FUNC(AC_lnT0), \
- FUNC(AC_lnrho0), \
- FUNC(AC_zeta), \
- FUNC(AC_trans),\
- /* Other */\
- FUNC(AC_bin_save_t), \
- /* Initial condition params */\
- FUNC(AC_ampl_lnrho), \
- FUNC(AC_ampl_uu), \
- FUNC(AC_angl_uu), \
- FUNC(AC_lnrho_edge),\
- 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 */\
- /* (deduced from other params do not set these directly!) */\
- FUNC(AC_G_CONST),\
- FUNC(AC_GM_star),\
- FUNC(AC_sq2GM_star),\
- FUNC(AC_cs2_sound), \
- FUNC(AC_inv_dsx), \
- FUNC(AC_inv_dsy), \
- FUNC(AC_inv_dsz),
-
-#define AC_FOR_USER_REAL3_PARAM_TYPES(FUNC)
-// clang-format on
-
-/*
- * =============================================================================
- * User-defined vertex buffers
- * =============================================================================
- */
-// clang-format off
-#if LENTROPY
-#define AC_FOR_VTXBUF_HANDLES(FUNC) \
- FUNC(VTXBUF_LNRHO), \
- FUNC(VTXBUF_UUX), \
- FUNC(VTXBUF_UUY), \
- FUNC(VTXBUF_UUZ), \
- FUNC(VTXBUF_AX), \
- FUNC(VTXBUF_AY), \
- FUNC(VTXBUF_AZ), \
- FUNC(VTXBUF_ENTROPY),
-#elif LMAGNETIC
-#define AC_FOR_VTXBUF_HANDLES(FUNC) \
- FUNC(VTXBUF_LNRHO), \
- FUNC(VTXBUF_UUX), \
- FUNC(VTXBUF_UUY), \
- FUNC(VTXBUF_UUZ), \
- FUNC(VTXBUF_AX), \
- FUNC(VTXBUF_AY), \
- FUNC(VTXBUF_AZ),
-#elif LHYDRO
-#define AC_FOR_VTXBUF_HANDLES(FUNC) \
- FUNC(VTXBUF_LNRHO), \
- FUNC(VTXBUF_UUX), \
- FUNC(VTXBUF_UUY), \
- FUNC(VTXBUF_UUZ),
-#else
-#define AC_FOR_VTXBUF_HANDLES(FUNC) \
- FUNC(VTXBUF_LNRHO),
-#endif
-// clang-format on
diff --git a/acc/mhd_solver/stencil_definition.sdh b/acc/mhd_solver/stencil_definition.sdh
new file mode 100644
index 0000000..31667fc
--- /dev/null
+++ b/acc/mhd_solver/stencil_definition.sdh
@@ -0,0 +1,117 @@
+#define LDENSITY (1)
+#define LHYDRO (1)
+#define LMAGNETIC (1)
+#define LENTROPY (1)
+#define LTEMPERATURE (0)
+#define LFORCING (1)
+#define LUPWD (1)
+
+#define AC_THERMAL_CONDUCTIVITY (AcReal(0.001)) // TODO: make an actual config parameter
+
+// Int params
+uniform int AC_max_steps;
+uniform int AC_save_steps;
+uniform int AC_bin_steps;
+uniform int AC_bc_type;
+
+// Real params
+// Spacing
+uniform Scalar AC_dsx;
+uniform Scalar AC_dsy;
+uniform Scalar AC_dsz;
+uniform Scalar AC_dsmin;
+// physical grid
+uniform Scalar AC_xlen;
+uniform Scalar AC_ylen;
+uniform Scalar AC_zlen;
+uniform Scalar AC_xorig;
+uniform Scalar AC_yorig;
+uniform Scalar AC_zorig;
+// Physical units
+uniform Scalar AC_unit_density;
+uniform Scalar AC_unit_velocity;
+uniform Scalar AC_unit_length;
+// properties of gravitating star
+uniform Scalar AC_star_pos_x;
+uniform Scalar AC_star_pos_y;
+uniform Scalar AC_star_pos_z;
+uniform Scalar AC_M_star;
+// Run params
+uniform Scalar AC_cdt;
+uniform Scalar AC_cdtv;
+uniform Scalar AC_cdts;
+uniform Scalar AC_nu_visc;
+uniform Scalar AC_cs_sound;
+uniform Scalar AC_eta;
+uniform Scalar AC_mu0;
+uniform Scalar AC_cp_sound;
+uniform Scalar AC_gamma;
+uniform Scalar AC_cv_sound;
+uniform Scalar AC_lnT0;
+uniform Scalar AC_lnrho0;
+uniform Scalar AC_zeta;
+uniform Scalar AC_trans;
+// Other
+uniform Scalar AC_bin_save_t;
+// Initial condition params
+uniform Scalar AC_ampl_lnrho;
+uniform Scalar AC_ampl_uu;
+uniform Scalar AC_angl_uu;
+uniform Scalar AC_lnrho_edge;
+uniform Scalar AC_lnrho_out;
+// Forcing parameters. User configured.
+uniform Scalar AC_forcing_magnitude;
+uniform Scalar AC_relhel;
+uniform Scalar AC_kmin;
+uniform Scalar AC_kmax;
+// Forcing parameters. Set by the generator.
+uniform Scalar AC_forcing_phase;
+uniform Scalar AC_k_forcex;
+uniform Scalar AC_k_forcey;
+uniform Scalar AC_k_forcez;
+uniform Scalar AC_kaver;
+uniform Scalar AC_ff_hel_rex;
+uniform Scalar AC_ff_hel_rey;
+uniform Scalar AC_ff_hel_rez;
+uniform Scalar AC_ff_hel_imx;
+uniform Scalar AC_ff_hel_imy;
+uniform Scalar AC_ff_hel_imz;
+// Additional helper params // (deduced from other params do not set these directly!)
+uniform Scalar AC_G_CONST;
+uniform Scalar AC_GM_star;
+uniform Scalar AC_sq2GM_star;
+uniform Scalar AC_cs2_sound;
+uniform Scalar AC_inv_dsx;
+uniform Scalar AC_inv_dsy;
+uniform Scalar AC_inv_dsz;
+
+/*
+ * =============================================================================
+ * User-defined vertex buffers
+ * =============================================================================
+ */
+#if LENTROPY
+uniform ScalarField VTXBUF_LNRHO;
+uniform ScalarField VTXBUF_UUX;
+uniform ScalarField VTXBUF_UUY;
+uniform ScalarField VTXBUF_UUZ;
+uniform ScalarField VTXBUF_AX;
+uniform ScalarField VTXBUF_AY;
+uniform ScalarField VTXBUF_AZ;
+uniform ScalarField VTXBUF_ENTROPY;
+#elif LMAGNETIC
+uniform ScalarField VTXBUF_LNRHO;
+uniform ScalarField VTXBUF_UUX;
+uniform ScalarField VTXBUF_UUY;
+uniform ScalarField VTXBUF_UUZ;
+uniform ScalarField VTXBUF_AX;
+uniform ScalarField VTXBUF_AY;
+uniform ScalarField VTXBUF_AZ;
+#elif LHYDRO
+uniform ScalarField VTXBUF_LNRHO;
+uniform ScalarField VTXBUF_UUX;
+uniform ScalarField VTXBUF_UUY;
+uniform ScalarField VTXBUF_UUZ;
+#else
+uniform ScalarField VTXBUF_LNRHO;
+#endif
diff --git a/acc/mhd_solver/stencil_process.sps b/acc/mhd_solver/stencil_process.sps
index ee8656e..e605b9f 100644
--- a/acc/mhd_solver/stencil_process.sps
+++ b/acc/mhd_solver/stencil_process.sps
@@ -1,28 +1,4 @@
-// Declare uniforms (i.e. device constants)
-uniform Scalar cs2_sound;
-uniform Scalar nu_visc;
-uniform Scalar cp_sound;
-uniform Scalar cv_sound;
-uniform Scalar mu0;
-uniform Scalar eta;
-uniform Scalar gamma;
-uniform Scalar zeta;
-
-uniform Scalar dsx;
-uniform Scalar dsy;
-uniform Scalar dsz;
-
-uniform Scalar lnT0;
-uniform Scalar lnrho0;
-
-uniform int nx_min;
-uniform int ny_min;
-uniform int nz_min;
-uniform int nx;
-uniform int ny;
-uniform int nz;
-
-
+#include "stencil_definition.sdh"
Vector
value(in VectorField uu)
@@ -37,7 +13,7 @@ upwd_der6(in VectorField uu, in ScalarField lnrho)
Scalar uux = fabs(value(uu).x);
Scalar uuy = fabs(value(uu).y);
Scalar uuz = fabs(value(uu).z);
- return (Scalar){uux*der6x_upwd(lnrho) + uuy*der6y_upwd(lnrho) + uuz*der6z_upwd(lnrho)};
+ return (Scalar){uux * der6x_upwd(lnrho) + uuy * der6y_upwd(lnrho) + uuz * der6z_upwd(lnrho)};
}
#endif
@@ -48,10 +24,11 @@ gradients(in VectorField uu)
}
Scalar
-continuity(in VectorField uu, in ScalarField lnrho) {
+continuity(in VectorField uu, in ScalarField lnrho)
+{
return -dot(value(uu), gradient(lnrho))
#if LUPWD
- //This is a corrective hyperdiffusion term for upwinding.
+ // This is a corrective hyperdiffusion term for upwinding.
+ upwd_der6(uu, lnrho)
#endif
- divergence(uu);
@@ -59,133 +36,136 @@ continuity(in VectorField uu, in ScalarField lnrho) {
#if LENTROPY
Vector
-momentum(in VectorField uu, in ScalarField lnrho, in ScalarField ss, in VectorField aa) {
- const Matrix S = stress_tensor(uu);
- 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
+momentum(in VectorField uu, in ScalarField lnrho, in ScalarField ss, in VectorField aa)
+{
+ const Matrix S = stress_tensor(uu);
+ const Scalar cs2 = AC_cs2_sound * exp(AC_gamma * value(ss) / AC_cp_sound +
+ (AC_gamma - 1) * (value(lnrho) - AC_lnrho0));
+ const Vector j = (Scalar(1.) / AC_mu0) *
+ (gradient_of_divergence(aa) - laplace_vec(aa)); // Current density
const Vector B = curl(aa);
- //TODO: DOES INTHERMAL VERSTION INCLUDE THE MAGNETIC FIELD?
+ // TODO: DOES INTHERMAL VERSTION INCLUDE THE MAGNETIC FIELD?
const Scalar inv_rho = Scalar(1.) / exp(value(lnrho));
// Regex replace CPU constants with get\(AC_([a-zA-Z_0-9]*)\)
// \1
- const Vector mom = - mul(gradients(uu), value(uu))
- - cs2 * ((Scalar(1.) / cp_sound) * gradient(ss) + gradient(lnrho))
- + inv_rho * cross(j, B)
- + nu_visc * (
- laplace_vec(uu)
- + Scalar(1. / 3.) * gradient_of_divergence(uu)
- + Scalar(2.) * mul(S, gradient(lnrho))
- )
- + zeta * gradient_of_divergence(uu);
+ const Vector mom = -mul(gradients(uu), value(uu)) -
+ cs2 * ((Scalar(1.) / AC_cp_sound) * gradient(ss) + gradient(lnrho)) +
+ inv_rho * cross(j, B) +
+ AC_nu_visc *
+ (laplace_vec(uu) + Scalar(1. / 3.) * gradient_of_divergence(uu) +
+ Scalar(2.) * mul(S, gradient(lnrho))) +
+ AC_zeta * gradient_of_divergence(uu);
return mom;
}
#elif LTEMPERATURE
Vector
-momentum(in VectorField uu, in ScalarField lnrho, in ScalarField tt) {
- Vector mom;
+momentum(in VectorField uu, in ScalarField lnrho, in ScalarField tt)
+{
+ Vector mom;
- const Matrix S = stress_tensor(uu);
+ const Matrix S = stress_tensor(uu);
- const Vector pressure_term = (cp_sound - cv_sound) * (gradient(tt) + value(tt) * gradient(lnrho));
+ const Vector pressure_term = (AC_cp_sound - AC_cv_sound) *
+ (gradient(tt) + value(tt) * gradient(lnrho));
- mom = -mul(gradients(uu), value(uu)) -
- pressure_term +
- nu_visc *
- (laplace_vec(uu) + Scalar(1. / 3.) * gradient_of_divergence(uu) +
- Scalar(2.) * mul(S, gradient(lnrho))) + zeta * gradient_of_divergence(uu);
+ mom = -mul(gradients(uu), value(uu)) - pressure_term +
+ AC_nu_visc * (laplace_vec(uu) + Scalar(1. / 3.) * gradient_of_divergence(uu) +
+ Scalar(2.) * mul(S, gradient(lnrho))) +
+ AC_zeta * gradient_of_divergence(uu);
- #if LGRAVITY
- mom = mom - (Vector){0, 0, -10.0};
- #endif
+#if LGRAVITY
+ mom = mom - (Vector){0, 0, -10.0};
+#endif
- return mom;
+ return mom;
}
#else
Vector
-momentum(in VectorField uu, in ScalarField lnrho) {
- Vector mom;
+momentum(in VectorField uu, in ScalarField lnrho)
+{
+ Vector mom;
- const Matrix S = stress_tensor(uu);
+ const Matrix S = stress_tensor(uu);
// Isothermal: we have constant speed of sound
- mom = -mul(gradients(uu), value(uu)) -
- cs2_sound * gradient(lnrho) +
- nu_visc *
- (laplace_vec(uu) + Scalar(1. / 3.) * gradient_of_divergence(uu) +
- Scalar(2.) * mul(S, gradient(lnrho))) + zeta * gradient_of_divergence(uu);
+ mom = -mul(gradients(uu), value(uu)) - AC_cs2_sound * gradient(lnrho) +
+ AC_nu_visc * (laplace_vec(uu) + Scalar(1. / 3.) * gradient_of_divergence(uu) +
+ Scalar(2.) * mul(S, gradient(lnrho))) +
+ AC_zeta * gradient_of_divergence(uu);
- #if LGRAVITY
- mom = mom - (Vector){0, 0, -10.0};
- #endif
+#if LGRAVITY
+ mom = mom - (Vector){0, 0, -10.0};
+#endif
- return mom;
+ return mom;
}
#endif
-
Vector
-induction(in VectorField uu, in VectorField aa) {
- // Note: We do (-nabla^2 A + nabla(nabla dot A)) instead of (nabla x (nabla
- // x A)) in order to avoid taking the first derivative twice (did the math,
- // yes this actually works. See pg.28 in arXiv:astro-ph/0109497)
- // u cross B - ETA * mu0 * (mu0^-1 * [- laplace A + grad div A ])
- const Vector B = curl(aa);
- const Vector grad_div = gradient_of_divergence(aa);
- const Vector lap = laplace_vec(aa);
+induction(in VectorField uu, in VectorField aa)
+{
+ // Note: We do (-nabla^2 A + nabla(nabla dot A)) instead of (nabla x (nabla
+ // x A)) in order to avoid taking the first derivative twice (did the math,
+ // yes this actually works. See pg.28 in arXiv:astro-ph/0109497)
+ // u cross B - AC_eta * AC_mu0 * (AC_mu0^-1 * [- laplace A + grad div A ])
+ const Vector B = curl(aa);
+ const Vector grad_div = gradient_of_divergence(aa);
+ const Vector lap = laplace_vec(aa);
- // Note, mu0 is cancelled out
- const Vector ind = cross(value(uu), B) - eta * (grad_div - lap);
+ // Note, AC_mu0 is cancelled out
+ const Vector ind = cross(value(uu), B) - AC_eta * (grad_div - lap);
- return ind;
+ return ind;
}
-
#if LENTROPY
Scalar
-lnT( in ScalarField ss, in ScalarField lnrho) {
- const Scalar lnT = lnT0 + gamma * value(ss) / cp_sound +
- (gamma - Scalar(1.)) * (value(lnrho) - lnrho0);
- return lnT;
+lnT(in ScalarField ss, in ScalarField lnrho)
+{
+ const Scalar lnT = AC_lnT0 + AC_gamma * value(ss) / AC_cp_sound +
+ (AC_gamma - Scalar(1.)) * (value(lnrho) - AC_lnrho0);
+ return lnT;
}
// Nabla dot (K nabla T) / (rho T)
Scalar
-heat_conduction( in ScalarField ss, in ScalarField lnrho) {
- const Scalar inv_cp_sound = AcReal(1.) / cp_sound;
+heat_conduction(in ScalarField ss, in ScalarField lnrho)
+{
+ const Scalar inv_AC_cp_sound = AcReal(1.) / AC_cp_sound;
- const Vector grad_ln_chi = - gradient(lnrho);
+ const Vector grad_ln_chi = -gradient(lnrho);
- const Scalar first_term = gamma * inv_cp_sound * laplace(ss) +
- (gamma - AcReal(1.)) * laplace(lnrho);
- const Vector second_term = gamma * inv_cp_sound * gradient(ss) +
- (gamma - AcReal(1.)) * gradient(lnrho);
- const Vector third_term = gamma * (inv_cp_sound * gradient(ss) +
- gradient(lnrho)) + grad_ln_chi;
+ const Scalar first_term = AC_gamma * inv_AC_cp_sound * laplace(ss) +
+ (AC_gamma - AcReal(1.)) * laplace(lnrho);
+ const Vector second_term = AC_gamma * inv_AC_cp_sound * gradient(ss) +
+ (AC_gamma - AcReal(1.)) * gradient(lnrho);
+ const Vector third_term = AC_gamma * (inv_AC_cp_sound * gradient(ss) + gradient(lnrho)) +
+ grad_ln_chi;
- const Scalar chi = AC_THERMAL_CONDUCTIVITY / (exp(value(lnrho)) * cp_sound);
- return cp_sound * chi * (first_term + dot(second_term, third_term));
+ const Scalar chi = AC_THERMAL_CONDUCTIVITY / (exp(value(lnrho)) * AC_cp_sound);
+ return AC_cp_sound * chi * (first_term + dot(second_term, third_term));
}
Scalar
-heating(const int i, const int j, const int k) {
- return 1;
+heating(const int i, const int j, const int k)
+{
+ return 1;
}
Scalar
-entropy(in ScalarField ss, in VectorField uu, in ScalarField lnrho, in VectorField aa) {
- const Matrix S = stress_tensor(uu);
+entropy(in ScalarField ss, in VectorField uu, in ScalarField lnrho, in VectorField aa)
+{
+ const Matrix S = stress_tensor(uu);
const Scalar inv_pT = Scalar(1.) / (exp(value(lnrho)) * exp(lnT(ss, lnrho)));
- const Vector j = (Scalar(1.) / mu0) * (gradient_of_divergence(aa) - laplace_vec(aa)); // Current density
- const Scalar RHS = H_CONST - C_CONST
- + eta * (mu0) * dot(j, j)
- + Scalar(2.) * exp(value(lnrho)) * nu_visc * contract(S)
- + zeta * exp(value(lnrho)) * divergence(uu) * divergence(uu);
+ const Vector j = (Scalar(1.) / AC_mu0) *
+ (gradient_of_divergence(aa) - laplace_vec(aa)); // Current density
+ const Scalar RHS = H_CONST - C_CONST + AC_eta * (AC_mu0)*dot(j, j) +
+ Scalar(2.) * exp(value(lnrho)) * AC_nu_visc * contract(S) +
+ AC_zeta * exp(value(lnrho)) * divergence(uu) * divergence(uu);
- return - dot(value(uu), gradient(ss))
- + inv_pT * RHS
- + heat_conduction(ss, lnrho);
+ return -dot(value(uu), gradient(ss)) + inv_pT * RHS + heat_conduction(ss, lnrho);
}
#endif
@@ -193,14 +173,15 @@ entropy(in ScalarField ss, in VectorField uu, in ScalarField lnrho, in VectorFie
Scalar
heat_transfer(in VectorField uu, in ScalarField lnrho, in ScalarField tt)
{
- const Matrix S = stress_tensor(uu);
- const Scalar heat_diffusivity_k = 0.0008; //8e-4;
- return -dot(value(uu), gradient(tt)) + heat_diffusivity_k * laplace(tt) + heat_diffusivity_k * dot(gradient(lnrho), gradient(tt)) + nu_visc * contract(S) * (Scalar(1.) / cv_sound) - (gamma - 1) * value(tt) * divergence(uu);
+ const Matrix S = stress_tensor(uu);
+ const Scalar heat_diffusivity_k = 0.0008; // 8e-4;
+ return -dot(value(uu), gradient(tt)) + heat_diffusivity_k * laplace(tt) +
+ heat_diffusivity_k * dot(gradient(lnrho), gradient(tt)) +
+ AC_nu_visc * contract(S) * (Scalar(1.) / AC_cv_sound) -
+ (AC_gamma - 1) * value(tt) * divergence(uu);
}
#endif
-
-
#if LFORCING
Vector
simple_vortex_forcing(Vector a, Vector b, Scalar magnitude)
@@ -214,13 +195,13 @@ simple_outward_flow_forcing(Vector a, Vector b, Scalar 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
+// The Pencil Code forcing_hel_noshear(), manual Eq. 222, inspired forcing function with adjustable
+// helicity
Vector
helical_forcing(Scalar magnitude, Vector k_force, Vector xx, Vector ff_re, Vector ff_im, Scalar phi)
{
// JP: This looks wrong:
- // 1) Should it be dsx * nx instead of dsx * ny?
+ // 1) Should it be AC_dsx * AC_nx instead of AC_dsx * AC_ny?
// 2) Should you also use globalGrid.n instead of the local n?
// MV: You are rigth. Made a quickfix. I did not see the error because multigpu is split
// in z direction not y direction.
@@ -229,24 +210,23 @@ helical_forcing(Scalar magnitude, Vector k_force, Vector xx, Vector ff_re, Vecto
// MV: Good idea. No an immediate priority.
// Fun related article:
// https://randomascii.wordpress.com/2014/10/09/intel-underestimates-error-bounds-by-1-3-quintillion/
- xx.x = xx.x*(2.0*M_PI/(dsx*globalGridN.x));
- xx.y = xx.y*(2.0*M_PI/(dsy*globalGridN.y));
- xx.z = xx.z*(2.0*M_PI/(dsz*globalGridN.z));
+ xx.x = xx.x * (2.0 * M_PI / (AC_dsx * globalGridN.x));
+ xx.y = xx.y * (2.0 * M_PI / (AC_dsy * globalGridN.y));
+ xx.z = xx.z * (2.0 * M_PI / (AC_dsz * globalGridN.z));
- Scalar cos_phi = cos(phi);
- Scalar sin_phi = sin(phi);
- Scalar cos_k_dot_x = cos(dot(k_force, xx));
- Scalar sin_k_dot_x = sin(dot(k_force, xx));
+ Scalar cos_phi = cos(phi);
+ Scalar sin_phi = sin(phi);
+ Scalar cos_k_dot_x = cos(dot(k_force, xx));
+ Scalar sin_k_dot_x = sin(dot(k_force, xx));
// Phase affect only the x-component
- //Scalar real_comp = cos_k_dot_x;
- //Scalar imag_comp = sin_k_dot_x;
- Scalar real_comp_phase = cos_k_dot_x*cos_phi - sin_k_dot_x*sin_phi;
- Scalar imag_comp_phase = cos_k_dot_x*sin_phi + sin_k_dot_x*cos_phi;
+ // Scalar real_comp = cos_k_dot_x;
+ // Scalar imag_comp = sin_k_dot_x;
+ Scalar real_comp_phase = cos_k_dot_x * cos_phi - sin_k_dot_x * sin_phi;
+ Scalar imag_comp_phase = cos_k_dot_x * sin_phi + sin_k_dot_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};
+ 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;
}
@@ -254,37 +234,39 @@ helical_forcing(Scalar magnitude, Vector k_force, Vector xx, Vector ff_re, Vecto
Vector
forcing(int3 globalVertexIdx, Scalar dt)
{
- Vector a = Scalar(.5) * (Vector){globalGridN.x * dsx,
- globalGridN.y * dsy,
- globalGridN.z * dsz}; // source (origin)
- Vector xx = (Vector){(globalVertexIdx.x - nx_min) * dsx,
- (globalVertexIdx.y - ny_min) * dsy,
- (globalVertexIdx.z - nz_min) * dsz}; // sink (current index)
- const Scalar cs2 = cs2_sound;
- const Scalar cs = sqrt(cs2);
+ Vector a = Scalar(.5) * (Vector){globalGridN.x * AC_dsx, globalGridN.y * AC_dsy,
+ globalGridN.z * AC_dsz}; // source (origin)
+ Vector xx = (Vector){(globalVertexIdx.x - DCONST(AC_nx_min)) * AC_dsx,
+ (globalVertexIdx.y - DCONST(AC_ny_min)) * AC_dsy,
+ (globalVertexIdx.z - DCONST(AC_nz_min)) * AC_dsz}; // sink (current index)
+ const Scalar cs2 = AC_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)};
+ // Placeholders until determined properly
+ Scalar magnitude = AC_forcing_magnitude;
+ Scalar phase = AC_forcing_phase;
+ Vector k_force = (Vector){AC_k_forcex, AC_k_forcey, AC_k_forcez};
+ Vector ff_re = (Vector){AC_ff_hel_rex, AC_ff_hel_rey, AC_ff_hel_rez};
+ Vector ff_im = (Vector){AC_ff_hel_imx, AC_ff_hel_imy, AC_ff_hel_imz};
+ // Determine that forcing funtion type at this point.
+ // Vector force = simple_vortex_forcing(a, xx, magnitude);
+ // Vector force = simple_outward_flow_forcing(a, xx, magnitude);
+ Vector force = helical_forcing(magnitude, k_force, xx, ff_re, ff_im, phase);
- //Determine that forcing funtion type at this point.
- //Vector force = simple_vortex_forcing(a, xx, magnitude);
- //Vector force = simple_outward_flow_forcing(a, xx, magnitude);
- Vector force = helical_forcing(magnitude, k_force, xx, ff_re,ff_im, phase);
+ // Scaling N = magnitude*cs*sqrt(k*cs/dt) * dt
+ const Scalar NN = cs * sqrt(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;
- //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}; }
+ if (is_valid(force)) {
+ return force;
+ }
+ else {
+ return (Vector){0, 0, 0};
+ }
}
#endif // LFORCING
@@ -294,12 +276,11 @@ in ScalarField lnrho(VTXBUF_LNRHO);
out ScalarField out_lnrho(VTXBUF_LNRHO);
in VectorField uu(VTXBUF_UUX, VTXBUF_UUY, VTXBUF_UUZ);
-out VectorField out_uu(VTXBUF_UUX,VTXBUF_UUY,VTXBUF_UUZ);
-
+out VectorField out_uu(VTXBUF_UUX, VTXBUF_UUY, VTXBUF_UUZ);
#if LMAGNETIC
-in VectorField aa(VTXBUF_AX,VTXBUF_AY,VTXBUF_AZ);
-out VectorField out_aa(VTXBUF_AX,VTXBUF_AY,VTXBUF_AZ);
+in VectorField aa(VTXBUF_AX, VTXBUF_AY, VTXBUF_AZ);
+out VectorField out_aa(VTXBUF_AX, VTXBUF_AY, VTXBUF_AZ);
#endif
#if LENTROPY
@@ -313,26 +294,27 @@ out ScalarField out_tt(VTXBUF_TEMPERATURE);
#endif
Kernel void
-solve(Scalar dt) {
+solve(Scalar dt)
+{
out_lnrho = rk3(out_lnrho, lnrho, continuity(uu, lnrho), dt);
- #if LMAGNETIC
+#if LMAGNETIC
out_aa = rk3(out_aa, aa, induction(uu, aa), dt);
- #endif
+#endif
- #if LENTROPY
- out_uu = rk3(out_uu, uu, momentum(uu, lnrho, ss, aa), dt);
- out_ss = rk3(out_ss, ss, entropy(ss, uu, lnrho, aa), dt);
- #elif LTEMPERATURE
- out_uu = rk3(out_uu, uu, momentum(uu, lnrho, tt), dt);
- out_tt = rk3(out_tt, tt, heat_transfer(uu, lnrho, tt), dt);
- #else
- out_uu = rk3(out_uu, uu, momentum(uu, lnrho), dt);
- #endif
+#if LENTROPY
+ out_uu = rk3(out_uu, uu, momentum(uu, lnrho, ss, aa), dt);
+ out_ss = rk3(out_ss, ss, entropy(ss, uu, lnrho, aa), dt);
+#elif LTEMPERATURE
+ out_uu = rk3(out_uu, uu, momentum(uu, lnrho, tt), dt);
+ out_tt = rk3(out_tt, tt, heat_transfer(uu, lnrho, tt), dt);
+#else
+ out_uu = rk3(out_uu, uu, momentum(uu, lnrho), dt);
+#endif
- #if LFORCING
+#if LFORCING
if (step_number == 2) {
out_uu = out_uu + forcing(globalVertexIdx, dt);
}
- #endif
+#endif
}
diff --git a/acc/src/acc.l b/acc/src/acc.l
index b180ae8..76104d8 100644
--- a/acc/src/acc.l
+++ b/acc/src/acc.l
@@ -15,7 +15,7 @@ L [a-zA-Z_]
"void" { return VOID; } /* Rest of the types inherited from C */
"int" { return INT; }
"int3" { return INT3; }
-"ScalarField" { return SCALAR; }
+"ScalarField" { return SCALARFIELD; }
"VectorField" { return VECTOR; }
"Kernel" { return KERNEL; } /* Function specifiers */
diff --git a/acc/src/acc.y b/acc/src/acc.y
index 7138b1b..0bd1d19 100644
--- a/acc/src/acc.y
+++ b/acc/src/acc.y
@@ -16,7 +16,7 @@ int yyget_lineno();
%token CONSTANT IN OUT UNIFORM
%token IDENTIFIER NUMBER
%token RETURN
-%token SCALAR VECTOR MATRIX
+%token SCALAR VECTOR MATRIX SCALARFIELD
%token VOID INT INT3
%token IF ELSE FOR WHILE ELIF
%token LEQU LAND LOR LLEQU
@@ -209,6 +209,7 @@ type_specifier: VOID
| SCALAR { $$ = astnode_create(NODE_TYPE_SPECIFIER, NULL, NULL); $$->token = SCALAR; }
| VECTOR { $$ = astnode_create(NODE_TYPE_SPECIFIER, NULL, NULL); $$->token = VECTOR; }
| MATRIX { $$ = astnode_create(NODE_TYPE_SPECIFIER, NULL, NULL); $$->token = MATRIX; }
+ | SCALARFIELD { $$ = astnode_create(NODE_TYPE_SPECIFIER, NULL, NULL); $$->token = SCALARFIELD; }
;
identifier: IDENTIFIER { $$ = astnode_create(NODE_IDENTIFIER, NULL, NULL); astnode_set_buffer(yytext, $$); }
diff --git a/acc/src/code_generator.c b/acc/src/code_generator.c
index 2de30d7..e560a70 100644
--- a/acc/src/code_generator.c
+++ b/acc/src/code_generator.c
@@ -54,12 +54,13 @@ static const char* translation_table[TRANSLATION_TABLE_SIZE] = {
[WHILE] = "while",
[FOR] = "for",
// Type specifiers
- [VOID] = "void",
- [INT] = "int",
- [INT3] = "int3",
- [SCALAR] = "AcReal",
- [VECTOR] = "AcReal3",
- [MATRIX] = "AcMatrix",
+ [VOID] = "void",
+ [INT] = "int",
+ [INT3] = "int3",
+ [SCALAR] = "AcReal",
+ [VECTOR] = "AcReal3",
+ [MATRIX] = "AcMatrix",
+ [SCALARFIELD] = "AcReal",
// Type qualifiers
[KERNEL] = "template static "
"__global__", //__launch_bounds__(RK_THREADBLOCK_SIZE,
@@ -228,14 +229,18 @@ translate_latest_symbol(void)
}
// UNIFORM
else if (symbol->type_qualifier == UNIFORM) {
+ // if (compilation_type != STENCIL_HEADER) {
+ // printf("ERROR: %s can only be used in stencil headers\n", translation_table[UNIFORM]);
+ //}
/* Do nothing */
}
// IN / OUT
else if (symbol->type != SYMBOLTYPE_FUNCTION_PARAMETER &&
(symbol->type_qualifier == IN || symbol->type_qualifier == OUT)) {
- printf("static __device__ const %s %s%s", symbol->type_specifier == SCALAR ? "int" : "int3",
- inout_name_prefix, symbol_table[handle].identifier);
+ printf("static __device__ const %s %s%s",
+ symbol->type_specifier == SCALARFIELD ? "int" : "int3", inout_name_prefix,
+ symbol_table[handle].identifier);
if (symbol->type_specifier == VECTOR)
printf(" = make_int3");
}
@@ -373,6 +378,8 @@ traverse(const ASTNode* node)
// printf("%s%s", inout_name_prefix, symbol->identifier);
//}
if (symbol->type_qualifier == UNIFORM) {
+ printf("DCONST(%s) ", symbol->identifier);
+ /*
if (symbol->type_specifier == SCALAR)
printf("DCONST_REAL(AC_%s) ", symbol->identifier);
else if (symbol->type_specifier == INT)
@@ -380,6 +387,7 @@ traverse(const ASTNode* node)
else
printf("INVALID UNIFORM type specifier %s with %s\n",
translate(symbol->type_specifier), symbol->identifier);
+ */
}
else {
// Do a regular translation
@@ -549,6 +557,68 @@ generate_preprocessed_structures(void)
");
}
+static void
+generate_header(void)
+{
+ printf("\n#pragma once\n");
+
+ // Int params
+ printf("#define AC_FOR_USER_INT_PARAM_TYPES(FUNC)");
+ for (int i = 0; i < num_symbols; ++i) {
+ if (symbol_table[i].type_specifier == INT) {
+ printf("\\\nFUNC(%s),", symbol_table[i].identifier);
+ }
+ }
+ printf("\n\n");
+
+ // Int3 params
+ printf("#define AC_FOR_USER_INT3_PARAM_TYPES(FUNC)");
+ for (int i = 0; i < num_symbols; ++i) {
+ if (symbol_table[i].type_specifier == INT3) {
+ printf("\\\nFUNC(%s),", symbol_table[i].identifier);
+ }
+ }
+ printf("\n\n");
+
+ // Scalar params
+ printf("#define AC_FOR_USER_REAL_PARAM_TYPES(FUNC)");
+ for (int i = 0; i < num_symbols; ++i) {
+ if (symbol_table[i].type_specifier == SCALAR) {
+ printf("\\\nFUNC(%s),", symbol_table[i].identifier);
+ }
+ }
+ printf("\n\n");
+
+ // Vector params
+ printf("#define AC_FOR_USER_REAL3_PARAM_TYPES(FUNC)");
+ for (int i = 0; i < num_symbols; ++i) {
+ if (symbol_table[i].type_specifier == VECTOR) {
+ printf("\\\nFUNC(%s),", symbol_table[i].identifier);
+ }
+ }
+ printf("\n\n");
+
+ // Scalar fields
+ printf("#define AC_FOR_VTXBUF_HANDLES(FUNC)");
+ for (int i = 0; i < num_symbols; ++i) {
+ if (symbol_table[i].type_specifier == SCALARFIELD) {
+ printf("\\\nFUNC(%s),", symbol_table[i].identifier);
+ }
+ }
+ printf("\n\n");
+
+ /*
+ printf("\n");
+ printf("typedef struct {\n");
+ for (int i = 0; i < num_symbols; ++i) {
+ if (symbol_table[i].type_qualifier == PREPROCESSED)
+ printf("%s %s;\n", translate(symbol_table[i].type_specifier),
+ symbol_table[i].identifier);
+ }
+ printf("} %sData;\n", translate(SCALAR));
+ */
+}
+
int
main(int argc, char** argv)
{
@@ -557,7 +627,7 @@ main(int argc, char** argv)
compilation_type = STENCIL_ASSEMBLY;
else if (!strcmp(argv[1], "-sps"))
compilation_type = STENCIL_PROCESS;
- else if (!strcmp(argv[1], "-hh"))
+ else if (!strcmp(argv[1], "-sdh"))
compilation_type = STENCIL_HEADER;
else
printf("Unknown flag %s. Generating stencil assembly.\n", argv[1]);
@@ -584,6 +654,8 @@ main(int argc, char** argv)
traverse(root);
if (compilation_type == STENCIL_ASSEMBLY)
generate_preprocessed_structures();
+ else if (compilation_type == STENCIL_HEADER)
+ generate_header();
// print_symbol_table();
diff --git a/include/astaroth_defines.h b/include/astaroth_defines.h
index 0d340a0..58525ba 100644
--- a/include/astaroth_defines.h
+++ b/include/astaroth_defines.h
@@ -45,6 +45,8 @@ typedef struct {
#endif // __CUDACC__
// Library flags
+#define STENCIL_ORDER (6)
+#define NGHOST (STENCIL_ORDER / 2)
#define VERBOSE_PRINTING (1)
// Built-in types and parameters
diff --git a/scripts/compile_acc.sh b/scripts/compile_acc.sh
index 8648079..7ebbd4c 100755
--- a/scripts/compile_acc.sh
+++ b/scripts/compile_acc.sh
@@ -8,17 +8,18 @@ fi
KERNEL_DIR=${AC_HOME}"/src/core/kernels"
ACC_DIR=${AC_HOME}"/acc"
-ACC_DEFAULT_HEADER="mhd_solver/stencil_defines.h"
ACC_DEFAULT_SAS="mhd_solver/stencil_assembly.sas"
ACC_DEFAULT_SPS="mhd_solver/stencil_process.sps"
+ACC_DEFAULT_HEADER="mhd_solver/stencil_definition.sdh"
+ACC_DEFAULT_INCLUDE_DIR="mhd_solver"
${ACC_DIR}/clean.sh
${ACC_DIR}/build_acc.sh
-
-ACC_HEADER=${ACC_DEFAULT_HEADER}
ACC_SAS=${ACC_DEFAULT_SAS}
ACC_SPS=${ACC_DEFAULT_SPS}
+ACC_HEADER=${ACC_DEFAULT_HEADER}
+ACC_INCLUDE_DIR=${ACC_DEFAULT_INCLUDE_DIR}
while [ "$#" -gt 0 ]
do
@@ -56,9 +57,17 @@ echo "Header file:" ${ACC_DIR}/${ACC_HEADER}
echo "Assembly file: ${ACC_DIR}/${ACC_SAS}"
echo "Process file: ${ACC_DIR}/${ACC_SPS}"
-cd ${KERNEL_DIR}
-${ACC_DIR}/compile.sh ${ACC_DIR}/${ACC_SAS} ${ACC_DIR}/${ACC_HEADER}
-${ACC_DIR}/compile.sh ${ACC_DIR}/${ACC_SPS} ${ACC_DIR}/${ACC_HEADER}
+cd ${ACC_DIR}/${ACC_INCLUDE_DIR}
+echo ${PWD}
+${ACC_DIR}/compile.sh ${ACC_DIR}/${ACC_SAS}
+${ACC_DIR}/compile.sh ${ACC_DIR}/${ACC_SPS}
+${ACC_DIR}/compile.sh ${ACC_DIR}/${ACC_HEADER}
-echo "Linking: " ${ACC_DIR}/${ACC_HEADER} " -> " ${AC_HOME}/include/stencil_defines.h
-ln -sf ${ACC_DIR}/${ACC_HEADER} ${AC_HOME}/include/stencil_defines.h
+echo "Moving stencil_assembly.cuh -> ${AC_HOME}/src/core/kernels"
+mv stencil_assembly.cuh ${AC_HOME}/src/core/kernels
+
+echo "Moving stencil_process.cuh -> ${AC_HOME}/src/core/kernels"
+mv stencil_process.cuh ${AC_HOME}/src/core/kernels
+
+echo "Moving stencil_defines.cuh -> ${AC_HOME}/include"
+mv stencil_defines.h ${AC_HOME}/include
diff --git a/src/core/device.cu b/src/core/device.cu
index 8b89a78..15f0c1e 100644
--- a/src/core/device.cu
+++ b/src/core/device.cu
@@ -40,26 +40,31 @@ typedef struct {
} VertexBufferArray;
__constant__ AcMeshInfo d_mesh_info;
-static inline int __device__
+static int __device__ __forceinline__
DCONST(const AcIntParam param)
{
return d_mesh_info.int_params[param];
}
-static inline int3 __device__
+static int3 __device__ __forceinline__
DCONST(const AcInt3Param param)
{
return d_mesh_info.int3_params[param];
}
-static inline AcReal __device__
+static AcReal __device__ __forceinline__
DCONST(const AcRealParam param)
{
return d_mesh_info.real_params[param];
}
-static inline AcReal3 __device__
+static AcReal3 __device__ __forceinline__
DCONST(const AcReal3Param param)
{
return d_mesh_info.real3_params[param];
}
+constexpr VertexBufferHandle
+DCONST(const VertexBufferHandle handle)
+{
+ return handle;
+}
#define DCONST_INT(x) DCONST(x)
#define DCONST_INT3(x) DCONST(x)
#define DCONST_REAL(x) DCONST(x)
@@ -103,7 +108,7 @@ struct device_s {
};
// clang-format off
-static __global__ void dummy_kernel(void) {}
+static __global__ void dummy_kernel(void) { DCONST((AcIntParam)0); DCONST((AcInt3Param)0); DCONST((AcRealParam)0); DCONST((AcReal3Param)0); }
// clang-format on
AcResult
diff --git a/src/mpitest/CMakeLists.txt b/src/mpitest/CMakeLists.txt
index c64105d..d6c5309 100644
--- a/src/mpitest/CMakeLists.txt
+++ b/src/mpitest/CMakeLists.txt
@@ -8,5 +8,5 @@ set(CMAKE_C_STANDARD_REQUIRED ON)
find_package(MPI REQUIRED)
add_executable(mpitest main.c)
-target_include_directories(mpitest PRIVATE ${MPI_C_INCLUDE_PATH})
-target_link_libraries(mpitest PRIVATE ${MPI_C_LIBRARIES} astaroth_core)
+target_include_directories(mpitest PRIVATE ${CMAKE_SOURCE_DIR}/src/standalone ${MPI_C_INCLUDE_PATH})
+target_link_libraries(mpitest astaroth_core astaroth_standalone ${MPI_C_LIBRARIES})
diff --git a/src/mpitest/main.c b/src/mpitest/main.c
index a07522e..2ca34d3 100644
--- a/src/mpitest/main.c
+++ b/src/mpitest/main.c
@@ -16,13 +16,120 @@
You should have received a copy of the GNU General Public License
along with Astaroth. If not, see .
*/
+/**
+ Running: mpirun -np
+*/
+#undef NDEBUG // Assert always
+#include
#include
#include
+#include
#include "astaroth.h"
+#include "autotest.h"
#include
+// From Astaroth Standalone
+#include "config_loader.h"
+#include "model/host_memory.h"
+
+static void
+distribute_mesh(const AcMesh* src, AcMesh* dst)
+{
+ MPI_Datatype datatype = MPI_FLOAT;
+ if (sizeof(AcReal) == 8)
+ datatype = MPI_DOUBLE;
+
+ int process_id, num_processes;
+ MPI_Comm_rank(MPI_COMM_WORLD, &process_id);
+ MPI_Comm_size(MPI_COMM_WORLD, &num_processes);
+
+ const size_t count = acVertexBufferSize(dst->info);
+ for (int i = 0; i < NUM_VTXBUF_HANDLES; ++i) {
+
+ // Communicate to self
+ if (process_id == 0) {
+ assert(src);
+ assert(dst);
+ memcpy(&dst->vertex_buffer[i][0], //
+ &src->vertex_buffer[i][0], //
+ count * sizeof(src->vertex_buffer[i][0]));
+ }
+ // Communicate to others
+ for (int j = 1; j < num_processes; ++j) {
+ if (process_id == 0) {
+ assert(src);
+
+ // Send
+ // TODO RECHECK THESE j INDICES
+ const size_t src_idx = j * dst->info.int_params[AC_mx] *
+ dst->info.int_params[AC_my] * src->info.int_params[AC_nz] /
+ num_processes;
+
+ MPI_Send(&src->vertex_buffer[i][src_idx], count, datatype, j, 0, MPI_COMM_WORLD);
+ }
+ else {
+ assert(dst);
+
+ // Recv
+ const size_t dst_idx = 0;
+ MPI_Status status;
+ MPI_Recv(&dst->vertex_buffer[i][dst_idx], count, datatype, 0, 0, MPI_COMM_WORLD,
+ &status);
+ }
+ }
+ }
+}
+
+static void
+gather_mesh(const AcMesh* src, AcMesh* dst)
+{
+ MPI_Datatype datatype = MPI_FLOAT;
+ if (sizeof(AcReal) == 8)
+ datatype = MPI_DOUBLE;
+
+ int process_id, num_processes;
+ MPI_Comm_rank(MPI_COMM_WORLD, &process_id);
+ MPI_Comm_size(MPI_COMM_WORLD, &num_processes);
+
+ size_t count = acVertexBufferSize(src->info);
+
+ for (int i = 0; i < NUM_VTXBUF_HANDLES; ++i) {
+ // Communicate to self
+ if (process_id == 0) {
+ assert(src);
+ assert(dst);
+ memcpy(&dst->vertex_buffer[i][0], //
+ &src->vertex_buffer[i][0], //
+ count * sizeof(AcReal));
+ }
+
+ // Communicate to others
+ for (int j = 1; j < num_processes; ++j) {
+ if (process_id == 0) {
+ // Recv
+ // const size_t dst_idx = j * acVertexBufferCompdomainSize(dst->info);
+ const size_t dst_idx = j * dst->info.int_params[AC_mx] *
+ dst->info.int_params[AC_my] * dst->info.int_params[AC_nz] /
+ num_processes;
+
+ assert(dst_idx + count <= acVertexBufferSize(dst->info));
+ MPI_Status status;
+ MPI_Recv(&dst->vertex_buffer[i][dst_idx], count, datatype, j, 0, MPI_COMM_WORLD,
+ &status);
+ }
+ else {
+ // Send
+ const size_t src_idx = 0;
+
+ assert(src_idx + count <= acVertexBufferSize(src->info));
+ MPI_Send(&src->vertex_buffer[i][src_idx], count, datatype, 0, 0, MPI_COMM_WORLD);
+ }
+ }
+ }
+}
+
int
main(void)
{
@@ -37,14 +144,39 @@ main(void)
MPI_Get_processor_name(processor_name, &name_len);
printf("Processor %s. Process %d of %d.\n", processor_name, process_id, num_processes);
- AcMeshInfo info = {
- .int_params[AC_nx] = 128,
- .int_params[AC_ny] = 64,
- .int_params[AC_nz] = 32,
- };
- acInit(info);
- acIntegrate(0.1f);
- acQuit();
+ AcMeshInfo mesh_info;
+ load_config(&mesh_info);
+ update_config(&mesh_info);
+
+ AcMesh* main_mesh = NULL;
+ ModelMesh* model_mesh = NULL;
+ if (process_id == 0) {
+ main_mesh = acmesh_create(mesh_info);
+ acmesh_init_to(INIT_TYPE_RANDOM, main_mesh);
+ model_mesh = modelmesh_create(mesh_info);
+ acmesh_to_modelmesh(*main_mesh, model_mesh);
+ }
+
+ AcMeshInfo submesh_info = mesh_info;
+ submesh_info.int_params[AC_nz] /= num_processes;
+ update_config(&submesh_info);
+
+ AcMesh* submesh = acmesh_create(submesh_info);
+
+ /////////////////////
+ distribute_mesh(main_mesh, submesh);
+ gather_mesh(submesh, main_mesh);
+ /////////////////////////
+ // Autotest
+ bool is_acceptable = verify_meshes(*model_mesh, *main_mesh);
+ /////
+
+ acmesh_destroy(submesh);
+
+ if (process_id == 0) {
+ modelmesh_destroy(model_mesh);
+ acmesh_destroy(main_mesh);
+ }
MPI_Finalize();
return EXIT_SUCCESS;
diff --git a/src/standalone/CMakeLists.txt b/src/standalone/CMakeLists.txt
index ea1d04c..0a61ede 100644
--- a/src/standalone/CMakeLists.txt
+++ b/src/standalone/CMakeLists.txt
@@ -25,10 +25,11 @@ add_compile_options(-pipe ${OpenMP_CXX_FLAGS})
add_compile_options(-Wall -Wextra -Werror -Wdouble-promotion -Wfloat-conversion)# -Wshadow)
## Compile and link
-add_library(astaroth_standalone ${SOURCES})
+add_library(astaroth_standalone STATIC ${SOURCES})
+target_link_libraries(astaroth_standalone PRIVATE astaroth_core "${OpenMP_CXX_FLAGS}" ${SDL2_LIBRARY})
add_executable(ac_run main.cc)
-target_link_libraries(ac_run PRIVATE astaroth_standalone astaroth_core "${OpenMP_CXX_FLAGS}" ${SDL2_LIBRARY})
+target_link_libraries(ac_run PRIVATE astaroth_standalone)
# Define the config directory
if (ALTER_CONF)
diff --git a/src/standalone/model/model_rk3.cc b/src/standalone/model/model_rk3.cc
index 756e3a7..5fab4b4 100644
--- a/src/standalone/model/model_rk3.cc
+++ b/src/standalone/model/model_rk3.cc
@@ -31,6 +31,16 @@
#include "host_memory.h"
#include "model_boundconds.h"
+// Standalone flags
+#define LDENSITY (1)
+#define LHYDRO (1)
+#define LMAGNETIC (1)
+#define LENTROPY (1)
+#define LTEMPERATURE (0)
+#define LFORCING (1)
+#define LUPWD (1)
+#define AC_THERMAL_CONDUCTIVITY (AcReal(0.001)) // TODO: make an actual config parameter
+
typedef struct {
ModelScalar x, y, z;
} ModelVector;