Merge branch 'master' into sink_20190723

Conflicts:
	acc/mhd_solver/stencil_process.sps

I've mannaully resolved the conflict, only that I'm leaving int3 globalVertexIdx as is, as instructed by Miikka.

acc/mhd_solver/stencil_assembly.sas

@@ -1,11 +1,11 @@
 Preprocessed Scalar
-value(in Scalar vertex)
+value(in ScalarField vertex)
 {
     return vertex[vertexIdx];
 }
 
 Preprocessed Vector
-gradient(in Scalar vertex)
+gradient(in ScalarField vertex)
 {
     return (Vector){derx(vertexIdx, vertex),
                     dery(vertexIdx, vertex),
@@ -15,7 +15,7 @@ gradient(in Scalar vertex)
 #if LUPWD
 
 Preprocessed Scalar
-der6x_upwd(in Scalar vertex)
+der6x_upwd(in ScalarField vertex)
 {
     Scalar inv_ds = DCONST_REAL(AC_inv_dsx);
 
@@ -30,7 +30,7 @@ der6x_upwd(in Scalar vertex)
 }
 
 Preprocessed Scalar
-der6y_upwd(in Scalar vertex)
+der6y_upwd(in ScalarField vertex)
 {
     Scalar inv_ds = DCONST_REAL(AC_inv_dsy);
 
@@ -45,7 +45,7 @@ der6y_upwd(in Scalar vertex)
 }
 
 Preprocessed Scalar
-der6z_upwd(in Scalar vertex)
+der6z_upwd(in ScalarField vertex)
 {
     Scalar inv_ds = DCONST_REAL(AC_inv_dsz);
 
@@ -62,7 +62,7 @@ der6z_upwd(in Scalar vertex)
 #endif
 
 Preprocessed Matrix
-hessian(in Scalar vertex)
+hessian(in ScalarField vertex)
 {
     Matrix hessian;
 

acc/mhd_solver/stencil_process.sps

@@ -23,14 +23,14 @@ uniform int ny;
 uniform int nz;
 
 Vector
-value(in Vector uu)
+value(in VectorField uu)
 {
     return (Vector){value(uu.x), value(uu.y), value(uu.z)};
 }
 
 #if LUPWD
 Scalar
-upwd_der6(in Vector uu, in Scalar lnrho)
+upwd_der6(in VectorField uu, in ScalarField lnrho)
 {
     Scalar uux = fabs(value(uu).x);
     Scalar uuy = fabs(value(uu).y);
@@ -40,13 +40,13 @@ upwd_der6(in Vector uu, in Scalar lnrho)
 #endif
 
 Matrix
-gradients(in Vector uu)
+gradients(in VectorField uu)
 {
     return (Matrix){gradient(uu.x), gradient(uu.y), gradient(uu.z)};
 }
 
 Scalar
-continuity(in Vector uu, in Scalar lnrho) {
+continuity(in VectorField uu, in ScalarField lnrho) {
     return -dot(value(uu), gradient(lnrho))
 #if LUPWD
            //This is a corrective hyperdiffusion term for upwinding.
@@ -124,7 +124,7 @@ accretion_profile(int3 globalVertexIdx, in Scalar lnrho){
 
 #if LENTROPY
 Vector
-momentum(int3 globalVertexIdx, in Vector uu, in Scalar lnrho, in Scalar ss, in Vector aa) {
+momentum(int3 globalVertexIdx, 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
@@ -152,13 +152,10 @@ momentum(int3 globalVertexIdx, in Vector uu, in Scalar lnrho, in Scalar ss, in V
 }
 #elif LTEMPERATURE
 Vector
-momentum(int3 globalVertexIdx, in Vector uu, in Scalar lnrho, in Scalar tt) {
+momentum(int3 globalVertexIdx, in VectorField uu, in ScalarField lnrho, in ScalarField tt) {
     Vector mom;
-
     const Matrix S = stress_tensor(uu);
-
     const Vector pressure_term = (cp_sound - cv_sound) * (gradient(tt) + value(tt) * gradient(lnrho));
-
     mom = -mul(gradients(uu), value(uu)) -
     pressure_term +
     nu_visc *
@@ -175,9 +172,8 @@ momentum(int3 globalVertexIdx, in Vector uu, in Scalar lnrho, in Scalar tt) {
 }
 #else
 Vector
-momentum(int3 globalVertexIdx, in Vector uu, in Scalar lnrho) {
+momentum(int3 globalVertexIdx, in VectorField uu, in ScalarField lnrho) {
     Vector mom;
-
     const Matrix S = stress_tensor(uu);
 
     // Isothermal: we have constant speed of sound
@@ -200,7 +196,7 @@ momentum(int3 globalVertexIdx, in Vector uu, in Scalar lnrho) {
 
 
 Vector
-induction(in Vector uu, in Vector 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)
@@ -218,7 +214,7 @@ induction(in Vector uu, in Vector aa) {
 
 #if LENTROPY
 Scalar
-lnT( in Scalar ss, in Scalar lnrho) {
+lnT( in ScalarField ss, in ScalarField lnrho) {
   const Scalar lnT = lnT0 + gamma * value(ss) / cp_sound +
     (gamma - Scalar(1.)) * (value(lnrho) - lnrho0);
   return lnT;
@@ -226,7 +222,7 @@ lnT( in Scalar ss, in Scalar lnrho) {
 
 // Nabla dot (K nabla T) / (rho T)
 Scalar
-heat_conduction( in Scalar ss, in Scalar lnrho) {
+heat_conduction( in ScalarField ss, in ScalarField lnrho) {
   const Scalar inv_cp_sound = AcReal(1.) / cp_sound;
 
   const Vector grad_ln_chi = - gradient(lnrho);
@@ -249,7 +245,7 @@ heating(const int i, const int j, const int k) {
 }
 
 Scalar
-entropy(in Scalar ss, in Vector uu, in Scalar lnrho, in Vector aa) {
+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
@@ -266,7 +262,7 @@ entropy(in Scalar ss, in Vector uu, in Scalar lnrho, in Vector aa) {
 
 #if LTEMPERATURE
 Scalar
-heat_transfer(in Vector uu, in Scalar lnrho, in Scalar tt)
+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;
@@ -365,26 +361,26 @@ forcing(int3 globalVertexIdx, Scalar dt)
 
 // Declare input and output arrays using locations specified in the
 // array enum in astaroth.h
-in Scalar lnrho = VTXBUF_LNRHO;
+in ScalarField lnrho(VTXBUF_LNRHO);
-out Scalar out_lnrho = VTXBUF_LNRHO;
+out ScalarField out_lnrho(VTXBUF_LNRHO);
 
-in Vector uu = (int3) {VTXBUF_UUX, VTXBUF_UUY, VTXBUF_UUZ};
+in VectorField uu(VTXBUF_UUX, VTXBUF_UUY, VTXBUF_UUZ);
-out Vector out_uu = (int3) {VTXBUF_UUX,VTXBUF_UUY,VTXBUF_UUZ};
+out VectorField out_uu(VTXBUF_UUX,VTXBUF_UUY,VTXBUF_UUZ);
 
 
 #if LMAGNETIC
-in Vector aa = (int3) {VTXBUF_AX,VTXBUF_AY,VTXBUF_AZ};
+in VectorField aa(VTXBUF_AX,VTXBUF_AY,VTXBUF_AZ);
-out Vector out_aa = (int3) {VTXBUF_AX,VTXBUF_AY,VTXBUF_AZ};
+out VectorField out_aa(VTXBUF_AX,VTXBUF_AY,VTXBUF_AZ);
 #endif
 
 #if LENTROPY
-in Scalar ss = VTXBUF_ENTROPY;
+in ScalarField ss(VTXBUF_ENTROPY);
-out Scalar out_ss = VTXBUF_ENTROPY;
+out ScalarField out_ss(VTXBUF_ENTROPY);
 #endif
 
 #if LTEMPERATURE
-in Scalar tt = VTXBUF_TEMPERATURE;
+in ScalarField tt(VTXBUF_TEMPERATURE);
-out Scalar out_tt = VTXBUF_TEMPERATURE;
+out ScalarField out_tt(VTXBUF_TEMPERATURE);
 #endif
 
 #if LSINK

acc/pseudodisk/stencil_process_gravx.sps

@@ -43,19 +43,19 @@ distance_x(Vector a, Vector b)
 }
 
 Vector
-value(in Vector uu)
+value(in VectorField uu)
 {
     return (Vector){value(uu.x), value(uu.y), value(uu.z)};
 }
 
 Matrix
-gradients(in Vector uu)
+gradients(in VectorField uu)
 {
     return (Matrix){gradient(uu.x), gradient(uu.y), gradient(uu.z)};
 }
 
 Scalar
-continuity(in Vector uu, in Scalar lnrho) {
+continuity(in VectorField uu, in ScalarField lnrho) {
     return -dot(value(uu), gradient(lnrho)) - divergence(uu);
 }
 
@@ -79,7 +79,7 @@ grav_force_line(const int3 vertexIdx)
 
 #if LENTROPY
 Vector
-momentum(in Vector uu, in Scalar lnrho, in Scalar ss, in Vector aa, const int3 vertexIdx) {
+momentum(in VectorField uu, in ScalarField lnrho, in ScalarField ss, in VectorField aa, const int3 vertexIdx) {
   Vector mom;
 
   const Matrix S = stress_tensor(uu);
@@ -104,7 +104,7 @@ momentum(in Vector uu, in Scalar lnrho, in Scalar ss, in Vector aa, const int3 v
 }
 #else
 Vector
-momentum(in Vector uu, in Scalar lnrho, const int3 vertexIdx) {
+momentum(in VectorField uu, in ScalarField lnrho, const int3 vertexIdx) {
   Vector mom;
 
   const Matrix S = stress_tensor(uu);
@@ -123,7 +123,7 @@ momentum(in Vector uu, in Scalar lnrho, const int3 vertexIdx) {
 
 
 Vector
-induction(in Vector uu, in Vector 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)
@@ -141,7 +141,7 @@ induction(in Vector uu, in Vector aa) {
 
 #if LENTROPY
 Scalar
-lnT( in Scalar ss, in Scalar lnrho) {
+lnT( in ScalarField ss, in ScalarField lnrho) {
   const Scalar lnT = LNT0 + value(ss) / cp_sound +
     (gamma - AcReal(1.)) * (value(lnrho) - LNRHO0);
   return lnT;
@@ -149,7 +149,7 @@ lnT( in Scalar ss, in Scalar lnrho) {
 
 // Nabla dot (K nabla T) / (rho T)
 Scalar
-heat_conduction( in Scalar ss, in Scalar lnrho) {
+heat_conduction( in ScalarField ss, in ScalarField lnrho) {
   const Scalar inv_cp_sound = AcReal(1.) / cp_sound;
 
   const Vector grad_ln_chi = (Vector) {
@@ -174,7 +174,7 @@ heating(const int i, const int j, const int k) {
 }
 
 Scalar
-entropy(in Scalar ss, in Vector uu, in Scalar lnrho, in Vector aa) {
+entropy(in ScalarField ss, in VectorField uu, in ScalarField lnrho, in VectorField aa) {
     const Matrix S = stress_tensor(uu);
 
     // nabla x nabla x A / mu0 = nabla(nabla dot A) - nabla^2(A)
@@ -193,21 +193,21 @@ entropy(in Scalar ss, in Vector uu, in Scalar lnrho, in Vector aa) {
 
 // Declare input and output arrays using locations specified in the
 // array enum in astaroth.h
-in Scalar lnrho = VTXBUF_LNRHO;
+in ScalarField lnrho(VTXBUF_LNRHO);
-out Scalar out_lnrho = VTXBUF_LNRHO;
+out ScalarField out_lnrho(VTXBUF_LNRHO);
 
-in Vector uu = (int3) {VTXBUF_UUX, VTXBUF_UUY, VTXBUF_UUZ};
+in VectorField uu(VTXBUF_UUX, VTXBUF_UUY, VTXBUF_UUZ);
-out Vector out_uu = (int3) {VTXBUF_UUX,VTXBUF_UUY,VTXBUF_UUZ};
+out VectorField out_uu(VTXBUF_UUX,VTXBUF_UUY,VTXBUF_UUZ);
 
 
 #if LMAGNETIC
-in Vector aa = (int3) {VTXBUF_AX,VTXBUF_AY,VTXBUF_AZ};
+in VectorField aa(VTXBUF_AX,VTXBUF_AY,VTXBUF_AZ);
-out Vector out_aa = (int3) {VTXBUF_AX,VTXBUF_AY,VTXBUF_AZ};
+out VectorField out_aa(VTXBUF_AX,VTXBUF_AY,VTXBUF_AZ);
 #endif
 
 #if LENTROPY
-in Scalar ss = VTXBUF_ENTROPY;
+in ScalarField ss(VTXBUF_ENTROPY);
-out Scalar out_ss = VTXBUF_ENTROPY;
+out ScalarField out_ss(VTXBUF_ENTROPY);
 #endif
 
 Kernel void

acc/pseudodisk/stencil_process_isotherm_gravx.sps

@@ -40,19 +40,19 @@ distance_x(Vector a, Vector b)
 }
 
 Vector
-value(in Vector uu)
+value(in VectorField uu)
 {
     return (Vector){value(uu.x), value(uu.y), value(uu.z)};
 }
 
 Matrix
-gradients(in Vector uu)
+gradients(in VectorField uu)
 {
     return (Matrix){gradient(uu.x), gradient(uu.y), gradient(uu.z)};
 }
 
 Scalar
-continuity(in Vector uu, in Scalar lnrho) {
+continuity(in VectorField uu, in ScalarField lnrho) {
     return -dot(value(uu), gradient(lnrho)) - divergence(uu);
 }
 
@@ -77,7 +77,7 @@ grav_force_line(const int3 vertexIdx)
 
 
 Vector
-momentum(in Vector uu, in Scalar lnrho, const int3 vertexIdx) {
+momentum(in VectorField uu, in ScalarField lnrho, const int3 vertexIdx) {
   Vector mom;
 
   const Matrix S = stress_tensor(uu);
@@ -94,7 +94,7 @@ momentum(in Vector uu, in Scalar lnrho, const int3 vertexIdx) {
 }
 
 Vector
-induction(in Vector uu, in Vector 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)
@@ -111,15 +111,16 @@ induction(in Vector uu, in Vector aa) {
 
 // Declare input and output arrays using locations specified in the
 // array enum in astaroth.h
-in Scalar lnrho = VTXBUF_LNRHO;
+in ScalarField lnrho(VTXBUF_LNRHO);
-out Scalar out_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);
 
-in Vector uu = (int3) {VTXBUF_UUX, VTXBUF_UUY, VTXBUF_UUZ};
-out Vector out_uu = (int3) {VTXBUF_UUX,VTXBUF_UUY,VTXBUF_UUZ};
 
 #if LMAGNETIC
-in Vector aa = (int3) {VTXBUF_AX,VTXBUF_AY,VTXBUF_AZ};
+in VectorField aa(VTXBUF_AX,VTXBUF_AY,VTXBUF_AZ);
-out Vector out_aa = (int3) {VTXBUF_AX,VTXBUF_AY,VTXBUF_AZ};
+out VectorField out_aa(VTXBUF_AX,VTXBUF_AY,VTXBUF_AZ);
 #endif
 
 Kernel void
@@ -132,38 +133,3 @@ solve(Scalar dt) {
 
   WRITE(out_uu, RK3(out_uu, uu, momentum(uu, lnrho, vertexIdx), dt));
 }
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-

acc/pseudodisk/stencil_process_isotherm_linegrav.sps

@@ -40,19 +40,19 @@ distance(Vector a, Vector b)
 }
 
 Vector
-value(in Vector uu)
+value(in VectorField uu)
 {
     return (Vector){value(uu.x), value(uu.y), value(uu.z)};
 }
 
 Matrix
-gradients(in Vector uu)
+gradients(in VectorField uu)
 {
     return (Matrix){gradient(uu.x), gradient(uu.y), gradient(uu.z)};
 }
 
 Scalar
-continuity(in Vector uu, in Scalar lnrho) {
+continuity(in VectorField uu, in ScalarField lnrho) {
     return -dot(value(uu), gradient(lnrho)) - divergence(uu);
 }
 
@@ -82,7 +82,7 @@ grav_force_line(const int3 vertexIdx)
 
 
 Vector
-momentum(in Vector uu, in Scalar lnrho, const int3 vertexIdx) {
+momentum(in VectorField uu, in ScalarField lnrho, const int3 vertexIdx) {
   Vector mom;
 
   const Matrix S = stress_tensor(uu);
@@ -99,7 +99,7 @@ momentum(in Vector uu, in Scalar lnrho, const int3 vertexIdx) {
 }
 
 Vector
-induction(in Vector uu, in Vector 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)
@@ -116,15 +116,15 @@ induction(in Vector uu, in Vector aa) {
 
 // Declare input and output arrays using locations specified in the
 // array enum in astaroth.h
-in Scalar lnrho = VTXBUF_LNRHO;
+in ScalarField lnrho(VTXBUF_LNRHO);
-out Scalar out_lnrho = VTXBUF_LNRHO;
+out ScalarField out_lnrho(VTXBUF_LNRHO);
 
-in Vector uu = (int3) {VTXBUF_UUX, VTXBUF_UUY, VTXBUF_UUZ};
+in VectorField uu(VTXBUF_UUX, VTXBUF_UUY, VTXBUF_UUZ);
-out Vector out_uu = (int3) {VTXBUF_UUX,VTXBUF_UUY,VTXBUF_UUZ};
+out VectorField out_uu(VTXBUF_UUX,VTXBUF_UUY,VTXBUF_UUZ);
 
 #if LMAGNETIC
-in Vector aa = (int3) {VTXBUF_AX,VTXBUF_AY,VTXBUF_AZ};
+in VectorField aa(VTXBUF_AX,VTXBUF_AY,VTXBUF_AZ);
-out Vector out_aa = (int3) {VTXBUF_AX,VTXBUF_AY,VTXBUF_AZ};
+out VectorField out_aa(VTXBUF_AX,VTXBUF_AY,VTXBUF_AZ);
 #endif
 
 Kernel void
@@ -137,38 +137,3 @@ solve(Scalar dt) {
 
   WRITE(out_uu, RK3(out_uu, uu, momentum(uu, lnrho, vertexIdx), dt));
 }
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-

acc/pseudodisk/stencil_process_linegrav.sps

@@ -43,19 +43,19 @@ distance_x(Vector a, Vector b)
 }
 
 Vector
-value(in Vector uu)
+value(in VectorField uu)
 {
     return (Vector){value(uu.x), value(uu.y), value(uu.z)};
 }
 
 Matrix
-gradients(in Vector uu)
+gradients(in VectorField uu)
 {
     return (Matrix){gradient(uu.x), gradient(uu.y), gradient(uu.z)};
 }
 
 Scalar
-continuity(in Vector uu, in Scalar lnrho) {
+continuity(in VectorField uu, in ScalarField lnrho) {
     return -dot(value(uu), gradient(lnrho)) - divergence(uu);
 }
 
@@ -84,7 +84,7 @@ grav_force_line(const int3 vertexIdx)
 
 #if LENTROPY
 Vector
-momentum(in Vector uu, in Scalar lnrho, in Scalar ss, in Vector aa, const int3 vertexIdx) {
+momentum(in VectorField uu, in ScalarField lnrho, in ScalarField ss, in VectorField aa, const int3 vertexIdx) {
   Vector mom;
 
   const Matrix S = stress_tensor(uu);
@@ -109,7 +109,7 @@ momentum(in Vector uu, in Scalar lnrho, in Scalar ss, in Vector aa, const int3 v
 }
 #else
 Vector
-momentum(in Vector uu, in Scalar lnrho, const int3 vertexIdx) {
+momentum(in VectorField uu, in ScalarField lnrho, const int3 vertexIdx) {
   Vector mom;
 
   const Matrix S = stress_tensor(uu);
@@ -128,7 +128,7 @@ momentum(in Vector uu, in Scalar lnrho, const int3 vertexIdx) {
 
 
 Vector
-induction(in Vector uu, in Vector 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)
@@ -146,7 +146,7 @@ induction(in Vector uu, in Vector aa) {
 
 #if LENTROPY
 Scalar
-lnT( in Scalar ss, in Scalar lnrho) {
+lnT( in ScalarField ss, in ScalarField lnrho) {
   const Scalar lnT = LNT0 + value(ss) / cp_sound +
     (gamma - AcReal(1.)) * (value(lnrho) - LNRHO0);
   return lnT;
@@ -154,7 +154,7 @@ lnT( in Scalar ss, in Scalar lnrho) {
 
 // Nabla dot (K nabla T) / (rho T)
 Scalar
-heat_conduction( in Scalar ss, in Scalar lnrho) {
+heat_conduction( in ScalarField ss, in ScalarField lnrho) {
   const Scalar inv_cp_sound = AcReal(1.) / cp_sound;
 
   const Vector grad_ln_chi = (Vector) {
@@ -179,7 +179,7 @@ heating(const int i, const int j, const int k) {
 }
 
 Scalar
-entropy(in Scalar ss, in Vector uu, in Scalar lnrho, in Vector aa) {
+entropy(in ScalarField ss, in VectorField uu, in ScalarField lnrho, in VectorField aa) {
     const Matrix S = stress_tensor(uu);
 
     // nabla x nabla x A / mu0 = nabla(nabla dot A) - nabla^2(A)
@@ -198,21 +198,20 @@ entropy(in Scalar ss, in Vector uu, in Scalar lnrho, in Vector aa) {
 
 // Declare input and output arrays using locations specified in the
 // array enum in astaroth.h
-in Scalar lnrho = VTXBUF_LNRHO;
+in ScalarField lnrho(VTXBUF_LNRHO);
-out Scalar out_lnrho = VTXBUF_LNRHO;
+out ScalarField out_lnrho(VTXBUF_LNRHO);
-
-in Vector uu = (int3) {VTXBUF_UUX, VTXBUF_UUY, VTXBUF_UUZ};
-out Vector out_uu = (int3) {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);
 
 #if LMAGNETIC
-in Vector aa = (int3) {VTXBUF_AX,VTXBUF_AY,VTXBUF_AZ};
+in VectorField aa(VTXBUF_AX,VTXBUF_AY,VTXBUF_AZ);
-out Vector out_aa = (int3) {VTXBUF_AX,VTXBUF_AY,VTXBUF_AZ};
+out VectorField out_aa(VTXBUF_AX,VTXBUF_AY,VTXBUF_AZ);
 #endif
 
 #if LENTROPY
-in Scalar ss = VTXBUF_ENTROPY;
+in ScalarField ss(VTXBUF_ENTROPY);
-out Scalar out_ss = VTXBUF_ENTROPY;
+out ScalarField out_ss(VTXBUF_ENTROPY);
 #endif
 
 Kernel void

acc/samples/common_header.h

@@ -1,422 +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/>.
-*/
-
-/**
- * @file
- * \brief Brief info.
- *
- * Provides an interface to Astaroth. Contains all the necessary configuration
- * structs and functions for running the code on multiple GPUs.
- *
- * All interface functions declared here (such as acInit()) operate all GPUs
- * available in the node under the hood, and the user does not need any
- * information about the decomposition, synchronization or such to use these
- * functions.
- *
- */
-#pragma once
-
-/* Prevent name mangling */
-#ifdef __cplusplus
-extern "C" {
-#endif
-
-#include <float.h>        // FLT_EPSILON, etc
-#include <stdlib.h>       // size_t
-#include <vector_types.h> // CUDA vector types (float4, etc)
-
-
-/*
- * =============================================================================
- * Flags for auto-optimization
- * =============================================================================
- */
-#define AUTO_OPTIMIZE (0) // DEPRECATED TODO remove
-#define BOUNDCONDS_OPTIMIZE (0)
-#define GENERATE_BENCHMARK_DATA (0)
-
-// Device info
-#define REGISTERS_PER_THREAD (255)
-#define MAX_REGISTERS_PER_BLOCK (65536)
-#define MAX_THREADS_PER_BLOCK (1024)
-#define MAX_TB_DIM (MAX_THREADS_PER_BLOCK)
-#define NUM_ITERATIONS (10)
-#define WARP_SIZE (32)
-
-
-/*
- * =============================================================================
- * Compile-time constants used during simulation (user definable)
- * =============================================================================
- */
-#define STENCIL_ORDER (6)
-
-///////////// PAD TEST
-// NOTE: works only with nx is divisible by 32
-//#define PAD_LEAD (32 - STENCIL_ORDER/2)
-//#define PAD_SIZE (32 - STENCIL_ORDER)
-///////////// PAD TEST
-
-// L-prefix inherited from the old Astaroth, no idea what it means
-// MV: L means a Logical switch variale, something having true of false value.
-#define LFORCING (0) // Note: forcing is disabled currently in the files generated by acc (compiler of our DSL)
-#define LMAGNETIC (1)
-#define LENTROPY (1)
-#define LTEMPERATURE (0)
-
-#define AC_THERMAL_CONDUCTIVITY (AcReal(0.001)) // TODO: make an actual config parameter
-
-/*
- * =============================================================================
- * Identifiers used to construct the parameter lists for AcMeshInfo
- * (IntParamType and RealParamType)
- * (user definable)
- * =============================================================================
- */
-// clang-format off
-#define AC_FOR_INT_PARAM_TYPES(FUNC)\
-        /* cparams */\
-        FUNC(AC_nx), \
-        FUNC(AC_ny), \
-        FUNC(AC_nz), \
-        FUNC(AC_mx), \
-        FUNC(AC_my), \
-        FUNC(AC_mz), \
-        FUNC(AC_nx_min), \
-        FUNC(AC_ny_min), \
-        FUNC(AC_nz_min), \
-        FUNC(AC_nx_max), \
-        FUNC(AC_ny_max), \
-        FUNC(AC_nz_max), \
-        /* Other */\
-        FUNC(AC_max_steps), \
-        FUNC(AC_save_steps), \
-        FUNC(AC_bin_steps), \
-        FUNC(AC_bc_type), \
-        /* Additional */\
-        FUNC(AC_mxy),\
-        FUNC(AC_nxy),\
-        FUNC(AC_nxyz)
-#define AC_FOR_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_relhel), \
-        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),\
-        /* 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)
-// clang-format on
-
-/*
- * =============================================================================
- * Identifiers for VertexBufferHandle
- * (i.e. the arrays used to construct AcMesh)
- * (user definable)
- * =============================================================================
- */
-// clang-format off
-#define AC_FOR_HYDRO_VTXBUF_HANDLES(FUNC)\
-        FUNC(VTXBUF_LNRHO), \
-        FUNC(VTXBUF_UUX), \
-        FUNC(VTXBUF_UUY), \
-        FUNC(VTXBUF_UUZ), \
-        // FUNC(VTXBUF_DYE),
-
-#if LMAGNETIC
-#define AC_FOR_MAGNETIC_VTXBUF_HANDLES(FUNC)\
-        FUNC(VTXBUF_AX), \
-        FUNC(VTXBUF_AY), \
-        FUNC(VTXBUF_AZ),
-#else
-#define AC_FOR_MAGNETIC_VTXBUF_HANDLES(FUNC)
-#endif
-
-#if LENTROPY
-#define AC_FOR_ENTROPY_VTXBUF_HANDLES(FUNC)\
-        FUNC(VTXBUF_ENTROPY),
-#else
-#define AC_FOR_ENTROPY_VTXBUF_HANDLES(FUNC)
-#endif
-
-#if LTEMPERATURE
-#define AC_FOR_TEMPERATURE_VTXBUF_HANDLES(FUNC)\
-        FUNC(VTXBUF_TEMPERATURE),
-#else
-#define AC_FOR_TEMPERATURE_VTXBUF_HANDLES(FUNC)
-#endif
-
-#define AC_FOR_VTXBUF_HANDLES(FUNC)\
-        AC_FOR_HYDRO_VTXBUF_HANDLES(FUNC)\
-        AC_FOR_MAGNETIC_VTXBUF_HANDLES(FUNC)\
-        AC_FOR_ENTROPY_VTXBUF_HANDLES(FUNC)\
-        AC_FOR_TEMPERATURE_VTXBUF_HANDLES(FUNC)
-// clang-format on
-
-/*
- * =============================================================================
- * Single/double precision switch
- * =============================================================================
- */
-#if AC_DOUBLE_PRECISION == 1
-typedef double AcReal;
-typedef double3 AcReal3;
-#define AC_REAL_MAX (DBL_MAX)
-#define AC_REAL_MIN (DBL_MIN)
-#define AC_REAL_EPSILON (DBL_EPSILON)
-#else
-typedef float AcReal;
-typedef float3 AcReal3;
-#define AC_REAL_MAX (FLT_MAX)
-#define AC_REAL_MIN (FLT_MIN)
-#define AC_REAL_EPSILON (FLT_EPSILON)
-#endif
-
-typedef struct {
-    AcReal3 row[3];
-} AcMatrix;
-
-/*
- * =============================================================================
- * Helper macros
- * =============================================================================
- */
-#define AC_GEN_ID(X) X
-#define AC_GEN_STR(X) #X
-
-/*
- * =============================================================================
- * Error codes
- * =============================================================================
- */
-typedef enum { AC_SUCCESS = 0, AC_FAILURE = 1 } AcResult;
-
-/*
- * =============================================================================
- * Reduction types
- * =============================================================================
- */
-typedef enum {
-    RTYPE_MAX,
-    RTYPE_MIN,
-    RTYPE_RMS,
-    RTYPE_RMS_EXP,
-    NUM_REDUCTION_TYPES
-} ReductionType;
-
-/*
- * =============================================================================
- * Definitions for the enums and structs for AcMeshInfo (DO NOT TOUCH)
- * =============================================================================
- */
-typedef enum {
-    AC_FOR_INT_PARAM_TYPES(AC_GEN_ID),
-    NUM_INT_PARAMS
-} AcIntParam;
-
-typedef enum {
-    AC_FOR_REAL_PARAM_TYPES(AC_GEN_ID),
-    NUM_REAL_PARAMS
-} AcRealParam;
-
-extern const char* intparam_names[];  // Defined in astaroth.cu
-extern const char* realparam_names[]; // Defined in astaroth.cu
-
-typedef struct {
-    int int_params[NUM_INT_PARAMS];
-    AcReal real_params[NUM_REAL_PARAMS];
-} AcMeshInfo;
-
-/*
- * =============================================================================
- * Definitions for the enums and structs for AcMesh (DO NOT TOUCH)
- * =============================================================================
- */
-typedef enum {
-    AC_FOR_VTXBUF_HANDLES(AC_GEN_ID) NUM_VTXBUF_HANDLES
-} VertexBufferHandle;
-
-extern const char* vtxbuf_names[]; // Defined in astaroth.cu
-
-/*
-typedef struct {
-    AcReal* data;
-} VertexBuffer;
-*/
-
-// NOTE: there's no particular benefit declaring AcMesh a class, since
-// a library user may already have allocated memory for the vertex_buffers.
-// But then we would allocate memory again when the user wants to start
-// filling the class with data. => Its better to consider AcMesh as a
-// payload-only struct
-typedef struct {
-    AcReal* vertex_buffer[NUM_VTXBUF_HANDLES];
-    AcMeshInfo info;
-} AcMesh;
-
-#define acVertexBufferSize(mesh_info)                                              \
-    ((size_t)(mesh_info.int_params[AC_mx] * mesh_info.int_params[AC_my] *      \
-              mesh_info.int_params[AC_mz]))
-
-#define acVertexBufferSizeBytes(mesh_info)                                        \
-    (sizeof(AcReal) * acVertexBufferSize(mesh_info))
-
-#define acVertexBufferCompdomainSize(mesh_info)                                   \
-    (mesh_info.int_params[AC_nx] * mesh_info.int_params[AC_ny] *               \
-     mesh_info.int_params[AC_nz])
-
-#define acVertexBufferCompdomainSizeBytes(mesh_info)                             \
-    (sizeof(AcReal) * acVertexBufferCompdomainSize(mesh_info))
-
-#define acVertexBufferIdx(i, j, k, mesh_info)                                      \
-    ((i) + (j)*mesh_info.int_params[AC_mx] +                                   \
-     (k)*mesh_info.int_params[AC_mx] * mesh_info.int_params[AC_my])
-
-/*
- * =============================================================================
- * Astaroth interface
- * =============================================================================
- */
-/** Starting point of all GPU computation. Handles the allocation and
-initialization of *all memory needed on all GPUs in the node*. In other words,
-setups everything GPU-side so that calling any other GPU interface function
-afterwards does not result in illegal memory accesses. */
-AcResult acInit(const AcMeshInfo& mesh_info);
-
-/** Splits the host_mesh and distributes it among the GPUs in the node */
-AcResult acLoad(const AcMesh& host_mesh);
-AcResult acLoadWithOffset(const AcMesh& host_mesh, const int3& start, const int num_vertices);
-
-/** Does all three steps of the RK3 integration and computes the boundary
-conditions when necessary. Note that the boundary conditions are not applied
-after the final integration step.
-The result can be fetched to CPU memory with acStore(). */
-AcResult acIntegrate(const AcReal& dt);
-
-/** Performs a single RK3 step without computing boundary conditions. */
-AcResult acIntegrateStep(const int& isubstep, const AcReal& dt);
-
-/** Applies boundary conditions on the GPU meshs and communicates the
- ghost zones among GPUs if necessary */
-AcResult acBoundcondStep(void);
-
-/** Performs a scalar reduction on all GPUs in the node and returns the result.
- */
-AcReal acReduceScal(const ReductionType& rtype, const VertexBufferHandle& a);
-
-/** Performs a vector reduction on all GPUs in the node and returns the result.
- */
-AcReal acReduceVec(const ReductionType& rtype, const VertexBufferHandle& a,
-                   const VertexBufferHandle& b, const VertexBufferHandle& c);
-
-/** Stores the mesh distributed among GPUs of the node back to a single host
- * mesh */
-AcResult acStore(AcMesh* host_mesh);
-AcResult acStoreWithOffset(const int3& start, const int num_vertices, AcMesh* host_mesh);
-
-/** Frees all GPU allocations and resets all devices in the node. Should be
- * called at exit. */
-AcResult acQuit(void);
-
-/** Synchronizes all devices. All calls to Astaroth are asynchronous by default
-    unless otherwise stated. */
-AcResult acSynchronize(void);
-
-/* End extern "C" */
-#ifdef __cplusplus
-}
-#endif
-
-/*
- * =============================================================================
- * Notes
- * =============================================================================
- */
-/*
-typedef enum {
-    VTX_BUF_LNRHO,
-    VTX_BUF_UUX,
-    VTX_BUF_UUY,
-    VTX_BUF_UUZ,
-    NUM_VERTEX_BUFFER_HANDLES
-} VertexBufferHandle
-
-// LNRHO etc
-typedef struct {
-    AcReal* data;
-} VertexBuffer;
-
-// Host
-typedef struct {
-    VertexBuffer vertex_buffers[NUM_VERTEX_BUFFER_HANDLES];
-    MeshInfo info;
-} Mesh;
-
-// Device
-typedef struct {
-    VertexBuffer in[NUM_VERTEX_BUFFER_HANDLES];
-    VertexBuffer out[NUM_VERTEX_BUFFER_HANDLES];
-} VertexBufferArray;
-*/

acc/samples/sample_stencil_assembly.sas

@@ -1,49 +0,0 @@
-// TODO comments and reformatting
-
-//Scalar
-//dostuff(in Scalar uux)
-//{
-//   return uux[vertexIdx.x, vertexIdx.y, vertexIdx.z];
-//}
-
-// stencil_assembly.in
-Preprocessed Scalar
-some_exotic_stencil_computation(in Scalar uux)
-{
-    //#if STENCIL_ORDER == 2
-    //    const Scalar coefficients[] = {1, 1, 1};
-    //#else if STENCIL_ORDER == 4
-    //    const Scalar coefficients[] = {....};
-    //#endif
-
-    int i = vertexIdx.x;
-    int j = vertexIdx.y;
-    int k = vertexIdx.z;
-    const Scalar coefficients[] = {1, 2, 3};
-
-    return coefficients[0] * uux[i-1, j, k] + 
-           coefficients[1] * uux[i, j, k] + 
-           coefficients[2] * uux[i+1, j, k];
-}
-
-// stencil_process.in
-//in Scalar uux_in = VTXBUF_UUX;
-//out Scalar uux_out = VTXBUF_UUX;
-
-
-//Kernel
-//solve(Scalar dt)
-//{
-//    uux_out = some_exotic_stencil(uux_in);
-//}
-
-
-
-
-
-
-
-
-
-
-

acc/samples/sample_stencil_process.sps

@@ -1,149 +0,0 @@
-// TODO comments and reformatting
-
-uniform Scalar dsx;
-uniform Scalar dsy;
-uniform Scalar dsz;
-
-uniform Scalar GM_star;
-// Other uniforms types than Scalar or int not yet supported
-
-// BUILTIN
-//Scalar dot(...){}
-
-// BUILTIN
-//Scalar distance(Vector a, Vector b) { return sqrt(dot(a, b)); }
-
-// BUILTIN
-// Scalar first_derivative(Scalar pencil[], Scalar inv_ds) { return pencil[3] * inv_ds; }
-
-Scalar first_derivative(Scalar pencil[], Scalar inv_ds)
-{
-    Scalar res = 0;
-    for (int i = 0; i < STENCIL_ORDER+1; ++i) {
-        res = res + pencil[i];
-    }
-    return inv_ds * res;
-}
-
-Scalar distance(Vector a, Vector b)
-{
-    return sqrt(a.x * b.x + a.y * b.y + a.z * b.z); 
-}
-
-Scalar
-gravity_potential(int i, int j, int k)
-{
-    Vector star_pos = (Vector){0, 0, 0};
-    Vector vertex_pos = (Vector){dsx * i, dsy * j, dsz * k};
-    return GM_star / distance(star_pos, vertex_pos);
-}
-
-Scalar
-gradx_gravity_potential(int i, int j, int k)
-{
-    Scalar pencil[STENCIL_ORDER + 1];
-    for (int offset = -STENCIL_ORDER; offset <= STENCIL_ORDER; ++offset) {
-        pencil[offset+STENCIL_ORDER] = gravity_potential(i + offset, j, k);
-    }
-
-    Scalar inv_ds = Scalar(1.) / dsx;
-    return first_derivative(pencil, inv_ds);
-}
-
-Scalar
-grady_gravity_potential(int i, int j, int k)
-{
-    Scalar pencil[STENCIL_ORDER + 1];
-    for (int offset = -STENCIL_ORDER; offset <= STENCIL_ORDER; ++offset) {
-        pencil[offset+STENCIL_ORDER] = gravity_potential(i, j + offset, k);
-    }
-
-    Scalar inv_ds = Scalar(1.) / dsy;
-    return first_derivative(pencil, inv_ds);
-}
-
-Scalar
-gradz_gravity_potential(int i, int j, int k)
-{
-    Scalar pencil[STENCIL_ORDER + 1];
-    for (int offset = -STENCIL_ORDER; offset <= STENCIL_ORDER; ++offset) {
-        pencil[offset+STENCIL_ORDER] = gravity_potential(i, j, k + offset);
-    }
-
-    Scalar inv_ds = Scalar(1.) / dsz;
-    return first_derivative(pencil, inv_ds);
-}
-
-Vector
-momentum(int i, int j, int k, in Vector uu)
-{
-
-    Vector gravity_potential = (Vector){gradx_gravity_potential(i, j, k),
-                                      grady_gravity_potential(i, j, k),
-                                      gradz_gravity_potential(i, j, k)};
-
-
-    return gravity_potential;
-}
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-

acc/src/acc.l

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

acc/src/acc.y

@@ -101,6 +101,7 @@ exec_statement: declaration
               ;
 
 assignment: declaration '=' expression                                                  { $$ = astnode_create(NODE_UNKNOWN, $1, $3); $$->infix = '='; }
+          | declaration '(' expression_list ')'                                         { $$ = astnode_create(NODE_UNKNOWN, $1, $3); $$->infix = '('; $$->postfix = ')'; } // C++ style initializer
           | expression '=' expression                                                   { $$ = astnode_create(NODE_UNKNOWN, $1, $3); $$->infix = '='; }
           ;
 

acc/src/code_generator.c

@@ -229,8 +229,11 @@ translate_latest_symbol(void)
     // IN / OUT
     else if (symbol->type != SYMBOLTYPE_FUNCTION_PARAMETER &&
              (symbol->type_qualifier == IN || symbol->type_qualifier == OUT)) {
-        const char* inout_type_qualifier = "static __device__ const auto";
+
-        printf("%s %s%s", inout_type_qualifier, inout_name_prefix, symbol_table[handle].identifier);
+        printf("static __device__ const %s %s%s", symbol->type_specifier == SCALAR ? "int" : "int3",
+               inout_name_prefix, symbol_table[handle].identifier);
+        if (symbol->type_specifier == VECTOR)
+            printf(" = make_int3");
     }
     // OTHER
     else {
@@ -335,8 +338,8 @@ traverse(const ASTNode* node)
     // Preprocessed parameter boilerplate
     if (node->type == NODE_TYPE_QUALIFIER && node->token == PREPROCESSED)
         inside_preprocessed = true;
-    static const char
+    static const char preprocessed_parameter_boilerplate
-        preprocessed_parameter_boilerplate[] = "const int3 vertexIdx, const int3 globalVertexIdx, ";
+        [] = "const int3& vertexIdx, const int3& globalVertexIdx, ";
     if (inside_preprocessed && node->type == NODE_FUNCTION_PARAMETER_DECLARATION)
         printf("%s ", preprocessed_parameter_boilerplate);
     // BOILERPLATE END////////////////////////////////////////////////////////
@@ -491,8 +494,8 @@ generate_preprocessed_structures(void)
 
     // FILLING THE DATA STRUCT
     printf("static __device__ __forceinline__ AcRealData\
-            read_data(const int3 vertexIdx,\
+            read_data(const int3& vertexIdx,\
-                const int3 globalVertexIdx,\
+                const int3& globalVertexIdx,\
             AcReal* __restrict__ buf[], const int handle)\
             {\n\
                 %sData data;\n",
@@ -527,8 +530,8 @@ generate_preprocessed_structures(void)
         } AcReal3Data;\
         \
         static __device__ __forceinline__ AcReal3Data\
-        read_data(const int3 vertexIdx,\
+        read_data(const int3& vertexIdx,\
-                  const int3 globalVertexIdx,\
+                  const int3& globalVertexIdx,\
                   AcReal* __restrict__ buf[], const int3& handle)\
         {\
             AcReal3Data data;\

include/astaroth.h

@@ -49,6 +49,10 @@ AcResult acInit(const AcMeshInfo mesh_info);
  * called at exit. */
 AcResult acQuit(void);
 
+/** Checks whether there are any CUDA devices available. Returns AC_SUCCESS if there is 1 or more,
+ * AC_FAILURE otherwise. */
+AcResult acCheckDeviceAvailability(void);
+
 /** Synchronizes a specific stream. All streams are synchronized if STREAM_ALL is passed as a
  * parameter*/
 AcResult acSynchronizeStream(const Stream stream);
@@ -66,8 +70,8 @@ AcResult acStore(AcMesh* host_mesh);
  * substep and the user is responsible for calling acBoundcondStep before reading the data. */
 AcResult acIntegrate(const AcReal dt);
 
-/** Applies periodic boundary conditions for the Mesh distributed among the devices visible to the
+/** Applies periodic boundary conditions for the Mesh distributed among the devices visible to
- * caller*/
+ * the caller*/
 AcResult acBoundcondStep(void);
 
 /** Does a scalar reduction with the data stored in some vertex buffer */
@@ -81,6 +85,14 @@ AcReal acReduceVec(const ReductionType rtype, const VertexBufferHandle a,
  */
 AcResult acStoreWithOffset(const int3 dst, const size_t num_vertices, AcMesh* host_mesh);
 
+/** Will potentially be deprecated in later versions. Added only to fix backwards compatibility with
+ * PC for now.*/
+AcResult acIntegrateStep(const int isubstep, const AcReal dt);
+AcResult acIntegrateStepWithOffset(const int isubstep, const AcReal dt, const int3 start,
+                                   const int3 end);
+AcResult acSynchronize(void);
+AcResult acLoadWithOffset(const AcMesh host_mesh, const int3 src, const int num_vertices);
+
 #ifdef __cplusplus
 } // extern "C"
 #endif

include/astaroth_defines.h

@@ -44,8 +44,6 @@ typedef struct {
 } double3;
 #endif // __CUDACC__
 
-#include "stencil_defines.h"
-
 // Library flags
 #define VERBOSE_PRINTING (1)
 
@@ -68,6 +66,8 @@ typedef struct {
     AcReal3 row[3];
 } AcMatrix;
 
+#include "stencil_defines.h" // User-defined header
+
 // clang-format off
 #define AC_FOR_BUILTIN_INT_PARAM_TYPES(FUNC)\
         FUNC(AC_nx), \

src/core/astaroth.cu

@@ -19,6 +19,9 @@
 // #include "astaroth_defines.h"
 #include "astaroth.h"
 
+#include "errchk.h"
+#include "math_utils.h" // int3 + int3
+
 #define AC_GEN_STR(X) #X
 const char* intparam_names[]   = {AC_FOR_BUILTIN_INT_PARAM_TYPES(AC_GEN_STR) //
                                 AC_FOR_USER_INT_PARAM_TYPES(AC_GEN_STR)};
@@ -33,10 +36,13 @@ const char* vtxbuf_names[]     = {AC_FOR_VTXBUF_HANDLES(AC_GEN_STR)};
 
 static const int num_nodes = 1;
 static Node nodes[num_nodes];
+static int3 grid_n;
 
 AcResult
 acInit(const AcMeshInfo mesh_info)
 {
+    grid_n = (int3){mesh_info.int_params[AC_nx], mesh_info.int_params[AC_ny],
+                    mesh_info.int_params[AC_nz]};
     return acNodeCreate(0, mesh_info, &nodes[0]);
 }
 
@@ -46,6 +52,23 @@ acQuit(void)
     return acNodeDestroy(nodes[0]);
 }
 
+AcResult
+acCheckDeviceAvailability(void)
+{
+    int device_count; // Separate from num_devices to avoid side effects
+    ERRCHK_CUDA_ALWAYS(cudaGetDeviceCount(&device_count));
+    if (device_count > 0)
+        return AC_SUCCESS;
+    else
+        return AC_FAILURE;
+}
+
+AcResult
+acSynchronize(void)
+{
+    return acNodeSynchronizeStream(nodes[0], STREAM_ALL);
+}
+
 AcResult
 acSynchronizeStream(const Stream stream)
 {
@@ -80,6 +103,20 @@ acIntegrate(const AcReal dt)
     return acNodeIntegrate(nodes[0], dt);
 }
 
+AcResult
+acIntegrateStep(const int isubstep, const AcReal dt)
+{
+    const int3 start = (int3){NGHOST, NGHOST, NGHOST};
+    const int3 end   = start + grid_n;
+    return acNodeIntegrateSubstep(nodes[0], STREAM_DEFAULT, isubstep, start, end, dt);
+}
+
+AcResult
+acIntegrateStepWithOffset(const int isubstep, const AcReal dt, const int3 start, const int3 end)
+{
+    return acNodeIntegrateSubstep(nodes[0], STREAM_DEFAULT, isubstep, start, end, dt);
+}
+
 AcResult
 acBoundcondStep(void)
 {
@@ -108,3 +145,9 @@ acStoreWithOffset(const int3 dst, const size_t num_vertices, AcMesh* host_mesh)
 {
     return acNodeStoreMeshWithOffset(nodes[0], STREAM_DEFAULT, dst, dst, num_vertices, host_mesh);
 }
+
+AcResult
+acLoadWithOffset(const AcMesh host_mesh, const int3 src, const int num_vertices)
+{
+    return acNodeLoadMeshWithOffset(nodes[0], STREAM_DEFAULT, host_mesh, src, src, num_vertices);
+}

src/core/kernels/integration.cuh

@@ -46,6 +46,9 @@ IDX(const int3 idx)
     return DEVICE_VTXBUF_IDX(idx.x, idx.y, idx.z);
 }
 
+#define make_int3(a, b, c)                                                                         \
+    (int3) { (int)a, (int)b, (int)c }
+
 static __forceinline__ AcMatrix
 create_rotz(const AcReal radians)
 {

src/standalone/model/model_rk3.cc

@@ -282,14 +282,11 @@ der6x_upwd(const int i, const int j, const int k, const ModelScalar* arr)
 {
     ModelScalar inv_ds = get(AC_inv_dsx);
 
-    return ModelScalar(1.0/60.0)*inv_ds* (
+    return ModelScalar(1.0 / 60.0) * inv_ds *
-          -ModelScalar(    20.0)* arr[IDX(i,   j,   k)]
+           (-ModelScalar(20.0) * arr[IDX(i, j, k)] +
-          +ModelScalar(    15.0)*(arr[IDX(i+1, j,   k)]
+            ModelScalar(15.0) * (arr[IDX(i + 1, j, k)] + arr[IDX(i - 1, j, k)]) -
-          +                       arr[IDX(i-1, j,   k)])
+            ModelScalar(6.0) * (arr[IDX(i + 2, j, k)] + arr[IDX(i - 2, j, k)]) +
-          -ModelScalar(     6.0)*(arr[IDX(i+2, j,   k)]
+            arr[IDX(i + 3, j, k)] + arr[IDX(i - 3, j, k)]);
-          +                       arr[IDX(i-2, j,   k)])
-          +                       arr[IDX(i+3, j,   k)]
-          +                       arr[IDX(i-3, j,   k)]);
 }
 
 static inline ModelScalar
@@ -297,14 +294,11 @@ der6y_upwd(const int i, const int j, const int k, const ModelScalar* arr)
 {
     ModelScalar inv_ds = get(AC_inv_dsy);
 
-    return ModelScalar(1.0/60.0)*inv_ds* (
+    return ModelScalar(1.0 / 60.0) * inv_ds *
-          -ModelScalar(    20.0)* arr[IDX(i,   j,   k)]
+           (-ModelScalar(20.0) * arr[IDX(i, j, k)] +
-          +ModelScalar(    15.0)*(arr[IDX(i,   j+1, k)]
+            ModelScalar(15.0) * (arr[IDX(i, j + 1, k)] + arr[IDX(i, j - 1, k)]) -
-          +                       arr[IDX(i,   j-1, k)])
+            ModelScalar(6.0) * (arr[IDX(i, j + 2, k)] + arr[IDX(i, j - 2, k)]) +
-          -ModelScalar(     6.0)*(arr[IDX(i,   j+2, k)]
+            arr[IDX(i, j + 3, k)] + arr[IDX(i, j - 3, k)]);
-          +                       arr[IDX(i,   j-2, k)])
-          +                       arr[IDX(i,   j+3, k)]
-          +                       arr[IDX(i,   j-3, k)]);
 }
 
 static inline ModelScalar
@@ -312,14 +306,11 @@ der6z_upwd(const int i, const int j, const int k, const ModelScalar* arr)
 {
     ModelScalar inv_ds = get(AC_inv_dsz);
 
-    return ModelScalar(1.0/60.0)*inv_ds* (
+    return ModelScalar(1.0 / 60.0) * inv_ds *
-          -ModelScalar(    20.0)* arr[IDX(i,   j,   k  )]
+           (-ModelScalar(20.0) * arr[IDX(i, j, k)] +
-          +ModelScalar(    15.0)*(arr[IDX(i,   j,   k+1)]
+            ModelScalar(15.0) * (arr[IDX(i, j, k + 1)] + arr[IDX(i, j, k - 1)]) -
-          +                       arr[IDX(i,   j,   k-1)])
+            ModelScalar(6.0) * (arr[IDX(i, j, k + 2)] + arr[IDX(i, j, k - 2)]) +
-          -ModelScalar(     6.0)*(arr[IDX(i,   j,   k+2)]
+            arr[IDX(i, j, k + 3)] + arr[IDX(i, j, k - 3)]);
-          +                       arr[IDX(i,   j,   k-2)])
-          +                       arr[IDX(i,   j,   k+3)]
-          +                       arr[IDX(i,   j,   k-3)]);
 }
 #endif
 
@@ -339,7 +330,8 @@ compute_gradient(const int i, const int j, const int k, const ModelScalar* arr)
 static inline ModelVector
 compute_upwind(const int i, const int j, const int k, const ModelScalar* arr)
 {
-    return (ModelVector){der6x_upwd(i, j, k, arr), der6y_upwd(i, j, k, arr), der6z_upwd(i, j, k, arr)};
+    return (ModelVector){der6x_upwd(i, j, k, arr), der6y_upwd(i, j, k, arr),
+                         der6z_upwd(i, j, k, arr)};
 }
 #endif
 
@@ -416,8 +408,6 @@ gradients(const ModelVectorData& data)
     return (ModelMatrix){gradient(data.x), gradient(data.y), gradient(data.z)};
 }
 
-
-
 /*
  * =============================================================================
  * Level 0.3 (Built-in functions available during the Stencil Processing Stage)
@@ -571,10 +561,10 @@ contract(const ModelMatrix& mat)
 ModelScalar
 upwd_der6(const ModelVectorData& uu, const ModelScalarData& lnrho)
 {
-    ModelScalar uux = fabs(value(uu).x);
+    ModelScalar uux = fabsl(value(uu).x);
-    ModelScalar uuy = fabs(value(uu).y);
+    ModelScalar uuy = fabsl(value(uu).y);
-    ModelScalar uuz = fabs(value(uu).z);
+    ModelScalar uuz = fabsl(value(uu).z);
-    return uux*lnrho.upwind.x + uuy*lnrho.upwind.y + uuz*lnrho.upwind.z;
+    return uux * lnrho.upwind.x + uuy * lnrho.upwind.y + uuz * lnrho.upwind.z;
 }
 #endif
 
@@ -583,7 +573,7 @@ continuity(const ModelVectorData& uu, const ModelScalarData& 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);
@@ -760,24 +750,23 @@ helical_forcing(ModelScalar magnitude, ModelVector k_force, ModelVector xx, Mode
                 ModelVector ff_im, ModelScalar phi)
 {
     (void)magnitude; // WARNING: unused
-    xx.x = xx.x*(2.0*M_PI/(get(AC_dsx)*get(AC_nx)));
+    xx.x = xx.x * (2.0 * M_PI / (get(AC_dsx) * get(AC_nx)));
-    xx.y = xx.y*(2.0*M_PI/(get(AC_dsy)*get(AC_ny)));
+    xx.y = xx.y * (2.0 * M_PI / (get(AC_dsy) * get(AC_ny)));
-    xx.z = xx.z*(2.0*M_PI/(get(AC_dsz)*get(AC_nz)));
+    xx.z = xx.z * (2.0 * M_PI / (get(AC_dsz) * get(AC_nz)));
 
     ModelScalar cos_phi     = cosl(phi);
     ModelScalar sin_phi     = sinl(phi);
     ModelScalar cos_k_dot_x = cosl(dot(k_force, xx));
     ModelScalar sin_k_dot_x = sinl(dot(k_force, xx));
     // 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;
-    ModelScalar real_comp_phase = cos_k_dot_x*cos_phi - sin_k_dot_x*sin_phi;
+    ModelScalar real_comp_phase = cos_k_dot_x * cos_phi - sin_k_dot_x * sin_phi;
-    ModelScalar imag_comp_phase = cos_k_dot_x*sin_phi + sin_k_dot_x*cos_phi;
+    ModelScalar imag_comp_phase = cos_k_dot_x * sin_phi + sin_k_dot_x * cos_phi;
 
-
+    ModelVector force = (ModelVector){ff_re.x * real_comp_phase - ff_im.x * imag_comp_phase,
-    ModelVector force = (ModelVector){ ff_re.x*real_comp_phase - ff_im.x*imag_comp_phase,
+                                      ff_re.y * real_comp_phase - ff_im.y * imag_comp_phase,
-                             ff_re.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;
 }