adding initial work on RANS tutorial

doc/source/tutorials/ktauTutorial/channel.oudf

@@ -0,0 +1,20 @@
+@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];
+    }
+  }
+}
+
+void scalarDirichletConditions(bcData *bc)
+{
+  bc->s = 0;
+  if (bc->scalarId == 0) {
+    bc->s = 0;
+  }
+}

doc/source/tutorials/ktauTutorial/channel.par

@@ -0,0 +1,43 @@
+[GENERAL]
+polynomialOrder = 7
+#startFrom = "r.fld"+time=0
+stopAt = endTime
+endTime = 50
+dt = 5e-02
+timeStepper = tombo2
+
+checkpointControl = simulationTime 
+checkpointInterval = 50
+
+constFlowRate = meanVelocity=1.0 + direction=X
+
+[PROBLEMTYPE]
+equation = navierStokes+variableViscosity
+
+[PRESSURE]
+residualTol = 1e-04
+
+[VELOCITY]
+boundaryTypeMap = wall, slipY
+residualTol = 1e-06
+density = 1.0
+viscosity = -43500. 
+
+[TEMPERATURE]
+#solver = none
+boundaryTypeMap = inlet, insulated 
+residualTol = 1e-06
+rhoCp = 1.0
+conductivity = -43500. 
+
+[SCALAR01] # k
+boundaryTypeMap = inlet, insulated 
+residualTol = 1e-08
+rho = 1.0
+diffusivity = -43500. 
+
+[SCALAR02] # tau
+boundaryTypeMap = inlet, insulated 
+residualTol = 1e-06
+rho = 1.0
+diffusivity = -43500. 

doc/source/tutorials/ktauTutorial/channel.udf

@@ -0,0 +1,128 @@
+#include "RANSktau.hpp"
+
+#include "ci.inc"
+
+static dfloat rho, mueLam;
+
+#ifdef __okl__
+
+#include "channel.oudf"
+
+#endif
+
+void extractLine(nrs_t *nrs, double time)
+{
+  const auto np = (platform->comm.mpiRank == 0) ? 200 : 0;
+  const auto offset = np;
+
+  static pointInterpolation_t *interpolator = nullptr;
+  static std::vector<dfloat> xp, yp, zp;
+  static deviceMemory<dfloat> o_Up;
+
+  if (!interpolator) {
+    auto mesh = nrs->meshV;
+    const auto yMin = platform->linAlg->min(mesh->Nlocal, mesh->o_y, platform->comm.mpiComm);
+    const auto yMax = platform->linAlg->max(mesh->Nlocal, mesh->o_y, platform->comm.mpiComm);
+
+    if (np) {
+      const auto x0 = 7.0;
+      const auto z0 = 0.5;
+
+      xp.push_back(x0);
+      yp.push_back(yMin);
+      zp.push_back(z0);
+
+      const auto betaY = 2.2;
+      const auto dy = (yMax - yMin) / (np - 1);
+      for (int i = 1; i < np - 1; i++) {
+        xp.push_back(x0);
+        yp.push_back(tanh(betaY * (i * dy - 1)) / tanh(betaY));
+        zp.push_back(z0);
+      }
+
+      xp.push_back(x0);
+      yp.push_back(yMax);
+      zp.push_back(z0);
+      o_Up.resize(nrs->NVfields * offset);
+    }
+
+    interpolator = new pointInterpolation_t(nrs->meshV);
+    interpolator->setPoints(np, xp.data(), yp.data(), zp.data());
+    interpolator->find();
+  }
+
+  interpolator->eval(nrs->NVfields, nrs->fieldOffset, nrs->o_U, offset, o_Up);
+
+  if (platform->comm.mpiRank == 0) {
+    std::vector<dfloat> Up(nrs->NVfields * np);
+    o_Up.copyTo(Up);
+
+    std::ofstream f("profile.dat", std::ios::app);
+    for (int i = 0; i < np; i++) {
+      f << std::scientific << time << " " << xp[i] << " " << yp[i] << " " << zp[i] << " "
+        << Up[i + 0 * offset] << " " << Up[i + 1 * offset] << " " << Up[i + 2 * offset] << std::endl;
+    }
+    f.close();
+  }
+}
+
+void userq(double time)
+{
+  RANSktau::updateSourceTerms();
+}
+
+void uservp(double time)
+{
+  auto mesh = nrs->meshV;
+  auto cds = nrs->cds;
+
+  RANSktau::updateProperties();
+
+  dfloat conductivity;
+  platform->options.getArgs("SCALAR00 DIFFUSIVITY", conductivity);
+  const dfloat Pr_t = 0.7;
+  auto o_mue_t = RANSktau::o_mue_t();
+  auto o_temp_mue = cds->o_diff + 0 * cds->fieldOffset[0];
+  scalarScaledAdd(mesh->Nlocal, conductivity, 1 / Pr_t, o_mue_t, o_temp_mue);
+}
+
+void UDF_LoadKernels(deviceKernelProperties& kernelInfo)
+{
+#if 0
+  {
+    auto props = kernelInfo;
+    props.define("p_sigma_k") = 0.6;
+    RANSktau::buildKernel(props);
+  }
+#endif
+}
+
+void UDF_Setup0(MPI_Comm comm, setupAide &options)
+{
+  options.getArgs("CI-MODE", ciMode);
+  if (ciMode) {
+    ciSetup(comm, options);
+  }
+}
+
+void UDF_Setup()
+{
+  nrs->userProperties = &uservp;
+  nrs->userScalarSource = &userq;
+
+  const int scalarFieldStart = 1;
+  platform->options.getArgs("VISCOSITY", mueLam);
+  platform->options.getArgs("DENSITY", rho);
+
+  RANSktau::setup(mueLam, rho, scalarFieldStart);
+}
+
+void UDF_ExecuteStep(double time, int tstep)
+{
+  if (ciMode) {
+    ciTestErrors(nrs, time, tstep);
+  }
+  if (nrs->isCheckpointStep) {
+    extractLine(nrs, time);
+  }
+}

doc/source/tutorials/ktauTutorial/channel.usr

@@ -0,0 +1,79 @@
+C
+C  USER SPECIFIED ROUTINES: 
+C
+C     - boundary conditions 
+C     - initial conditions  
+C     - variable properties 
+C     - forcing function for fluid (f)
+C     - forcing function for passive scalar (q)
+C     - general purpose routine for checking errors etc.        
+C
+c-----------------------------------------------------------------------
+      subroutine useric (ix,iy,iz,ieg)
+      include 'SIZE'
+      include 'TOTAL'
+      include 'NEKUSE'
+
+      ux   = 1.0  
+      uy   = 0.0 
+      uz   = 0.0
+
+      if (ifield.eq.2) then
+         temp = 1.0
+      elseif (ifield.eq.3) then
+         temp = 0.01
+      elseif (ifield.eq.4) then
+         temp = 0.1
+      endif
+
+      return
+      end
+c-----------------------------------------------------------------------
+      subroutine userchk
+      include 'SIZE'
+      include 'TOTAL'
+
+      return
+      end
+c-----------------------------------------------------------------------
+      subroutine usrdat   ! This routine to modify element vertices
+      include 'SIZE'      ! _before_ mesh is generated, which 
+      include 'TOTAL'     ! guarantees GLL mapping of mesh.
+
+      return
+      end
+c-----------------------------------------------------------------------
+      subroutine usrdat2()  ! This routine to modify mesh coordinates
+      include 'SIZE'
+      include 'TOTAL'
+
+      parameter(BETAM = 2.8)
+
+      call rescale_x(xm1, 0.0,8.0)
+      call rescale_x(ym1,-1.0,0.0)
+      call rescale_x(zm1, 0.0,1.0)
+
+      ntot = nx1*ny1*nz1*nelt
+
+      do i=1,ntot
+         ym1(i,1,1,1) = tanh(BETAM*ym1(i,1,1,1))/tanh(BETAM)
+      enddo
+
+      do iel=1,nelt
+      do ifc=1,2*ndim
+        boundaryID(ifc,iel) = 0 
+         if (cbc(ifc,iel,1) .eq. 'W  ') boundaryID(ifc,iel) = 1
+         if (cbc(ifc,iel,1) .eq. 'SYM') boundaryID(ifc,iel) = 2
+      enddo
+      enddo
+
+      return
+      end
+c-----------------------------------------------------------------------
+      subroutine usrdat3
+      include 'SIZE'
+      include 'TOTAL'
+
+      return
+      end
+c-----------------------------------------------------------------------

doc/source/tutorials/ktauTutorial/ci.inc

@@ -0,0 +1,74 @@
+#include <math.h>
+#include <vector>
+
+static int ciMode = 0;
+
+#define EPS 4e-3
+
+void ciSetup(MPI_Comm comm, setupAide &options)
+{
+  options.setArgs("POLYNOMIAL DEGREE", std::string("7"));
+  options.setArgs("CUBATURE POLYNOMIAL DEGREE", std::string("10"));
+  options.setArgs("CONSTANT FLOW RATE", "TRUE");
+  options.setArgs("CONSTANT FLOW DIRECTION", "X");
+  options.setArgs("RESTART FILE NAME", std::string("r.fld+time=0"));
+  options.setArgs("SOLUTION OUTPUT INTERVAL", "-1");
+  options.setArgs("VISCOSITY", std::to_string(1.0 / 43500.0));
+  options.setArgs("DENSITY", std::string("1.0"));
+  options.setArgs("NUMBER TIMESTEPS", std::string("-1"));
+  options.setArgs("END TIME", std::string("10"));
+  options.setArgs("TIME INTEGRATOR", "TOMBO2");
+  options.setArgs("ADVECTION TYPE", "CONVECTIVE+CUBATURE");
+  options.setArgs("DT", std::string("5e-2"));
+  options.setArgs("SUBCYCLING STEPS", std::string("0"));
+
+  options.setArgs("PRESSURE KRYLOV SOLVER", "PGMRES+FLEXIBLE");
+  options.setArgs("PRESSURE SOLVER TOLERANCE", std::string("1e-4"));
+  options.setArgs("PRESSURE PRECONDITIONER", "MULTIGRID");
+  options.setArgs("PRESSURE MULTIGRID COARSE SOLVE", "TRUE");
+  options.setArgs("PRESSURE MULTIGRID SMOOTHER", "FOURTHOPTCHEBYSHEV+ASM");
+  options.setArgs("PRESSURE MULTIGRID CHEBYSHEV DEGREE", std::string("3"));
+  options.setArgs("PRESSURE MULTIGRID CHEBYSHEV MAX EIGENVALUE BOUND FACTOR", std::string("1.1"));
+  options.setArgs("PRESSURE INITIAL GUESS", "PROJECTION-ACONJ");
+
+  options.setArgs("VELOCITY SOLVER TOLERANCE", std::string("1e-6"));
+  options.setArgs("VELOCITY BLOCK SOLVER", "TRUE");
+  options.setArgs("VELOCITY INITIAL GUESS", "EXTRAPOLATION");
+
+  options.setArgs("SCALAR00 SOLVER TOLERANCE", std::string("1e-6"));
+  options.setArgs("SCALAR00 INITIAL GUESS", "EXTRAPOLATION");
+
+  options.setArgs("SCALAR01 SOLVER TOLERANCE", std::string("1e-8"));
+  options.setArgs("SCALAR01 INITIAL GUESS", "EXTRAPOLATION");
+
+  options.setArgs("SCALAR02 SOLVER TOLERANCE", std::string("1e-6"));
+  options.setArgs("SCALAR02 INITIAL GUESS", "EXTRAPOLATION");
+}
+
+void ciTestErrors(nrs_t *nrs, double time, int tstep)
+{
+  if (!nrs->lastStep) {
+    return;
+  }
+
+  mesh_t *mesh = nrs->meshV;
+
+  std::vector<int> bidWall = {1};
+  occa::memory o_bidWall = platform->device.malloc<int>(bidWall.size(), bidWall.data());
+
+  auto o_Sij = nrs->strainRate();
+
+  auto forces = nrs->aeroForces(bidWall.size(), o_bidWall, o_Sij);
+  o_Sij.free();
+
+  auto o_tmp = platform->o_memPool.reserve<dfloat>(mesh->Nlocal);
+  platform->linAlg->fill(mesh->Nlocal, 1.0, o_tmp);
+  const auto areaWall = mesh->surfaceAreaMultiplyIntegrate(bidWall.size(), o_bidWall, o_tmp);
+
+  // https://turbulence.oden.utexas.edu/channel2015/data/LM_Channel_2000_mean_prof.dat
+  const auto utauRef = 4.58794e-02;
+  const auto utau = sqrt(std::abs(forces->forceViscous()[0]) / areaWall.at(0));
+  const auto utauErr = std::abs((utau - utauRef) / utauRef);
+
+  CiEvalTest("tau", fabs(utauErr) < EPS, "with relative LinfErr:" + to_string_f(utauErr));
+}

doc/source/tutorials/ktauTutorial/input.box

@@ -0,0 +1,22 @@
+base.rea
+-3                     spatial dimension  ( < 0 --> generate .rea/.re2 pair)
+1                      number of fields
+#=======================================================================
+#
+#    Example of .box file for Taylor-Green
+#
+#    If nelx (y or z) < 0, then genbox automatically generates the
+#                          grid spacing in the x (y or z) direction
+#                          with a geometric ratio given by "ratio".
+#                          ( ratio=1 implies uniform spacing )
+#
+#    Note that the character bcs _must_ have 3 spaces.
+#
+#=======================================================================
+#
+Box
+-5  -12  -4                                          nelx,nely,nelz for Box
+0 8 1.                                               x0,x1,gain  (rescaled in usrdat)
+0 1 1.                                               y0,y1,gain  (rescaled in usrdat)
+0 1 1.                                               z0,z1,gain
+P  ,P  ,W  ,SYM,P  ,P                                bc's  (3 chars each!)