Commit dead6b19 authored by ylan's avatar ylan
Browse files

[peb mesh] add annulus 3344 case

parent 863986da
## 3344 (non-touching) pebbles in am annulus cylinder (Last updated: Jun. 12, 2020)
cyl3344 N, Version 3.0, E1121214
Running with dt=5e-4 for time > 26
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-flat (?)
- Sphere radius = 1.0 (Wall)
- inner Cylinder radius = 10.468 (W01), (converted to Wall in usr file)
- outer Cylinder radius = 27.743 (W02), (converted to Wall in usr file)
- Zbot of the cylinder = 12.942 (W03), (converted to Outflow in usr file)
- Ztop of the cylinder = -10.022 (W03), (converted to Inflow in usr file)
- 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 -27.743 -27.743 -10.022
xyz max 27.743 27.744 12.942
mesh metrics:
GLL grid spacing min/max : 3.42E-04 4.30E-01
scaled Jacobian min/max/avg: 1.66E-02 9.99E-01 3.88E-01
aspect ratio min/max/avg: 1.06E+00 5.89E+01 1.19E+01
```
- Manage to run to time=26 with dt=5e-4 in NekRS (version untracked)
// 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 = r3.fld
stopAt = numSteps
numSteps = 200000
dt = 5.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 = -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"
#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();
}
}
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
c if (mod(istep,500).ne.0) then
c call compute_cfl(cfl,vx,vy,vz,1.0)
c call outpost(cflf,vy,vz,pr,t,'cfl')
c 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 ' ! Dong's outflow
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-----------------------------------------------------------------------
This source diff could not be displayed because it is too large. You can view the blob instead.
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