[peb-mesh] add peb11k mesh

- This is not stable for Re10k, but still good enough for scaling tests.

11145_pebble_yuhsiang_E3p5M/README.md

+## 11145 (non-touching) pebbles in a cylinder (Last updated: Jun. 12, 2020)
+cyl11145 N, Version 3.0, E3,575,076
+
+This mesh is not very stable for Re10k. There is (at least) a known location having high distorted Cell.  
+Run dt=5e-4 is fine for the first 14000 steps.   
+After that dt=4e-4 will have CFL too high stop at 14000+5996th step (time=9.99800000e+00)   
+We never manage to run over an convection unit  
+
+This is a better mesh using improved techniques developed within last few months.  
+Please note that the improvement has not finished.  
+
+- New features:
+  - A more robust and cleaner code with lots of user intreface improvement
+  - Add Taubin smoothing before edge collapse
+  - Switch to a passive tessellation if the first attempt didn't do well. 
+  - Add local optimization for better performance.
+  - dump co2 directly  
+  - DEM results from project chrono on nurburg gpu
+  - Fix some bugs as well  
+
+- Geometry and BC:
+  - Sphere distribution: flat-free
+  - Sphere radius = 1.0 (Wall)
+  - Cylinder radius = 40.672 (SYM), (convert SYM to Wall in usr file)
+  - Zbot of the cylinder =  14.655 (Outflow)
+  - Ztop of the cylinder = -12.627 (Inflow)
+  - We use dimless scale to run Re10000 with vz=1.0
+
+- Case files:
+    - .re2, a curved (projected at lx1=3) hex20 re2 mesh, converted by a hex20 rea: xn3.rea+reatore2, R\_peb=1.0
+      ```
+        xyz min    -40.672      -40.672      -12.627    
+        xyz max     40.672       40.672       14.655
+      
+        mesh metrics:
+        GLL grid spacing min/max    : 3.67E-04 3.32E-01
+        scaled Jacobian  min/max/avg: 1.52E-02 9.97E-01 3.65E-01
+        aspect ratio     min/max/avg: 1.06E+00 6.06E+01 1.27E+01
+      ```
+    - Manage to run to time=7.0 with dt=5e-4 and the follow commit of NekRS
+      ```
+        commit c26b570af4a85e8ef8425b7e7057f0926f5aaace
+        Author: Stefan Kerkemeier <stefanke@jlselogin2.ftm.alcf.anl.gov>
+        Date:   Fri Aug 7 12:06:07 2020 +0000
+        
+            add missing file and fix gsh for scalars
+      ```
+    - peb.11k.dat, locations of the pebble centers under R\_peb=1.0  

11145_pebble_yuhsiang_E3p5M/peb11k.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 = 1.0;
+}
+
+void insVelocityNeumannConditions3D(bcData *bc)
+{
+}
+
+void cdsDirichletConditions3D(bcData *bc)
+{
+  bc->sP = 0.0;
+}
+void cdsNeumannConditions3D(bcData *bc)
+{
+  bc->sF = 1.0;
+}
+
+// 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;
+}

11145_pebble_yuhsiang_E3p5M/peb11k.par

+[OCCA]
+backend = CUDA
+deviceNumber = 0
+
+[GENERAL]
+#verbose = true
+polynomialOrder = 7
+startFrom = r7.fld
+stopAt = numSteps
+numSteps = 20000
+
+dt = 4.0e-4
+timeStepper = tombo2
+extrapolation = subCycling
+subCyclingSteps = 2
+
+writeControl = TIMESTEP
+writeInterval = 1000 
+
+filtering = hpfrt
+filterWeight = 200
+filterModes = 2
+
+[PRESSURE]
+#preconditioner = semg_amg
+#amgSolver = parAlmond
+#galerkinCoarseOperator = yes
+residualTol = 1e-04
+
+[VELOCITY]
+#solver = pcg+block
+boundaryTypeMap = inlet, outlet, wall, wall
+density = 1.0
+viscosity = -5000.0
+residualTol = 1e-06
+
+#[TEMPERATURE]
+#solver = none
+#boundaryTypeMap = inlet, outlet, insulated, flux
+#rhoCp = 1.0
+#conductivity = -100
+#residualTol = 1e-06
+
+#[SCALAR01]
+#solver = none
+
+[BOOMERAMG]
+#coarsenType = 10
+#interpolationType = 6
+#smootherType = -1
+#iterations = 2 # testing
+#strongThreshold = 0.25
+#nonGalerkinTol = 0.1

11145_pebble_yuhsiang_E3p5M/peb11k.udf

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

11145_pebble_yuhsiang_E3p5M/peb11k.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   = 1.0
+
+      if (cbu.eq.'o  ') then
+         U0 = 1.0                  ! characteristic velocity
+         delta = 0.1               ! small positive constant
+         pa = dongOutflow(i,j,k,e,f,U0,delta)
+      endif
+
+      flux = 0.0
+      temp = 0.0
+c     if (cbc(f,e,1).eq.'W  '.AND.cbc(f,e,2).eq.'t  ') temp=1.0
+      if (cbc(f,e,1).eq.'W  '.AND.cbc(f,e,2).eq.'f  ') flux=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   = 1.0
+      temp = 0.0
+
+      return
+      end
+c-----------------------------------------------------------------------
+      subroutine userchk
+      include 'SIZE'
+      include 'TOTAL'
+
+      n=lx1*ly1*lz1*nelt
+
+      if (istep.eq.0) then
+        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
+
+c      if (mod(istep,500).eq.0) then
+c         umin=glmin(vx,n)
+c         umax=glmax(vx,n)
+c         vmin=glmin(vy,n)
+c         vmax=glmax(vy,n)
+c         wmin=glmin(vz,n)
+c         wmax=glmax(vz,n)
+c         tmin=glmin(t ,n)
+c         tmax=glmax(t ,n)
+c         if (nio.eq.0) write(6,1) 
+c     $    istep,time,umin,umax,vmin,vmax,wmin,wmax,tmin,tmax
+c    1    format(i9,1p9e11.3,' tmax')
+c
+c         call compute_cfl(cfl,vx,vy,vz,dt)
+c         call outpost(cflf,vy,vz,pr,t,'cfl')
+c      endif
+
+c      if (istep.gt.0) then
+c        ifield = 1 ! essential for lambda2 right now
+c        call lambda2(t(1,1,1,1,2))
+c          tmin=glmin(t(1,1,1,1,2) ,n)
+c          tmax=glmax(t(1,1,1,1,2) ,n)
+c          if (nio.eq.0) write(6,2)istep,time,tmin,tmax
+c    2     format(i9,1p3e11.3,' tmax2') 
+c
+c        ifxyo     = .true.
+c        ifvo      = .true.
+c        ifpro     = .true.
+c        ifto      = .true.
+c        ifpsco(1) = .true.
+cc       call outpost(vx,vy,vz,pr,t,'   ')
+c      endif
+c      call comp_Reynolds
+      if (mod(istep,500).ne.0) then
+        call compute_cfl(cfl,vx,vy,vz,1.0)
+        call outpost(cflf,vy,vz,pr,t,'cfl')
+      endif
+
+      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
+
+      ifheat = .false.
+      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.'W  ')then  ! pebbles
+          boundaryID(f,e) = 4
+          cbc(f,e,1)='W  '
+          cbc(f,e,2)='f  '
+        elseif(cbc(f,e,1).eq.'W01')then  ! side of cylinder
+          boundaryID(f,e) = 3
+          cbc(f,e,1)='W  '
+          cbc(f,e,2)='I  '
+        elseif(cbc(f,e,1).eq.'W02')then  ! z_bottom
+          boundaryID(f,e) = 1
+          cbc(f,e,1)='v  '
+          cbc(f,e,2)='t  '
+        elseif(cbc(f,e,1).eq.'W03')then  ! z_top
+          boundaryID(f,e) = 2
+          cbc(f,e,1)='o  ' ! Dong's outflow
+          cbc(f,e,2)='O  '
+        else
+          write(*,*)'bc incorrect',e,f,cbc(f,e,1)
+        endif
+      enddo
+      enddo
+
+      param(18) = 1 ! use algebraic residual norm
+
+      return
+      end
+c-----------------------------------------------------------------------
+      function dongOutflow(ix,iy,iz,iel,iside,u0,delta)
+
+      include 'SIZE'
+      include 'SOLN'
+      include 'GEOM'
+
+      real sn(3)
+
+      ux = vx(ix,iy,iz,iel)
+      uy = vy(ix,iy,iz,iel)
+      uz = vz(ix,iy,iz,iel)
+
+      call getSnormal(sn,ix,iy,iz,iside,iel)
+      vn = ux*sn(1) + uy*sn(2) + uz*sn(3)
+      S0 = 0.5*(1.0 - tanh(vn/u0/delta))
+
+      dongOutflow = -0.5*(ux*ux+uy*uy+uz*uz)*S0
+
+      return
+      end
+c-----------------------------------------------------------------------
+      subroutine comp_Reynolds
+      implicit none
+      include 'SIZE'
+      include 'TOTAL'
+
+      integer i,n,npeb
+      real Rcyl,Rpeb
+      real z0,z1,wd,fun_z,wt(lx1*ly1*lz1*lelv)
+      real zz,dis,glsum,glsc2
+      real wght_vol,area_peb,area_cyl1,area_cyl2,vz_bar,Dh1,Dh2
+
+      npeb = 1568
+      Rpeb = 1.0
+      Rcyl = 13.858 
+
+      ! min/max of pebble centers are -1.05733229e+01 1.15912238e+01
+      ! f ------\
+      !          \
+      !           \-----
+      !    o   ) . 
+      !    P 1 | wd |  
+      wd =  1.5
+      z0 = -10.573 - 1.0 - wd*0.5
+      z1 =  11.159 + 1.0 + wd*0.5
+
+
+      n=lx1*ly1*lz1*nelv
+      do i=1,n
+        zz=zm1(i,1,1,1)
+
+        if (zz>z1) then
+          dis = -abs(zz-z0)
+        elseif (zz>z0) then
+          dis = min(abs(z1-zz),abs(zz-z0))
+        else
+          dis = -abs(z0-zz)
+        endif
+        fun_z = tanh(dis*wd)
+
+        wt(i) = fun_z * bm1(i,1,1,1) 
+      enddo
+
+      wght_vol = glsum(wt,n)
+
+      area_peb = npeb*4.0*pi*Rpeb*Rpeb
+      area_cyl1 = 2.0*Rcyl*pi*(z1-z0)
+      area_cyl2 = 2.0*Rcyl*pi*(z1-z0+wd)
+
+      vz_bar=glsc2(wt,vz,n)/wght_vol
+
+      Dh1 = 4.0*wght_vol/(area_peb+area_cyl1)
+      Dh2 = 4.0*wght_vol/(area_peb+area_cyl2)
+
+      if(nid.eq.0)then
+        write(*,1)z0,z1,wd,wght_vol
+        write(*,2)area_peb,area_cyl1,area_cyl2
+        write(*,3)time,vz_bar,Dh1,Dh2
+      endif
+    1 format('CompRe, domain',1p4e11.3)
+    2 format('CompRe, s_area',1p3e11.3)
+    3 format('CompRe, vz/Dh ',1p4e11.3)
+
+      return
+      end
+c-----------------------------------------------------------------------

3260_pebble_yuhsiang_E1p1M/README.md

@@ -13,9 +13,11 @@ Please note that the improvement has not finished.
   - Switch to a passive tessellation if the first attempt didn't do well. 
   - Add local optimization for better performance.
   - dump co2 directly  
+  - DEM results from project chrono on nurburg gpu
   - Fix some bugs as well  
 
 - Geometry and BC:
+  - Sphere distribution: flat-free
   - Sphere radius = 1.0 (Wall)
   - Cylinder radius = 26.776 (SYM), (convert SYM to Wall in usr file)
   - Zbot of the cylinder =  12.117 (Outflow)