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

Merged in dsl_parameter_overhaul_2019-08-19 (pull request #6)

Browse Source 
DSL user parameter overhaul
This commit is contained in:
jpekkila
2019-08-30 02:13:35 +00:00
committed by Miikka Väisälä
parent bd549f5d28 9fa39b66fc
commit be9ec6293b
16 changed files with 589 additions and 411 deletions
Download Patch File Download Diff File Expand all files Collapse all files

12
acc/compile.sh
View File

@@ -8,17 +8,21 @@ FILENAME="${FULL_NAME%.*}"
EXTENSION="${FULL_NAME##*.}" EXTENSION="${FULL_NAME##*.}"
if [ "${EXTENSION}" = "sas" ]; then if [ "${EXTENSION}" = "sas" ]; then
echo "Generating stencil assembly stage ${FILENAME}.sas -> stencil_assembly.cuh"
COMPILE_FLAGS="-sas" # Generate stencil assembly stage COMPILE_FLAGS="-sas" # Generate stencil assembly stage
CUH_FILENAME="stencil_assembly.cuh" CUH_FILENAME="stencil_assembly.cuh"
echo "Generating stencil assembly stage ${FILENAME}.sas -> ${CUH_FILENAME}"
elif [ "${EXTENSION}" = "sps" ]; then elif [ "${EXTENSION}" = "sps" ]; then
echo "Generating stencil processing stage: ${FILENAME}.sps -> stencil_process.cuh"
COMPILE_FLAGS="-sps" # Generate stencil processing stage COMPILE_FLAGS="-sps" # Generate stencil processing stage
CUH_FILENAME="stencil_process.cuh" 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 else
echo "Error: unknown extension" ${EXTENSION} "of file" ${FULL_NAME} 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 exit
fi 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}

5
acc/mhd_solver/.gitignore vendored Normal file
View File

@@ -0,0 +1,5 @@
build
testbin
# Except this file
!.gitignore

67
acc/mhd_solver/stencil_assembly.sas
View File

@@ -1,3 +1,5 @@
#include "stencil_definition.sdh"
Preprocessed Scalar Preprocessed Scalar
value(in ScalarField vertex) value(in ScalarField vertex)
{ {
@@ -7,9 +9,7 @@ value(in ScalarField vertex)
Preprocessed Vector Preprocessed Vector
gradient(in ScalarField vertex) gradient(in ScalarField vertex)
{ {
return (Vector){derx(vertexIdx, vertex), return (Vector){derx(vertexIdx, vertex), dery(vertexIdx, vertex), derz(vertexIdx, vertex)};
dery(vertexIdx, vertex),
derz(vertexIdx, vertex)};
} }
#if LUPWD #if LUPWD
@@ -17,46 +17,46 @@ gradient(in ScalarField vertex)
Preprocessed Scalar Preprocessed Scalar
der6x_upwd(in ScalarField vertex) 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* ( return (Scalar){Scalar(1.0 / 60.0) * inv_ds *
- Scalar(20.0)* vertex[vertexIdx.x, vertexIdx.y, vertexIdx.z] (-Scalar(20.0) * vertex[vertexIdx.x, vertexIdx.y, vertexIdx.z] +
+ Scalar(15.0)*(vertex[vertexIdx.x+1, vertexIdx.y, vertexIdx.z] Scalar(15.0) * (vertex[vertexIdx.x + 1, vertexIdx.y, vertexIdx.z] +
+ 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] Scalar(6.0) * (vertex[vertexIdx.x + 2, vertexIdx.y, vertexIdx.z] +
+ 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] +
+ vertex[vertexIdx.x-3, vertexIdx.y, vertexIdx.z])}; vertex[vertexIdx.x - 3, vertexIdx.y, vertexIdx.z])};
} }
Preprocessed Scalar Preprocessed Scalar
der6y_upwd(in ScalarField vertex) 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* ( return (Scalar){Scalar(1.0 / 60.0) * inv_ds *
-Scalar( 20.0)* vertex[vertexIdx.x, vertexIdx.y, vertexIdx.z] (-Scalar(20.0) * vertex[vertexIdx.x, vertexIdx.y, vertexIdx.z] +
+Scalar( 15.0)*(vertex[vertexIdx.x, vertexIdx.y+1, vertexIdx.z] Scalar(15.0) * (vertex[vertexIdx.x, vertexIdx.y + 1, vertexIdx.z] +
+ 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] Scalar(6.0) * (vertex[vertexIdx.x, vertexIdx.y + 2, vertexIdx.z] +
+ 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] +
+ vertex[vertexIdx.x, vertexIdx.y-3, vertexIdx.z])}; vertex[vertexIdx.x, vertexIdx.y - 3, vertexIdx.z])};
} }
Preprocessed Scalar Preprocessed Scalar
der6z_upwd(in ScalarField vertex) 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* ( return (Scalar){Scalar(1.0 / 60.0) * inv_ds *
-Scalar( 20.0)* vertex[vertexIdx.x, vertexIdx.y, vertexIdx.z] (-Scalar(20.0) * vertex[vertexIdx.x, vertexIdx.y, vertexIdx.z] +
+Scalar( 15.0)*(vertex[vertexIdx.x, vertexIdx.y, vertexIdx.z+1] Scalar(15.0) * (vertex[vertexIdx.x, vertexIdx.y, vertexIdx.z + 1] +
+ 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] 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 - 2]) +
+ vertex[vertexIdx.x, vertexIdx.y, vertexIdx.z+3] vertex[vertexIdx.x, vertexIdx.y, vertexIdx.z + 3] +
+ vertex[vertexIdx.x, vertexIdx.y, vertexIdx.z-3])}; vertex[vertexIdx.x, vertexIdx.y, vertexIdx.z - 3])};
} }
#endif #endif
@@ -66,9 +66,10 @@ hessian(in ScalarField vertex)
{ {
Matrix hessian; Matrix hessian;
hessian.row[0] = (Vector){derxx(vertexIdx, vertex), derxy(vertexIdx, vertex), derxz(vertexIdx, vertex)}; hessian.row[0] = (Vector){derxx(vertexIdx, vertex), derxy(vertexIdx, vertex),
hessian.row[1] = (Vector){hessian.row[0].y, deryy(vertexIdx, vertex), deryz(vertexIdx, vertex)}; derxz(vertexIdx, vertex)};
hessian.row[2] = (Vector){hessian.row[0].z, hessian.row[1].z, derzz(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; return hessian;
} }

163
acc/mhd_solver/stencil_defines.h
View File

@@ -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 <http://www.gnu.org/licenses/>.
*/
#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

117
acc/mhd_solver/stencil_definition.sdh Normal file
View File

@@ -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

332
acc/mhd_solver/stencil_process.sps
View File

@@ -1,28 +1,4 @@
// Declare uniforms (i.e. device constants) #include "stencil_definition.sdh"
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;
Vector Vector
value(in VectorField uu) value(in VectorField uu)
@@ -37,7 +13,7 @@ upwd_der6(in VectorField uu, in ScalarField lnrho)
Scalar uux = fabs(value(uu).x); Scalar uux = fabs(value(uu).x);
Scalar uuy = fabs(value(uu).y); Scalar uuy = fabs(value(uu).y);
Scalar uuz = fabs(value(uu).z); 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 #endif
@@ -48,10 +24,11 @@ gradients(in VectorField uu)
} }
Scalar Scalar
continuity(in VectorField uu, in ScalarField lnrho) { continuity(in VectorField uu, in ScalarField lnrho)
{
return -dot(value(uu), gradient(lnrho)) return -dot(value(uu), gradient(lnrho))
#if LUPWD #if LUPWD
//This is a corrective hyperdiffusion term for upwinding. // This is a corrective hyperdiffusion term for upwinding.
+ upwd_der6(uu, lnrho) + upwd_der6(uu, lnrho)
#endif #endif
- divergence(uu); - divergence(uu);
@@ -59,133 +36,136 @@ continuity(in VectorField uu, in ScalarField lnrho) {
#if LENTROPY #if LENTROPY
Vector Vector
momentum(in VectorField uu, in ScalarField lnrho, in ScalarField ss, in VectorField aa) { 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 Matrix S = stress_tensor(uu);
const Vector j = (Scalar(1.) / mu0) * (gradient_of_divergence(aa) - laplace_vec(aa)); // Current density 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); 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)); 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]*)\)
// \1 // \1
const Vector mom = - mul(gradients(uu), value(uu)) const Vector mom = -mul(gradients(uu), value(uu)) -
- cs2 * ((Scalar(1.) / cp_sound) * gradient(ss) + gradient(lnrho)) cs2 * ((Scalar(1.) / AC_cp_sound) * gradient(ss) + gradient(lnrho)) +
+ inv_rho * cross(j, B) inv_rho * cross(j, B) +
+ nu_visc * ( AC_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))) +
+ Scalar(2.) * mul(S, gradient(lnrho)) AC_zeta * gradient_of_divergence(uu);
)
+ zeta * gradient_of_divergence(uu);
return mom; return mom;
} }
#elif LTEMPERATURE #elif LTEMPERATURE
Vector Vector
momentum(in VectorField uu, in ScalarField lnrho, in ScalarField tt) { momentum(in VectorField uu, in ScalarField lnrho, in ScalarField tt)
Vector mom; {
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)) - mom = -mul(gradients(uu), value(uu)) - pressure_term +
pressure_term + AC_nu_visc * (laplace_vec(uu) + Scalar(1. / 3.) * gradient_of_divergence(uu) +
nu_visc * Scalar(2.) * mul(S, gradient(lnrho))) +
(laplace_vec(uu) + Scalar(1. / 3.) * gradient_of_divergence(uu) + AC_zeta * gradient_of_divergence(uu);
Scalar(2.) * mul(S, gradient(lnrho))) + zeta * gradient_of_divergence(uu);
#if LGRAVITY #if LGRAVITY
mom = mom - (Vector){0, 0, -10.0}; mom = mom - (Vector){0, 0, -10.0};
#endif #endif
return mom; return mom;
} }
#else #else
Vector Vector
momentum(in VectorField uu, in ScalarField lnrho) { momentum(in VectorField uu, in ScalarField lnrho)
Vector mom; {
Vector mom;
const Matrix S = stress_tensor(uu); const Matrix S = stress_tensor(uu);
// Isothermal: we have constant speed of sound // Isothermal: we have constant speed of sound
mom = -mul(gradients(uu), value(uu)) - mom = -mul(gradients(uu), value(uu)) - AC_cs2_sound * gradient(lnrho) +
cs2_sound * gradient(lnrho) + AC_nu_visc * (laplace_vec(uu) + Scalar(1. / 3.) * gradient_of_divergence(uu) +
nu_visc * Scalar(2.) * mul(S, gradient(lnrho))) +
(laplace_vec(uu) + Scalar(1. / 3.) * gradient_of_divergence(uu) + AC_zeta * gradient_of_divergence(uu);
Scalar(2.) * mul(S, gradient(lnrho))) + zeta * gradient_of_divergence(uu);
#if LGRAVITY #if LGRAVITY
mom = mom - (Vector){0, 0, -10.0}; mom = mom - (Vector){0, 0, -10.0};
#endif #endif
return mom; return mom;
} }
#endif #endif
Vector Vector
induction(in VectorField uu, in VectorField 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, // Note: We do (-nabla^2 A + nabla(nabla dot A)) instead of (nabla x (nabla
// yes this actually works. See pg.28 in arXiv:astro-ph/0109497) // x A)) in order to avoid taking the first derivative twice (did the math,
// u cross B - ETA * mu0 * (mu0^-1 * [- laplace A + grad div A ]) // yes this actually works. See pg.28 in arXiv:astro-ph/0109497)
const Vector B = curl(aa); // u cross B - AC_eta * AC_mu0 * (AC_mu0^-1 * [- laplace A + grad div A ])
const Vector grad_div = gradient_of_divergence(aa); const Vector B = curl(aa);
const Vector lap = laplace_vec(aa); const Vector grad_div = gradient_of_divergence(aa);
const Vector lap = laplace_vec(aa);
// Note, mu0 is cancelled out // Note, AC_mu0 is cancelled out
const Vector ind = cross(value(uu), B) - eta * (grad_div - lap); const Vector ind = cross(value(uu), B) - AC_eta * (grad_div - lap);
return ind; return ind;
} }
#if LENTROPY #if LENTROPY
Scalar Scalar
lnT( in ScalarField ss, in ScalarField lnrho) { lnT(in ScalarField ss, in ScalarField lnrho)
const Scalar lnT = lnT0 + gamma * value(ss) / cp_sound + {
(gamma - Scalar(1.)) * (value(lnrho) - lnrho0); const Scalar lnT = AC_lnT0 + AC_gamma * value(ss) / AC_cp_sound +
return lnT; (AC_gamma - Scalar(1.)) * (value(lnrho) - AC_lnrho0);
return lnT;
} }
// Nabla dot (K nabla T) / (rho T) // Nabla dot (K nabla T) / (rho T)
Scalar Scalar
heat_conduction( in ScalarField ss, in ScalarField lnrho) { heat_conduction(in ScalarField ss, in ScalarField lnrho)
const Scalar inv_cp_sound = AcReal(1.) / cp_sound; {
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) + const Scalar first_term = AC_gamma * inv_AC_cp_sound * laplace(ss) +
(gamma - AcReal(1.)) * laplace(lnrho); (AC_gamma - AcReal(1.)) * laplace(lnrho);
const Vector second_term = gamma * inv_cp_sound * gradient(ss) + const Vector second_term = AC_gamma * inv_AC_cp_sound * gradient(ss) +
(gamma - AcReal(1.)) * gradient(lnrho); (AC_gamma - AcReal(1.)) * gradient(lnrho);
const Vector third_term = gamma * (inv_cp_sound * gradient(ss) + const Vector third_term = AC_gamma * (inv_AC_cp_sound * gradient(ss) + gradient(lnrho)) +
gradient(lnrho)) + grad_ln_chi; grad_ln_chi;
const Scalar chi = AC_THERMAL_CONDUCTIVITY / (exp(value(lnrho)) * cp_sound); const Scalar chi = AC_THERMAL_CONDUCTIVITY / (exp(value(lnrho)) * AC_cp_sound);
return cp_sound * chi * (first_term + dot(second_term, third_term)); return AC_cp_sound * chi * (first_term + dot(second_term, third_term));
} }
Scalar Scalar
heating(const int i, const int j, const int k) { heating(const int i, const int j, const int k)
return 1; {
return 1;
} }
Scalar Scalar
entropy(in ScalarField ss, in VectorField uu, in ScalarField lnrho, in VectorField aa) { entropy(in ScalarField ss, in VectorField uu, in ScalarField lnrho, in VectorField aa)
const Matrix S = stress_tensor(uu); {
const Matrix S = stress_tensor(uu);
const Scalar inv_pT = Scalar(1.) / (exp(value(lnrho)) * exp(lnT(ss, lnrho))); 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 Vector j = (Scalar(1.) / AC_mu0) *
const Scalar RHS = H_CONST - C_CONST (gradient_of_divergence(aa) - laplace_vec(aa)); // Current density
+ eta * (mu0) * dot(j, j) const Scalar RHS = H_CONST - C_CONST + AC_eta * (AC_mu0)*dot(j, j) +
+ Scalar(2.) * exp(value(lnrho)) * nu_visc * contract(S) Scalar(2.) * exp(value(lnrho)) * AC_nu_visc * contract(S) +
+ zeta * exp(value(lnrho)) * divergence(uu) * divergence(uu); AC_zeta * exp(value(lnrho)) * divergence(uu) * divergence(uu);
return - dot(value(uu), gradient(ss)) return -dot(value(uu), gradient(ss)) + inv_pT * RHS + heat_conduction(ss, lnrho);
+ inv_pT * RHS
+ heat_conduction(ss, lnrho);
} }
#endif #endif
@@ -193,14 +173,15 @@ entropy(in ScalarField ss, in VectorField uu, in ScalarField lnrho, in VectorFie
Scalar Scalar
heat_transfer(in VectorField uu, in ScalarField lnrho, in ScalarField tt) heat_transfer(in VectorField uu, in ScalarField lnrho, in ScalarField tt)
{ {
const Matrix S = stress_tensor(uu); const Matrix S = stress_tensor(uu);
const Scalar heat_diffusivity_k = 0.0008; //8e-4; 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); 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 #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)
@@ -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 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
// The Pencil Code forcing_hel_noshear(), manual Eq. 222, inspired forcing function with adjustable helicity // helicity
Vector Vector
helical_forcing(Scalar magnitude, Vector k_force, Vector xx, Vector ff_re, Vector ff_im, Scalar phi) helical_forcing(Scalar magnitude, Vector k_force, Vector xx, Vector ff_re, Vector ff_im, Scalar phi)
{ {
// JP: This looks wrong: // 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? // 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 // MV: You are rigth. Made a quickfix. I did not see the error because multigpu is split
// in z direction not y direction. // 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. // MV: Good idea. No an immediate priority.
// Fun related article: // Fun related article:
// https://randomascii.wordpress.com/2014/10/09/intel-underestimates-error-bounds-by-1-3-quintillion/ // 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.x = xx.x * (2.0 * M_PI / (AC_dsx * globalGridN.x));
xx.y = xx.y*(2.0*M_PI/(dsy*globalGridN.y)); xx.y = xx.y * (2.0 * M_PI / (AC_dsy * globalGridN.y));
xx.z = xx.z*(2.0*M_PI/(dsz*globalGridN.z)); xx.z = xx.z * (2.0 * M_PI / (AC_dsz * globalGridN.z));
Scalar cos_phi = cos(phi); Scalar cos_phi = cos(phi);
Scalar sin_phi = sin(phi); Scalar sin_phi = sin(phi);
Scalar cos_k_dot_x = cos(dot(k_force, xx)); Scalar cos_k_dot_x = cos(dot(k_force, xx));
Scalar sin_k_dot_x = sin(dot(k_force, xx)); Scalar sin_k_dot_x = sin(dot(k_force, xx));
// Phase affect only the x-component // Phase affect only the x-component
//Scalar real_comp = cos_k_dot_x; // Scalar real_comp = cos_k_dot_x;
//Scalar imag_comp = sin_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 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 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,
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.y*real_comp_phase - ff_im.y*imag_comp_phase, ff_re.z * real_comp_phase - ff_im.z * imag_comp_phase};
ff_re.z*real_comp_phase - ff_im.z*imag_comp_phase};
return force; return force;
} }
@@ -254,37 +234,39 @@ helical_forcing(Scalar magnitude, Vector k_force, Vector xx, Vector ff_re, Vecto
Vector Vector
forcing(int3 globalVertexIdx, Scalar dt) forcing(int3 globalVertexIdx, Scalar dt)
{ {
Vector a = Scalar(.5) * (Vector){globalGridN.x * dsx, Vector a = Scalar(.5) * (Vector){globalGridN.x * AC_dsx, globalGridN.y * AC_dsy,
globalGridN.y * dsy, globalGridN.z * AC_dsz}; // source (origin)
globalGridN.z * dsz}; // source (origin) Vector xx = (Vector){(globalVertexIdx.x - DCONST(AC_nx_min)) * AC_dsx,
Vector xx = (Vector){(globalVertexIdx.x - nx_min) * dsx, (globalVertexIdx.y - DCONST(AC_ny_min)) * AC_dsy,
(globalVertexIdx.y - ny_min) * dsy, (globalVertexIdx.z - DCONST(AC_nz_min)) * AC_dsz}; // sink (current index)
(globalVertexIdx.z - nz_min) * dsz}; // sink (current index) const Scalar cs2 = AC_cs2_sound;
const Scalar cs2 = cs2_sound; const Scalar cs = sqrt(cs2);
const Scalar cs = sqrt(cs2);
//Placeholders until determined properly // Placeholders until determined properly
Scalar magnitude = DCONST_REAL(AC_forcing_magnitude); Scalar magnitude = AC_forcing_magnitude;
Scalar phase = DCONST_REAL(AC_forcing_phase); Scalar phase = AC_forcing_phase;
Vector k_force = (Vector){ DCONST_REAL(AC_k_forcex), DCONST_REAL(AC_k_forcey), DCONST_REAL(AC_k_forcez)}; Vector k_force = (Vector){AC_k_forcex, AC_k_forcey, 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_re = (Vector){AC_ff_hel_rex, AC_ff_hel_rey, 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)}; 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. // Scaling N = magnitude*cs*sqrt(k*cs/dt) * dt
//Vector force = simple_vortex_forcing(a, xx, magnitude); const Scalar NN = cs * sqrt(AC_kaver * cs);
//Vector force = simple_outward_flow_forcing(a, xx, magnitude); // MV: Like in the Pencil Code. I don't understandf the logic here.
Vector force = helical_forcing(magnitude, k_force, xx, ff_re,ff_im, phase); 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 if (is_valid(force)) {
const Scalar NN = cs*sqrt(DCONST_REAL(AC_kaver)*cs); return force;
//MV: Like in the Pencil Code. I don't understandf the logic here. }
force.x = sqrt(dt)*NN*force.x; else {
force.y = sqrt(dt)*NN*force.y; return (Vector){0, 0, 0};
force.z = sqrt(dt)*NN*force.z; }
if (is_valid(force)) { return force; }
else { return (Vector){0, 0, 0}; }
} }
#endif // LFORCING #endif // LFORCING
@@ -294,12 +276,11 @@ in ScalarField lnrho(VTXBUF_LNRHO);
out ScalarField out_lnrho(VTXBUF_LNRHO); out ScalarField out_lnrho(VTXBUF_LNRHO);
in VectorField uu(VTXBUF_UUX, VTXBUF_UUY, VTXBUF_UUZ); 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 #if LMAGNETIC
in VectorField 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); out VectorField out_aa(VTXBUF_AX, VTXBUF_AY, VTXBUF_AZ);
#endif #endif
#if LENTROPY #if LENTROPY
@@ -313,26 +294,27 @@ out ScalarField out_tt(VTXBUF_TEMPERATURE);
#endif #endif
Kernel void Kernel void
solve(Scalar dt) { solve(Scalar dt)
{
out_lnrho = rk3(out_lnrho, lnrho, continuity(uu, lnrho), 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); out_aa = rk3(out_aa, aa, induction(uu, aa), dt);
#endif #endif
#if LENTROPY #if LENTROPY
out_uu = rk3(out_uu, uu, momentum(uu, lnrho, ss, aa), dt); out_uu = rk3(out_uu, uu, momentum(uu, lnrho, ss, aa), dt);
out_ss = rk3(out_ss, ss, entropy(ss, uu, lnrho, aa), dt); out_ss = rk3(out_ss, ss, entropy(ss, uu, lnrho, aa), dt);
#elif LTEMPERATURE #elif LTEMPERATURE
out_uu = rk3(out_uu, uu, momentum(uu, lnrho, tt), dt); out_uu = rk3(out_uu, uu, momentum(uu, lnrho, tt), dt);
out_tt = rk3(out_tt, tt, heat_transfer(uu, lnrho, tt), dt); out_tt = rk3(out_tt, tt, heat_transfer(uu, lnrho, tt), dt);
#else #else
out_uu = rk3(out_uu, uu, momentum(uu, lnrho), dt); out_uu = rk3(out_uu, uu, momentum(uu, lnrho), dt);
#endif #endif
#if LFORCING #if LFORCING
if (step_number == 2) { if (step_number == 2) {
out_uu = out_uu + forcing(globalVertexIdx, dt); out_uu = out_uu + forcing(globalVertexIdx, dt);
} }
#endif #endif
} }

2
acc/src/acc.l
View File

@@ -15,7 +15,7 @@ L [a-zA-Z_]
"void" { return VOID; } /* Rest of the types inherited from C */ "void" { return VOID; } /* Rest of the types inherited from C */
"int" { return INT; } "int" { return INT; }
"int3" { return INT3; } "int3" { return INT3; }
"ScalarField" { return SCALAR; } "ScalarField" { return SCALARFIELD; }
"VectorField" { return VECTOR; } "VectorField" { return VECTOR; }
"Kernel" { return KERNEL; } /* Function specifiers */ "Kernel" { return KERNEL; } /* Function specifiers */

3
acc/src/acc.y
View File

@@ -16,7 +16,7 @@ int yyget_lineno();
%token CONSTANT IN OUT UNIFORM %token CONSTANT IN OUT UNIFORM
%token IDENTIFIER NUMBER %token IDENTIFIER NUMBER
%token RETURN %token RETURN
%token SCALAR VECTOR MATRIX %token SCALAR VECTOR MATRIX SCALARFIELD
%token VOID INT INT3 %token VOID INT INT3
%token IF ELSE FOR WHILE ELIF %token IF ELSE FOR WHILE ELIF
%token LEQU LAND LOR LLEQU %token LEQU LAND LOR LLEQU
@@ -209,6 +209,7 @@ type_specifier: VOID
| SCALAR { $$ = astnode_create(NODE_TYPE_SPECIFIER, NULL, NULL); $$->token = SCALAR; } | SCALAR { $$ = astnode_create(NODE_TYPE_SPECIFIER, NULL, NULL); $$->token = SCALAR; }
| VECTOR { $$ = astnode_create(NODE_TYPE_SPECIFIER, NULL, NULL); $$->token = VECTOR; } | VECTOR { $$ = astnode_create(NODE_TYPE_SPECIFIER, NULL, NULL); $$->token = VECTOR; }
| MATRIX { $$ = astnode_create(NODE_TYPE_SPECIFIER, NULL, NULL); $$->token = MATRIX; } | 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, $$); } identifier: IDENTIFIER { $$ = astnode_create(NODE_IDENTIFIER, NULL, NULL); astnode_set_buffer(yytext, $$); }

90
acc/src/code_generator.c
View File

@@ -54,12 +54,13 @@ static const char* translation_table[TRANSLATION_TABLE_SIZE] = {
[WHILE] = "while", [WHILE] = "while",
[FOR] = "for", [FOR] = "for",
// Type specifiers // Type specifiers
[VOID] = "void", [VOID] = "void",
[INT] = "int", [INT] = "int",
[INT3] = "int3", [INT3] = "int3",
[SCALAR] = "AcReal", [SCALAR] = "AcReal",
[VECTOR] = "AcReal3", [VECTOR] = "AcReal3",
[MATRIX] = "AcMatrix", [MATRIX] = "AcMatrix",
[SCALARFIELD] = "AcReal",
// Type qualifiers // Type qualifiers
[KERNEL] = "template <int step_number> static " [KERNEL] = "template <int step_number> static "
"__global__", //__launch_bounds__(RK_THREADBLOCK_SIZE, "__global__", //__launch_bounds__(RK_THREADBLOCK_SIZE,
@@ -228,14 +229,18 @@ translate_latest_symbol(void)
} }
// UNIFORM // UNIFORM
else if (symbol->type_qualifier == 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 */ /* Do nothing */
} }
// IN / OUT // IN / OUT
else if (symbol->type != SYMBOLTYPE_FUNCTION_PARAMETER && else if (symbol->type != SYMBOLTYPE_FUNCTION_PARAMETER &&
(symbol->type_qualifier == IN || symbol->type_qualifier == OUT)) { (symbol->type_qualifier == IN || symbol->type_qualifier == OUT)) {
printf("static __device__ const %s %s%s", symbol->type_specifier == SCALAR ? "int" : "int3", printf("static __device__ const %s %s%s",
inout_name_prefix, symbol_table[handle].identifier); symbol->type_specifier == SCALARFIELD ? "int" : "int3", inout_name_prefix,
symbol_table[handle].identifier);
if (symbol->type_specifier == VECTOR) if (symbol->type_specifier == VECTOR)
printf(" = make_int3"); printf(" = make_int3");
} }
@@ -373,6 +378,8 @@ traverse(const ASTNode* node)
// printf("%s%s", inout_name_prefix, symbol->identifier); // printf("%s%s", inout_name_prefix, symbol->identifier);
//} //}
if (symbol->type_qualifier == UNIFORM) { if (symbol->type_qualifier == UNIFORM) {
printf("DCONST(%s) ", symbol->identifier);
/*
if (symbol->type_specifier == SCALAR) if (symbol->type_specifier == SCALAR)
printf("DCONST_REAL(AC_%s) ", symbol->identifier); printf("DCONST_REAL(AC_%s) ", symbol->identifier);
else if (symbol->type_specifier == INT) else if (symbol->type_specifier == INT)
@@ -380,6 +387,7 @@ traverse(const ASTNode* node)
else else
printf("INVALID UNIFORM type specifier %s with %s\n", printf("INVALID UNIFORM type specifier %s with %s\n",
translate(symbol->type_specifier), symbol->identifier); translate(symbol->type_specifier), symbol->identifier);
*/
} }
else { else {
// Do a regular translation // 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 int
main(int argc, char** argv) main(int argc, char** argv)
{ {
@@ -557,7 +627,7 @@ main(int argc, char** argv)
compilation_type = STENCIL_ASSEMBLY; compilation_type = STENCIL_ASSEMBLY;
else if (!strcmp(argv[1], "-sps")) else if (!strcmp(argv[1], "-sps"))
compilation_type = STENCIL_PROCESS; compilation_type = STENCIL_PROCESS;
else if (!strcmp(argv[1], "-hh")) else if (!strcmp(argv[1], "-sdh"))
compilation_type = STENCIL_HEADER; compilation_type = STENCIL_HEADER;
else else
printf("Unknown flag %s. Generating stencil assembly.\n", argv[1]); printf("Unknown flag %s. Generating stencil assembly.\n", argv[1]);
@@ -584,6 +654,8 @@ main(int argc, char** argv)
traverse(root); traverse(root);
if (compilation_type == STENCIL_ASSEMBLY) if (compilation_type == STENCIL_ASSEMBLY)
generate_preprocessed_structures(); generate_preprocessed_structures();
else if (compilation_type == STENCIL_HEADER)
generate_header();
// print_symbol_table(); // print_symbol_table();

2
include/astaroth_defines.h
View File

@@ -45,6 +45,8 @@ typedef struct {
#endif // __CUDACC__ #endif // __CUDACC__
// Library flags // Library flags
#define STENCIL_ORDER (6)
#define NGHOST (STENCIL_ORDER / 2)
#define VERBOSE_PRINTING (1) #define VERBOSE_PRINTING (1)
// Built-in types and parameters // Built-in types and parameters

25
scripts/compile_acc.sh
View File

@@ -8,17 +8,18 @@ fi
KERNEL_DIR=${AC_HOME}"/src/core/kernels" KERNEL_DIR=${AC_HOME}"/src/core/kernels"
ACC_DIR=${AC_HOME}"/acc" ACC_DIR=${AC_HOME}"/acc"
ACC_DEFAULT_HEADER="mhd_solver/stencil_defines.h"
ACC_DEFAULT_SAS="mhd_solver/stencil_assembly.sas" ACC_DEFAULT_SAS="mhd_solver/stencil_assembly.sas"
ACC_DEFAULT_SPS="mhd_solver/stencil_process.sps" 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}/clean.sh
${ACC_DIR}/build_acc.sh ${ACC_DIR}/build_acc.sh
ACC_HEADER=${ACC_DEFAULT_HEADER}
ACC_SAS=${ACC_DEFAULT_SAS} ACC_SAS=${ACC_DEFAULT_SAS}
ACC_SPS=${ACC_DEFAULT_SPS} ACC_SPS=${ACC_DEFAULT_SPS}
ACC_HEADER=${ACC_DEFAULT_HEADER}
ACC_INCLUDE_DIR=${ACC_DEFAULT_INCLUDE_DIR}
while [ "$#" -gt 0 ] while [ "$#" -gt 0 ]
do do
@@ -56,9 +57,17 @@ echo "Header file:" ${ACC_DIR}/${ACC_HEADER}
echo "Assembly file: ${ACC_DIR}/${ACC_SAS}" echo "Assembly file: ${ACC_DIR}/${ACC_SAS}"
echo "Process file: ${ACC_DIR}/${ACC_SPS}" echo "Process file: ${ACC_DIR}/${ACC_SPS}"
cd ${KERNEL_DIR} cd ${ACC_DIR}/${ACC_INCLUDE_DIR}
${ACC_DIR}/compile.sh ${ACC_DIR}/${ACC_SAS} ${ACC_DIR}/${ACC_HEADER} echo ${PWD}
${ACC_DIR}/compile.sh ${ACC_DIR}/${ACC_SPS} ${ACC_DIR}/${ACC_HEADER} ${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 echo "Moving stencil_assembly.cuh -> ${AC_HOME}/src/core/kernels"
ln -sf ${ACC_DIR}/${ACC_HEADER} ${AC_HOME}/include/stencil_defines.h 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

15
src/core/device.cu
View File

@@ -40,26 +40,31 @@ typedef struct {
} VertexBufferArray; } VertexBufferArray;
__constant__ AcMeshInfo d_mesh_info; __constant__ AcMeshInfo d_mesh_info;
static inline int __device__ static int __device__ __forceinline__
DCONST(const AcIntParam param) DCONST(const AcIntParam param)
{ {
return d_mesh_info.int_params[param]; return d_mesh_info.int_params[param];
} }
static inline int3 __device__ static int3 __device__ __forceinline__
DCONST(const AcInt3Param param) DCONST(const AcInt3Param param)
{ {
return d_mesh_info.int3_params[param]; return d_mesh_info.int3_params[param];
} }
static inline AcReal __device__ static AcReal __device__ __forceinline__
DCONST(const AcRealParam param) DCONST(const AcRealParam param)
{ {
return d_mesh_info.real_params[param]; return d_mesh_info.real_params[param];
} }
static inline AcReal3 __device__ static AcReal3 __device__ __forceinline__
DCONST(const AcReal3Param param) DCONST(const AcReal3Param param)
{ {
return d_mesh_info.real3_params[param]; return d_mesh_info.real3_params[param];
} }
constexpr VertexBufferHandle
DCONST(const VertexBufferHandle handle)
{
return handle;
}
#define DCONST_INT(x) DCONST(x) #define DCONST_INT(x) DCONST(x)
#define DCONST_INT3(x) DCONST(x) #define DCONST_INT3(x) DCONST(x)
#define DCONST_REAL(x) DCONST(x) #define DCONST_REAL(x) DCONST(x)
@@ -103,7 +108,7 @@ struct device_s {
}; };
// clang-format off // 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 // clang-format on
AcResult AcResult

4
src/mpitest/CMakeLists.txt
View File

@@ -8,5 +8,5 @@ set(CMAKE_C_STANDARD_REQUIRED ON)
find_package(MPI REQUIRED) find_package(MPI REQUIRED)
add_executable(mpitest main.c) add_executable(mpitest main.c)
target_include_directories(mpitest PRIVATE ${MPI_C_INCLUDE_PATH}) target_include_directories(mpitest PRIVATE ${CMAKE_SOURCE_DIR}/src/standalone ${MPI_C_INCLUDE_PATH})
target_link_libraries(mpitest PRIVATE ${MPI_C_LIBRARIES} astaroth_core) target_link_libraries(mpitest astaroth_core astaroth_standalone ${MPI_C_LIBRARIES})

148
src/mpitest/main.c
View File

@@ -16,13 +16,120 @@
You should have received a copy of the GNU General Public License You should have received a copy of the GNU General Public License
along with Astaroth. If not, see <http://www.gnu.org/licenses/>. along with Astaroth. If not, see <http://www.gnu.org/licenses/>.
*/ */
/**
Running: mpirun -np <num processes> <executable>
*/
#undef NDEBUG // Assert always
#include <assert.h>
#include <stdio.h> #include <stdio.h>
#include <stdlib.h> #include <stdlib.h>
#include <string.h>
#include "astaroth.h" #include "astaroth.h"
#include "autotest.h"
#include <mpi.h> #include <mpi.h>
// 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 int
main(void) main(void)
{ {
@@ -37,14 +144,39 @@ main(void)
MPI_Get_processor_name(processor_name, &name_len); MPI_Get_processor_name(processor_name, &name_len);
printf("Processor %s. Process %d of %d.\n", processor_name, process_id, num_processes); printf("Processor %s. Process %d of %d.\n", processor_name, process_id, num_processes);
AcMeshInfo info = { AcMeshInfo mesh_info;
.int_params[AC_nx] = 128, load_config(&mesh_info);
.int_params[AC_ny] = 64, update_config(&mesh_info);
.int_params[AC_nz] = 32,
}; AcMesh* main_mesh = NULL;
acInit(info); ModelMesh* model_mesh = NULL;
acIntegrate(0.1f); if (process_id == 0) {
acQuit(); 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(); MPI_Finalize();
return EXIT_SUCCESS; return EXIT_SUCCESS;

5
src/standalone/CMakeLists.txt
View File

@@ -25,10 +25,11 @@ add_compile_options(-pipe ${OpenMP_CXX_FLAGS})
add_compile_options(-Wall -Wextra -Werror -Wdouble-promotion -Wfloat-conversion)# -Wshadow) add_compile_options(-Wall -Wextra -Werror -Wdouble-promotion -Wfloat-conversion)# -Wshadow)
## Compile and link ## 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) 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 # Define the config directory
if (ALTER_CONF) if (ALTER_CONF)

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

@@ -31,6 +31,16 @@
#include "host_memory.h" #include "host_memory.h"
#include "model_boundconds.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 { typedef struct {
ModelScalar x, y, z; ModelScalar x, y, z;
} ModelVector; } ModelVector;