From d0b27a0347f43cd1de2bb9e4078f3f051b575555 Mon Sep 17 00:00:00 2001
From: Miikka Vaisala <mvaisala@asiaa.sinica.edu.tw>
Date: Mon, 24 Jun 2019 16:32:21 +0800
Subject: [PATCH] Upwinding terms now compile.

Not tested yet.
---
 acc/mhd_solver/stencil_assembly.sas | 45 +++++++++++++++--------------
 acc/mhd_solver/stencil_process.sps  | 22 ++++++++------
 2 files changed, 37 insertions(+), 30 deletions(-)

diff --git a/acc/mhd_solver/stencil_assembly.sas b/acc/mhd_solver/stencil_assembly.sas
index b9f30a6..380bf93 100644
--- a/acc/mhd_solver/stencil_assembly.sas
+++ b/acc/mhd_solver/stencil_assembly.sas
@@ -14,36 +14,39 @@ gradient(in Scalar vertex)
 }
 
 Preprocessed Scalar
-der6x_upwd(in Scalar vertex), 
+der6x_upwd(in Scalar vertex) 
 {
-    return (Scalar){(1.0/60.0)*inv_ds* (
-                    - 20.0* vertex[vertexIdx.x,   vertexIdx.y, vertexIdx.z] 
-                    + 15.0*(vertex[vertexIdx.x, vertexIdx.y+1, vertexIdx.z] + vertex[vertexIdx.x, vertexIdx.y-1, vertexIdx.z]) 
-                    -  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];
-                    )}
+    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, vertexIdx.y+1, vertexIdx.z] + vertex[vertexIdx.x, vertexIdx.y-1, vertexIdx.z]) 
+                    - Scalar( 6.0)*(vertex[vertexIdx.x, vertexIdx.y+2, vertexIdx.z] + vertex[vertexIdx.x, vertexIdx.y-2, vertexIdx.z])
+                    +       vertex[vertexIdx.x, vertexIdx.y+3, vertexIdx.z] + vertex[vertexIdx.x, vertexIdx.y-3, vertexIdx.z])};
 }
 
 Preprocessed Scalar
-der6y_upwd(in Scalar vertex), 
+der6y_upwd(in Scalar vertex) 
 {
-    return (Scalar){(1.0/60.0)*inv_ds* (
-                    - 20.0* vertex[vertexIdx.x,   vertexIdx.y, vertexIdx.z] 
-                    + 15.0*(vertex[vertexIdx.x, vertexIdx.y+1, vertexIdx.z] + vertex[vertexIdx.x, vertexIdx.y-1, vertexIdx.z]) 
-                    -  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];
-                    )}
+    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] + vertex[vertexIdx.x, vertexIdx.y-1, vertexIdx.z]) 
+                    -Scalar(  6.0)*(vertex[vertexIdx.x, vertexIdx.y+2, vertexIdx.z] + vertex[vertexIdx.x, vertexIdx.y-2, vertexIdx.z])
+                    +       vertex[vertexIdx.x, vertexIdx.y+3, vertexIdx.z] + vertex[vertexIdx.x, vertexIdx.y-3, vertexIdx.z])};
 }
 
 Preprocessed Scalar
-der6z_upwd(in Scalar vertex), 
+der6z_upwd(in Scalar vertex) 
 {
-    return (Scalar){(1.0/60.0)*inv_ds* (
-                    - 20.0* vertex[vertexIdx.x,   vertexIdx.y, vertexIdx.z] 
-                    + 15.0*(vertex[vertexIdx.x, vertexIdx.y, vertexIdx.z+1] + vertex[vertexIdx.x, vertexIdx.y, vertexIdx.z-1]) 
-                    -  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];
-                    )}
+    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] + vertex[vertexIdx.x, vertexIdx.y, vertexIdx.z-1]) 
+                    -Scalar(  6.0)*(vertex[vertexIdx.x, vertexIdx.y, vertexIdx.z+2] + vertex[vertexIdx.x, vertexIdx.y, vertexIdx.z-2])
+                    +       vertex[vertexIdx.x, vertexIdx.y, vertexIdx.z+3] + vertex[vertexIdx.x, vertexIdx.y, vertexIdx.z-3])};
 }
 
 Preprocessed Matrix
diff --git a/acc/mhd_solver/stencil_process.sps b/acc/mhd_solver/stencil_process.sps
index 5128f8c..213d879 100644
--- a/acc/mhd_solver/stencil_process.sps
+++ b/acc/mhd_solver/stencil_process.sps
@@ -3,7 +3,7 @@
 #define LTEMPERATURE (0)
 #define LGRAVITY (0)
 #define LFORCING (1)
-#define LUPWD (0)
+#define LUPWD (1)
 
 
 // Declare uniforms (i.e. device constants)
@@ -28,13 +28,6 @@ uniform int ny;
 uniform int nz;
 
 
-Scalar
-upwd_der6(in Vector uu, in Scalar lnrho)
-{
-    return (Scalar){uu.x*der6x_upwd(lnrho) + uu.y*der6y_upwd(lnrho) + uu.z*der6z_upwd(lnrho)}
-}
-
-
 
 Vector
 value(in Vector uu)
@@ -42,6 +35,17 @@ value(in Vector uu)
     return (Vector){value(uu.x), value(uu.y), value(uu.z)};
 }
 
+#if LUPWD
+Scalar
+upwd_der6(in Vector uu, in Scalar lnrho)
+{
+    Scalar uux = value(uu).x;
+    Scalar uuy = value(uu).y;
+    Scalar uuz = value(uu).z;
+    return (Scalar){uux*der6x_upwd(lnrho) + uuy*der6y_upwd(lnrho) + uuz*der6z_upwd(lnrho)};
+}
+#endif
+
 Matrix
 gradients(in Vector uu)
 {
@@ -54,7 +58,7 @@ continuity(in Vector uu, in Scalar lnrho) {
 #if LUPWD
            //This is a corrective hyperdiffusion term for upwinding. 
            + upwd_der6(uu, lnrho)
-#ENDIF
+#endif
            - divergence(uu);
 }