Commit b441a7bc authored by ylan's avatar ylan
Browse files

[peb1568] updateto mesh t2

- This one gives larger dt
- I've run this to time=25
parent 2a0b657c
## 1568 (non-touching) pebbles in a cylinder (Last updated: Mar. 02, 2020)
## 1568 (non-touching) pebbles in a cylinder (Last updated: Mar. 21, 2020)
Peb1568 N, Version 2.0 test 1, E524386
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.
- slightly adjust dxmin to have larger dt
- Geometry and BC:
......@@ -25,7 +22,16 @@ Please note that the improvement has not finished.
xyz max 13.858 13.858 17.691
mesh metrics:
GLL grid spacing min/max : 8.26E-05 3.02E-01
scaled Jacobian min/max/avg: 5.87E-03 9.91E-01 3.71E-01
aspect ratio min/max/avg: 1.12E+00 4.10E+02 1.26E+01
GLL grid spacing min/max : 2.21E-04 3.02E-01
scaled Jacobian min/max/avg: 2.59E-02 9.91E-01 3.71E-01
aspect ratio min/max/avg: 1.12E+00 1.08E+02 1.26E+01
```
- Manage to run to time=25 with this setup and the follow commit of NekRS
```
commit 25ff78464ddae123c6a30b28f2cd6d3d5e03afb7
Author: Stefan Kerkemeier <stefanke@jlselogin2.ftm.alcf.anl.gov>
Date: Mon Apr 20 13:33:18 2020 +0000
change writeInterval for turbPipe example
```
This diff is collapsed.
module load cmake gcc cuda
. ~/.bashrc
export NEKRS_HOME=$HOME/.local/nekrs_repo0204
export NEKRS_HOME=$HOME/.local/nekrs_v20_dbg
PATH=${NEKRS_HOME}/bin:${PATH}
export OCCA_CACHE_DIR=/gpfs/alpine/proj-shared/csc262/ylan/cache_repo0204
export OCCA_CACHE_DIR=/gpfs/alpine/proj-shared/csc262/ylan/cache_v20_dbg
which nekrs
echo $OCCA_CACHE_DIR
****** PARAMETERS *****
2.60999990 NEKTON VERSION
3 DIMENSIONAL RUN
103 PARAMETERS FOLLOW
1.00000 p001 DENSITY
-1000.000 p002 VISCOS
0.00000 p003
0.00000 p004
0.00000 p005
0.00000 p006
1.00000 p007 RHOCP
-1000.00000 p008 CONDUCT
0.00000 p009
0.00000 p010 FINTIME
2000.00 p011 NSTEPS
-4.00000E-04 p012 DT
0.00000 p013 IOCOMM
0.00000 p014 IOTIME
0.00000 p015 IOSTEP
0.00000 p016 p16 PSSOLVER
0.00000 p017
0.00000 p018
0.00000 p019
0.00000 p020
0.100000E-05 p021 DIVERGENCE
0.100000E-07 p022 HELMHOLTZ
0.00000 p023 NPSCAL
0.100000E-09 p024 TOLREL
0.100000E-09 p025 TOLABS
0.250000 p026 COURANT
2.00000 p027 TORDER
0.00000 p028 TORDER: mesh velocity (0: p28=p27)
0.00000 p029 = magnetic visc if > 0, = -1/Rm if < 0
0.00000 p030 > 0 properties set in uservp()
0.00000 p031 NPERT: #perturbation modes
0.00000 p032 #BCs in re2 file, if > 0
0.00000 p033
0.00000 p034
0.00000 p035
0.00000 p036
0.00000 p037
0.00000 p038
0.00000 p039
0.00000 p040
0.00000 p041
0.00000 p042 p42 0=gmres/1=pcg
0.00000 p043 p43 0=semg/1=schwarz
0.00000 p044 p44 0=E-based/1=A-based prec.
0.00000 p045 p45 Relaxation factor for DTFS
0.00000 p046
0.00000 p047 viscosity for mesh elasticity solver
0.00000 p048
0.00000 p049
0.00000 p050
0.00000 p051
0.00000 p052 p52 IOHIS
0.00000 p053
0.00000 p054 fixed flow rate dir: |p54|=1,2,3=x,y,z
0.00000 p055 vol.flow rate (p54>0) or Ubar (p54<0)
0.00000 p056
0.00000 p057
0.00000 p058
0.00000 p059
0.00000 p060 !=0 --> init. velocity to small nonzero
0.00000 p061
0.00000 p062 >0 --> force byte_swap for output
8 0.00000 p063 =8 --> force 8-byte output
0.00000 p064 =1 --> perturbation restart
0.00000 p065 #iofiles (eg, 0 or 64); <0 --> sep. dirs
6.00000 p066 p66 write fmt:ONLY postx uses rea v
6.00000 p067
0.00000 p068 iastep: freq for avg_all (0=iostep)
0.00000 p069
0.00000 p070
0.00000 p071
0.00000 p072
0.00000 p073
0.00000 p074 p74 verbose Helmholtz
0.00000 p075
0.00000 p076
0.00000 p077
0.00000 p078
0.00000 p079
0.00000 p080
0.00000 p081
0.00000 p082
0.00000 p083
0.00000 p084 !=0 --> sets initial timestep if p12>0
0.00000 p085 dt ratio if p84 !=0, for timesteps>0
0.00000 p086 p86 reserved
0.00000 p087
0.00000 p088
0.00000 p089 p89 reserved
0.00000 p090
0.00000 p091
0.00000 p092
20.00000 p093 Number of previous pressure solns saved
5.00000 p094 start projecting velocity after p94 step
5.00000 p095 start projecting pressure after p95 step
0.00000 p096
0.00000 p097
0.00000 p098
3.00000 p099 dealiasing: <0 off
0.00000 p100
0.00000 p101 Number of additional modes to filter
0.00000 p102
0.500000E-01 p103 weight of stabilizing filter (.01)
4 Lines of passive scalar data follows2 CONDUCT; 2RHOCP
1.00000 1.00000 1.00000 1.00000 1.00000
1.00000 1.00000 1.00000 1.00000
1.00000 1.00000 1.00000 1.00000 1.00000
1.00000 1.00000 1.00000 1.00000
13 LOGICAL SWITCHES FOLLOW
T IFFLOW
T IFHEAT
T IFTRAN
T T T T T F F F F F F IFNAV & IFADVC (convection in P.S. fields)
F F F F F F T T T T T T IFTMSH (IF mesh for this field is T mesh)
F IFAXIS
T IFSTRS
F IFSPLIT
F IFMGRID
F IFMODEL
F IFKEPS
F IFMVBD
F IFCHAR
8.00000 8.00000 -4.00000 -4.00000 XFAC,YFAC,XZERO,YZERO
**MESH DATA** 6 lines are X,Y,Z;X,Y,Z. Columns corners 1-4;5-8
-524386 3 524386 NELT,NDIM,NELV
0 PRESOLVE/RESTART OPTIONS *****
7 INITIAL CONDITIONS *****
C Default
C Default
C Default
C Default
C Default
C Default
C Default
***** DRIVE FORCE DATA ***** BODY FORCE, FLOW, Q
4 Lines of Drive force data follow
C
C
C
C
***** Variable Property Data ***** Overrrides Parameter data.
1 Lines follow.
0 PACKETS OF DATA FOLLOW
***** HISTORY AND INTEGRAL DATA *****
0 POINTS. Hcode, I,J,H,IEL
***** OUTPUT FIELD SPECIFICATION *****
6 SPECIFICATIONS FOLLOW
T COORDINATES
T VELOCITY
T PRESSURE
T TEMPERATURE
F TEMPERATURE GRADIENT
0 PASSIVE SCALARS
***** OBJECT SPECIFICATION *****
0 Surface Objects
0 Volume Objects
0 Edge Objects
0 Point Objects
......@@ -11,6 +11,15 @@ 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)
{
......
[OCCA]
backend = CUDA
deviceNumber = 0 #LOCAL-RANK
deviceNumber = 0
[GENERAL]
#verbose = true
polynomialOrder = 7
startFrom = r17.fld
stopAt = endTime
endTime = 200
startFrom = r5.fld
stopAt = numSteps
numSteps = 20000
dealiasing = yes
dt = 1.e-4 #1.0e-3
dt = 5.0e-4
timeStepper = tombo2
extrapolation = subCycling
subCyclingSteps = 2
targetCFL = 3.8
writeControl = TIMESTEP
writeInterval = 1000
writeInterval = 1000
filtering = hpfrt
filterWeight = 250
filterCutoffRatio = 0.9
filterWeight = 200
filterModes = 2
[PRESSURE]
preconditioner = semg_amg
#preconditioner = semg_amg
#amgSolver = parAlmond
galerkincoarseoperator = yes
#galerkinCoarseOperator = yes
residualTol = 1e-04
residualProj = yes
[VELOCITY]
boundaryTypeMap = inlet, outlet, wall
#solver = pcg+block
boundaryTypeMap = inlet, outlet, wall, wall
density = 1.0
viscosity = -5000.0 #0.25 0.000405
viscosity = -10000.0
residualTol = 1e-06
#[TEMPERATURE]
#solver=none
#boundaryTypeMap = t, O, I, t
#rhoCp = 10
#conductivity = 5000
#residualTol = 1e-08
#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
......@@ -5,6 +5,8 @@
#include <math.h>
#include "udf.hpp"
int my_iostep=100;
/* UDF Functions */
void UDF_LoadKernels(ins_t *ins)
......@@ -14,10 +16,13 @@ void UDF_LoadKernels(ins_t *ins)
void UDF_Setup(ins_t *ins)
{
// get IC from nek
if (!ins->readRestartFile) nek_copyTo(ins, ins->startTime);
if (!ins->readRestartFile) nek_copyTo(ins->startTime);
}
void UDF_ExecuteStep(ins_t *ins, dfloat time, int tstep)
{
if (ins->isOutputStep) nek_userchk();
if (ins->isOutputStep) {
nek_ocopyFrom(time, tstep);
nek_userchk();
}
}
......@@ -61,9 +61,16 @@ c NOTE ::: This routine MAY NOT be called by every process
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
if (cbc(f,e,1).eq.'W '.AND.cbc(f,e,2).eq.'t ') temp=1.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
......@@ -86,28 +93,9 @@ c-----------------------------------------------------------------------
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)
......@@ -126,28 +114,38 @@ c call outpost(vx,vy,vz,t,pr,' ')
endif
c 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')
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
return
end
c-----------------------------------------------------------------------
subroutine userqtl ! Set thermal divergence
call userqtl_scig
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
return
end
......@@ -171,7 +169,7 @@ c-----------------------------------------------------------------------
include 'TOTAL'
integer e,f
c ifheat = .true.
ifheat = .false.
nface = 2*ldim
do e=1,nelv
......@@ -180,28 +178,112 @@ c ifheat = .true.
cbc(f,e,1)='E '
cbc(f,e,2)='E '
elseif(cbc(f,e,1).eq.'v ')then ! z_bottom
boundaryID(f,e) = 1
cbc(f,e,2)='t '
elseif(cbc(f,e,1).eq.'O ')then ! z_top
boundaryID(f,e) = 2
cbc(f,e,1)='o ' ! Dong's outflow
cbc(f,e,2)='O '
elseif(cbc(f,e,1).eq.'W ')then ! pebbles
boundaryID(f,e) = 4
cbc(f,e,2)='f '
elseif(cbc(f,e,1).eq.'SYM')then ! side of cylinder
boundaryID(f,e) = 3
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
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
......
#!/bin/bash
#/bin/bash
bin=nekrs
project_id=nfi114
project_id=csc262
time="00:30"
NEKRS_HYPRE_NUM_THREADS=1
#source $NEKRS_INSTALL_DIR/bin/nekrs.bashrc
export OCCA_VERBOSE
: ${CPUONLY:=0}
export NEKRS_HOME=$HOME/.local/nekrs_v20_dbg
#export OCCA_CACHE_DIR=/gpfs/alpine/proj-shared/${project_id}/${USER}/.occa
export OCCA_CACHE_DIR=/gpfs/alpine/proj-shared/csc262/ylan/cache_v20_dbg2
export NEKRS_HYPRE_NUM_THREADS
export OCCA_CXX="/sw/summit/xl/16.1.1-3/xlC/16.1.1/bin/xlc"
export OCCA_CXXFLAGS="-O3 -qarch=pwr9 -qhot -DUSE_OCCA_MEM_BYTE_ALIGN=64"
export OCCA_LDFLAGS="/sw/summit/xl/16.1.1-3/xlC/16.1.1/lib/libibmc++.a"
#export OMPI_LD_PRELOAD_POSTPEND=$OLCF_SPECTRUM_MPI_ROOT/lib/libmpitrace.so
if [ $# -ne 3 ]; then
echo "usage: [CPUONLY=1] $0 <casename> <number of nodes> <hh:mm>"
exit 0
fi
module load gcc
bin=$NEKRS_HOME/bin/nekrs
case=$1
nodes=$2
gpu_per_node=6
cores_per_socket=21
let nn=$nodes*$gpu_per_node
let ntasks=nn
time=$3
#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
# Bin Min Nodes Max Nodes Max Walltime (Hours) Aging Boost (Days)
# 1 2,765 4,608 24.0 15
# 2 922 2,764 24.0 10
# 3 92 921 12.0 0
# 4 46 91 6.0 0
# 5 1 45 2.0 0
if [ $CPUONLY -eq 1 ]; then
let nn=2*$nodes
let ntasks=$nn*$cores_per_socket
export NEKRS_BACKEND="CPU"
fi
if [ ! -f $bin ]; then
echo "Cannot find" $bin
exit 1
fi