Commit 8fae897a authored by ylan's avatar ylan
Browse files

add 127k case

parent d63a569b
......@@ -9,9 +9,9 @@ Pebble bed, outer geometries:
- ann = annulus
- box = box
|cases| #peb | E | E/peb | p\_ratio1 | p\_ratio2 |
|:---|---:|---:|---:|---:|---:|
| cyl146 | 146 | 62,132 | 425.56 | 91.95 |
|cases| #peb | E | E/peb | p\_ratio1 | p\_ratio2 | note |
|:---|---:|---:|---:|---:|---:|---:|
| cyl146 | 146 | 62,132 | 425.56 | 72.83 | | Rs=1.5 |
| box1053 | 1503 | 376,828 | 250.72 | 34.90 |
| cyl1568 | 1568 | 524,386 | 334.43 | 55.32 |
| cyl3260 | 3260 | 1,121,214 | 343.93 | 54.85 |
......@@ -19,6 +19,7 @@ Pebble bed, outer geometries:
| cyl11k | 11145 | 3,575,076 | 320.78 | 53.41 |
| cyl44k | 44257 | 13,032,440 | 294.47 | 51.32 |
| cyl49k | 49967 | 14,864,766 | 297.49 | 51.28 |
| cyl49k | 127602 | 39,240,190 | 307.52 | 51.39 |
**p_ratio1** is an upper estimate of the prosperity ratio computed by `(1 - Nsph*4/3*pi*r^3 / vol_geom)` where `vol_geom` is comptued by the `grep xyz` with zlen = (zmax-zmin-8) for 8=extruded length. The more pebbles, the more accurate this estimation is.
**p_ratio2** (WIP) is computed directly inside user file.
......
## 127602 (non-touching) pebbles in an annulus cylinder (Last updated: Jan. 27, 2021)
ann127602 N, Version 3.2, E39,240,190
- New features (after Aug, 2020):
- fixed edge insert
- mesh smoother+optimizer in Nek5000
- tol adjustment in edge collapse
- Notes
- DEM data from Yiqi (StarCCM++?)
- smaller edge collapse tol depends on geom + few local locations
- Geometries, pebble locations are stored at `ann127k.dat`
- Geometry and BC in re2 file:
- linear mesh, no curved sides
- Sphere distribution: free-free (selected from some range in Z, free = non-flat)
- Sphere radius = 1.0 (Wall)
- Cylinder (inner) radius = 24.256 (W01, convert to Wall in production run)
- Cylinder (outer) radius = 87.895 (W02, convert to Wall in production run)
- Zbot of the cylinder = -27.010 (W03, convert to Inflow in production run)
- Ztop of the cylinder = -29.031 (W04, convert to Outflow in production run)
```
xyz min -87.895 -87.895 -27.010
xyz max 87.895 87.895 29.031
mesh metrics: (linear mesh)
GLL grid spacing min/max : 6.77E-04 3.52E-01
scaled Jacobian min/max/avg: 1.70E-02 9.99E-01 3.89E-01
aspect ratio min/max/avg: 1.05E+00 9.23E+01 1.58E+01
```
- Mesh file (`mshmesh_retry2.fld+X`)
- lx1=3 GLL points projected onto curved sides
- no issues for interpolation onto lx1=8 and mid-multigrid levels at lx1=4
- no issues for lx1=2 corase grid
- Production run
- Use re2 for boundary conditions, but use restart file for the curved mesh
- Time is NOT set to 0 when reading the mesh file, physical time = simulation time - 20.0
- We use dimless scale to run Re10000 with vz=1.0
- NekRS version (old)
```
commit 636fe6281c1db7e45ef0a223e16745bd0fb95a16
Author: Malachi <malachi2@illinois.edu>
Date: Mon Nov 2 09:08:19 2020 -0600
Time all reduce and copy to device in ellipticResidualProjection (#168)
```
- Runs (check `dat_logfiles`):
```
- run1
step= 1 t= 2.00000100e+01 dt=1.0e-05 C= 0.07
step= 1170 t= 2.00117000e+01 dt=1.0e-05 C= 0.07 U: 2 V: 2 W: 2 P: 15 eTime= 7.32e-01, 8.82609e+02 s
- run2
step= 1 t= 2.00100500e+01 dt=5.0e-05 C= 0.34
step= 6590 t= 2.03395000e+01 dt=5.0e-05 C= 0.35 U: 2 V: 2 W: 2 P: 12 eTime= 7.06e-01, 4.36998e+03 s
- run3
step= 1 t= 2.03101000e+01 dt=1.0e-04 C= 0.72
step= 3960 t= 2.07060000e+01 dt=1.0e-04 C= 0.84 U: 2 V: 2 W: 2 P: 45 eTime= 2.04e+00, 5.60561e+03 s
- run4
step= 1 t= 2.06103000e+01 dt=3.0e-04 C= 2.59
step= 4090 t= 2.18370000e+01 dt=3.0e-04 C= 3.66 U: 3 V: 3 W: 3 P: 120 eTime= 4.87e+00, 1.33436e+04 s
- run5 start using Pr-proj=F
step= 1 t= 2.18102000e+01 dt=2.0e-04 C= 1.91
step= 5740 t= 2.29580000e+01 dt=2.0e-04 C= 1.92 U: 2 V: 2 W: 3 P: 56 eTime= 2.29e+00, 1.31354e+04 s
- run6
step= 1 t= 2.28102000e+01 dt=2.0e-04 C= 1.89
step= 5580 t= 2.39260000e+01 dt=2.0e-04 C= 2.84 U: 3 V: 2 W: 3 P: 54 eTime= 2.45e+00, 1.31585e+04 s
- run7
step= 1 t= 2.38101000e+01 dt=1.0e-04 C= 1.21
step= 6430 t= 2.44530000e+01 dt=1.0e-04 C= 0.95 U: 2 V: 2 W: 2 P: 43 eTime= 2.07e+00, 1.30440e+04 s
- run8
step= 1 t= 2.44102000e+01 dt=2.0e-04 C= 2.10
step= 5090 t= 2.54280000e+01 dt=2.0e-04 C= 1.89 U: 2 V: 2 W: 3 P: 55 eTime= 2.76e+00, 1.31923e+04 s
```
// 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;
}
[OCCA]
backend = CUDA
deviceNumber = 0
[GENERAL]
#verbose = true
polynomialOrder = 7
#startFrom = mshmesh_retry2.fld+X
startFrom = r8.fld
stopAt = numSteps
numSteps = 200000
dt = 2.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
residualProj = false
[VELOCITY]
#solver = pcg+block
boundaryTypeMap = inlet, outlet, wall, wall
density = 1.0
viscosity = -10000.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
//
// 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->startTime);
}
void UDF_ExecuteStep(ins_t *ins, dfloat time, int tstep)
{
if (ins->isOutputStep) {
nek_ocopyFrom(time, tstep);
nek_userchk();
}
}
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'
integer icalld
save icalld
data icalld /0/
integer elist(lelt)
common /imy_fld/ elist
real ca,cb,cc,cd,tmp(lx1,ly1,lz1),tmp2(lx1,ly1,lz1),plmax,plmin
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.
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 inner
boundaryID(f,e) = 3
cbc(f,e,1)='W '
cbc(f,e,2)='I '
elseif(cbc(f,e,1).eq.'W02')then ! side of cylinder outer
boundaryID(f,e) = 3
cbc(f,e,1)='W '
cbc(f,e,2)='I '
elseif(cbc(f,e,1).eq.'W03')then ! z_bottom
boundaryID(f,e) = 1
cbc(f,e,1)='v '
cbc(f,e,2)='t '
elseif(cbc(f,e,1).eq.'W04')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-----------------------------------------------------------------------
__ ____ _____
____ ___ / /__ / __ \/ ___/
/ __ \ / _ \ / //_// /_/ /\__ \
/ / / // __// ,< / _, _/___/ /
/_/ /_/ \___//_/|_|/_/ |_|/____/ v20.1
COPYRIGHT (c) 2019-2020 UCHICAGO ARGONNE, LLC
MPI tasks: 3600
using OCCA_DIR: /ccs/home/ylan/.local/nekrs_repo1105/
using OCCA_CACHE_DIR: /mnt/bb/ylan/
Initializing device
active occa mode: CUDA
loading udf ... done (7.76773s)
loading nek ... done (6.99713s)
Reading /gpfs/alpine/csc262/proj-shared/ylan/pebble_scale/peb127k_nekrs/ann127k.re2
WARNING: lgmres might be too low!
nelgt/nelgv/lelt: 39240190 39240190 10903
lx1/lx2/lx3/lxd: 8 8 8 1
partioning elements to MPI ranks
reading mesh
reading bc for ifld 1
done :: read .re2 file 3.3 sec
reading /gpfs/alpine/csc262/proj-shared/ylan/pebble_scale/peb127k_nekrs/ann127k.co2
running RCB ... finished in 0.203548 s
running RSB WARNING: Lanczos failed to converge during partitioning! finished in 101.462312s
done :: partioning 107.34 sec
reading mesh
reading bc for ifld 1
done :: read .re2 file 1.7 sec
setup mesh topology
Right-handed check complete for 39240190 elements. OK.
setvert3d: 8 5298992176 13774873216 5298992176 5298992176
gs_setup: 332319256 unique labels shared
pairwise times (avg, min, max): 0.000634839 0.000492354 0.00091067
crystal router : 0.00985698 0.00967818 0.0102903
used all_to_all method: pairwise
handle bytes (avg, min, max): 3.82833e+07 37178172 41774396
buffer bytes (avg, min, max): 3.02819e+06 1864128 6963120
setupds time 5.9509E+00 seconds 0 8 5298992176 39240190
nElements max/min/bal: 10901 10900 1.00
nMessages max/min/avg: 47 5 15.63
msgSize max/min/avg: 75415 1 12623.64
msgSizeSum max/min/avg: 435195 116508 189261.96
max multiplicity 58
done :: setup mesh topology
call usrdat
done :: usrdat
generate geometry data
NOTE: All elements deformed , param(59) ^=0
done :: generate geometry data
call usrdat2
done :: usrdat2
4.7502E-14 4.7474E-14 2.0580E-14 4.7502E-14 4.7474E-14 2.0580E-14 xyz repair 1
7.2452E-14 8.3957E-14 1.7860E-14 1.3863E-13 1.4979E-13 3.3787E-14 xyz repair 2
3.8452E-14 4.2654E-14 1.7385E-14 5.6655E-14 5.8446E-14 1.8954E-14 xyz repair 3
3.7945E-14 3.8022E-14 9.5320E-15 3.7945E-14 3.8022E-14 9.5320E-15 xyz repair 4
regenerate geometry data 1
NOTE: All elements deformed , param(59) ^=0
done :: regenerate geometry data 1
regenerate geometry data 1
NOTE: All elements deformed , param(59) ^=0
done :: regenerate geometry data 1
verify mesh topology
-87.895485019939414 87.895485250339206 Xrange
-87.895485299360658 87.895485292005461 Yrange
-27.010040198416132 29.031337995374976 Zrange
done :: verify mesh topology
mesh metrics:
GLL grid spacing min/max : 6.77E-04 3.52E-01
scaled Jacobian min/max/avg: 1.70E-02 9.99E-01 3.89E-01
aspect ratio min/max/avg: 1.05E+00 9.23E+01 1.58E+01
call usrdat3
done :: usrdat3
gridpoints unique/tot: 13774873216 20090977280
dofs vel/pr: 13138052054 13770971052
set initial conditions
Checking restart options: mshmesh_retry2.fld x
nekuic (1) for ifld 1
Reading checkpoint data
FILE:/gpfs/alpine/csc262/proj-shared/ylan/pebble_scale/peb127k_nekrs/mshmesh_retry2.fld
Mapping restart data from Nold= 3 to Nnew= 8.
0 2.0000E+01 done :: Read checkpoint data
avg data-throughput = 3700.7MBps
io-nodes = 3600
call nekuic for vel
xyz min -87.895 -87.895 -27.010
uvwpt min 0.0000 0.0000 1.0000 0.0000 0.0000
xyz max 87.895 87.895 29.031
uvwpt max 0.0000 0.0000 1.0000 0.0000 0.0000
Restart: recompute geom. factors.
regenerate geometry data 1
NOTE: All elements deformed , param(59) ^=0
done :: regenerate geometry data 1
done :: set initial conditions
calling nek_userchk
xyz min -87.895 -87.895 -27.010
xyz max 87.895 87.895 29.031
0 2.0000E+01 Write checkpoint
my_mfo_outfld 6.0000000000000000 T
my_mfo_outfld 0 0 10901 4144 39240190 13508633
FILE:/gpfs/alpine/csc262/proj-shared/ylan/pebble_scale/peb127k_nekrs/selann127k0.f00001
min/max: -87.895 87.895 -87.895 87.895 -27.010 29.031
min/max: 0.0000 0.0000 0.0000 0.0000 1.0000 1.0000
min/max: 0.0000 0.0000
0 2.0000E+01 done :: Write checkpoint
file size = 185.E+03MB
avg data-throughput = 1012.8MB/s
io-nodes = 3600
NboundaryIDs: 4
NboundaryFaces: 13070566
Nq: 8 cubNq: 11
J in range [1.19791e-06,0.215898] and max Skew = 749.237
setvert3d: 8 5298992176 13774873216 5298992176 5298992176
min NinternalElements: 6599 (ratio: 0.61)
loading udf kernels ... done
loading gs kernels ... done (44.1238s)
timing oogs modes: 0.0123243s 0.00329238s 0.00356857s used mode: 2
==================VELOCITY SETUP=========================
bID 1 -> bcType fixedValue
bID 2 -> bcType zeroGradient
bID 3 -> bcType zeroValue