[peb146] add a better mesh for peb146 case, run with NekRS

- E62k

146_pebble_yuhsiang_E62k/README.md

+## 146 (non-touching) pebbles in a cylinder (Last updated: Feb. 16, 2020)
+Peb146 N, Version 2.0 test 1, E62132
+
+This is a better mesh using improve technique developped within last few months.  
+Please note that the improvement has not finished.
+
+- New features:
+  - the resolution is refined, close to the octant refinement, but the elements is only 4x
+  - remove/merge small facet/edge redistribute the points 
+  - no negative Jacobain in NekRS multigrids 
+  - better convergence in pressure solver, so NekRS can start from constant initial flow.
+
+
+- Geometry and BC:
+  - Sphere radius = 1.0 (Wall)
+  - Cylinder radius = 7.0701 (SYM), (convert SYM to Wall in usr file)
+  - Zbot of the cylinder = 10.7524 (Outflow)
+  - Ztop of the cylinder = -8.9079 (Inflow)
+  - We scale coordinates xml=scale\*xm1 for scale=1.5/R_peb for production run
+
+
+- Relative location
+```
+ sphere                                             facets                                        oppisite sphere
+  (BC: Wall)| skin boundary layer | layer 1 | layer 2 ||  layer 2 | layer 1 | skin boundary layer |(BC: Wall)
+
+ sphere                                             facets                   Side of the cylinder
+  (BC: Wall)| skin boundary layer | layer 1 | layer 2 || skin boundary layers |(BC: SYM)
+
+ sphere                                             facets                          Bottom
+  (BC: Wall)| skin boundary layer | layer 1 | layer 2 || 3 layers of inflow region  |(BC: Inflow)
+
+ sphere                                             facets                          Top
+  (BC: Wall)| skin boundary layer | layer 1 | layer 2 || 7 layers of outfloe region |(BC: Outflow)
+```
+
+- Case files:
+    - xn3b.re2, a curved (projected at lx1=3) hex20 re2 mesh, converted by a hex20 rea: xn3.rea+reatore2, R_peb=1.5
+      ```
+        xyz min    -10.607      -10.606      -13.362    
+        xyz max     10.606       10.606       16.129
+
+        mesh metrics:
+        GLL grid spacing min/max    : 1.51E-03 4.80E-01 4.13E-01
+        scaled Jacobian  min/max/avg: 4.31E-02 9.77E-01 4.19E-01
+        aspect ratio     min/max/avg: 1.07E+00 5.69E+01 7.14E+00
+      ```

146_pebble_yuhsiang_E62k/note

+export NEKRS_HOME=$HOME/.local/nekrs_repo0204
+PATH=${NEKRS_HOME}/bin:${PATH}
+export OCCA_CACHE_DIR=/gpfs/alpine/proj-shared/csc262/ylan/cache_repo0204
+
+which nekrs
+echo $OCCA_CACHE_DIR

146_pebble_yuhsiang_E62k/submit.sh

+#!/bin/bash
+
+bin=nekrs
+project_id=nfi114
+project_id=csc262
+time="00:30"
+
+#source $NEKRS_INSTALL_DIR/bin/nekrs.bashrc
+export OCCA_VERBOSE
+
+case=$1
+
+#if [ $# -ne 3 ]; then
+#  echo "usage: $0 <setup file> <number of nodes> <GPUs per node> <hh:mm>"
+#  echo "usage: $0 <setup file> <number of nodes> <GPUs per node>"
+#  exit 0
+#fi
+if [ $# -eq 4 ]; then
+  time=$4
+fi
+
+if [ ! -f $case.par ]; then
+  echo "Cannot find" $case.par
+  exit 1
+fi
+
+mkdir -p $OCCA_CACHE_DIR 2>/dev/null
+
+nodes=$2
+gpu_per_node=$3
+let nn=$nodes*$gpu_per_node
+
+cmd="bsub -nnodes $nodes -W $time -P $project_id -J nekRS jsrun -n$nn -r$gpu_per_node -a1 -c1 -g1 -b packed:1 -d packed $bin --setup $case"
+echo $cmd
+echo $cmd > qq
+$cmd

146_pebble_yuhsiang_E62k/xn3b.oudf

+// Boundary conditions
+/* wall 1, inflow 2, outflow 3, x-slip 4, y-slip 5, z-slip 6 */
+void insVelocityDirichletConditions3D(bcData *bc)
+{                                                                        
+  bc->uP = 0.0;
+  bc->vP = 0.0;
+  bc->wP = 17.0;
+}
+
+void insVelocityNeumannConditions3D(bcData *bc)
+{
+}
+
+// Stabilized outflow (Dong et al)
+void insPressureDirichletConditions3D(bcData *bc)
+{
+  const dfloat iU0delta = 10.0;
+  const dfloat un = bc->uM*bc->nx + bc->vM*bc->ny + bc->wM*bc->nz;
+  const dfloat s0 = 0.5 * (1.0 - tanh(un*iU0delta)); 
+  bc->pP = -0.5 * (bc->uM*bc->uM + bc->vM*bc->vM + bc->wM*bc->wM) * s0;
+}

146_pebble_yuhsiang_E62k/xn3b.par

+[OCCA]
+backend = CUDA
+deviceNumber = 0 #LOCAL-RANK
+
+[GENERAL]
+polynomialOrder = 7
+#startFrom = w2.fld
+stopAt = endTime
+endTime = 200
+
+dealiasing = yes
+
+dt = 2.0e-4
+#timeStepper = bdf2
+timeStepper = tombo2
+extrapolation = subCycling
+subCyclingSteps = 2
+#targetCFL = 3.8
+
+writeControl = TIMESTEP
+writeInterval = 1000
+
+filtering = hpfrt
+filterWeight = 250
+filterCutoffRatio = 0.9 
+
+[PRESSURE]
+preconditioner = semg_amg
+#amgSolver = parAlmond
+#galerkincoarsematrix = yes
+#galerkincoarsegrid = yes
+galerkincoarseoperator = yes
+residualTol = 1e-04
+residualProj = yes
+
+
+[VELOCITY]
+boundaryTypeMap = inlet, outlet, wall
+density = 1.0
+viscosity = 0.25 #0.000405
+residualTol = 1e-06
+
+#[TEMPERATURE]
+#solver=none
+#boundaryTypeMap = t, O, I, t
+#rhoCp = 10
+#conductivity = 5000
+#residualTol = 1e-08

146_pebble_yuhsiang_E62k/xn3b.udf

+//
+// nekRS User Defined File
+//
+
+#include <math.h>
+#include "udf.hpp"
+
+/* UDF Functions */                                                      
+
+void UDF_LoadKernels(ins_t *ins)
+{
+}
+
+void UDF_Setup(ins_t *ins)
+{
+  // get IC from nek
+  if (!ins->readRestartFile) nek_copyTo(ins, ins->startTime);
+}
+
+void UDF_ExecuteStep(ins_t *ins, dfloat time, int tstep)
+{
+  if (ins->isOutputStep) nek_userchk();
+}

146_pebble_yuhsiang_E62k/xn3b.usr

+c-----------------------------------------------------------------------
+      subroutine uservp(i,j,k,eg) ! set variable properties
+      include 'SIZE'
+      include 'TOTAL'
+      include 'NEKUSE'
+
+      integer i,j,k,e,eg
+     
+c     e = gllel(eg)
+
+      udiff  = 0.0
+      utrans = 0.0
+
+      return
+      end
+c-----------------------------------------------------------------------
+      subroutine userf(i,j,k,eg) ! set acceleration term
+c
+c     Note: this is an acceleration term, NOT a force!
+c     Thus, ffx will subsequently be multiplied by rho(x,t).
+c
+      include 'SIZE'
+      include 'TOTAL'
+      include 'NEKUSE'
+
+      integer i,j,k,e,eg
+
+c     e = gllel(eg)
+
+      ffx = 0.0
+      ffy = 0.0
+      ffz = 0.0
+
+      return
+      end
+c-----------------------------------------------------------------------
+      subroutine userq(i,j,k,eg) ! set source term
+      include 'SIZE'
+      include 'TOTAL'
+      include 'NEKUSE'
+
+      integer i,j,k,e,eg
+
+c     e = gllel(eg)
+
+      qvol   = 0.0
+
+      return
+      end
+c-----------------------------------------------------------------------
+      subroutine userbc(i,j,k,f,eg) ! set up boundary conditions
+c     NOTE ::: This routine MAY NOT be called by every process
+      include 'SIZE'
+      include 'TOTAL'
+      include 'NEKUSE'
+
+      integer i,j,k,f,e,eg
+      e=gllel(eg)
+
+      ux   = 0.0
+      uy   = 0.0
+      uz   = 17.0
+
+      flux = 0.0
+      temp = 0.0
+      if (cbc(f,e,1).eq.'W  '.AND.cbc(f,e,2).eq.'t  ') temp=1.0
+
+      return
+      end
+c-----------------------------------------------------------------------
+      subroutine useric(i,j,k,eg) ! set up initial conditions
+      include 'SIZE'
+      include 'TOTAL'
+      include 'NEKUSE'
+      integer i,j,k,e,eg
+
+      ux   = 0.0
+      uy   = 0.0
+      uz   = 17.0
+      temp = 0.0
+
+      return
+      end
+c-----------------------------------------------------------------------
+      subroutine userchk
+      include 'SIZE'
+      include 'TOTAL'
+
+c      common /myoutflow/ d(lx1,ly1,lz1,lelt),m1(lx1*ly1*lz1,lelt)
+c      real m1
+c      rq  = 2.
+c      uin = 0.
+c      call turb_outflow(d,m1,rq,uin)
+
+      n=lx1*ly1*lz1*nelt
+
+      if (istep.eq.0) then
+c        time=0.0
+c
+c        Rfld = 0.9   ! radii of peb from restart file
+c        Rtar = 1.5   ! target radii matching Re/density
+c        scale = 1.5/0.9
+c        do i=1,n
+c          xm1(i,1,1,1)=xm1(i,1,1,1)*scale
+c          ym1(i,1,1,1)=ym1(i,1,1,1)*scale
+c          zm1(i,1,1,1)=zm1(i,1,1,1)*scale
+c        enddo
+c        call fix_geom
+c        call geom_reset(1)
+
+        xxmax = glmax(xm1,n)
+        yymax = glmax(ym1,n)
+        zzmax = glmax(zm1,n)
+        xxmin = glmin(xm1,n)
+        yymin = glmin(ym1,n)
+        zzmin = glmin(zm1,n)
+        if (nio.eq.0) then
+          write(6,18) xxmin,yymin,zzmin
+          write(6,19) xxmax,yymax,zzmax
+   18     format(' xyz min  ',5g13.5)
+   19     format(' xyz max  ',5g13.5)
+        endif
+
+        ifxyo=.true.
+c        call outpost(vx,vy,vz,t,pr,'   ')
+      endif
+
+      if (mod(istep,500).eq.0) then
+         umin=glmin(vx,n)
+         umax=glmax(vx,n)
+         vmin=glmin(vy,n)
+         vmax=glmax(vy,n)
+         wmin=glmin(vz,n)
+         wmax=glmax(vz,n)
+         tmin=glmin(t ,n)
+         tmax=glmax(t ,n)
+         if (nio.eq.0) write(6,1) 
+     $    istep,time,umin,umax,vmin,vmax,wmin,wmax,tmin,tmax
+    1    format(i9,1p9e11.3,' tmax')
+      endif
+
+      return
+      end
+c-----------------------------------------------------------------------
+      subroutine userqtl ! Set thermal divergence
+
+      call userqtl_scig
+
+      return
+      end
+c-----------------------------------------------------------------------
+      subroutine usrdat   ! This routine to modify element vertices
+      include 'SIZE'
+      include 'TOTAL'
+
+      return
+      end
+c-----------------------------------------------------------------------
+      subroutine usrdat3
+      include 'SIZE'
+      include 'TOTAL'
+
+      return
+      end
+c-----------------------------------------------------------------------
+      subroutine usrdat2
+      include 'SIZE'
+      include 'TOTAL'
+      integer e,f
+
+c      ifheat = .true.
+      nface = 2*ldim
+
+      do e=1,nelv
+      do f=1,nface
+        if(cbc(f,e,1).eq.'E  '.or.cbc(f,e,1).eq.'   ')then
+          cbc(f,e,1)='E  '
+          cbc(f,e,2)='E  '
+        elseif(cbc(f,e,1).eq.'v  ')then  ! z_bottom
+          cbc(f,e,2)='t  '
+        elseif(cbc(f,e,1).eq.'O  ')then  ! z_top
+          cbc(f,e,2)='O  '
+        elseif(cbc(f,e,1).eq.'SYM')then  ! side of cylinder
+          cbc(f,e,1)='W  '
+          cbc(f,e,2)='I  '
+        elseif(cbc(f,e,1).eq.'W  ')then  ! pebbles
+          cbc(f,e,2)='t  '
+        else
+          write(*,*)'bc incorrect',e,f,cbc(f,e,1)
+        endif
+      enddo
+      enddo
+
+      do iel=1,nelt
+      do ifc=1,2*ndim
+         if (cbc(ifc,iel,1) .eq. 'v  ') boundaryID(ifc,iel) = 1
+         if (cbc(ifc,iel,1) .eq. 'O  ') boundaryID(ifc,iel) = 2
+         if (cbc(ifc,iel,1) .eq. 'W  ') boundaryID(ifc,iel) = 3
+      enddo
+      enddo
+
+
+      return
+      end
+c-----------------------------------------------------------------------