From 1e6740f9996b59f4cd90376b4e7b41cdd178c55e Mon Sep 17 00:00:00 2001
From: jpekkila <johannes.pekkila@aalto.fi>
Date: Mon, 1 Jul 2019 18:56:13 +0300
Subject: [PATCH] Added the equations for hydro only for both CPU and GPU.
 NOTE: NOT RIGOROUSLY CHECKED FOR CORRECTNESS. I just took the equations used
 with entropy and removed the terms which included entropy and magnetic fields

---
 acc/mhd_solver/stencil_process.sps | 32 ++++++++++++------------
 src/standalone/model/model_rk3.cc  | 39 ++++++++----------------------
 2 files changed, 25 insertions(+), 46 deletions(-)

diff --git a/acc/mhd_solver/stencil_process.sps b/acc/mhd_solver/stencil_process.sps
index 52e9402..c1a2214 100644
--- a/acc/mhd_solver/stencil_process.sps
+++ b/acc/mhd_solver/stencil_process.sps
@@ -100,23 +100,21 @@ momentum(in Vector uu, in Scalar lnrho, in Scalar tt) {
 #else
 Vector
 momentum(in Vector uu, in Scalar lnrho) {
-  Vector mom;
-
-  const Matrix S = stress_tensor(uu);
-
-    // Isothermal: we have constant speed of sound
-
-  mom = -mul(gradients(uu), value(uu)) -
-    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);
-
-  #if LGRAVITY
-  mom = mom - (Vector){0, 0, -10.0};
-  #endif
-
-  return mom;
+    // !!!!!!!!!!!!!!!!%JP: NOTE TODO IMPORTANT!!!!!!!!!!!!!!!!!!!!!!!!
+    // NOT CHECKED FOR CORRECTNESS: USE AT YOUR OWN RISK
+    const Matrix S = stress_tensor(uu);
+    const Scalar cs2 = cs2_sound * exp((gamma - 1) * (value(lnrho) - LNRHO0));
+    // Regex replace CPU constants with get\(AC_([a-zA-Z_0-9]*)\)
+    // \1
+    const Vector mom = - mul(gradients(uu), value(uu))
+                                                       - cs2 * 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);
+    return mom;
 }
 #endif
 
diff --git a/src/standalone/model/model_rk3.cc b/src/standalone/model/model_rk3.cc
index dd04bcf..a18f1b8 100644
--- a/src/standalone/model/model_rk3.cc
+++ b/src/standalone/model/model_rk3.cc
@@ -587,38 +587,19 @@ momentum(const ModelVectorData& uu, const ModelScalarData& lnrho
                                                ModelScalar(2.) * mul(S, gradient(lnrho))) +
                             get(AC_zeta) * gradient_of_divergence(uu);
     return mom;
-#endif
+#else
+    // !!!!!!!!!!!!!!!!%JP: NOTE TODO IMPORTANT!!!!!!!!!!!!!!!!!!!!!!!!
+    // NOT CHECKED FOR CORRECTNESS: USE AT YOUR OWN RISK
+    const ModelMatrix S   = stress_tensor(uu);
+    const ModelScalar cs2 = get(AC_cs2_sound) * expl((get(AC_gamma) - 1) * (value(lnrho) - LNRHO0));
 
-#if 0
-    const ModelMatrix S = stress_tensor(uu);
-
-    //#if LENTROPY
-    //const ModelScalar lnrho0 = 1; // TODO correct lnrho0
-    const ModelScalar cs02 = get(AC_cs2_sound); // TODO better naming
-    const ModelScalar cs2 = cs02;// * expl(get(AC_gamma) * value(ss) / get(AC_cp_sound) + (get(AC_gamma)-ModelScalar(1.l)) * (value(lnrho) - lnrho0));
-
-    mom = -mul(gradients(uu), value(uu)) -
-    cs2 * ((ModelScalar(1.) / get(AC_cp_sound)) * gradient(ss) + gradient(lnrho)) +
-    get(AC_nu_visc) *
-    (laplace_vec(uu) + ModelScalar(1.l / 3.l) * gradient_of_divergence(uu) +
-      ModelScalar(2.l) * mul(S, gradient(lnrho))) + get(AC_zeta) * gradient_of_divergence(uu);
-
-    const ModelVector grad_div = gradient_of_divergence(aa);
-    const ModelVector lap = laplace_vec(aa);
-    const ModelVector j = (ModelScalar(1.l) / get(AC_mu0)) * (grad_div - lap);
-    const ModelVector B = curl(aa);
-    mom = mom + (ModelScalar(1.l) / expl(value(lnrho))) * cross(j, B);
-    //#else // Basic hydro
-        const ModelScalar cs02 = get(AC_cs2_sound);
-        mom = -mul(gradients(uu), value(uu)) -
-          cs02 * gradient(lnrho) +
-          get(AC_nu_visc) *
-              (laplace_vec(uu) + ModelScalar(1. / 3.) * gradient_of_divergence(uu) +
-               ModelScalar(2.) * mul(S, gradient(lnrho))) + get(AC_zeta) * gradient_of_divergence(uu);
-    //#endif
+    const ModelVector mom = -mul(gradients(uu), value(uu)) - cs2 * gradient(lnrho) +
+                            get(AC_nu_visc) * (laplace_vec(uu) +
+                                               ModelScalar(1. / 3.) * gradient_of_divergence(uu) +
+                                               ModelScalar(2.) * mul(S, gradient(lnrho))) +
+                            get(AC_zeta) * gradient_of_divergence(uu);
     return mom;
 #endif
-    return (ModelVector){NAN, NAN, NAN}; // TODO HYDRO ONLY MODEL SOLUTION
 }
 
 static inline ModelVector