Add Persson's AVM for scalars (#295)

CMakeLists.txt

@@ -199,7 +199,8 @@ set(SRC
     src/core/utils/tinyexpr.c
     src/core/setupAide.cpp
     src/core/cfl.cpp
-    src/core/filter.cpp
+    src/regularization/filter.cpp
+    src/regularization/avm.cpp
     src/core/bcMap.cpp
     src/core/setup.cpp
     src/timeStepper/timeStepper.cpp
@@ -252,6 +253,7 @@ target_include_directories(nekrs-lib
   src/lib
   src/io
   src/udf
+  src/regularization
   src/linAlg
   src/timeStepper
   src/lns

examples/double_shear/double_shear.oudf

@@ -0,0 +1,12 @@
+// Boundary conditions
+void velocityDirichletConditions(bcData *bc)
+{                                                                        
+  bc->u = 0.0;
+  bc->v = 0.0;
+  bc->w = 0.0;
+}
+
+void scalarDirichletConditions(bcData *bc)
+{                                                                        
+  bc->s = 0.0;
+}

examples/double_shear/double_shear.par

@@ -0,0 +1,39 @@
+[OCCA]
+backend = CUDA
+deviceNumber = LOCAL-RANK
+
+[GENERAL] 
+numSteps = 2000
+polynomialOrder = 5
+dt = 1e-03
+timestepper = bdf3
+
+[PRESSURE]
+residualTol = 1e-08
+
+[VELOCITY]
+residualTol = 1e-8
+density = 1.
+viscosity = -100000
+
+[TEMPERATURE]
+residualProj = true
+residualTol = 1e-8
+rhocp = 1
+conductivity = -100000
+regularization=none
+
+[SCALAR01]
+residualProj = true
+residualTol = 1e-8
+rho = 1
+diffusivity = -100000
+regularization=none
+
+[SCALAR02]
+dt = 1e-03
+residualProj = true
+residualTol = 1e-8
+rho = 1
+diffusivity = -100000
+regularization=avm+c0

examples/double_shear/double_shear.udf

@@ -0,0 +1,25 @@
+//
+// nekRS User Defined File
+//
+
+#include "udf.hpp"
+
+void UDF_LoadKernels(nrs_t *nrs)
+{
+}
+
+void UDF_Setup0(MPI_Comm comm, setupAide &options)
+{
+}
+
+void UDF_Setup(nrs_t *nrs)
+{
+}
+
+void UDF_ExecuteStep(nrs_t *nrs, dfloat time, int tstep)
+{
+  if (nrs->isOutputStep) {
+    nek::ocopyToNek(time, tstep);
+    nek::userchk();
+  }
+}

examples/double_shear/double_shear.usr

@@ -0,0 +1,184 @@
+c-----------------------------------------------------------------------
+c    2D Double Shear layer with doubly periodic BCs.
+c    This a test for the high-pass filtered relaxation-term forcing 
+c    (hprtf) stabilization method. 
+c    The test case used is from: 
+c    https://doi.org/10.1016/S0764-4442(00)01763-8
+c    with a reduced resolution to create a smaller test case.
+c    Without over-integration the simulation is unstable.
+c    The following stabilization parameters are used for the test case:
+c    
+c    The short test checks the total kinetic energy in the flow at the 
+c    end of simulation as a quantitative test.
+c    The simulation undergoes numerical instbility if no stabilization is used.
+c-----------------------------------------------------------------------
+      subroutine uservp (ix,iy,iz,ieg)
+      include 'SIZE'
+      include 'NEKUSE'
+
+      return
+      end
+c-----------------------------------------------------------------------
+      subroutine userf  (ix,iy,iz,ieg)
+
+      include 'SIZE'
+      include 'NEKUSE'
+      include 'PARALLEL'
+
+      integer ie,ix,iy,iz,ieg
+
+      ie     = gllel(ieg)
+
+      ffx = 0.
+      ffy = 0.
+      ffz = 0.         
+
+      return
+      end
+c-----------------------------------------------------------------------
+      subroutine userq  (ix,iy,iz,ieg)
+      include 'SIZE'
+      include 'TOTAL'
+      include 'NEKUSE'
+
+      qvol   = 0.0
+      source = 0.0
+
+      return
+      end
+c-----------------------------------------------------------------------
+      subroutine userchk
+
+      include 'SIZE'
+      include 'SOLN'
+      include 'INPUT'
+      include 'TSTEP'
+
+      real vort,w1,w2
+
+      integer lt,nt
+      parameter (lt=lx1*ly1*lz1*lelv)
+      common /scrns/ vort(lt,3), w1(lt), w2(lt)
+
+      real en, enx, eny
+      real gl2norm
+
+      ifto =.true.
+
+      call comp_vort3(vort,w1,w2,vx,vy,vz)
+      call copy(T,vort,lt)
+
+      nt=nx1*ny1*nz1*nelv
+
+      if (istep.eq.nsteps.or.lastep.eq.1) then
+        enx = gl2norm(vx,nt)
+        eny = gl2norm(vy,nt)
+        en = enx**2 + eny**2    
+        if (nid.eq.0) write(6,*) 'L2 norm: ', en
+      endif
+
+      return
+      end
+c-----------------------------------------------------------------------
+      subroutine userbc (ix,iy,iz,iside,ieg)
+
+      include 'SIZE'
+      include 'NEKUSE'
+
+      integer ix,iy,iz,iside,ieg
+
+      ux=0.0
+      uy=0.0
+      uz=0.0
+
+      temp=0.0
+
+      return
+      end
+c-----------------------------------------------------------------------
+      subroutine useric (ix,iy,iz,ieg)
+
+      include 'SIZE'
+      include 'NEKUSE'
+      include 'SOLN'
+
+      integer ix,iy,iz,ieg
+
+      real pi
+      real one
+
+      one = 1.0
+
+      pi = 4.*atan(one)
+
+      if (y.le.0.5) then
+          ux = tanh(30.*(y-0.25))
+      else
+          ux = tanh(30.*(0.75-y))
+      endif
+
+      uy = 0.05*sin(2.*pi*x)
+      uz = 0.0
+
+      temp = ux
+
+      return
+      end
+c-----------------------------------------------------------------------
+      subroutine usrdat   ! This routine to modify element vertices
+
+      include 'SIZE'      
+
+      return
+      end
+c-----------------------------------------------------------------------
+      subroutine usrdat2   ! This routine to modify mesh coordinates
+
+      include 'SIZE'
+      include 'TOTAL'
+      logical ifaln,ifnorx,ifnory,ifnorz
+      do iel=1,nelt
+      do ifc=5,6
+         cbc(ifc,iel,1) = 'P  ' ! required for kludge (see below)
+      enddo
+      enddo
+
+
+      do iel=1,nelt
+      do ifc=1,2*ndim
+         call chknord(ifaln,ifnorx,ifnory,ifnorz,ifc,iel)
+         if (cbc(ifc,iel,1) .eq. 'v  ' .and. ifnorx) then 
+           boundaryID(ifc,iel) = 1
+         endif
+         if (cbc(ifc,iel,1) .eq. 'v  ') cbc(ifc,iel,2) = 't  ' 
+         if (cbc(ifc,iel,1) .eq. 'P  ') cbc(ifc,iel,2) = 'P  ' 
+      enddo
+      enddo
+
+      return
+      end
+
+c-----------------------------------------------------------------------
+      subroutine usrdat3
+
+      return
+      end
+c-----------------------------------------------------------------------
+      subroutine usrsetvert(glo_num,nel,nx,ny,nz) ! to modify glo_num
+            integer*8 glo_num(1)
+
+      include 'SIZE'
+
+      ! 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/double_shear/input.box

@@ -0,0 +1,13 @@
+one3d.rea
+-3                     spatial dimension
+1                     number of fields
+#========================================================
+#
+#========================================================
+#
+Box
+-18 -18 -1            nelx,nely,nelz for Box)
+0 1 1               x0 x1 ratio  !! Rescaled in usrdat2().
+0 1 1               y0 y1 ratio 
+0 1 1               z0 z1 ratio 
+P  ,P  ,P  ,P  ,P  ,P   BCs (3 characters each!)