Commit fe6a2def authored by Patrick Shriwise's avatar Patrick Shriwise
Browse files

Merge remote-tracking branch 'origin/1k_openmc_model' into 1k_openmc_model

parents 31d82c60 01e8efc5
[Mesh]
type = NekMesh
[]
[Variables]
[dummy]
[]
[]
[AuxVariables]
[avg_flux]
[]
[]
[Problem]
type = NekProblem
[]
[Executioner]
type = Transient
num_steps = 1000
[]
[Outputs]
[exo]
type = Exodus
output_dimension = 3
interval = 100
execute_on = 'timestep_end'
[]
[]
[Postprocessors]
[total_flux]
type = Receiver
default = 0
[]
[]
[Mesh]
type = FileMesh
file = spheres_generator_in.e
[]
[Kernels]
[hc]
type = HeatConduction
variable = temp
[]
[heat]
type = BodyForce
#value = 35.37
value = 70
variable = temp
[]
[]
[BCs]
inactive = 'outside'
[outside]
type = DirichletBC
variable = temp
boundary = '1'
value = 300
[]
[match_nek]
type = MatchedValueBC
variable = temp
boundary = '1'
v = 'nek_temp'
[]
[]
[Materials]
[hc]
type = GenericConstantMaterial
prop_values = '0.2' # 20 W/mK -> 0.2 W/cmK
prop_names = 'thermal_conductivity'
[]
[]
[Executioner]
type = Transient
petsc_options_iname = '-pc_type -pc_hypre_type'
num_steps = 10000000
petsc_options_value = 'hypre boomeramg'
dt = 1e-4
nl_rel_tol = 1e-5
[]
[Variables]
[temp]
[]
[]
[Outputs]
inteval = 100
exodus = true
[]
[MultiApps]
[nek]
type = TransientMultiApp
app_type = NekApp
input_files = 'nek.i'
[]
[]
[Transfers]
[nek_temp]
type = MultiAppNearestNodeTransfer
source_variable = temp
direction = from_multiapp
multi_app = nek
variable = nek_temp
fixed_meshes = true
[]
[avg_flux]
type = MultiAppNearestNodeTransfer
source_variable = avg_flux
direction = to_multiapp
multi_app = nek
variable = avg_flux
fixed_meshes = true
[]
[total_flux_to_nek]
type = MultiAppPostprocessorTransfer
to_postprocessor = total_flux
direction = to_multiapp
from_postprocessor = total_flux
multi_app = nek
[]
[]
[AuxVariables]
[nek_temp]
[]
[avg_flux]
family = MONOMIAL
order = CONSTANT
initial_condition = 0
[]
[]
[AuxKernels]
[avg_flux]
type = FluxAverageAux
coupled = 'temp'
diffusivity = thermal_conductivity
variable = avg_flux
boundary = '1'
[]
[]
[Postprocessors]
[total_flux_check]
# Should add up to the same thing as total_flux
type = SideIntegralVariablePostprocessor
variable = avg_flux
boundary = 1
[]
[total_flux]
type = SideFluxIntegral
diffusivity = thermal_conductivity
variable = 'temp'
boundary = '1'
[]
[average_flux]
type = SideFluxAverage
diffusivity = thermal_conductivity
variable = 'temp'
boundary = '1'
[]
[]
// Boundary conditions
void insVelocityDirichletConditions3D(bcData *bc)
{
bc->uP = 0.0;
bc->vP = 0.0;
bc->wP = 0.17;
}
void cdsDirichletConditions3D(bcData *bc)
{
bc->sP = 573.0;
}
void cdsNeumannConditions3D(bcData *bc)
{
bc->sF = bc->wrk[bc->idM];
}
// 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 = r5.fld
stopAt = numSteps
numSteps = 20000
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 <iostream>
/* 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);
ins->o_usrwrk.free(); free(ins->usrwrk);
ins->usrwrk = (dfloat*) calloc(ins->mesh->Nelements*ins->mesh->Np, sizeof(dfloat));
ins->o_usrwrk = ins->mesh->device.malloc(ins->mesh->Nelements*ins->mesh->Np*sizeof(dfloat), ins->usrwrk);
ins->cds->o_usrwrk = ins->o_usrwrk;
// printf ("> Test %d \n", ins->mesh->Nelements*ins->mesh->Np*sizeof(dfloat));
}
void UDF_ExecuteStep(ins_t *ins, dfloat time, int tstep)
{
nek_userchk();
//std::cout << "Finished userchk" << std::endl;
int nwbdt=nekData.cbscnrs[0];
auto *wrk1 = nekData.cbscnrs;
wrk1 = wrk1 + nwbdt*5*4 + 2;
ins->o_usrwrk.copyFrom(wrk1,ins->mesh->Nelements*ins->mesh->Np*sizeof(dfloat));
//auto *ptr=nekData.cbscnrs;
//ins->o_usrwrk.copyFrom(wrk, mesh->Nelements*mesh->Np*sizeof(dfloat));
//std::cout << "Finished UDF_ExecuteStep" << std::endl;
}
include 'nek_moose.f'
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 = 573.0
return
end
c-----------------------------------------------------------------------
subroutine userchk
include 'SIZE'
include 'TOTAL'
include 'NEKMOOSE'
c implicit none
COMMON /SCNRS/ SC_NRS(LX1*LY1*LZ1*LELT*7)
real SC_NRS
common /flux_reconstructed/flux_recon(lx1,ly1,lz1,lelt)
integer n_count
save n_count
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
n_count=n_count+1
! The first call is at time step zero
if ((1+1+nw_bdt*4*5+n).gt.(7*n)) then
if (nid.eq.0) write(6,*) "> Insufficient scratch memory"
call exitt()
endif
if (n_count.gt.2) then
if (mod(n_count,2).eq.1) then
if (nid.eq.0) write(6,*) "> DOWNLOADING FLUX"
do i=1,nw_bdt*4
pc_f(i)=sc_nrs(1+i+nw_bdt*4*4)*1.5*1.5
enddo
flux_moose=sc_nrs(1+1+nw_bdt*4*5)
call flux_reconstruction()
do i=1,n
sc_nrs(1+1+nw_bdt*4*5+i)=flux_recon(i,1,1,1)
enddo
endif
if (mod(n_count,2).eq.0) then
if (nid.eq.0) write(6,*) "> LOADING TEMPERATURE"
tmin=glmin(t,n)
tmax=glmax(t,n)
if (nid.eq.0) then
write(6,*) "> Temperature: ", tmin," - ", tmax
endif
call nek_interpolation()
do i=1,nw_bdt*4
sc_nrs(1+i+nw_bdt*4*3)=pc_t(i)
enddo
endif
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'
include 'NEKMOOSE'
COMMON /SCNRS/ SC_NRS(LX1*LY1*LZ1*LELT*7)
real SC_NRS
common /flux_reconstructed/flux_recon(lx1,ly1,lz1,lelt)
integer f,e
n=lx1*ly1*lz1*nelt
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.'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 '
else
write(*,*)'bc incorrect',e,f,cbc(f,e,1)
endif
enddo
enddo
param(18) = 1 ! use algebraic residual norm
call setup_topo
call rzero(flux_recon,n)
call nek_pointscloud()
c storing information in common
c parameter (lsurf_m=60000)
c common /point_cloudx/ pc_x(lsurf_m*4)
c common /point_cloudy/ pc_y(lsurf_m*4)
c common /point_cloudz/ pc_z(lsurf_m*4)
c common /point_cloudt/ pc_t(lsurf_m*4)
c common /point_cloudf/ pc_f(lsurf_m*4)
c common /ref_element/ pc_flag(lsurf_m*4)
c common /tot_surf/ nw_bdt
! Units in cm, radius is 1.5 cm
sc_nrs(1)=nw_bdt
do i=1,nw_bdt*4
sc_nrs(1+i)=pc_x(i)*1.5
sc_nrs(1+i+nw_bdt*4)=pc_y(i)*1.5
sc_nrs(1+i+nw_bdt*4*2)=pc_z(i)*1.5
enddo
if (nid.eq.0) then
do i=1,nw_bdt*4
write (6,*)i,pc_x(i),pc_y(i),pc_z(i)
enddo
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