diff --git a/acc/mhd_solver/stencil_assembly.sas b/acc/mhd_solver/stencil_assembly.sas
index 53b5a45..5c4f14a 100644
--- a/acc/mhd_solver/stencil_assembly.sas
+++ b/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,54 +15,54 @@ 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);
 
     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] 
+                    + 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] 
+                    - 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]
                     +               vertex[vertexIdx.x-3, vertexIdx.y, vertexIdx.z])};
 }
 
 Preprocessed Scalar
-der6y_upwd(in Scalar vertex)
+der6y_upwd(in ScalarField vertex)
 {
     Scalar inv_ds = DCONST_REAL(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] 
+                    +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] 
+                    -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]
                     +               vertex[vertexIdx.x, vertexIdx.y-3, vertexIdx.z])};
 }
 
 Preprocessed Scalar
-der6z_upwd(in Scalar vertex)
+der6z_upwd(in ScalarField vertex)
 {
     Scalar inv_ds = DCONST_REAL(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] 
+                    +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] 
+                    -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]
                     +               vertex[vertexIdx.x, vertexIdx.y, vertexIdx.z-3])};
 }
 
 #endif
 
 Preprocessed Matrix
-hessian(in Scalar vertex)
+hessian(in ScalarField vertex)
 {
     Matrix hessian;
 
diff --git a/acc/mhd_solver/stencil_process.sps b/acc/mhd_solver/stencil_process.sps
index 813f51a..ee8656e 100644
--- a/acc/mhd_solver/stencil_process.sps
+++ b/acc/mhd_solver/stencil_process.sps
@@ -25,14 +25,14 @@ 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);
@@ -42,13 +42,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.
@@ -59,7 +59,7 @@ continuity(in Vector uu, in Scalar lnrho) {
 
 #if LENTROPY
 Vector
-momentum(in Vector uu, in Scalar lnrho, in Scalar ss, in Vector aa) {
+momentum(in VectorField uu, in ScalarField lnrho, in ScalarField ss, in VectorField aa) {
     const Matrix S = stress_tensor(uu);
     const Scalar cs2 = cs2_sound * exp(gamma * value(ss) / cp_sound + (gamma - 1) * (value(lnrho) - lnrho0));
     const Vector  j = (Scalar(1.) / mu0) * (gradient_of_divergence(aa) - laplace_vec(aa)); // Current density
@@ -82,7 +82,7 @@ momentum(in Vector uu, in Scalar lnrho, in Scalar ss, in Vector aa) {
 }
 #elif LTEMPERATURE
 Vector
-momentum(in Vector uu, in Scalar lnrho, in Scalar tt) {
+momentum(in VectorField uu, in ScalarField lnrho, in ScalarField tt) {
   Vector mom;
 
   const Matrix S = stress_tensor(uu);
@@ -103,7 +103,7 @@ momentum(in Vector uu, in Scalar lnrho, in Scalar tt) {
 }
 #else
 Vector
-momentum(in Vector uu, in Scalar lnrho) {
+momentum(in VectorField uu, in ScalarField lnrho) {
   Vector mom;
 
   const Matrix S = stress_tensor(uu);
@@ -126,7 +126,7 @@ momentum(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)
@@ -144,7 +144,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;
@@ -152,7 +152,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);
@@ -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);
     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
@@ -191,7 +191,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;
@@ -290,26 +290,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;
-out Scalar out_lnrho = VTXBUF_LNRHO;
+in ScalarField 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};
-out Vector out_aa = (int3) {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
-in Scalar ss = VTXBUF_ENTROPY;
-out Scalar out_ss = VTXBUF_ENTROPY;
+in ScalarField ss(VTXBUF_ENTROPY);
+out ScalarField out_ss(VTXBUF_ENTROPY);
 #endif
 
 #if LTEMPERATURE
-in Scalar tt = VTXBUF_TEMPERATURE;
-out Scalar out_tt = VTXBUF_TEMPERATURE;
+in ScalarField tt(VTXBUF_TEMPERATURE);
+out ScalarField out_tt(VTXBUF_TEMPERATURE);
 #endif
 
 Kernel void
diff --git a/acc/pseudodisk/stencil_process_gravx.sps b/acc/pseudodisk/stencil_process_gravx.sps
index 0ccb4f9..3473a0b 100644
--- a/acc/pseudodisk/stencil_process_gravx.sps
+++ b/acc/pseudodisk/stencil_process_gravx.sps
@@ -36,31 +36,31 @@ uniform Scalar inv_dsx;
 uniform Scalar inv_dsy;
 uniform Scalar inv_dsz;
 
-Scalar 
-distance_x(Vector a, Vector b) 
-{ 
-    return sqrt(dot(a-b, a-b)); 
+Scalar
+distance_x(Vector a, Vector b)
+{
+    return sqrt(dot(a-b, a-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);
 }
 
-// Gravitation for in negative x-direction. 
-Vector 
+// Gravitation for in negative x-direction.
+Vector
 grav_force_line(const int3 vertexIdx)
 {
     Vector vertex_pos = (Vector){dsx * vertexIdx.x - xorig, dsy * vertexIdx.y - yorig, dsz * vertexIdx.z - zorig};
@@ -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;
-out Scalar out_lnrho = VTXBUF_LNRHO;
+in ScalarField 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};
-out Vector out_aa = (int3) {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
-in Scalar ss = VTXBUF_ENTROPY;
-out Scalar out_ss = VTXBUF_ENTROPY;
+in ScalarField ss(VTXBUF_ENTROPY);
+out ScalarField out_ss(VTXBUF_ENTROPY);
 #endif
 
 Kernel void
diff --git a/acc/pseudodisk/stencil_process_isotherm_gravx.sps b/acc/pseudodisk/stencil_process_isotherm_gravx.sps
index 9584774..8241d18 100644
--- a/acc/pseudodisk/stencil_process_isotherm_gravx.sps
+++ b/acc/pseudodisk/stencil_process_isotherm_gravx.sps
@@ -33,32 +33,32 @@ uniform Scalar inv_dsx;
 uniform Scalar inv_dsy;
 uniform Scalar inv_dsz;
 
-Scalar 
-distance_x(Vector a, Vector b) 
-{ 
-    return sqrt(dot(a-b, a-b)); 
+Scalar
+distance_x(Vector a, Vector b)
+{
+    return sqrt(dot(a-b, a-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);
 }
 
 
 // "Line-like" gravity with no y-component
-Vector 
+Vector
 grav_force_line(const int3 vertexIdx)
 {
     Vector vertex_pos = (Vector){dsx * vertexIdx.x - xorig, dsy * vertexIdx.y - yorig, dsz * vertexIdx.z - zorig};
@@ -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);
@@ -86,15 +86,15 @@ momentum(in Vector uu, in Scalar lnrho, const int3 vertexIdx) {
     cs2_sound * gradient(lnrho) +
     nu_visc *
     (laplace_vec(uu) + Scalar(1. / 3.) * gradient_of_divergence(uu) +
-      Scalar(2.) * mul(S, gradient(lnrho))) + zeta * gradient_of_divergence(uu) 
+      Scalar(2.) * mul(S, gradient(lnrho))) + zeta * gradient_of_divergence(uu)
       + grav_force_line(vertexIdx);
-  
+
 
   return mom;
 }
 
 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;
-out Scalar out_lnrho = VTXBUF_LNRHO;
+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);
 
-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};
-out Vector out_aa = (int3) {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
 
 Kernel void
@@ -132,38 +133,3 @@ solve(Scalar dt) {
 
   WRITE(out_uu, RK3(out_uu, uu, momentum(uu, lnrho, vertexIdx), dt));
 }
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
diff --git a/acc/pseudodisk/stencil_process_isotherm_linegrav.sps b/acc/pseudodisk/stencil_process_isotherm_linegrav.sps
index 9f90e7c..94207fe 100644
--- a/acc/pseudodisk/stencil_process_isotherm_linegrav.sps
+++ b/acc/pseudodisk/stencil_process_isotherm_linegrav.sps
@@ -33,32 +33,32 @@ uniform Scalar inv_dsx;
 uniform Scalar inv_dsy;
 uniform Scalar inv_dsz;
 
-Scalar 
-distance(Vector a, Vector b) 
-{ 
-    return sqrt(dot(a-b, a-b)); 
+Scalar
+distance(Vector a, Vector b)
+{
+    return sqrt(dot(a-b, a-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);
 }
 
 
 // "Line-like" gravity with no y-component
-Vector 
+Vector
 grav_force_line(const int3 vertexIdx)
 {
     Vector vertex_pos = (Vector){dsx * vertexIdx.x - xorig, dsy * vertexIdx.y - yorig, dsz * vertexIdx.z - zorig};
@@ -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);
@@ -91,15 +91,15 @@ momentum(in Vector uu, in Scalar lnrho, const int3 vertexIdx) {
     cs2_sound * gradient(lnrho) +
     nu_visc *
     (laplace_vec(uu) + Scalar(1. / 3.) * gradient_of_divergence(uu) +
-      Scalar(2.) * mul(S, gradient(lnrho))) + zeta * gradient_of_divergence(uu) 
+      Scalar(2.) * mul(S, gradient(lnrho))) + zeta * gradient_of_divergence(uu)
       + grav_force_line(vertexIdx);
-  
+
 
   return mom;
 }
 
 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;
-out Scalar out_lnrho = VTXBUF_LNRHO;
+in ScalarField 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};
-out Vector out_aa = (int3) {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
 
 Kernel void
@@ -137,38 +137,3 @@ solve(Scalar dt) {
 
   WRITE(out_uu, RK3(out_uu, uu, momentum(uu, lnrho, vertexIdx), dt));
 }
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
diff --git a/acc/pseudodisk/stencil_process_linegrav.sps b/acc/pseudodisk/stencil_process_linegrav.sps
index e42e680..a9d732f 100644
--- a/acc/pseudodisk/stencil_process_linegrav.sps
+++ b/acc/pseudodisk/stencil_process_linegrav.sps
@@ -36,31 +36,31 @@ uniform Scalar inv_dsx;
 uniform Scalar inv_dsy;
 uniform Scalar inv_dsz;
 
-Scalar 
-distance_x(Vector a, Vector b) 
-{ 
-    return sqrt(dot(a-b, a-b)); 
+Scalar
+distance_x(Vector a, Vector b)
+{
+    return sqrt(dot(a-b, a-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);
 }
 
 // "Line-like" gravity with no y-component
-Vector 
+Vector
 grav_force_line(const int3 vertexIdx)
 {
     Vector vertex_pos = (Vector){dsx * vertexIdx.x - xorig, dsy * vertexIdx.y - yorig, dsz * vertexIdx.z - zorig};
@@ -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;
-out Scalar 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 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);
 
 #if LMAGNETIC
-in Vector aa = (int3) {VTXBUF_AX,VTXBUF_AY,VTXBUF_AZ};
-out Vector out_aa = (int3) {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
-in Scalar ss = VTXBUF_ENTROPY;
-out Scalar out_ss = VTXBUF_ENTROPY;
+in ScalarField ss(VTXBUF_ENTROPY);
+out ScalarField out_ss(VTXBUF_ENTROPY);
 #endif
 
 Kernel void
diff --git a/acc/samples/common_header.h b/acc/samples/common_header.h
deleted file mode 100644
index 4701348..0000000
--- a/acc/samples/common_header.h
+++ /dev/null
@@ -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;
-*/
diff --git a/acc/samples/sample_stencil_assembly.sas b/acc/samples/sample_stencil_assembly.sas
deleted file mode 100644
index 4ddd64c..0000000
--- a/acc/samples/sample_stencil_assembly.sas
+++ /dev/null
@@ -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);
-//}
-
-
-
-
-
-
-
-
-
-
-
diff --git a/acc/samples/sample_stencil_process.sps b/acc/samples/sample_stencil_process.sps
deleted file mode 100644
index 219e40e..0000000
--- a/acc/samples/sample_stencil_process.sps
+++ /dev/null
@@ -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;
-}
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
diff --git a/acc/src/acc.l b/acc/src/acc.l
index e68fe8b..b180ae8 100644
--- a/acc/src/acc.l
+++ b/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; }
diff --git a/acc/src/acc.y b/acc/src/acc.y
index db49225..7138b1b 100644
--- a/acc/src/acc.y
+++ b/acc/src/acc.y
@@ -20,7 +20,7 @@ int yyget_lineno();
 %token VOID INT INT3
 %token IF ELSE FOR WHILE ELIF
 %token LEQU LAND LOR LLEQU
-%token KERNEL PREPROCESSED 
+%token KERNEL PREPROCESSED
 %token INPLACE_INC INPLACE_DEC
 
 %%
@@ -72,11 +72,11 @@ selection_statement: IF expression else_selection_statement
                    ;
 
 else_selection_statement: compound_statement                                            { $$ = astnode_create(NODE_UNKNOWN, $1, NULL); }
-                        | compound_statement elif_selection_statement                   { $$ = astnode_create(NODE_UNKNOWN, $1, $2); }                        
+                        | compound_statement elif_selection_statement                   { $$ = astnode_create(NODE_UNKNOWN, $1, $2); }
                         | compound_statement ELSE compound_statement                    { $$ = astnode_create(NODE_UNKNOWN, $1, $3); $$->infix = ELSE; }
                         ;
 
-elif_selection_statement: ELIF expression else_selection_statement                      { $$ = astnode_create(NODE_UNKNOWN, $2, $3); $$->prefix = ELIF; }   
+elif_selection_statement: ELIF expression else_selection_statement                      { $$ = astnode_create(NODE_UNKNOWN, $2, $3); $$->prefix = ELIF; }
                         ;
 
 iteration_statement: WHILE expression compound_statement                                { $$ = astnode_create(NODE_UNKNOWN, $2, $3); $$->prefix = WHILE; }
@@ -101,8 +101,9 @@ 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 = '='; }
-          ; 
+          ;
 
 return_statement: /* Empty */                                                           { $$ = astnode_create(NODE_UNKNOWN, NULL, NULL); }
                 | expression                                                            { $$ = astnode_create(NODE_UNKNOWN, $1, NULL); }
@@ -126,7 +127,7 @@ array_declaration: identifier '[' ']'
                  | identifier '[' expression ']'                                        { $$ = astnode_create(NODE_UNKNOWN, $1, $3); $$->infix = '['; $$->postfix = ']'; }
                  ;
 
-type_declaration: type_specifier                                                        { $$ = astnode_create(NODE_UNKNOWN, $1, NULL); }                 
+type_declaration: type_specifier                                                        { $$ = astnode_create(NODE_UNKNOWN, $1, NULL); }
                 | type_qualifier type_specifier                                         { $$ = astnode_create(NODE_UNKNOWN, $1, $2); }
                 ;
 
@@ -181,7 +182,7 @@ binary_operator: '+'
                | '/'                                                                    { $$ = astnode_create(NODE_UNKNOWN, NULL, NULL); $$->infix = yytext[0]; }
                | '*'                                                                    { $$ = astnode_create(NODE_UNKNOWN, NULL, NULL); $$->infix = yytext[0]; }
                | '<'                                                                    { $$ = astnode_create(NODE_UNKNOWN, NULL, NULL); $$->infix = yytext[0]; }
-               | '>'                                                                    { $$ = astnode_create(NODE_UNKNOWN, NULL, NULL); $$->infix = yytext[0]; }           
+               | '>'                                                                    { $$ = astnode_create(NODE_UNKNOWN, NULL, NULL); $$->infix = yytext[0]; }
                | LEQU                                                                   { $$ = astnode_create(NODE_UNKNOWN, NULL, NULL); astnode_set_buffer(yytext, $$); }
                | LAND                                                                   { $$ = astnode_create(NODE_UNKNOWN, NULL, NULL); astnode_set_buffer(yytext, $$); }
                | LOR                                                                    { $$ = astnode_create(NODE_UNKNOWN, NULL, NULL); astnode_set_buffer(yytext, $$); }
diff --git a/acc/src/code_generator.c b/acc/src/code_generator.c
index 93ca42f..bd05ea0 100644
--- a/acc/src/code_generator.c
+++ b/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
-        preprocessed_parameter_boilerplate[] = "const int3 vertexIdx, const int3 globalVertexIdx, ";
+    static const char preprocessed_parameter_boilerplate
+        [] = "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,\
-                const int3 globalVertexIdx,\
+            read_data(const int3& vertexIdx,\
+                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,\
-                  const int3 globalVertexIdx,\
+        read_data(const int3& vertexIdx,\
+                  const int3& globalVertexIdx,\
                   AcReal* __restrict__ buf[], const int3& handle)\
         {\
             AcReal3Data data;\
diff --git a/src/core/kernels/integration.cuh b/src/core/kernels/integration.cuh
index fec3b95..50db9ae 100644
--- a/src/core/kernels/integration.cuh
+++ b/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)
 {