[FLASH-BUGS] FLASH 2.2 bugs (fwd)
Sean Matt
matt at physics.mcmaster.ca
Fri Jan 17 09:22:20 CST 2003
I was asked to repost this email, since the post size limit was previously
too small to allow the attachments.
---------- Forwarded message ----------
Date: Thu, 16 Jan 2003 11:55:57 -0500 (EST)
From: Sean Matt <matt at physics.mcmaster.ca>
To: Jonathan Dursi <ljdursi at eugenia.asci.uchicago.edu>
Cc: flash-bugs at flash.uchicago.edu, Robi Banerjee <banerjee at physics.mcmaster.ca>
Subject: Re: [FLASH-BUGS] FLASH 2.2 bugs
Jonathan,
I'm glad to hear about the other bug fixes. The setup we've been
using is one of our own, called "evrard." I've attached to this email the
"Config" "flash.par" "init_block.F90" and "ref_marking.F90" files that we
have in our "evrard" setups directory. Please let me know if you have any
problems receiving them.
-Sean
-------------- next part --------------
# Configuration file for spherical dust collapse problem
# (Colgate, S. A., & White, R. H. 1966, ApJ, 143, 626)
REQUIRES driver/time_dep
REQUIRES hydro
REQUIRES gravity/poisson/multigrid
REQUIRES materials/eos/gamma
REQUIRES mesh/amr/paramesh2.0/second_order
# R_init Initial radius of cloud
# rho_0 Initial density of cloud
# T_ambient Initial ambient temperature (everywhere)
# x/y/zctr Coordinates of the center of the cloud
# delta_ref Refine a block if the max density contrast is
# greater than this
# delta_deref Derefine a block if the max density contrast is
# less than this
PARAMETER T_ambient REAL 1.
PARAMETER R_init REAL 0.05
PARAMETER rho_0 REAL 1.
PARAMETER xctr REAL 0.5
PARAMETER yctr REAL 0.5
PARAMETER zctr REAL 0.5
PARAMETER jeans_ref_sqr REAL 4.
PARAMETER jeans_dref_sqr REAL 64.
PARAMETER delta_ref REAL 0.1
PARAMETER delta_deref REAL 0.01
PARAMETER reference_density REAL 1.
# For testing purposes
# VARIABLE grvx NORENORM # x-component of grav. acceleration
# VARIABLE grvy NORENORM # y-component of grav. acceleration
# VARIABLE grvz NORENORM # z-component of grav. acceleration
-------------- next part --------------
# Runtime parameters for the gas collapse problem.
print_tstep_loc = .true.
reference_density = 1e6
#jeans_ref_sqr = 160000
#jeans_dref_sqr = 64e30
eint_switch = 1e-4
# Parameters for initial model
rho_0 = 9947
R_init = 6E8
T_ambient = 1000.
xctr = 1e5
yctr = 1e5
zctr = 1e5
# Gas ratio of specific heats
gamma = 1.6666667
# Computational volume parameters
# Grid dimensionality and geometry
igeomx = 0
igeomy = 0
igeomz = 0
# Size of computational volume
xmin = 0.
xmax = 2E5
ymin = 0.
ymax = 2E5
zmin = 0.
zmax = 2E5
Nblockx = 1
Nblocky = 1
Nblockz = 1
# Boundary conditions
xl_boundary_type = "reflecting"
xr_boundary_type = "reflecting"
yl_boundary_type = "reflecting"
yr_boundary_type = "reflecting"
zl_boundary_type = "reflecting"
zr_boundary_type = "reflecting"
grav_boundary_type = "periodic"
# Simulation (grid, time, I/O) parameters
cfl = 0.8
lrefine_min = 2
lrefine_max = 5
refine_cutoff = .8
derefine_cutoff = 0.2
basenm = "evrard"
restart = .false.
trstrt = 2621
nrstrt = 2500
tplot = .436
nend = 9999
tmax = 7
refine_var_1 = "dens"
#refine_var_2 = "pres"
mpole_lmax = 5
igrav = 1
plot_var_1 = "dens"
plot_var_2 = "pres"
plot_var_3 = "velx"
plot_var_4 = "temp"
run_comment = "evrard probled"
log_file = "evrard.log"
smlrho = 1.E-10
smallp = 1.E-10
smalle = 1.E-10
smallt = 1.E-12
dtini = 1.E-4
dtmin = 1.E-40
#dtmax = .002621
nriem = 100
cvisc = 0.
#ppm_modifystates = .true.
conserved_var = .true.
monotone = .true.
mgrid_max_residual_norm = 1e-5
mgrid_solve_max_iter = 10000
#mgrid_nsmooth = 20
mgrid_max_vcycles = 1000
-------------- next part --------------
!*******************************************************************************
! Subroutine: init_block()
! Description: Initializes fluid data (density, pressure, velocity, etc.) for
! a specified block. This version sets up the dust cloud
! collapse problem.
! References: Colgate, S. A., & White, R. H. 1966, ApJ, 143, 626
! Monchmeyer, R., & Muller, E. 1989, A&A, 217, 351
! Parameters: block_no The number of the block to initialize
! Runtime: R_init Initial radius of cloud
! rho_0 Initial density of cloud
! T_ambient Initial ambient temperature (everywhere)
! x/y/zctr Coordinates of the center of the cloud
subroutine init_block (block_no)
!===============================================================================
use logfile
use physical_constants
use runtime_parameters
use dBase, ONLY: k2d, k3d, &
& ionmax, nguard, &
& dBasePropertyInteger, &
& dBasePropertyReal, dBaseKeyNumber, &
& dBaseBlockSize, dBaseBlockCoord, &
& dBaseGetDataPtrAllBlocks, &
& dBaseSpecies
implicit none
real, pointer, dimension(:,:,:,:,:) :: data
integer, PARAMETER :: MAXDIMS = 3
logical, save :: first_call = .true.
integer, save :: nxb, nyb, nzb
integer, save :: ivelx, ively, ivelz, itemp
integer, save :: ipres, iener, idens, igame, igamc
integer, save :: inuc_begin
integer, save :: ndim, mype, masterpe
real, save :: xmin, xmax, ymin, ymax, zmin, zmax
real, save :: smallp, r_init, gamma, t_ambient
real, save :: rho_0, smlrho, xctr, yctr, zctr
real :: size(MAXDIMS), coord(MAXDIMS)
integer block_no, i, j, k, n, imax, jmax, kmax, jlo
integer Nint, ii, jj, kk
real delx, xx, dely, yy, delz, zz, velocity, distinv
real vfrac, xdist, ydist, zdist, dist, vctr
real Nintinv, sum_rho, sum_p, sum_vx, sum_vy, sum_vz, &
& Nintinv1
real xxmin, xxmax, yymin, yymax, zzmin, zzmax
real pi, gascon, Newton, m_p, boltz
real sndspd, presfac, ek
integer N_ext
integer N_prof
parameter (N_prof = 10000)
real r_prof(N_prof), rho_prof(N_prof), p_prof(N_prof)
real v_prof(N_prof), dr_prof
save r_prof, rho_prof, p_prof, v_prof
save sndspd, presfac
real :: l, dldr, dmdr, r, m, hdr, const, mconst, deltar, lpred, mpred, B
logical :: iso_set = .true.
integer :: nsubstep
!===============================================================================
! Initialize scalar quantities we will need.
data => dBaseGetDataPtrAllBlocks() ! get soln vector
call get_constant_from_db ("pi", pi)
call get_constant_from_db ("ideal gas constant", gascon)
call get_constant_from_db ("Newton", Newton)
call get_constant_from_db ("proton mass", m_p)
call get_constant_from_db ("Boltzmann", boltz)
if (first_call) then
call get_parm_from_context("xmin", xmin)
call get_parm_from_context("xmax", xmax)
call get_parm_from_context("ymin", ymin)
call get_parm_from_context("ymax", ymax)
call get_parm_from_context("zmax", zmax)
call get_parm_from_context("zmin", zmin)
call get_parm_from_context("smallp", smallp)
call get_parm_from_context("r_init", r_init)
call get_parm_from_context("gamma", gamma)
call get_parm_from_context("t_ambient", t_ambient)
call get_parm_from_context("rho_0", rho_0)
call get_parm_from_context("smlrho", smlrho)
call get_parm_from_context("xctr", xctr)
call get_parm_from_context("yctr", yctr)
call get_parm_from_context("zctr", zctr)
nxb = dBasePropertyInteger("xBlockSize")
nyb = dBasePropertyInteger("yBlockSize")
nzb = dBasePropertyInteger("zBlockSize")
myPE = dBasePropertyInteger("MyProcessor")
masterPE = dBasePropertyInteger("MasterProcessor")
ndim = dBasePropertyInteger("Dimensionality")
ivelx = dBaseKeyNumber("velx")
ively = dBaseKeyNumber("vely")
ivelz = dBaseKeyNumber("velz")
idens = dBaseKeyNumber("dens")
ipres = dBaseKeyNumber("pres")
itemp = dBaseKeyNumber("temp")
igame = dBaseKeyNumber("game")
igamc = dBaseKeyNumber("gamc")
iener = dBaseKeyNumber("ener")
inuc_begin= dBaseSpecies(1)
! first call set to false later ...
end if
imax = 2*nguard + nxb ! Maximum values of the cell index ranges
jmax = 2*nguard*k2d + nyb
kmax = 2*nguard*k3d + nzb
! Coordinates of the edges of the block
size = dBaseBlockSize(block_no)
coord = dBaseBlockCoord(block_no)
xxmax = coord(1) + 0.5*size(1)
xxmin = coord(1) - 0.5*size(1)
yymax = coord(2) + 0.5*size(2)
yymin = coord(2) - 0.5*size(2)
zzmax = coord(3) + 0.5*size(3)
zzmin = coord(3) - 0.5*size(3)
! Cell size
delx = size(1) / real(nxb)
dely = size(2) / real(nyb)
delz = size(3) / real(nzb)
!-------------------------------------------------------------------------------
! Write a message to stdout describing the problem setup.
if (first_call .and. MyPE .eq. MasterPE) then
write (*,*)
call write_logfile &
& ("flash: initializing for gas collapse problem.")
write (*,*) &
& 'flash: initializing for gas collapse problem.'
write (*,*)
1 format (1X, 1P, 4(A13, E12.6, :, 1X))
2 format (1X, 1P, A13, I12)
endif
! Construct the radial samples needed for the initialization.
if (first_call) then
first_call = .false.
! rho2 = fltarr(n+1)
! rad2 = fltarr(n+1)
! mr2 = fltarr(n+1)
deltar = sqrt(xctr*xctr*2)/N_prof
const = 4*pi*Newton*(rho_0)**2/3.3e13/(gamma*(1)**2)
mconst = 4*pi*rho_0*(xctr*.7)**3
nsubstep = 4
l = 0.0
dldr = 0.0
dmdr = 0
r = 0.0
m = 0.0
hdr = deltar/nsubstep*.5
! rho_prof(0) = rho_0*exp(l)
! r_prof(0) = r
!print*,'calculating table'
! do i=1,N_prof
! do j=1,nsubstep
! r = r + hdr
! mpred = m + hdr*dmdr
! lpred = l + hdr*dldr
! dmdr = r**2*exp(lpred)
! dldr = -const*mpred*exp(lpred*(1-gamma))/(r**2)
! r = r + hdr
! m = m + 2*hdr*dmdr
! l = l + 2*hdr*dldr
! dldr = -const*m*exp(l*(1-gamma))/(r**2)
! dmdr = r**2*exp(l)
! enddo
! rho_prof(i) = rho_0*exp(l)
! r_prof(i) = r
! enddo
r= 1
r_prof(1) = 0
rho_prof(1) = 19.5
m = 4.0/3.0*pi*(1)*19.5
hdr=hdr*2
dldr = 0
l=rho_0
do i=1,N_prof
do j=1,nsubstep
if (r < xctr*.4) then
dldr = -1/r**2!-m*newton*(l**(2-gamma))/(gamma*(3.3e13/(rho_0**gamma))*r**2)! *(l**(gamma-2)))
l = 1e10/r
else
!isothermal solution
if (iso_set) then
B = -l*newton*m/(dldr*r**2)
l = l*1e-4
! l=1e-1
print*,'switched to isothermal, B = ',B
iso_set = .false.
endif
dldr =-l/(B)*newton*m/r**2
l = hdr*dldr+l
! print*,'dldr = ',dldr
endif
m = m + 4./3.*pi*((r)**3-(r-hdr)**3)*l
r = hdr+r
!!print*,i,' ',j,' dldr ',dldr,' m ',m,' hdr ',hdr,' l ',l
enddo
!!print*,l,' ',r
!!call abort
!! print*,'dPdr ', gamma*(3.3e13/(rho_0**gamma))*l**(gamma-1)*dldr
!! print*,'other side ',l*newton*m/r**2
rho_prof(i) = l
r_prof(i) = r
enddo
endif
!!
!!-------------------------------------------------------------------------------
! Loop over cells in the block. For each, compute the physical
! position of its left and right edge and its center as well as
! its physical width. Then decide whether it is inside the
! initial radius or outside and initialize the hydro variables
! appropriately.
Nint = 7
Nintinv = 1./float(Nint)
Nintinv1= 1./(float(Nint)-1.)
data(idens,:,:,:,block_no) = smlrho
do k = 1, kmax
do j = 1, jmax
do i = 1, imax
sum_rho = 1e-10
sum_p = 0
sum_vx = 0.
sum_vy = 0.
sum_vz = 0.
do kk = 0, (Nint-1)*k3d
zz = zzmin + delz*(-.5+real(k-nguard-1)+kk*Nintinv1)
zdist = (zz - zctr) * k3d
do jj = 0, (Nint-1)*k2d
yy = yymin + dely*(-.5+real(j-nguard-1)+jj*Nintinv1)
ydist = (yy - yctr) * k2d
do ii = 0, Nint-1
xx = xxmin + delx*(-.5+real(i-nguard-1)+ii*Nintinv1)
xdist = xx - xctr
dist = sqrt( xdist**2 + ydist**2 + zdist**2 )
distinv = 1. / max( dist, 1.E-10 )
! if (dist < (xctr*.8)) then
! sum_rho = 9947*distinv
call find(r_prof,N_prof,dist,jlo)
if(jlo>0) then
sum_rho = sum_rho+rho_prof(jlo)*Nintinv*Nintinv*Nintinv
else
sum_rho = sum_rho+rho_prof(1)*4*Nintinv*Nintinv*Nintinv
endif
! data(itemp,i,j,k,block_no) = 1.5e7
! else
! data(itemp,i,j,k,block_no) = 1e-30
! sum_rho = 1e-5
! endif
! if(distinv == 1/1e-10) then
! sum_rho = 8*Nintinv*Nintinv*Nintinv
! endif
enddo
enddo
enddo
if (dist < (xctr*.4)) then
! data(itemp,i,j,k,block_no) = 1.7e9
! data(ipres,i,j,k,block_no) = (data(idens,i,j,k,block_no)*gascon)*data(itemp,i,j,k,block_no)
data(ipres,i,j,k,block_no) = sum_rho*(.01884)*(gamma-1)
data(itemp,i,j,k,block_no) = data(ipres,i,j,k,block_no) / &
& (data(idens,i,j,k,block_no)*gascon)
else
!isotherma9l9
call find(r_prof,N_prof,xctr*.3999,jlo)
data(ipres,i,j,k,block_no) = rho_prof(jlo)*(.01884)*(gamma-1) !1e-10!B* sum_rho
data(itemp,i,j,k,block_no) = data(ipres,i,j,k,block_no) / &
& (data(idens,i,j,k,block_no)*gascon)
endif
data(itemp,i,j,k,block_no) = 1e-8
data(idens,i,j,k,block_no) = max(sum_rho, &
& smlrho)
! data(ipres,i,j,k,block_no) = max(sum_p, &
! & smallp)
data(ipres,i,j,k,block_no) = (data(idens,i,j,k,block_no)*gascon)*data(itemp,i,j,k,block_no)
!************************** for constant pressure test
! data(ipres,i,j,k,block_no) = smallp * 100.
! data(itemp,i,j,k,block_no) = data(ipres,i,j,k,block_no) / &
! & (data(idens,i,j,k,block_no)*gascon)
!**************************
data(ivelx,i,j,k,block_no) = sum_vx
data(ively,i,j,k,block_no) = sum_vy
data(ivelz,i,j,k,block_no) = sum_vz
!************************** for constant velocity gradient/uniform density test
! data(idens,i,j,k,block_no) = rho_0
! data(ipres,i,j,k,block_no) = rho_0*presfac
! data(ivelx,i,j,k,block_no) = 3.4E9 *
! & (1.-(xxmin+delx*(i-nguard-0.5))/xmax)
! data(ively,i,j,k,block_no) = 0.
! data(ivelz,i,j,k,block_no) = 0.
!**************************
enddo
enddo
enddo
!-------------------------------------------------------------------------------
! Initialize the nuclear abundances. These are not of interest
! in this problem, so we set them to 1. everywhere.
do n = 1, ionmax
do k = 1, kmax
do j = 1, jmax
do i = 1, imax
data(inuc_begin+n-1,i,j,k,block_no) = 1.
enddo
enddo
enddo
enddo
! Compute the gas energy and set the gamma-values needed for
! the equation of state.
do k = 1, kmax
do j = 1, jmax
do i = 1, imax
data(igame,i,j,k,block_no) = gamma
data(igamc,i,j,k,block_no) = gamma
ek = 0.5 * (data(ivelx,i,j,k,block_no)**2 + &
& data(ively,i,j,k,block_no)**2 + &
& data(ivelz,i,j,k,block_no)**2)
data(iener,i,j,k,block_no) = data(ipres,i,j,k,block_no) / &
& (data(igame,i,j,k,block_no)-1.)
data(iener,i,j,k,block_no) = data(iener,i,j,k,block_no) / &
& data(idens,i,j,k,block_no)
data(iener,i,j,k,block_no) = data(iener,i,j,k,block_no) + ek
! data(iener,i,j,k,block_no) = max(data(iener,i,j,k,block_no),smallp)
enddo
enddo
enddo
!===============================================================================
end
!*******************************************************************************
! Routine: find()
! Description: Given a monotonically increasing table x(N) and a test value
! x0, return the index i of the largest table value less than
! or equal to x0 (or 0 if x0 < x(1)). Use binary search.
subroutine find (x, N, x0, i)
integer N, i, il, ir, im
real x(N), x0
if (x0 .lt. x(1)) then
i = 0
elseif (x0 .gt. x(N)) then
i = N
else
il = 1
ir = N
10 if (ir .eq. il+1) goto 20
im = (il + ir) / 2
if (x(im) .gt. x0) then
ir = im
else
il = im
endif
goto 10
20 i = il
endif
return
end
-------------- next part --------------
!******************************************************************************
! Routine: amr_test_refinement()
! Description: Mesh refinement/derefinement testing routine. This routine
! serves as a wrapper for the ref_marking() routine, which uses
! a second-derivative estimator on a chosen variable to
! determine whether each block should be refined or derefined.
! This version also refines on the magnitude of the density.
! Parameters: kmype - My processor #. Needed??
! mrefine_max - Maximum refinement level.
! mrefine_min - Minimum refinement level.
! iref - Variable ID (index to unk()) to refine on.
subroutine amr_test_refinement (kmype, mrefine_max, mrefine_min, &
& iref)
!==============================================================================
use runtime_parameters
implicit none
include 'mpif.h'
integer kmype, mrefine_max, mrefine_min, iref
real refine_cutoff, derefine_cutoff, refine_filter
integer ierr
logical, save :: first_call = .true.
!==============================================================================
if (first_call) then
call get_parm_from_context("refine_cutoff", refine_cutoff)
call get_parm_from_context("derefine_cutoff", derefine_cutoff)
call get_parm_from_context("refine_filter", refine_filter)
first_call = .false.
end if
if (mrefine_max.eq.1) then
#ifdef USEBARS
call MPI_Barrier (MPI_COMM_WORLD, ierr)
#endif
return
end if
call ref_marking (refine_cutoff, derefine_cutoff, refine_filter, &
& mrefine_max, iref)
!==============================================================================
return
end
!******************************************************************************
! Routine: ref_marking()
! Description: Compute the refinement criterion on all blocks and use it to
! determine whether each block should be refined or derefined.
! Written: K. Olson 8/97, based on Lohner (1987), Comput. Meth. in Appl.
! Mech. and Eng., vol. 61, pp. 323-338, and on notes and code
! written by R. deFainchtein.
! Inputs: ctore - Refinement cutoff. Refine only if
! estimator is greater than this value.
! Should be between 0 and 1.
! ctode - Derefinement cutoff. Derefine only if
! estimator is smaller than this value.
! Should be between 0 and 1.
! epsil - Filter coefficient used in computing the
! second-derivative estimator. See the
! User's Guide for details.
! lrefine_max - Maximum allowed level of refinement.
! ivar - Index into unk() of variable to refine on.
! Outputs: Blocks are marked appropriately (using the global variables
! refine() and derefine()). Refinements and derefinements are
! carried out by amr_refine_derefine().
subroutine ref_marking (ctore, ctode, epsil, lrefine_max, iref)
!==============================================================================
use runtime_parameters
use dBase, ONLY: dBaseGetDataPtrAllBlocks, &
& dBaseKeyNumber, &
& il_bnd, iu_bnd, jl_bnd, ju_bnd, &
& kl_bnd, ku_bnd, ndim, maxblocks, &
& nguard, k2d, k3d, &
& nxb, nyb, nzb
use dBaseIncludes, ONLY: lnblocks, nodetype, &
& bsize, neigh, parent, nchild, &
& child, refine, derefine, stay, &
& lrefine
implicit none
include 'mpif.h'
integer ndim2
parameter (ndim2=ndim*ndim)
logical first
real delx,dely,delz
real delu(ndim,il_bnd:iu_bnd,jl_bnd:ju_bnd,kl_bnd:ku_bnd)
real delua(ndim,il_bnd:iu_bnd,jl_bnd:ju_bnd,kl_bnd:ku_bnd)
real delu2(ndim2), delu3(ndim2), delu4(ndim2)
real num,denom,epsil,error(maxblocks),error_par(maxblocks), &
& errort
real ctore,ctode,J_sq,xx,dd
real error_max
real vx,vy,vz,bv,bx,by,bz,bx1,by1,bz1,rhbi,rhi,bk,et
real px,py,pz,rho,b2,v2,gf,vv,gr,dx,vt,press
real times,time_exe
integer lb,i,j,k,kk,ivar,lrefine_max,iref,nref
integer mype
integer ineigh,ineigh_proc,ierr
integer kstart,kend,jstart,jend,istart,iend
integer nsend,nrecv
integer reqr(2*maxblocks),reqs(2*maxblocks)
integer statr(MPI_STATUS_SIZE,maxblocks)
integer stats(MPI_STATUS_SIZE,maxblocks)
real :: delta_max(maxblocks), block_mass
real :: delta_max_par(maxblocks)
logical, save :: first_call = .true.
real, save :: delta_ref, delta_deref, avg_den
integer, save :: idens
real, pointer, dimension(:,:,:,:,:) :: unk
!==============================================================================
unk => dBaseGetDataPtrAllBlocks()
if (first_call) then
call get_parm_from_context &
& (global_parm_context, "delta_ref", delta_ref)
call get_parm_from_context &
& (global_parm_context, "delta_deref", delta_deref)
call get_parm_from_context &
& (global_parm_context, "reference_density", avg_den)
idens = dBaseKeyNumber("dens")
first_call = .false.
endif
call MPI_Comm_Rank (MPI_Comm_World, mype, ierr)
do lb = 1,lnblocks
error(lb) = 0.
delta_max(lb) = -1.E99
if (nodetype(lb).eq.1.or.nodetype(lb).eq.2) then
delx = 0.5e0*float(nxb)/bsize(1,lb)
#if N_DIM >= 2
dely = 0.5e0*float(nyb)/bsize(2,lb)
#endif
#if N_DIM == 3
delz = 0.5e0*float(nzb)/bsize(3,lb)
#endif
! Compute max density contrast relative to reference_density
do k = 1+k3d*nguard,nzb+k3d*nguard
do j = 1+k2d*nguard,nyb+k2d*nguard
do i = 1+nguard,nxb+nguard
delta_max(lb) = max(delta_max(lb), &
& unk(idens,i,j,k,lb)/avg_den-1.)
enddo
enddo
enddo
! Compute first derivatives
do k = 1+k3d,nzb+(k3d*((2*nguard)-1))
do j = 1+k2d,nyb+(k2d*((2*nguard)-1))
do i = 2,nxb+(2*nguard)-1
! d/dx
delu(1,i,j,k) = unk(iref,i+1,j,k,lb) - &
& unk(iref,i-1,j,k,lb)
delu(1,i,j,k) = delu(1,i,j,k)*delx
delua(1,i,j,k) = abs(unk(iref,i+1,j,k,lb)) + &
& abs(unk(iref,i-1,j,k,lb))
delua(1,i,j,k) = delua(1,i,j,k)*delx
#if N_DIM >= 2
! d/dy
delu(2,i,j,k) = unk(iref,i,j+1,k,lb) - &
& unk(iref,i,j-1,k,lb)
delu(2,i,j,k) = delu(2,i,j,k)*dely
delua(2,i,j,k) = abs(unk(iref,i,j+1,k,lb)) + &
& abs(unk(iref,i,j-1,k,lb))
delua(2,i,j,k) = delua(2,i,j,k)*dely
#endif
#if N_DIM == 3
! d/dz
delu(3,i,j,k) = unk(iref,i,j,k+1,lb) - &
& unk(iref,i,j,k-1,lb)
delu(3,i,j,k) = delu(3,i,j,k)*delz
delua(3,i,j,k) = abs(unk(iref,i,j,k+1,lb)) + &
& abs(unk(iref,i,j,k-1,lb))
delua(3,i,j,k) = delua(3,i,j,k)*delz
#endif
end do
end do
end do
! Compute second derivatives
! Test if at a block boundary
! Two guardcells
kstart = 2*k3d+1
kend = nzb+(k3d*((2*nguard)-2))
jstart = 2*k2d+1
jend = nyb+(k2d*((2*nguard)-2))
istart = 3
iend = nxb+(2*nguard)-2
! One guardcell
! kstart = 2*k3d+1+k3d
! kend = nzb+(k3d*((2*nguard)-2))-k3d
! jstart = 2*k2d+1+k2d
! jend = nyb+(k2d*((2*nguard)-2))-k2d
! istart = nguard
! iend = nxb+(2*nguard)-3
! No guardcells
! kstart = k3d*nguard+1
! kend = nzb+k3d*nguard
! jstart = k2d*nguard+1
! jend = nyb+k2d*nguard
! istart = nguard+1
! iend = nxb+nguard
if (neigh(1,1,lb).le.-20) istart = nguard+1
if (neigh(1,2,lb).le.-20) iend = nguard+nxb
#if N_DIM >= 2
if (neigh(1,3,lb).le.-20) jstart = nguard*k2d+1
if (neigh(1,4,lb).le.-20) jend = nguard*k2d+nyb
#endif
#if N_DIM == 3
if (neigh(1,5,lb).le.-20) kstart = nguard*k3d+1
if (neigh(1,6,lb).le.-20) kend = nguard*k3d+nzb
#endif
do k = kstart,kend
do j = jstart,jend
do i = istart,iend
! d/dxdx
delu2(1) = delu(1,i+1,j,k) - delu(1,i-1,j,k)
delu2(1) = delu2(1)*delx
delu3(1) = abs(delu(1,i+1,j,k)) + &
& abs(delu(1,i-1,j,k))
delu3(1) = delu3(1)*delx
delu4(1) = delua(1,i+1,j,k) + delua(1,i-1,j,k)
delu4(1) = delu4(1)*delx
#if N_DIM >= 2
! d/dydx
delu2(2) = delu(1,i,j+1,k) - delu(1,i,j-1,k)
delu2(2) = delu2(2)*dely
delu3(2) = abs(delu(1,i,j+1,k)) + &
& abs(delu(1,i,j-1,k))
delu3(2) = delu3(2)*dely
delu4(2) = delua(1,i,j+1,k) + delua(1,i,j-1,k)
delu4(2) = delu4(2)*dely
! d/dxdy
delu2(3) = delu(2,i+1,j,k) - delu(2,i-1,j,k)
delu2(3) = delu2(3)*delx
delu3(3) = abs(delu(2,i+1,j,k)) + &
& abs(delu(2,i-1,j,k))
delu3(3) = delu3(3)*delx
delu4(3) = delua(2,i+1,j,k) + delua(2,i-1,j,k)
delu4(3) = delu4(3)*delx
! d/dydy
delu2(4) = delu(2,i,j+1,k) - delu(2,i,j-1,k)
delu2(4) = delu2(4)*dely
delu3(4) = abs(delu(2,i,j+1,k)) + &
& abs(delu(2,i,j-1,k))
delu3(4) = delu3(4)*dely
delu4(4) = delua(2,i,j+1,k) + delua(2,i,j-1,k)
delu4(4) = delu4(4)*dely
#endif
#if N_DIM == 3
! d/dzdx
delu2(5) = delu(1,i,j,k+1) - delu(1,i,j,k-1)
delu2(5) = delu2(5)*delz
delu3(5) = abs(delu(1,i,j,k+1)) + &
& abs(delu(1,i,j,k-1))
delu3(5) = delu3(5)*delz
delu4(5) = delua(1,i,j,k+1) + delua(1,i,j,k-1)
delu4(5) = delu4(5)*delz
! d/dzdy
delu2(6) = delu(2,i,j,k+1) - delu(2,i,j,k-1)
delu2(6) = delu2(6)*delz
delu3(6) = abs(delu(2,i,j,k+1)) + &
& abs(delu(2,i,j,k-1))
delu3(6) = delu3(6)*delz
delu4(6) = delua(2,i,j,k+1) + delua(2,i,j,k-1)
delu4(6) = delu4(6)*delz
! d/dxdz
delu2(7) = delu(3,i+1,j,k) - delu(3,i-1,j,k)
delu2(7) = delu2(7)*delx
delu3(7) = abs(delu(3,i+1,j,k)) + &
& abs(delu(3,i-1,j,k))
delu3(7) = delu3(7)*delx
delu4(7) = delua(3,i+1,j,k) + delua(3,i-1,j,k)
delu4(7) = delu4(7)*delx
! d/dydz
delu2(8) = delu(3,i,j+1,k) - delu(3,i,j-1,k)
delu2(8) = delu2(8)*dely
delu3(8) = abs(delu(3,i,j+1,k)) + &
& abs(delu(3,i,j-1,k))
delu3(8) = delu3(8)*dely
delu4(8) = delua(3,i,j+1,k) + delua(3,i,j-1,k)
delu4(8) = delu4(8)*dely
! d/dzdz
delu2(9) = delu(3,i,j,k+1) - delu(3,i,j,k-1)
delu2(9) = delu2(9)*delz
delu3(9) = abs(delu(3,i,j,k+1)) + &
& abs(delu(3,i,j,k-1))
delu3(9) = delu3(9)*delz
delu4(9) = delua(3,i,j,k+1) + delua(3,i,j,k-1)
delu4(9) = delu4(9)*delz
#endif
! compute the error
num = 0.
denom = 0.
do kk = 1, ndim2
num = num + delu2(kk)**2
denom = denom + (delu3(kk) + &
& (epsil*delu4(kk)+1.e-20))**2
end do
! mz -- compare the square of the error
error(lb) = max(error(lb),num/denom)
end do
end do
end do
! store the maximum error for the current block
error(lb) = sqrt(error(lb))
end if
end do
! MARK FOR REFINEMENT OR DEREFINEMENT
! first communicate error of parent to its children
! parents collect messages from children
error_par(1:lnblocks) = 0.
delta_max_par(1:lnblocks) = 0.
nrecv = 0
do lb = 1,lnblocks
if(parent(1,lb).gt.-1) then
if (parent(2,lb).ne.mype) then
nrecv = nrecv + 1
call MPI_IRecv(error_par(lb), &
& 1, &
& MPI_DOUBLE_PRECISION, &
& parent(2,lb), &
& lb, &
& MPI_COMM_WORLD, &
& reqr(nrecv), &
& ierr)
nrecv = nrecv + 1
call MPI_IRecv(delta_max_par(lb), &
& 1, &
& MPI_DOUBLE_PRECISION, &
& parent(2,lb), &
& lb, &
& MPI_COMM_WORLD, &
& reqr(nrecv), &
& ierr)
else
error_par(lb) = error(parent(1,lb))
delta_max_par(lb) = delta_max(parent(1,lb))
end if
end if
end do
! parents send error to children
nsend = 0
do lb = 1,lnblocks
do j = 1,nchild
if(child(1,j,lb).gt.-1) then
if (child(2,j,lb).ne.mype) then
nsend = nsend + 1
call MPI_ISend(error(lb), &
& 1, &
& MPI_DOUBLE_PRECISION, &
& child(2,j,lb), & ! PE TO SEND TO
& child(1,j,lb), & ! THIS IS THE TAG
& MPI_COMM_WORLD, &
& reqs(nsend), &
& ierr)
nsend = nsend + 1
call MPI_ISend(delta_max(lb), &
& 1, &
& MPI_DOUBLE_PRECISION, &
& child(2,j,lb), & ! PE TO SEND TO
& child(1,j,lb), & ! THIS IS THE TAG
& MPI_COMM_WORLD, &
& reqs(nsend), &
& ierr)
end if
end if
end do
end do
if (nsend.gt.0) then
call MPI_Waitall (nsend, reqs, stats, ierr)
end if
if (nrecv.gt.0) then
call MPI_Waitall (nrecv, reqr, statr, ierr)
end if
do lb = 1,lnblocks
if (nodetype(lb).eq.1) then
! test for derefinement
if (.not.refine(lb).and..not.stay(lb) &
& .and.error(lb).le.ctode &
& .and.error_par(lb).le.ctode &
& .and.delta_max(lb).le.delta_deref &
& .and.delta_max_par(lb).le.delta_deref) then
derefine(lb) = .TRUE.
else
derefine(lb) = .FALSE.
end if
! test for refinement
if (error(lb).gt.ctore &
& .or.delta_max(lb).gt.delta_ref) then
derefine(lb) = .FALSE.
refine(lb) = .TRUE.
end if
if (error(lb).gt.ctode.or.error_par(lb).gt.ctode &
& .or.delta_max(lb).gt.delta_deref &
& .or.delta_max_par(lb).gt.delta_deref) &
& stay(lb) = .TRUE.
if (lrefine(lb).ge.lrefine_max) &
& refine(lb) = .FALSE.
end if
end do
!==============================================================================
return
end
More information about the flash-bugs
mailing list