From 57881e5961fa77237c317970271e643b4a5073ad Mon Sep 17 00:00:00 2001
From: Miikka Vaisala <mvaisala@asiaa.sinica.edu.tw>
Date: Fri, 21 Jun 2019 17:29:06 +0800
Subject: [PATCH] Implementetion for hyperdiffusion to be tested later.

---
 acc/mhd_solver/stencil_assembly.sas | 15 ++-------------
 acc/mhd_solver/stencil_process.sps  | 16 +++++++++++++++-
 2 files changed, 17 insertions(+), 14 deletions(-)

diff --git a/acc/mhd_solver/stencil_assembly.sas b/acc/mhd_solver/stencil_assembly.sas
index e095d08..b9f30a6 100644
--- a/acc/mhd_solver/stencil_assembly.sas
+++ b/acc/mhd_solver/stencil_assembly.sas
@@ -25,7 +25,7 @@ der6x_upwd(in Scalar vertex),
 }
 
 Preprocessed Scalar
-der6x_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] 
@@ -36,7 +36,7 @@ der6x_upwd(in Scalar vertex),
 }
 
 Preprocessed Scalar
-der6x_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] 
@@ -46,17 +46,6 @@ der6x_upwd(in Scalar vertex),
                     )}
 }
 
-Preprocessed Scalar
-der6_upwd(in Scalar vertex)
-{
-    der6 =  
-         
-         
-         
-
-    return der6;
-}
-
 Preprocessed Matrix
 hessian(in Scalar vertex)
 {
diff --git a/acc/mhd_solver/stencil_process.sps b/acc/mhd_solver/stencil_process.sps
index b87e8ed..5128f8c 100644
--- a/acc/mhd_solver/stencil_process.sps
+++ b/acc/mhd_solver/stencil_process.sps
@@ -3,6 +3,7 @@
 #define LTEMPERATURE (0)
 #define LGRAVITY (0)
 #define LFORCING (1)
+#define LUPWD (0)
 
 
 // Declare uniforms (i.e. device constants)
@@ -27,6 +28,14 @@ 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)
 {
@@ -41,7 +50,12 @@ gradients(in Vector uu)
 
 Scalar
 continuity(in Vector uu, in Scalar lnrho) {
-    return -dot(value(uu), gradient(lnrho)) - divergence(uu);
+    return -dot(value(uu), gradient(lnrho)) 
+#if LUPWD
+           //This is a corrective hyperdiffusion term for upwinding. 
+           + upwd_der6(uu, lnrho)
+#ENDIF
+           - divergence(uu);
 }
 
 #if LENTROPY