Commit 02c20cbf authored by ylan's avatar ylan
Browse files

add ann350k

parent 07e614d6
......@@ -20,6 +20,7 @@ Pebble bed, outer geometries:
| cyl44k | 44257 | 13,032,440 | 294.47 | 51.32 |
| cyl49k | 49967 | 14,864,766 | 297.49 | 51.28 |
| ann127k | 127602 | 39,240,190 | 307.52 | 51.39 |
| ann350k | 352625 | 98,782,067 | 280.13 | 60.08 | | Rs=0.91945
**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.
......
ann350k.co2 filter=lfs diff=lfs merge=lfs -text
ann350k.re2 filter=lfs diff=lfs merge=lfs -text
ann350k_t1r11.fld filter=lfs diff=lfs merge=lfs -text
## 352,625 (non-touching) pebbles in an annulus cylinder (Last updated: Feb. 15, 2021)
ann352625 N, Version 3.3, E98,782,067
- New features (after Jan, 2021):
- mesh smoother+optimizer in Nek5000
- tol adjustment in edge collapse
- new scale (dsep = 2.0 = k(2.5), Rs0=0.8388, Rs1=1.0)
- new boundary configureation (extrusion instead of projection)
- Fix edge insert, edge collapse, Eadj, and Qadj
- Notes
- DEM data from Yiqi (StarCCM++?)
- smaller edge collapse tol depends on geom + few local locations
- Geometries, pebble locations are stored at `ann350k.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 = 0.91945 (Wall)
- Cylinder (inner) radius = 23.101 (W01, convert to Wall in production run)
- Cylinder (outer) radius = 83.300 (W02, convert to Wall in production run)
- Zbot of the cylinder = -64.275 (W03, convert to Inflow in production run)
- Ztop of the cylinder = 66.437 (W04, convert to Outflow in production run)
```
xyz min -83.301 -83.301 -64.276
xyz max 83.301 83.301 66.438
mesh metrics: (linear mesh)
GLL grid spacing min/max : 1.18E-03 3.22E-01
scaled Jacobian min/max/avg: 4.33E-02 1.00E+00 4.09E-01
aspect ratio min/max/avg: 1.05E+00 4.60E+01 6.79E+00
mesh metrics: (curved mesh)
GLL grid spacing min/max : 4.41E-04 3.18E-01
scaled Jacobian min/max/avg: 2.66E-02 9.99E-01 3.84E-01
aspect ratio min/max/avg: 1.05E+00 7.72E+01 1.07E+01
```
- Mesh file (`ann350k_t1r11.fld+X+time=0.0`)
- 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
- We use dimless scale to run Re5000 with vz=1.0 (Re10k seems ok, but no time/allocation left)
- NekRS version (after udf changes)
```
commit 82a68b6bbf550c84e3757ffcf06b783578a133d3
Author: Stefan <stgeke@gmail.com>
Date: Wed Jan 20 21:02:30 2021 +0100
Fix zero usr file generation
```
- Runs (check `dat_logfiles`) (1600 nodes on Summit):
```
- run1, Re50
step= 500 t= 1.00000000e-01 dt=2.0e-04 C= 0.30 UVW: 3 P: 16 eTime= 8.79e-01, 7.41745e+02 s
- run2, Re5k
step= 5000 t= 3.50000000e-01 dt=5.0e-05 C= 0.10 UVW: 1 P: 7 eTime= 4.02e-01, 1.75611e+03 s
- run3, Re10k
step= 3000 t= 9.50000000e-01 dt=2.0e-04 C= 1.19 UVW: 1 P: 25 eTime= 1.09e+00, 2.44549e+03 s
- run4, Re5k
step= 2490 t= 1.24800000e+00 dt=2.0e-04 C= 0.95 UVW: 2 P: 35 eTime= 1.49e+00, 2.37206e+03 s
- run5, Re5k
step= 5000 t= 2.15000000e+00 dt=2.0e-04 C= 0.77 UVW: 2 P: 42 eTime= 1.82e+00, 5.22005e+03 s
- run6, Re5k
step= 10370 t= 4.02400000e+00 dt=2.0e-04 C= 0.83 UVW: 2 P: 23 eTime= 1.11e+00, 1.04075e+04 s
- run7, Re5
step= 11000 t= 6.15000000e+00 dt=2.0e-04 C= 0.96 UVW: 2 P: 28 eTime= 1.26e+00, 1.04812e+04 s
```
This diff is collapsed.
// Boundary conditions
void velocityDirichletConditions(bcData *bc)
{
bc->u = 0.0;
bc->v = 0.0;
bc->w = 1.0;
}
void scalarDirichletConditions(bcData *bc)
{
bc->s = 0.0;
}
void scalarNeumannConditions(bcData *bc)
{
bc->flux = 1.0;
}
// Stabilized outflow (Dong et al)
void pressureDirichletConditions(bcData *bc)
{
const dfloat iU0delta = 10.0;
const dfloat un = bc->u*bc->nx + bc->v*bc->ny + bc->w*bc->nz;
const dfloat s0 = 0.5 * (1.0 - tanh(un*iU0delta));
bc->p = -0.5 * (bc->u*bc->u + bc->v*bc->v + bc->w*bc->w) * s0;
}
[OCCA]
backend = CUDA
deviceNumber = 0
[GENERAL]
#verbose = true
polynomialOrder = 7
#startFrom = ann350k_t1r11.fld+X+time=0.0
startFrom = r11okr7.fld
stopAt = numSteps
numSteps = 50000
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
#residualProj = false
[VELOCITY]
#solver = pcg+block
boundaryTypeMap = wall, wall, wall, inlet, outlet
density = 1.0
viscosity = -5000.0
residualTol = 1e-06
#[TEMPERATURE]
#solver = none
#boundaryTypeMap = insulated, flux, flux, inlet, outlet
#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(nrs_t *nrs)
{
}
void UDF_Setup0(MPI_Comm comm, setupAide &options)
{
// force verbose output
options.setArgs("VERBOSE", "TRUE");
}
void UDF_Setup(nrs_t *nrs)
{
}
void UDF_ExecuteStep(nrs_t *nrs, dfloat time, int tstep)
{
mesh_t *mesh = nrs->mesh;
cds_t *cds = nrs->cds;
// For the first terminatVerboseOutput steps, allow verbose output
// from all solvers.
// After terminateVerboseOutput steps, allow verbose output
// only from solution projection.
const int terminateVerboseOutput = 20;
if(tstep == (terminateVerboseOutput+1)){
nrs->pSolver->options.setArgs("VERBOSE", "FALSE");
if(nrs->uvwSolver){
nrs->uvwSolver->options.setArgs("VERBOSE", "FALSE");
} else {
nrs->uSolver->options.setArgs("VERBOSE", "FALSE");
nrs->vSolver->options.setArgs("VERBOSE", "FALSE");
nrs->wSolver->options.setArgs("VERBOSE", "FALSE");
}
for(int scalarField = 0; scalarField < cds->NSfields; ++scalarField){
cds->solver[scalarField]->options.setArgs("VERBOSE", "FALSE");
}
}
if (nrs->isOutputStep) {
nek_ocopyFrom(time, tstep);
nek_userchk();
}
}
include 'my_errchk.f' ! my checking routines
include 'my_post.f' ! selected outpost
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.
call my_mesh_metric
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
call set_myboundaryID_dummy
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) = 1
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) = 2
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) = 4
cbc(f,e,1)='v '
cbc(f,e,2)='t '
elseif(cbc(f,e,1).eq.'W04')then ! z_top
boundaryID(f,e) = 5
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-----------------------------------------------------------------------
subroutine set_myboundaryID_dummy
implicit none
include 'SIZE'
include 'TOTAL'
integer lgeo
parameter(lgeo=7)
character*3 CB
character*3 cbc_msh(lgeo)
integer e,f,ig, ncnt(lgeo),iglsum,cbc_id
cbc_msh(1)='W '
cbc_msh(2)='W01'
cbc_msh(3)='W02'
cbc_msh(4)='W03'
cbc_msh(5)='W04'
cbc_msh(6)='W05'
cbc_msh(7)='W06'
call izero(ncnt,lgeo)
do e=1,nelv
do f=1,2*ldim
CB =cbc(f,e,1)
! bc(slot,face,e,ifield)
c peb_id=int(bc(5,f,e,1))
c opp_id=int(bc(4,f,e,1))
cbc_id=int(bc(3,f,e,1))
if (cbc_id.gt.lgeo) then
write(*,*)'Err, lgeo is too small',lglel(e),f,cbc_id
call exitt
endif
if(cbc_id.gt.0) then
cbc(f,e,1)=cbc_msh(cbc_id)
ncnt(cbc_id)=ncnt(cbc_id)+1
endif
enddo
enddo
do ig=1,lgeo
ncnt(ig)=iglsum(ncnt(ig),1)
enddo
if(nid.eq.0)then
write(*,*)'bc cnt: ',(ncnt(ig),ig=1,lgeo)
endif
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-----------------------------------------------------------------------
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment