From 0208d55e4ebc50973ea7fc5b991edb1f46f935bf Mon Sep 17 00:00:00 2001 From: jpekkila Date: Mon, 19 Aug 2019 16:40:47 +0300 Subject: [PATCH 01/11] Moved STENCIL_ORDER and NGHOST out of user-defined parameter as these are actually internal defines used to configure the built-in functions. Additionally, renamed all explicitly declared uniforms from dsx -> AC_dsx in the DSL in preparation for having clear connection between DSL uniforms and the library parameter handles created by the user (AcRealParam etc) --- acc/mhd_solver/.gitignore | 5 ++ acc/mhd_solver/stencil_process.sps | 105 +++++++++++------------------ acc/src/code_generator.c | 6 ++ include/astaroth_defines.h | 2 + src/core/kernels/stencil_header.hh | 33 +++++++++ 5 files changed, 87 insertions(+), 64 deletions(-) create mode 100644 acc/mhd_solver/.gitignore create mode 100644 src/core/kernels/stencil_header.hh 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_process.sps b/acc/mhd_solver/stencil_process.sps index ee8656e..18141f7 100644 --- a/acc/mhd_solver/stencil_process.sps +++ b/acc/mhd_solver/stencil_process.sps @@ -1,27 +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_header.hh" Vector @@ -61,8 +38,8 @@ continuity(in VectorField uu, in ScalarField lnrho) { 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 + 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? const Scalar inv_rho = Scalar(1.) / exp(value(lnrho)); @@ -70,14 +47,14 @@ momentum(in VectorField uu, in ScalarField lnrho, in ScalarField ss, in VectorFi // 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)) + - cs2 * ((Scalar(1.) / AC_cp_sound) * gradient(ss) + gradient(lnrho)) + inv_rho * cross(j, B) - + nu_visc * ( + + AC_nu_visc * ( laplace_vec(uu) + Scalar(1. / 3.) * gradient_of_divergence(uu) + Scalar(2.) * mul(S, gradient(lnrho)) ) - + zeta * gradient_of_divergence(uu); + + AC_zeta * gradient_of_divergence(uu); return mom; } #elif LTEMPERATURE @@ -87,13 +64,13 @@ momentum(in VectorField uu, in ScalarField lnrho, in ScalarField tt) { 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 * + AC_nu_visc * (laplace_vec(uu) + Scalar(1. / 3.) * gradient_of_divergence(uu) + - Scalar(2.) * mul(S, gradient(lnrho))) + zeta * gradient_of_divergence(uu); + Scalar(2.) * mul(S, gradient(lnrho))) + AC_zeta * gradient_of_divergence(uu); #if LGRAVITY mom = mom - (Vector){0, 0, -10.0}; @@ -111,10 +88,10 @@ momentum(in VectorField uu, in ScalarField lnrho) { // Isothermal: we have constant speed of sound mom = -mul(gradients(uu), value(uu)) - - cs2_sound * gradient(lnrho) + - nu_visc * + AC_cs2_sound * gradient(lnrho) + + AC_nu_visc * (laplace_vec(uu) + Scalar(1. / 3.) * gradient_of_divergence(uu) + - Scalar(2.) * mul(S, gradient(lnrho))) + zeta * gradient_of_divergence(uu); + Scalar(2.) * mul(S, gradient(lnrho))) + AC_zeta * gradient_of_divergence(uu); #if LGRAVITY mom = mom - (Vector){0, 0, -10.0}; @@ -130,13 +107,13 @@ 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 ]) + // 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; } @@ -145,27 +122,27 @@ induction(in VectorField uu, in VectorField aa) { #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); + 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; + const Scalar inv_AC_cp_sound = AcReal(1.) / AC_cp_sound; 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) + + 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 @@ -177,11 +154,11 @@ Scalar 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 Vector j = (Scalar(1.) / AC_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); + + 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 @@ -195,7 +172,7 @@ 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); + 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 @@ -220,7 +197,7 @@ 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,9 +206,9 @@ 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); @@ -254,13 +231,13 @@ 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; + Vector a = Scalar(.5) * (Vector){globalGridN.x * AC_dsx, + globalGridN.y * AC_dsy, + globalGridN.z * AC_dsz}; // source (origin) + Vector xx = (Vector){(globalVertexIdx.x - AC_nx_min) * AC_dsx, + (globalVertexIdx.y - AC_ny_min) * AC_dsy, + (globalVertexIdx.z - AC_nz_min) * AC_dsz}; // sink (current index) + const Scalar cs2 = AC_cs2_sound; const Scalar cs = sqrt(cs2); //Placeholders until determined properly diff --git a/acc/src/code_generator.c b/acc/src/code_generator.c index 2de30d7..01cf1d1 100644 --- a/acc/src/code_generator.c +++ b/acc/src/code_generator.c @@ -228,6 +228,9 @@ 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 @@ -373,6 +376,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 +385,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 diff --git a/include/astaroth_defines.h b/include/astaroth_defines.h index ee5d70c..96ab3fd 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/src/core/kernels/stencil_header.hh b/src/core/kernels/stencil_header.hh new file mode 100644 index 0000000..798e3cc --- /dev/null +++ b/src/core/kernels/stencil_header.hh @@ -0,0 +1,33 @@ +#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 + +// Declare uniforms (i.e. device constants) +uniform Scalar AC_cs2_sound; +uniform Scalar AC_nu_visc; +uniform Scalar AC_cp_sound; +uniform Scalar AC_cv_sound; +uniform Scalar AC_mu0; +uniform Scalar AC_eta; +uniform Scalar AC_gamma; +uniform Scalar AC_zeta; + +uniform Scalar AC_dsx; +uniform Scalar AC_dsy; +uniform Scalar AC_dsz; + +uniform Scalar AC_lnT0; +uniform Scalar AC_lnrho0; + +uniform int AC_nx_min; +uniform int AC_ny_min; +uniform int AC_nz_min; +uniform int AC_nx; +uniform int AC_ny; +uniform int AC_nz; From bcdd827a4f7de3d77ecc8798d81f28b208529658 Mon Sep 17 00:00:00 2001 From: jpekkila Date: Mon, 19 Aug 2019 17:05:56 +0300 Subject: [PATCH 02/11] Added a proper declarations for all user-specified uniform. Note: built-in uniforms are not correctly translated into CUDA --- acc/mhd_solver/stencil_process.sps | 18 ++--- src/core/kernels/stencil_header.hh | 122 +++++++++++++++++++++++++---- 2 files changed, 115 insertions(+), 25 deletions(-) diff --git a/acc/mhd_solver/stencil_process.sps b/acc/mhd_solver/stencil_process.sps index 18141f7..a2219b6 100644 --- a/acc/mhd_solver/stencil_process.sps +++ b/acc/mhd_solver/stencil_process.sps @@ -234,18 +234,18 @@ forcing(int3 globalVertexIdx, Scalar dt) Vector a = Scalar(.5) * (Vector){globalGridN.x * AC_dsx, globalGridN.y * AC_dsy, globalGridN.z * AC_dsz}; // source (origin) - Vector xx = (Vector){(globalVertexIdx.x - AC_nx_min) * AC_dsx, - (globalVertexIdx.y - AC_ny_min) * AC_dsy, - (globalVertexIdx.z - AC_nz_min) * AC_dsz}; // sink (current index) + 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)}; + 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. @@ -254,7 +254,7 @@ forcing(int3 globalVertexIdx, Scalar dt) 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(DCONST_REAL(AC_kaver)*cs); + 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; diff --git a/src/core/kernels/stencil_header.hh b/src/core/kernels/stencil_header.hh index 798e3cc..f7de1ba 100644 --- a/src/core/kernels/stencil_header.hh +++ b/src/core/kernels/stencil_header.hh @@ -8,26 +8,116 @@ #define AC_THERMAL_CONDUCTIVITY (AcReal(0.001)) // TODO: make an actual config parameter -// Declare uniforms (i.e. device constants) -uniform Scalar AC_cs2_sound; -uniform Scalar AC_nu_visc; -uniform Scalar AC_cp_sound; -uniform Scalar AC_cv_sound; -uniform Scalar AC_mu0; -uniform Scalar AC_eta; -uniform Scalar AC_gamma; -uniform Scalar AC_zeta; +// 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; -uniform int AC_nx_min; -uniform int AC_ny_min; -uniform int AC_nz_min; -uniform int AC_nx; -uniform int AC_ny; -uniform int AC_nz; +/* + * ============================================================================= + * 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 From d801ebdd41e5c893f96f5513f65f25df4194cbb0 Mon Sep 17 00:00:00 2001 From: jpekkila Date: Mon, 19 Aug 2019 17:35:03 +0300 Subject: [PATCH 03/11] Now parameters and vertexbuffers (fields) can be declared with the DSL only. TODO: translation from the DSL header to C --- src/core/device.cu | 5 ++++ src/core/kernels/stencil_header.hh | 46 +++++++++++++----------------- 2 files changed, 25 insertions(+), 26 deletions(-) diff --git a/src/core/device.cu b/src/core/device.cu index 8b89a78..a3db98a 100644 --- a/src/core/device.cu +++ b/src/core/device.cu @@ -60,6 +60,11 @@ 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) diff --git a/src/core/kernels/stencil_header.hh b/src/core/kernels/stencil_header.hh index f7de1ba..31667fc 100644 --- a/src/core/kernels/stencil_header.hh +++ b/src/core/kernels/stencil_header.hh @@ -90,34 +90,28 @@ uniform Scalar AC_inv_dsz; * 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), +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 -#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), +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 -#define AC_FOR_VTXBUF_HANDLES(FUNC) \ - FUNC(VTXBUF_LNRHO), \ - FUNC(VTXBUF_UUX), \ - FUNC(VTXBUF_UUY), \ - FUNC(VTXBUF_UUZ), +uniform ScalarField VTXBUF_LNRHO; +uniform ScalarField VTXBUF_UUX; +uniform ScalarField VTXBUF_UUY; +uniform ScalarField VTXBUF_UUZ; #else -#define AC_FOR_VTXBUF_HANDLES(FUNC) \ - FUNC(VTXBUF_LNRHO), +uniform ScalarField VTXBUF_LNRHO; #endif -// clang-format on From 51cf1f1068466d865e69228d2602a5099c6d6734 Mon Sep 17 00:00:00 2001 From: jpekkila Date: Mon, 19 Aug 2019 18:19:28 +0300 Subject: [PATCH 04/11] The C header is now generated from the DSL, stashing the changes just to be sure since I might overwrite something when updating the compilation scripts to work with this new scheme --- acc/compile.sh | 12 +- acc/mhd_solver/stencil_defines.h | 163 ------------------ .../mhd_solver/stencil_definition.sdh | 0 acc/mhd_solver/stencil_process.sps | 2 +- acc/src/acc.l | 2 +- acc/src/acc.y | 3 +- acc/src/code_generator.c | 79 ++++++++- config/astaroth.conf | 2 +- scripts/compile_acc.sh | 22 ++- 9 files changed, 99 insertions(+), 186 deletions(-) delete mode 100644 acc/mhd_solver/stencil_defines.h rename src/core/kernels/stencil_header.hh => acc/mhd_solver/stencil_definition.sdh (100%) 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/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/src/core/kernels/stencil_header.hh b/acc/mhd_solver/stencil_definition.sdh similarity index 100% rename from src/core/kernels/stencil_header.hh rename to acc/mhd_solver/stencil_definition.sdh diff --git a/acc/mhd_solver/stencil_process.sps b/acc/mhd_solver/stencil_process.sps index a2219b6..420193e 100644 --- a/acc/mhd_solver/stencil_process.sps +++ b/acc/mhd_solver/stencil_process.sps @@ -1,4 +1,4 @@ -#include "stencil_header.hh" +#include "stencil_definition.sdh" Vector 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 01cf1d1..c31ad86 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, @@ -555,6 +556,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) { @@ -563,7 +626,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]); @@ -590,6 +653,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/config/astaroth.conf b/config/astaroth.conf index 32f50a3..944364b 100644 --- a/config/astaroth.conf +++ b/config/astaroth.conf @@ -7,7 +7,7 @@ */ AC_nx = 128 AC_ny = 128 -AC_nz = 128 +AC_nz = 6 AC_dsx = 0.04908738521 AC_dsy = 0.04908738521 diff --git a/scripts/compile_acc.sh b/scripts/compile_acc.sh index 8648079..4649bf5 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=${ACC_DEFAULT_INCLUDE_DIR} while [ "$#" -gt 0 ] do @@ -56,9 +57,14 @@ 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} +${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 +#mv ${ACC_SAS} ${AC_HOME}/src/core/kernels +#mv ${ACC_SPS} ${AC_HOME}/src/core/kernels +#mv ${ACC_HEADER} ${AC_HOME}/include + +#echo "Linking: " ${ACC_DIR}/${ACC_HEADER} " -> " ${AC_HOME}/include/stencil_defines.h +#ln -sf ${ACC_DIR}/${ACC_HEADER} ${AC_HOME}/include/stencil_defines.h From 5b7408eb55d90f11998ff9d8b9cbe04dec228241 Mon Sep 17 00:00:00 2001 From: jpekkila Date: Mon, 19 Aug 2019 18:43:16 +0300 Subject: [PATCH 05/11] User config param overhaul complete, works. If I haven't missed anything, all fields and user parameters, and everything related to simulation can now be declared with the DSL. The only thing that you need to do is to fill the declared symbols with data, like with OpenGL and GLSL. --- acc/mhd_solver/stencil_assembly.sas | 8 +++++--- acc/src/code_generator.c | 5 +++-- scripts/compile_acc.sh | 15 +++++++++------ src/standalone/model/model_rk3.cc | 10 ++++++++++ 4 files changed, 27 insertions(+), 11 deletions(-) diff --git a/acc/mhd_solver/stencil_assembly.sas b/acc/mhd_solver/stencil_assembly.sas index 5c4f14a..e30d500 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) { @@ -17,7 +19,7 @@ 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] @@ -32,7 +34,7 @@ der6x_upwd(in ScalarField vertex) 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] @@ -47,7 +49,7 @@ der6y_upwd(in ScalarField vertex) 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] diff --git a/acc/src/code_generator.c b/acc/src/code_generator.c index c31ad86..e560a70 100644 --- a/acc/src/code_generator.c +++ b/acc/src/code_generator.c @@ -238,8 +238,9 @@ translate_latest_symbol(void) 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"); } diff --git a/scripts/compile_acc.sh b/scripts/compile_acc.sh index 4649bf5..7ebbd4c 100755 --- a/scripts/compile_acc.sh +++ b/scripts/compile_acc.sh @@ -19,7 +19,7 @@ ${ACC_DIR}/build_acc.sh ACC_SAS=${ACC_DEFAULT_SAS} ACC_SPS=${ACC_DEFAULT_SPS} ACC_HEADER=${ACC_DEFAULT_HEADER} -ACC_INCLUDE=${ACC_DEFAULT_INCLUDE_DIR} +ACC_INCLUDE_DIR=${ACC_DEFAULT_INCLUDE_DIR} while [ "$#" -gt 0 ] do @@ -58,13 +58,16 @@ echo "Assembly file: ${ACC_DIR}/${ACC_SAS}" echo "Process file: ${ACC_DIR}/${ACC_SPS}" 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} -#mv ${ACC_SAS} ${AC_HOME}/src/core/kernels -#mv ${ACC_SPS} ${AC_HOME}/src/core/kernels -#mv ${ACC_HEADER} ${AC_HOME}/include +echo "Moving stencil_assembly.cuh -> ${AC_HOME}/src/core/kernels" +mv stencil_assembly.cuh ${AC_HOME}/src/core/kernels -#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_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/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; From 73d393e4190d8802fe03295cfd0adad655bf3f92 Mon Sep 17 00:00:00 2001 From: jpekkila Date: Tue, 20 Aug 2019 18:40:38 +0300 Subject: [PATCH 06/11] Changed order for linking the MPI library to work around cmake error on CMP0004 --- src/mpitest/CMakeLists.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/mpitest/CMakeLists.txt b/src/mpitest/CMakeLists.txt index c64105d..ed19b4d 100644 --- a/src/mpitest/CMakeLists.txt +++ b/src/mpitest/CMakeLists.txt @@ -9,4 +9,4 @@ 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_link_libraries(mpitest PRIVATE astaroth_core ${MPI_C_LIBRARIES}) From 5d2b658fb0680531e885fbb42ee8786af3f06c99 Mon Sep 17 00:00:00 2001 From: jpekkila Date: Tue, 20 Aug 2019 18:41:26 +0300 Subject: [PATCH 07/11] Autoformatted the DSL files --- acc/mhd_solver/stencil_assembly.sas | 59 +++--- acc/mhd_solver/stencil_process.sps | 293 ++++++++++++++-------------- 2 files changed, 178 insertions(+), 174 deletions(-) diff --git a/acc/mhd_solver/stencil_assembly.sas b/acc/mhd_solver/stencil_assembly.sas index e30d500..85886e6 100644 --- a/acc/mhd_solver/stencil_assembly.sas +++ b/acc/mhd_solver/stencil_assembly.sas @@ -9,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 @@ -21,14 +19,14 @@ der6x_upwd(in ScalarField vertex) { 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 @@ -36,14 +34,14 @@ der6y_upwd(in ScalarField vertex) { 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 @@ -51,14 +49,14 @@ der6z_upwd(in ScalarField vertex) { 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 @@ -68,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_process.sps b/acc/mhd_solver/stencil_process.sps index 420193e..e605b9f 100644 --- a/acc/mhd_solver/stencil_process.sps +++ b/acc/mhd_solver/stencil_process.sps @@ -1,6 +1,5 @@ #include "stencil_definition.sdh" - Vector value(in VectorField uu) { @@ -14,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 @@ -25,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); @@ -36,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 = 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 +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.) / 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); + 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 = (AC_cp_sound - AC_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 + - 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); + 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)) - - 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); + 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 - 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); +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, AC_mu0 is cancelled out - const Vector ind = cross(value(uu), B) - AC_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 = AC_lnT0 + AC_gamma * value(ss) / AC_cp_sound + - (AC_gamma - Scalar(1.)) * (value(lnrho) - AC_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_AC_cp_sound = AcReal(1.) / AC_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 = 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 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)) * AC_cp_sound); - return AC_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.) / 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); + 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 @@ -170,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)) + AC_nu_visc * contract(S) * (Scalar(1.) / AC_cv_sound) - (AC_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) @@ -191,8 +195,8 @@ 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) { @@ -206,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/(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)); + 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; } @@ -231,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 * AC_dsx, - globalGridN.y * AC_dsy, + 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) + 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); + const Scalar cs = sqrt(cs2); - //Placeholders until determined properly + // 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 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(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 @@ -271,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 @@ -290,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 } From 39dcda4a040f33797aeff2f28a62a3eff9ef4096 Mon Sep 17 00:00:00 2001 From: jpekkila Date: Wed, 21 Aug 2019 14:28:46 +0300 Subject: [PATCH 08/11] Made warnings about unused functions go away (this is intended functionality and not all programs will use all types of device constants, thus unnecessary warning) --- src/core/device.cu | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/core/device.cu b/src/core/device.cu index a3db98a..15f0c1e 100644 --- a/src/core/device.cu +++ b/src/core/device.cu @@ -40,22 +40,22 @@ 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]; @@ -108,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 From 5867ff4b3e47a99863969a0a52264bd3c6f8c06d Mon Sep 17 00:00:00 2001 From: jpekkila Date: Wed, 21 Aug 2019 16:16:12 +0300 Subject: [PATCH 09/11] Stashing MPItest changes --- src/mpitest/CMakeLists.txt | 4 +- src/mpitest/main.c | 148 +++++++++++++++++++++++++++++++++++-- 2 files changed, 142 insertions(+), 10 deletions(-) diff --git a/src/mpitest/CMakeLists.txt b/src/mpitest/CMakeLists.txt index ed19b4d..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 astaroth_core ${MPI_C_LIBRARIES}) +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; From d52e002e5d1756317927ff10d9774ce25bce561e Mon Sep 17 00:00:00 2001 From: jpekkila Date: Wed, 21 Aug 2019 16:18:48 +0300 Subject: [PATCH 10/11] Made Astaroth Standalone a library component (still works as before but can be included in other projects which need f.ex. autotesting) --- src/standalone/CMakeLists.txt | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) 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) From 9fa39b66fc142da695caa33562e1f6972c9456a1 Mon Sep 17 00:00:00 2001 From: jpekkila Date: Mon, 26 Aug 2019 18:33:50 +0300 Subject: [PATCH 11/11] Reverted an accidentally modified astaroth.conf to the same version as in master --- config/astaroth.conf | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/config/astaroth.conf b/config/astaroth.conf index 944364b..32f50a3 100644 --- a/config/astaroth.conf +++ b/config/astaroth.conf @@ -7,7 +7,7 @@ */ AC_nx = 128 AC_ny = 128 -AC_nz = 6 +AC_nz = 128 AC_dsx = 0.04908738521 AC_dsy = 0.04908738521