From 7e72a47460510617c5bf3ced2b8dabe5a2bf5355 Mon Sep 17 00:00:00 2001
From: stgeke <stgeke@gmail.com>
Date: Thu, 26 Mar 2020 10:01:37 +0100
Subject: [PATCH] Fix multiscalar bug (#83)

---
 okl/core/cdsAdvectionHex3D.okl                |  30 +-
 okl/core/cdsSumMakefHex3D.okl                 |   4 +-
 okl/core/math.okl                             |  13 +
 okl/plugins/RANSktau.okl                      | 317 ++++++++++++++++++
 okl/{avg => plugins}/avg.okl                  |   0
 .../velRecycling.okl                          |   0
 src/core/cds.cpp                              |   1 -
 src/core/cds.h                                |   2 +-
 src/core/ins.h                                |   1 +
 src/core/insSetup.cpp                         |  14 +-
 src/core/parReader.cpp                        |  82 ++---
 src/core/runTime.cpp                          |  38 ++-
 src/nekInterface/nekInterface.f               |  11 +-
 src/{core => }/nrs.hpp                        |   0
 src/udf/udf.hpp                               |   4 -
 15 files changed, 433 insertions(+), 84 deletions(-)
 create mode 100644 okl/plugins/RANSktau.okl
 rename okl/{avg => plugins}/avg.okl (100%)
 rename okl/{velRecycling => plugins}/velRecycling.okl (100%)
 rename src/{core => }/nrs.hpp (100%)

diff --git a/okl/core/cdsAdvectionHex3D.okl b/okl/core/cdsAdvectionHex3D.okl
index ea29d84ba..8596e4052 100644
--- a/okl/core/cdsAdvectionHex3D.okl
+++ b/okl/core/cdsAdvectionHex3D.okl
@@ -40,16 +40,16 @@ SOFTWARE.
 
     for(int k=0;k<p_Nq;++k;@inner(2)){
       for(int j=0;j<p_Nq;++j;@inner(1)){
-  for(int i=0;i<p_Nq;++i;@inner(0)){
-    dlong id = i+j*p_Nq;
-    if(k==0)
-      s_D[0][id] = D[id];
-
-    id = e*p_Np + k*p_Nq*p_Nq + j*p_Nq + i;
-
-    const dfloat Sn = S[id+0*soffset];
-    s_S[k][j][i] = Sn;
-  }
+        for(int i=0;i<p_Nq;++i;@inner(0)){
+          dlong id = i+j*p_Nq;
+          if(k==0)
+            s_D[0][id] = D[id];
+      
+          id = e*p_Np + k*p_Nq*p_Nq + j*p_Nq + i;
+      
+          const dfloat Sn = S[id+soffset];
+          s_S[k][j][i] = Sn;
+        }
       }
     }
 
@@ -91,9 +91,9 @@ SOFTWARE.
     const dfloat u = U[id+0*voffset];
     const dfloat v = U[id+1*voffset];
     const dfloat w = U[id+2*voffset];
-    const dfloat rhoM = RHO[id];
+    const dfloat rhoM = RHO[id+soffset];
 
-    NS[id+0*soffset] =  rhoM*(u*dSdx + v*dSdy + w*dSdz);
+    NS[id] =  rhoM*(u*dSdx + v*dSdy + w*dSdz);
   }
       }
     }
@@ -166,7 +166,7 @@ SOFTWARE.
       s_U[b][a] = U[id+0*voffset]; 
       s_V[b][a] = U[id+1*voffset];
       s_W[b][a] = U[id+2*voffset];      
-      s_S[b][a] = S[id+0*soffset];      
+      s_S[b][a] = S[id+soffset];      
     }
   }
       }
@@ -329,9 +329,9 @@ SOFTWARE.
             const dlong id = e*p_Np + c*p_Nq*p_Nq + b*p_Nq + a;
             
             const dfloat IJW   = vgeo[gid + p_IJWID*p_Np];
-            const dfloat rhoM = RHO[id];
+            const dfloat rhoM = RHO[id+soffset];
             
-            NS[id+0*soffset] = rhoM * IJW*rhsS;
+            NS[id] = rhoM * IJW*rhsS;
           }
         }
       }
diff --git a/okl/core/cdsSumMakefHex3D.okl b/okl/core/cdsSumMakefHex3D.okl
index 02d3f733d..f1c3fc73d 100644
--- a/okl/core/cdsSumMakefHex3D.okl
+++ b/okl/core/cdsSumMakefHex3D.okl
@@ -59,13 +59,13 @@ SOFTWARE.
                 sum += rhoM * idt*NSm;
               }else {
                 for (int s=0;s<p_Nstages;s++) {
-                  const dfloat Sm = S[id+offset+s*sOffset];
+                  const dfloat Sm = S[id+offset + s*sOffset];
                   sum += rhoM*idt*extbdfB[s]*Sm;
                 }
               }
               // already multiplied by rho
               for (int s=0;s<p_Nstages;s++) {
-                const dfloat FSm = FS[id+offset+s*sOffset];
+                const dfloat FSm = FS[id+offset + s*sOffset];
                 sum += extbdfA[s]*FSm;
               }
               BF[id+offset] = JW*sum;
diff --git a/okl/core/math.okl b/okl/core/math.okl
index b6ce2f66a..f51af8b2d 100644
--- a/okl/core/math.okl
+++ b/okl/core/math.okl
@@ -37,3 +37,16 @@
     for(int t=0;t<p_blockSize;++t;@inner(0)) if(t<  1) wxx[b] = (s_wxx[0]>s_wxx[1]) ? s_wxx[0]:s_wxx[1];
   }
 }
+
+@kernel void scalarScaledAdd(const dlong N,
+             const dfloat a,
+             const dfloat b,
+             @restrict const dfloat * X,
+             @restrict dfloat * Y){
+
+  for(dlong n=0;n<N;++n;@tile(256,@outer,@inner)){
+    if(n<N){
+      Y[n] = a + b*X[n]; 
+    }
+  }
+}
diff --git a/okl/plugins/RANSktau.okl b/okl/plugins/RANSktau.okl
new file mode 100644
index 000000000..934decf0a
--- /dev/null
+++ b/okl/plugins/RANSktau.okl
@@ -0,0 +1,317 @@
+@kernel void computeHex3D(const dlong Nelements,
+		          const dlong offset,
+		          const dfloat rho,
+		          const dfloat mue,
+		          @restrict const  dfloat * vgeo,
+		          @restrict const  dfloat * D,
+		          @restrict const  dfloat * K,
+		          @restrict const  dfloat * TAU,
+		          @restrict const  dfloat * STMAG2,
+                          @restrict const  dfloat * OIOJSK,
+                          @restrict dfloat * SRCDIAG,
+		          @restrict dfloat * SRC){
+
+  for(dlong e=0; e<Nelements; ++e; @outer(0)) {
+    @shared dfloat s_k[p_Nq][p_Nq][p_Nq];
+    @shared dfloat s_tau[p_Nq][p_Nq][p_Nq];
+    @shared dfloat s_tauSqrt[p_Nq][p_Nq][p_Nq];
+
+    @shared dfloat s_D[p_Nq][p_Nq];
+
+    for(int k=0;k<p_Nq;++k;@inner(2)){
+      for(int j=0;j<p_Nq;++j;@inner(1)){
+	for(int i=0;i<p_Nq;++i;@inner(0)){
+	  dlong id = i+j*p_Nq;
+	  if(k==0) s_D[0][id] = D[id];
+
+	  id = e*p_Np + k*p_Nq*p_Nq + j*p_Nq + i;
+
+	  const dfloat kn   = K[id];
+	  const dfloat taun = TAU[id];
+
+	  s_k[k][j][i]       = kn;
+	  s_tau[k][j][i]     = taun;
+	  s_tauSqrt[k][j][i] = sqrt(taun);
+	}
+      }
+    }
+
+    @barrier("local");
+
+    for(int k=0;k<p_Nq;++k;@inner(2)){
+      for(int j=0;j<p_Nq;++j;@inner(1)){
+	for(int i=0;i<p_Nq;++i;@inner(0)){
+
+	  const dlong gid = e*p_Np*p_Nvgeo+ k*p_Nq*p_Nq + j*p_Nq +i;
+	  const dfloat drdx = vgeo[gid + p_RXID*p_Np];
+	  const dfloat drdy = vgeo[gid + p_RYID*p_Np];
+	  const dfloat drdz = vgeo[gid + p_RZID*p_Np];
+	  const dfloat dsdx = vgeo[gid + p_SXID*p_Np];
+	  const dfloat dsdy = vgeo[gid + p_SYID*p_Np];
+	  const dfloat dsdz = vgeo[gid + p_SZID*p_Np];
+	  const dfloat dtdx = vgeo[gid + p_TXID*p_Np];
+	  const dfloat dtdy = vgeo[gid + p_TYID*p_Np];
+	  const dfloat dtdz = vgeo[gid + p_TZID*p_Np];
+          const dfloat JW   = vgeo[gid + p_JWID*p_Np];
+
+	  dfloat dkdr = 0, dkds = 0, dkdt = 0;
+	  dfloat dtaudr = 0, dtauds = 0, dtaudt = 0;
+	  dfloat dtauSqrtdr = 0, dtauSqrtds = 0, dtauSqrtdt = 0;
+
+#pragma unroll p_Nq
+	  for (int n=0;n<p_Nq;n++) {
+	    const dfloat Dr = s_D[i][n];
+	    const dfloat Ds = s_D[j][n];
+	    const dfloat Dt = s_D[k][n];
+
+	    dkdr += Dr*s_k[k][j][n];
+	    dkds += Ds*s_k[k][n][i];
+	    dkdt += Dt*s_k[n][j][i];
+
+	    dtaudr += Dr*s_tau[k][j][n];
+	    dtauds += Ds*s_tau[k][n][i];
+	    dtaudt += Dt*s_tau[n][j][i];
+
+	    dtauSqrtdr += Dr*s_tauSqrt[k][j][n];
+	    dtauSqrtds += Ds*s_tauSqrt[k][n][i];
+	    dtauSqrtdt += Dt*s_tauSqrt[n][j][i];
+	  }
+
+	  const dfloat dkdx = drdx*dkdr + dsdx*dkds + dtdx*dkdt;
+	  const dfloat dkdy = drdy*dkdr + dsdy*dkds + dtdy*dkdt;
+	  const dfloat dkdz = drdz*dkdr + dsdz*dkds + dtdz*dkdt;
+
+	  const dfloat dtaudx = drdx*dtaudr + dsdx*dtauds + dtdx*dtaudt;
+	  const dfloat dtaudy = drdy*dtaudr + dsdy*dtauds + dtdy*dtaudt;
+	  const dfloat dtaudz = drdz*dtaudr + dsdz*dtauds + dtdz*dtaudt;
+
+	  const dfloat dtauSqrtdx = drdx*dtauSqrtdr + dsdx*dtauSqrtds + dtdx*dtauSqrtdt;
+	  const dfloat dtauSqrtdy = drdy*dtauSqrtdr + dsdy*dtauSqrtds + dtdy*dtauSqrtdt;
+	  const dfloat dtauSqrtdz = drdz*dtauSqrtdr + dsdz*dtauSqrtds + dtdz*dtauSqrtdt;
+
+          const dfloat xk  = -(dkdx*dtaudx + dkdy*dtaudy + dkdz*dtaudz);
+          const dfloat xt  = dtaudx*dtaudx + dtaudy*dtaudy + dtaudz*dtaudz;
+          const dfloat xtq = dtauSqrtdx*dtauSqrtdx + dtauSqrtdy*dtauSqrtdy + dtauSqrtdz*dtauSqrtdz;
+	
+	  const dlong id = e*p_Np + k*p_Nq*p_Nq + j*p_Nq + i;
+
+	  const dfloat stMag2 = STMAG2[id];
+          const dfloat OiOjSk = OIOJSK[id]; 
+
+	  const dfloat kk  = K[id];
+	  const dfloat tau = TAU[id];
+
+          const dfloat mu_t  = rho*p_alpinf_str*kk*tau;
+          dfloat itau = 0;
+          if(tau > 0) itau = 1/tau;
+
+          dfloat sigd = p_sigd_min;
+          dfloat f_beta_str = 1.0;
+          if (xk > 0) {
+            const dfloat xk3 = xk*xk*tau*tau;
+            sigd = p_sigd_max;
+            f_beta_str = (1.0 + p_fb_c1st*xk3)/(1.0 + p_fb_c2st*xk3);
+          }
+        
+          // compute source term for k 
+          const dfloat Y_k = rho*p_betainf_str*f_beta_str*itau;
+          const dfloat kSrc  = mu_t*stMag2;
+          const dfloat kDiag = Y_k; 
+
+          // compute rource term for omega
+          const dfloat x_w = abs(OiOjSk)*(tau*tau*tau*p_ibetainf_str3);
+          const dfloat f_b = (1.0 + p_fb_c1*x_w)/(1.0 + p_fb_c2*x_w);
+          dfloat tauSrc  = rho*(p_beta0*f_b - sigd*xk*tau);
+          dfloat tauDiag = rho*tau*p_alp_inf*stMag2 + 
+                           8.0*rho*p_alpinf_str*kk*xtq*p_sigma_tau;
+
+          // apply correction
+          const dfloat S_tau = 8.0*mue*xtq;
+          if (tau <= p_tiny) 
+            tauSrc  -= S_tau;
+          else
+            tauDiag += S_tau*itau;
+
+          SRC[id + 0*offset]     = kSrc;
+          SRC[id + 1*offset]     = tauSrc;
+          SRCDIAG[id + 0*offset] = kDiag;
+          SRCDIAG[id + 1*offset] = tauDiag;
+	}
+      }
+    }
+  }
+}
+
+@kernel void SijOijHex3D(const dlong Nelements,
+		         const dlong offset,
+		         @restrict const  dfloat * vgeo,
+		         @restrict const  dfloat * D,
+		         @restrict const  dfloat * U,
+		         @restrict dfloat * SO){
+
+  for(dlong e=0; e<Nelements; ++e; @outer(0)) {
+    @shared dfloat s_U[p_Nq][p_Nq][p_Nq];
+    @shared dfloat s_V[p_Nq][p_Nq][p_Nq];
+    @shared dfloat s_W[p_Nq][p_Nq][p_Nq];
+
+    @shared dfloat s_D[p_Nq][p_Nq];
+
+    for(int k=0;k<p_Nq;++k;@inner(2)){
+      for(int j=0;j<p_Nq;++j;@inner(1)){
+	for(int i=0;i<p_Nq;++i;@inner(0)){
+	  dlong id = i+j*p_Nq;
+	  if(k==0) s_D[0][id] = D[id];
+
+	  id = e*p_Np + k*p_Nq*p_Nq + j*p_Nq + i;
+
+          const dfloat Un = U[id+0*offset];
+          const dfloat Vn = U[id+1*offset];
+          const dfloat Wn = U[id+2*offset];
+
+          s_U[k][j][i] = Un;
+          s_V[k][j][i] = Vn;
+          s_W[k][j][i] = Wn;
+	}
+      }
+    }
+
+    @barrier("local");
+
+    for(int k=0;k<p_Nq;++k;@inner(2)){
+      for(int j=0;j<p_Nq;++j;@inner(1)){
+	for(int i=0;i<p_Nq;++i;@inner(0)){
+
+	  const dlong gid = e*p_Np*p_Nvgeo+ k*p_Nq*p_Nq + j*p_Nq +i;
+	  const dfloat drdx = vgeo[gid + p_RXID*p_Np];
+	  const dfloat drdy = vgeo[gid + p_RYID*p_Np];
+	  const dfloat drdz = vgeo[gid + p_RZID*p_Np];
+	  const dfloat dsdx = vgeo[gid + p_SXID*p_Np];
+	  const dfloat dsdy = vgeo[gid + p_SYID*p_Np];
+	  const dfloat dsdz = vgeo[gid + p_SZID*p_Np];
+	  const dfloat dtdx = vgeo[gid + p_TXID*p_Np];
+	  const dfloat dtdy = vgeo[gid + p_TYID*p_Np];
+	  const dfloat dtdz = vgeo[gid + p_TZID*p_Np];
+          const dfloat JW   = vgeo[gid + p_JWID*p_Np];
+
+	  dfloat dudr = 0, duds = 0, dudt = 0;
+	  dfloat dvdr = 0, dvds = 0, dvdt = 0;
+	  dfloat dwdr = 0, dwds = 0, dwdt = 0;
+
+#pragma unroll p_Nq
+          for (int n=0;n<p_Nq;n++) {
+            const dfloat Dr = s_D[i][n];
+            const dfloat Ds = s_D[j][n];
+            const dfloat Dt = s_D[k][n];
+            dudr += Dr*s_U[k][j][n];
+            duds += Ds*s_U[k][n][i];
+            dudt += Dt*s_U[n][j][i];
+
+            dvdr += Dr*s_V[k][j][n];
+            dvds += Ds*s_V[k][n][i];
+            dvdt += Dt*s_V[n][j][i];
+
+            dwdr += Dr*s_W[k][j][n];
+            dwds += Ds*s_W[k][n][i];
+            dwdt += Dt*s_W[n][j][i];
+          }
+
+          const dfloat dudx = drdx*dudr + dsdx*duds + dtdx*dudt;
+          const dfloat dudy = drdy*dudr + dsdy*duds + dtdy*dudt;
+          const dfloat dudz = drdz*dudr + dsdz*duds + dtdz*dudt;
+
+          const dfloat dvdx = drdx*dvdr + dsdx*dvds + dtdx*dvdt;
+          const dfloat dvdy = drdy*dvdr + dsdy*dvds + dtdy*dvdt;
+          const dfloat dvdz = drdz*dvdr + dsdz*dvds + dtdz*dvdt;
+
+          const dfloat dwdx = drdx*dwdr + dsdx*dwds + dtdx*dwdt;
+          const dfloat dwdy = drdy*dwdr + dsdy*dwds + dtdy*dwdt;
+          const dfloat dwdz = drdz*dwdr + dsdz*dwds + dtdz*dwdt;
+
+          const dlong id = e*p_Np + k*p_Nq*p_Nq + j*p_Nq + i;
+
+          SO[id+0*offset] = JW*2*dudx;
+          SO[id+1*offset] = JW*2*dvdy;
+          SO[id+2*offset] = JW*2*dwdz;
+          SO[id+3*offset] = JW*(dudy + dvdx);
+          SO[id+4*offset] = JW*(dvdz + dwdy);
+          SO[id+5*offset] = JW*(dudz + dwdx);
+
+          SO[id+6*offset] = JW*(dwdy - dvdz);
+          SO[id+7*offset] = JW*(dudz - dwdx);
+          SO[id+8*offset] = JW*(dvdx - dudy);
+	}
+      }
+    }
+  }
+}
+
+@kernel void SijOijMag2(const dlong N,
+		        const dlong offset,
+		        @restrict const  dfloat * SOIJ,
+		        @restrict dfloat * OIOJSK,
+		        @restrict dfloat * MAGSIJ){
+
+  for(dlong n=0;n<N;++n;@tile(256,@outer,@inner)){
+    if(n<N){
+      const dfloat s1 = SOIJ[n+0*offset]; 
+      const dfloat s2 = SOIJ[n+1*offset]; 
+      const dfloat s3 = SOIJ[n+2*offset]; 
+      const dfloat s4 = SOIJ[n+3*offset]; 
+      const dfloat s5 = SOIJ[n+4*offset]; 
+      const dfloat s6 = SOIJ[n+5*offset];
+
+      const dfloat o1 = SOIJ[n+6*offset];
+      const dfloat o2 = SOIJ[n+7*offset];
+      const dfloat o3 = SOIJ[n+8*offset];
+
+      const dfloat magSij = s1*s1 + s2*s2 + s3*s3 + 2*(s4*s4 + s5*s5 + s6*s6);
+      const dfloat OiOjSk = s1*(o2*o2 + o3*o3) + s2*(o1*o1 + o3*o3) + s3*(o1*o1 + o2*o2)
+                            + 2*(o1*o2*s4 + o2*o3*s5 - o1*o3*s6); 
+
+      MAGSIJ[n] = 0.5*magSij;
+      OIOJSK[n] = OiOjSk;
+    }
+  }
+}
+
+@kernel void limit(const dlong N,
+		   @restrict dfloat * K,
+		   @restrict dfloat * TAU){
+
+  for(dlong n=0;n<N;++n;@tile(256,@outer,@inner)){
+    if(n<N){ 
+      dfloat k   = K[n];
+      dfloat tau = TAU[n];
+      if(k < 0) k = 0.01*abs(k);
+      if(tau < 0) tau = 0.01*abs(tau);
+
+      K[n] = k;
+      TAU[n] = tau;
+    }
+  }
+}
+
+@kernel void mue(const dlong N,
+		 const dlong offset,
+                 const dfloat rho,
+                 const dfloat mueLam,
+	         @restrict const dfloat * K,
+	         @restrict const dfloat * TAU,
+	         @restrict dfloat * MUET,
+	         @restrict dfloat * MUE,
+	         @restrict dfloat * DIFF){
+
+  for(dlong n=0;n<N;++n;@tile(256,@outer,@inner)){
+    if(n<N){ 
+      const dfloat k   = K[n];
+      const dfloat tau = TAU[n];
+      const dfloat mut = rho*p_alpinf_str*k*tau;
+
+      MUET[n] = mut;
+
+      MUE[n] = mueLam + mut;
+      DIFF[n + 0*offset] = mueLam + p_sigma_k*mut;
+      DIFF[n + 1*offset] = mueLam + p_sigma_tau*mut;
+    }
+  }
+}
diff --git a/okl/avg/avg.okl b/okl/plugins/avg.okl
similarity index 100%
rename from okl/avg/avg.okl
rename to okl/plugins/avg.okl
diff --git a/okl/velRecycling/velRecycling.okl b/okl/plugins/velRecycling.okl
similarity index 100%
rename from okl/velRecycling/velRecycling.okl
rename to okl/plugins/velRecycling.okl
diff --git a/src/core/cds.cpp b/src/core/cds.cpp
index 385bf2dd4..7fa4367b9 100644
--- a/src/core/cds.cpp
+++ b/src/core/cds.cpp
@@ -36,7 +36,6 @@ occa::memory cdsSolve(const int is, cds_t *cds, dfloat time)
   cds->o_wrk0.copyFrom(cds->o_S, cds->Ntotal*sizeof(dfloat), 0, is*cds->fieldOffset*sizeof(dfloat)); 
   if (solver->Nmasked) mesh->maskKernel(solver->Nmasked, solver->o_maskIds, cds->o_wrk0);
 
-  const dfloat lambda = 1; // dummy value not used if coeff is variable
   cds->Niter[is] = ellipticSolve(solver, cds->TOL, cds->o_wrk1, cds->o_wrk0);
 
   cds->helmholtzAddBCKernel(mesh->Nelements,
diff --git a/src/core/cds.h b/src/core/cds.h
index 2a2d74656..6db2cf0e3 100644
--- a/src/core/cds.h
+++ b/src/core/cds.h
@@ -100,7 +100,7 @@ typedef struct {
   occa::memory o_rho, o_diff;
 
   dfloat *cU, *cSd, *cS, *FS, *BF; 
-  occa::memory o_cU, o_cSd, o_cS, o_FS, o_BF;
+  occa::memory o_cU, o_cSd, o_cS, o_FS, o_BF, o_BFDiag;
 
   occa::memory o_wrk0, o_wrk1, o_wrk2, o_wrk3, o_wrk4, o_wrk5, o_wrk6;
 
diff --git a/src/core/ins.h b/src/core/ins.h
index 4c28f7c2d..dfa9c50b6 100644
--- a/src/core/ins.h
+++ b/src/core/ins.h
@@ -126,6 +126,7 @@ typedef struct {
   occa::kernel pqKernel;
   occa::kernel ncKernel;
 
+  occa::kernel scalarScaledAddKernel;
   occa::kernel scaledAddKernel;
   occa::kernel subCycleVolumeKernel,  subCycleCubatureVolumeKernel ;
   occa::kernel subCycleSurfaceKernel, subCycleCubatureSurfaceKernel;;
diff --git a/src/core/insSetup.cpp b/src/core/insSetup.cpp
index 7a3e59c15..dbaf16c91 100644
--- a/src/core/insSetup.cpp
+++ b/src/core/insSetup.cpp
@@ -634,6 +634,10 @@ ins_t *insSetup(MPI_Comm comm, setupAide &options, int buildOnly)
       ins->maxKernel =  
         mesh->device.buildKernel(fileName.c_str(), kernelName.c_str(), kernelInfo);
 
+      kernelName = "scalarScaledAdd";
+      ins->scalarScaledAddKernel =  
+        mesh->device.buildKernel(fileName.c_str(), kernelName.c_str(), kernelInfo);
+
       // ===========================================================================
 
       fileName = oklpath + "insFilterRT" + suffix + ".okl";
@@ -974,12 +978,12 @@ cds_t *cdsSetup(ins_t *ins, mesh_t *mesh, setupAide options, occa::properties &k
       cds->filterRTKernel =
         mesh->device.buildKernel(fileName.c_str(), kernelName.c_str(), kernelInfo);
 
-      if(cds->Nsubsteps){
-        fileName = install_dir + "/libparanumal/okl/scaledAdd.okl";
-        kernelName = "scaledAddwOffset";
-        cds->scaledAddKernel = 
-          mesh->device.buildKernel(fileName.c_str(), kernelName.c_str(), kernelInfo);
+      fileName = install_dir + "/libparanumal/okl/scaledAdd.okl";
+      kernelName = "scaledAddwOffset";
+      cds->scaledAddKernel = 
+        mesh->device.buildKernel(fileName.c_str(), kernelName.c_str(), kernelInfo);
 
+      if(cds->Nsubsteps){
         fileName = oklpath + "cdsSubCycle" + suffix + ".okl"; 
         kernelName = "cdsSubCycleStrongCubatureVolume" + suffix;
         cds->subCycleStrongCubatureVolumeKernel =  mesh->device.buildKernel(fileName, kernelName, kernelInfo);
diff --git a/src/core/parReader.cpp b/src/core/parReader.cpp
index 70cecae89..88c72eddc 100644
--- a/src/core/parReader.cpp
+++ b/src/core/parReader.cpp
@@ -378,47 +378,53 @@ libParanumal::setupAide parRead(std::string &setupFile, MPI_Comm comm)
     nscal++;
     isStart++;
 
-    options.setArgs("SCALAR00 PRECONDITIONER", "JACOBI");
-
-    double s_residualTol;
-    if(ini.extract("temperature", "residualtol", s_residualTol)) {
-      options.setArgs("SCALAR00 SOLVER TOLERANCE", to_string_f(s_residualTol));
-    }
-
-    if(ini.extract("temperature", "conductivity", sbuf)) {
-      int err = 0;
-      double diffusivity = te_interp(sbuf.c_str(), &err);
-      if(err) ABORT("Invalid expression for conductivity!");
-      if(diffusivity < 0) diffusivity = fabs(1/diffusivity);
-      options.setArgs("SCALAR00 DIFFUSIVITY", to_string_f(diffusivity));
-    } else {
-      if(!variableProperties)
-        ABORT("Cannot find mandatory parameter TEMPERATURE::conductivity!"); 
-    }
-
-    if(ini.extract("temperature", "rhocp", sbuf)) {
-      int err = 0;
-      double rhoCp = te_interp(sbuf.c_str(), &err);
-      if(err) ABORT("Invalid expression for rhoCp!");
-      options.setArgs("SCALAR00 DENSITY", to_string_f(rhoCp));
+    string solver;
+    ini.extract("temperature", "solver", solver);
+    if(solver == "none") {
+      options.setArgs("SCALAR00 SOLVER", "NONE");
     } else {
-      if(!variableProperties)
-        ABORT("Cannot find mandatory parameter TEMPERATURE::rhoCp!"); 
+      options.setArgs("SCALAR00 PRECONDITIONER", "JACOBI");
+ 
+      double s_residualTol;
+      if(ini.extract("temperature", "residualtol", s_residualTol)) {
+        options.setArgs("SCALAR00 SOLVER TOLERANCE", to_string_f(s_residualTol));
+      }
+ 
+      if(ini.extract("temperature", "conductivity", sbuf)) {
+        int err = 0;
+        double diffusivity = te_interp(sbuf.c_str(), &err);
+        if(err) ABORT("Invalid expression for conductivity!");
+        if(diffusivity < 0) diffusivity = fabs(1/diffusivity);
+        options.setArgs("SCALAR00 DIFFUSIVITY", to_string_f(diffusivity));
+      } else {
+        if(!variableProperties)
+          ABORT("Cannot find mandatory parameter TEMPERATURE::conductivity!"); 
+      }
+ 
+      if(ini.extract("temperature", "rhocp", sbuf)) {
+        int err = 0;
+        double rhoCp = te_interp(sbuf.c_str(), &err);
+        if(err) ABORT("Invalid expression for rhoCp!");
+        options.setArgs("SCALAR00 DENSITY", to_string_f(rhoCp));
+      } else {
+        if(!variableProperties)
+          ABORT("Cannot find mandatory parameter TEMPERATURE::rhoCp!"); 
+      }
+ 
+      string s_bcMap;
+      if(ini.extract("temperature", "boundarytypemap", s_bcMap)) {
+        std::vector<std::string> sList;
+        sList = serializeString(s_bcMap);
+        bcMap::setup(sList, "scalar00");
+      } else {
+        ABORT("Cannot find mandatory parameter TEMPERATURE::boundaryTypeMap!");
+      } 
     }
-
-    string s_bcMap;
-    if(ini.extract("temperature", "boundarytypemap", s_bcMap)) {
-      std::vector<std::string> sList;
-      sList = serializeString(s_bcMap);
-      bcMap::setup(sList, "scalar00");
-    } else {
-      ABORT("Cannot find mandatory parameter TEMPERATURE::boundaryTypeMap!");
-    } 
-
-  } else {
-    if(equation == "lowmachns") 
-      ABORT("PROBLEMTYPE::equation = lowMachNS requires solving for temperature!");
   }
+ 
+  if(equation == "lowmachns" && ini.sections.count("temperature") == 0) 
+    ABORT("PROBLEMTYPE::equation = lowMachNS requires solving for temperature!");
+    
   //
   for (auto & sec : ini.sections) {
     string key = sec.first;
diff --git a/src/core/runTime.cpp b/src/core/runTime.cpp
index aed2bf62c..3261bebf5 100644
--- a/src/core/runTime.cpp
+++ b/src/core/runTime.cpp
@@ -175,22 +175,27 @@ void makeq(ins_t *ins, dfloat time, occa::memory o_BF)
 
     mesh_t *mesh;
     (is) ? mesh = cds->meshV : mesh = cds->mesh;
+    const dlong sOffset = is*cds->fieldOffset;
+    occa::memory o_adv = cds->o_wrk0;
 
     if(cds->options.compareArgs("FILTER STABILIZATION", "RELAXATION"))
       cds->filterRTKernel(
            cds->meshV->Nelements,
            ins->o_filterMT,
            ins->filterS,
-           is*cds->fieldOffset,
+           sOffset,
            cds->o_rho,
            cds->o_S,
            cds->o_FS);
 
-    occa::memory o_adv = cds->o_wrk0;
     if(cds->Nsubsteps) {
+
       o_adv = scalarStrongSubCycle(cds, time, is, cds->o_U, cds->o_S);
+
     } else {
-      if(cds->options.compareArgs("ADVECTION TYPE", "CUBATURE")) {
+
+     //ins->setScalarKernel(cds->fieldOffset, 1.0, cds->o_rho);
+     if(cds->options.compareArgs("ADVECTION TYPE", "CUBATURE")) {
         cds->advectionStrongCubatureVolumeKernel(
              cds->meshV->Nelements,
              mesh->o_vgeo,
@@ -199,7 +204,7 @@ void makeq(ins_t *ins, dfloat time, occa::memory o_BF)
              mesh->o_cubInterpT,
              mesh->o_cubProjectT,
              cds->vFieldOffset,
-             is*cds->fieldOffset,
+             sOffset,
              cds->o_Ue,
              cds->o_S,
              cds->o_rho,
@@ -210,20 +215,22 @@ void makeq(ins_t *ins, dfloat time, occa::memory o_BF)
              mesh->o_vgeo,
              mesh->o_Dmatrices,
              cds->vFieldOffset,
-             is*cds->fieldOffset,
+             sOffset,
              cds->o_Ue,
              cds->o_S,
              cds->o_rho,
              cds->o_wrk0);
       }
+
       ins->scaledAddKernel(
            cds->meshV->Nelements*cds->meshV->Np, 
            -1.0, 
            0*cds->fieldOffset, 
            cds->o_wrk0, 
            1.0, 
-           is*cds->fieldOffset, 
+           sOffset, 
            cds->o_FS);
+
     }
 
     cds->sumMakefKernel(
@@ -235,7 +242,7 @@ void makeq(ins_t *ins, dfloat time, occa::memory o_BF)
          cds->o_extbdfB,
          cds->o_extbdfC,
          cds->fieldOffset*cds->NSfields,
-         is*cds->fieldOffset,
+         sOffset,
          cds->o_S,
          o_adv,
          cds->o_FS,
@@ -276,9 +283,15 @@ void scalarSolve(ins_t *ins, dfloat time, dfloat dt, occa::memory o_S)
          cds->o_rho,
          cds->o_ellipticCoeff);
 
-    if(udf.sEqnAddToLhs) {
-      occa::memory o_LHS = cds->o_ellipticCoeff + cds->fieldOffset*sizeof(dfloat);
-      udf.sEqnAddToLhs(ins, is, time+dt, o_LHS);
+    if(cds->o_BFDiag.ptr()) {
+      cds->scaledAddKernel(
+           cds->Nlocal, 
+           1.0, 
+           is*cds->fieldOffset, 
+           cds->o_BFDiag, 
+           1.0, 
+           cds->fieldOffset, 
+           cds->o_ellipticCoeff); 
     }
 
     occa::memory o_Snew = cdsSolve(is, cds, time+dt);
@@ -394,10 +407,7 @@ void fluidSolve(ins_t *ins, dfloat time, dfloat dt, occa::memory o_U)
     ins->o_mue,
     ins->o_rho,
     ins->o_ellipticCoeff);
-  if(udf.uEqnAddToLhs) {
-    occa::memory o_LHS = ins->o_ellipticCoeff + ins->fieldOffset*sizeof(dfloat);
-    udf.uEqnAddToLhs(ins, time+dt, o_LHS);
-  }
+
   occa::memory o_Unew = tombo::velocitySolve(ins, time+dt);
   for (int s=ins->Nstages;s>1;s--) {
     o_U.copyFrom(
diff --git a/src/nekInterface/nekInterface.f b/src/nekInterface/nekInterface.f
index 273faaf45..604a0c1b2 100644
--- a/src/nekInterface/nekInterface.f
+++ b/src/nekInterface/nekInterface.f
@@ -256,8 +256,9 @@ subroutine nekf_resetio()
       real ts
       integer npscals, p63s
       logical ifxyos, ifvos, ifpos, iftos, ifpscos(ldimt1)
-      common /ifos/ fxyos, ifvos, ifpos, iftos, ifpscos, npscals, p63s,
-     $              ts
+      common /ros/  ts
+      common /ios/  npscals, p63s
+      common /ifos/ ifxyos, ifvos, ifpos, iftos, ifpscos
 
       time = ts
 
@@ -287,8 +288,10 @@ subroutine nekf_setio(ttime, xo, vo, po, so, ns, fp64)
       real ts
       integer npscals, p63s
       logical ifxyos, ifvos, ifpos, iftos, ifpscos(ldimt1)
-      common /ifos/ fxyos, ifvos, ifpos, iftos, ifpscos, npscals, p63s,
-     $              ts
+
+      common /ros/  ts
+      common /ios/  npscals, p63s
+      common /ifos/ ifxyos, ifvos, ifpos, iftos, ifpscos
 
       if(ns.gt.ldimt) call exitti('nekf_setifo: ns > ldimt$',ns) 
 
diff --git a/src/core/nrs.hpp b/src/nrs.hpp
similarity index 100%
rename from src/core/nrs.hpp
rename to src/nrs.hpp
diff --git a/src/udf/udf.hpp b/src/udf/udf.hpp
index 2eb84bf04..3f7eafae7 100644
--- a/src/udf/udf.hpp
+++ b/src/udf/udf.hpp
@@ -17,8 +17,6 @@ typedef void (*udfsetup)(ins_t *ins);
 typedef void (*udfloadKernels)(ins_t *ins);
 typedef void (*udfexecuteStep)(ins_t *ins, dfloat time, int tstep);
 
-typedef void (*udfuEqnAddToLhs)(ins_t *ins, dfloat time, occa::memory o_LHS);
-typedef void (*udfsEqnAddToLhs)(ins_t *ins, int i, dfloat time, occa::memory o_LHS);
 typedef void (*udfuEqnSource)(ins_t *ins, dfloat time, occa::memory o_U, occa::memory o_FU);
 typedef void (*udfsEqnSource)(ins_t *ins, dfloat time, occa::memory o_S, occa::memory o_SU);
 typedef void (*udfproperties)(ins_t *ins, dfloat time, occa::memory o_U, 
@@ -32,8 +30,6 @@ typedef struct
   udfsetup setup;
   udfloadKernels loadKernels;
   udfexecuteStep executeStep;
-  udfuEqnAddToLhs uEqnAddToLhs;
-  udfsEqnAddToLhs sEqnAddToLhs;
   udfuEqnSource uEqnSource;
   udfsEqnSource sEqnSource;
   udfproperties properties;