Importat latest changes (#86)

- use globalIds from nek
- add multiplicity cross-check
- add ktauChannel example
- add ktau RANS model
- add variable viscocity term to pressure rhs (but no stress ellipticOperator)

3rd_party/libparanumal/libparanumal.makefile

@@ -79,7 +79,6 @@ $(HDRDIR)/src/meshLoadReferenceNodesTet3D.o \
 $(HDRDIR)/src/meshOccaSetup2D.o \
 $(HDRDIR)/src/meshOccaSetup3D.o \
 $(HDRDIR)/src/meshOccaSetupQuad3D.o \
-$(HDRDIR)/src/meshParallelConnectNodes.o \
 $(HDRDIR)/src/meshParallelConnectOpt.o \
 $(HDRDIR)/src/meshParallelConsecutiveGlobalNumbering.o\
 $(HDRDIR)/src/meshParallelGatherScatterSetup.o \
@@ -123,7 +122,8 @@ $(HDRDIR)/src/hash.o\
 $(HDRDIR)/src/setupAide.o \
 $(HDRDIR)/src/readArray.o\
 $(HDRDIR)/src/occaHostMallocPinned.o \
-$(HDRDIR)/src/timer.o
+$(HDRDIR)/src/timer.o \
+#$(HDRDIR)/src/meshParallelConnectNodes.o
 
 ifeq ($(OS),Windows_NT)
     detected_OS := Windows

CMakeLists.txt

@@ -322,6 +322,7 @@ set(SRC
     src/mesh/meshSetup.cpp
     src/mesh/meshNekReader.cpp
     src/mesh/meshPhysicalNodesHex3D.cpp
+    src/mesh/meshParallelConnectNodes.cpp
     src/core/occaDeviceConfig.cpp
     src/nekInterface/nekInterfaceAdapter.cpp
     src/core/parReader.cpp
@@ -330,6 +331,7 @@ set(SRC
     src/core/timer.cpp
     src/plugins/avg.cpp
     src/plugins/velRecycling.cpp
+    src/plugins/RANSktau.cpp
     ## To get fortran flags
     src/dummy.f
 )

examples/ktauChannel/channel.oudf

@@ -0,0 +1,18 @@
+@kernel void cfill(const dlong Ntotal,
+                    const dlong offset,
+                    const dfloat CONST,
+                    @restrict dfloat * NU) 
+{
+  for(dlong n=0;n<Ntotal;++n;@tile(256,@outer,@inner)){
+    if(n<Ntotal){
+      const dlong id = n + offset;
+      NU[id] = CONST;
+    }
+  }
+}
+
+void cdsDirichletConditions3D(bcData *bc)
+{
+  bc->sP = 0;
+  if(bc->scalarId == 0) bc->sP = 0; 
+}

examples/ktauChannel/channel.par

@@ -0,0 +1,43 @@
+[OCCA]
+backend = CUDA
+deviceNumber = LOCAL-RANK
+
+[GENERAL]
+polynomialOrder = 7
+#startFrom = r.fld+time=0
+stopAt = endTime
+endTime = 500
+dt = 2e-02
+timeStepper = tombo2
+
+writeControl = runTime
+writeInterval = 50
+
+[PRESSURE]
+preconditioner = semg_amg
+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-06
+rho = 1.0
+diffusivity = -43500. 
+
+[SCALAR02] # tau
+boundaryTypeMap = inlet, insulated 
+residualTol = 1e-06
+rho = 1.0
+diffusivity = -43500. 

examples/ktauChannel/channel.udf

@@ -0,0 +1,75 @@
+//
+// nekRS User Defined File
+//
+#include <math.h>
+#include "udf.hpp"
+#include "plugins/RANSktau.hpp"
+
+/* User Functions */
+
+static dfloat rho, mueLam;
+occa::kernel userfKernel; 
+
+void userf(ins_t *ins, dfloat time, occa::memory o_U, occa::memory o_FU)
+{
+  const dfloat Re_tau = 2000.0; 
+  const dfloat Re_b   = rho/mueLam; 
+  const dfloat DPDX   = (Re_tau/Re_b)*(Re_tau/Re_b);
+  userfKernel(ins->Nlocal, 0*ins->fieldOffset, DPDX, o_FU);
+}
+
+void userq(ins_t *ins, dfloat time, occa::memory o_S, occa::memory o_FS)
+{
+  mesh_t *mesh = ins->mesh;
+  cds_t *cds   = ins->cds;
+
+  RANSktau::updateSourceTerms();
+}
+
+void uservp(ins_t *ins, dfloat time, occa::memory o_U, occa::memory o_S,
+            occa::memory o_UProp, occa::memory o_SProp)
+{
+  mesh_t *mesh = ins->mesh;
+  cds_t *cds   = ins->cds;
+
+  RANSktau::updateProperties();
+
+  dfloat conductivity;
+  ins->options.getArgs("SCALAR00 DIFFUSIVITY", conductivity);
+  const dfloat Pr_t = 0.7;
+  occa::memory o_mue_t = RANSktau::o_mue_t();
+  occa::memory o_temp_mue = cds->o_diff + 0*cds->fieldOffset*sizeof(dfloat);
+  ins->scalarScaledAddKernel(ins->Nlocal, mueLam, 1/Pr_t, o_mue_t, o_temp_mue); 
+}
+
+void UDF_LoadKernels(ins_t *ins)
+{
+  userfKernel = udfBuildKernel(ins, "cfill");
+  RANSktau::buildKernel(ins);
+}
+
+void UDF_Setup(ins_t *ins)
+{
+  mesh_t *mesh = ins->mesh;
+  cds_t *cds = ins->cds;
+
+  // get IC from nek
+  if (!ins->readRestartFile) nek_copyTo(ins->startTime);
+
+  udf.properties = &uservp;
+  udf.uEqnSource = &userf;
+  udf.sEqnSource = &userq;
+
+  const int scalarFieldStart = 1;
+  ins->options.getArgs("VISCOSITY", mueLam); 
+  ins->options.getArgs("DENSITY", rho); 
+  RANSktau::setup(ins, mueLam, rho, scalarFieldStart);
+}
+
+void UDF_ExecuteStep(ins_t *ins, dfloat time, int tstep)
+{
+  if (ins->isOutputStep) {
+    nek_ocopyFrom(time, tstep);
+    nek_userchk();
+  }
+}

examples/ktauChannel/channel.usr

@@ -0,0 +1,115 @@
+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'
+
+      real kmax
+
+      Re    = 1/param(2)
+      kmax  = 3.9201E-03
+      omax  = 0.5
+
+      ux   = 1.0  
+      uy   = 0.0 
+      uz   = 0.0
+
+      yd = 1 + y
+      tau = 0.0
+      if(yd.ne.0) then
+         omeg = omax + 6./Re/0.075/yd**2
+         tau = 1/omeg
+      endif
+
+      if (ifield.eq.2) then
+         temp = 1.0
+      elseif (ifield.eq.3) then
+         temp = kmax
+      elseif (ifield.eq.4) then
+         temp = tau
+      endif
+
+      return
+      end
+c-----------------------------------------------------------------------
+      subroutine userchk
+      include 'SIZE'
+      include 'TOTAL'
+
+c      ubar  = glsc2(vx,bm1,nx1*ny1*nz1*nelt)/volvm1
+c      if (nid.eq.0) write(6,*) 'ubar=', ubar
+
+      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 = 1.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
+         cbc(5,iel,1) = 'P  '
+         cbc(6,iel,1) = 'P  '
+      do ifc=1,2*ndim
+         cbc(ifc,iel,2)  = cbc(ifc,iel,1)
+         if (cbc(ifc,iel,1) .eq. 'W  ') cbc(ifc,iel,2) = 't  ' 
+         if (cbc(ifc,iel,1) .eq. 'SYM') cbc(ifc,iel,2) = 'I  ' 
+         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-----------------------------------------------------------------------
+      subroutine usrsetvert(glo_num,nel,nx,ny,nz) ! to modify glo_num
+      integer*8 glo_num(1)
+
+      ! kludge for periodic bc in z
+      nxy  = nx*ny
+      nxyz = nx*ny*nz
+      do iel = 1,nel
+         ioff = nxyz*(iel-1)
+         do ixy = 1,nxy
+            glo_num(ioff + nxy*(nz-1) + ixy) = glo_num(ioff + ixy)
+         enddo
+      enddo
+
+      return
+      end
+c-----------------------------------------------------------------------

examples/ktauChannel/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  -1                                          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,E  ,E                                bc's  (3 chars each!)

okl/core/insCurlHex3D.okl

@@ -28,7 +28,6 @@ SOFTWARE.
 @kernel void insCurlHex3D(const dlong Nelements, 
                            @restrict const  dfloat *  vgeo,
                            @restrict const  dfloat *  const D,
-                           @restrict const  dfloat *  const M,
                            const dlong offset,
                            @restrict const  dfloat *  U,
                            @restrict dfloat *  W){

okl/core/insDivergenceHex3D.okl

@@ -165,9 +165,7 @@
 @kernel void insDivergenceSurfaceTOMBOHex3D(const dlong Nelements,
 					    @restrict const  dfloat *  vgeo,
 					    @restrict const  dfloat *  sgeo,
-					    @restrict const  dfloat *  LIFTT,
 					    @restrict const  dlong  *  vmapM,
-					    @restrict const  dlong  *  mapB,
 					    @restrict const  int    *  EToBM,
 					    @restrict const  int    *  EToB,
 					    const dfloat time,

okl/core/insAxHex3D.okl → okl/core/insPressureAxHex3D.okl

@@ -25,12 +25,11 @@
 */
 
 // Computes local [lap(u) + lambda*u] = [-(grad(u), grad(phi)) + lambda*u] operation
-@kernel void insAxHex3D(const dlong Nelements,
+@kernel void insPressureAxHex3D(const dlong Nelements,
                         const dlong offset,
 		        @restrict const  dfloat *  ggeo,
 		        @restrict const  dfloat *  D,
 		        @restrict const  dfloat *  S,
-		        @restrict const  dfloat *  MM,
 		        @restrict const  dfloat *  q,
 		        @restrict const  dfloat *  lambda,
 		        @restrict dfloat *  Aq){