From 028a4e4316573b543c1af8403703925f8951d4d7 Mon Sep 17 00:00:00 2001 From: Marsha Berger Date: Thu, 20 Aug 2015 17:13:05 -0700 Subject: [PATCH 01/53] fixed indexing error --- src/2d/filpatch.f90 | 12 +++++++++++- 1 file changed, 11 insertions(+), 1 deletion(-) diff --git a/src/2d/filpatch.f90 b/src/2d/filpatch.f90 index d6206803b..0c5d73b9a 100644 --- a/src/2d/filpatch.f90 +++ b/src/2d/filpatch.f90 @@ -38,7 +38,8 @@ recursive subroutine filrecur(level,nvar,valbig,aux,naux,t,mitot,mjtot, & integer :: refinement_ratio_x, refinement_ratio_y integer :: unset_indices(4) real(kind=8) :: dx_fine, dy_fine, dx_coarse, dy_coarse - real(kind=8) :: xlow_coarse,ylow_coarse, xlow_fine, ylow_fine, xhi_fine,yhi_fine + real(kind=8) :: xlow_coarse,ylow_coarse, xlow_fine, ylow_fine, xhi_fine,yhi_fine + real(kind=8) :: xcent_fine, xcent_coarse, ycent_fine, ycent_coarse ! Interpolation variables real(kind=8) :: eta1, eta2, valp10, valm10, valc, valp01, valm01, dupc, dumc @@ -183,11 +184,20 @@ recursive subroutine filrecur(level,nvar,valbig,aux,naux,t,mitot,mjtot, & eta1 = (-0.5d0 + real(mod(i_fine - 1, refinement_ratio_x),kind=8)) & / real(refinement_ratio_x,kind=8) + + xcent_coarse = xlow_coarse + (i_coarse-.5d0)*dx_coarse + xcent_fine = xlower + (i_fine-1+ilo + .5d0)*dx_fine + eta1 = (xcent_fine-xcent_coarse) + do j_fine = 1,mjtot_patch j_coarse = 2 + (j_fine - (unset_indices(3) - jlo) - 1) / refinement_ratio_y eta2 = (-0.5d0 + real(mod(j_fine - 1, refinement_ratio_y),kind=8)) & / real(refinement_ratio_y,kind=8) + ycent_coarse = ylow_coarse + (j_coarse-.5d0)*dy_coarse + ycent_fine = ylower + (j_fine-1+jlo + .5d0)*dy_fine + eta2 = (ycent_fine-ycent_coarse) + if (flaguse(i_fine,j_fine) == 0) then do n=1,nvar From 1e2d77c86bd09d9b8f7972e2bdc47add8d44d2a4 Mon Sep 17 00:00:00 2001 From: Marsha Berger Date: Sat, 22 Aug 2015 15:29:25 -0700 Subject: [PATCH 02/53] larger format for ascii otuput --- src/2d/valout.f | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/2d/valout.f b/src/2d/valout.f index 2e514907c..68c3797ab 100644 --- a/src/2d/valout.f +++ b/src/2d/valout.f @@ -109,12 +109,12 @@ subroutine valout (lst, lend, time, nvar, naux) write(matunit1,1004) & xlow,hxposs(level) endif - 1002 format(e18.8,' xlow', /, - & e18.8,' ylow', /, - & e18.8,' dx', /, - & e18.8,' dy',/) - 1004 format(e18.8,' xlow', /, - & e18.8,' dx', /) + 1002 format(e26.16,' xlow', /, + & e26.16,' ylow', /, + & e26.16,' dx', /, + & e26.16,' dy',/) + 1004 format(e26.16,' xlow', /, + & e26.16,' dx', /) if (output_format == 1) then From 8c6dce43d954cc91b2637a554e735835fa700274 Mon Sep 17 00:00:00 2001 From: Marsha Berger Date: Sat, 22 Aug 2015 15:37:35 -0700 Subject: [PATCH 03/53] new test dir - should get zero error for linear solution --- dev/advection_2d_linearTest/Makefile | 159 +++++++++ dev/advection_2d_linearTest/bc2amr.f | 351 +++++++++++++++++++ dev/advection_2d_linearTest/flag2refine2.f90 | 107 ++++++ dev/advection_2d_linearTest/qinit.f | 25 ++ dev/advection_2d_linearTest/qtrue.f | 30 ++ dev/advection_2d_linearTest/setplot.py | 169 +++++++++ dev/advection_2d_linearTest/setprob.f | 22 ++ dev/advection_2d_linearTest/setrun.py | 348 ++++++++++++++++++ 8 files changed, 1211 insertions(+) create mode 100644 dev/advection_2d_linearTest/Makefile create mode 100644 dev/advection_2d_linearTest/bc2amr.f create mode 100644 dev/advection_2d_linearTest/flag2refine2.f90 create mode 100644 dev/advection_2d_linearTest/qinit.f create mode 100644 dev/advection_2d_linearTest/qtrue.f create mode 100644 dev/advection_2d_linearTest/setplot.py create mode 100644 dev/advection_2d_linearTest/setprob.f create mode 100644 dev/advection_2d_linearTest/setrun.py diff --git a/dev/advection_2d_linearTest/Makefile b/dev/advection_2d_linearTest/Makefile new file mode 100644 index 000000000..ac83816c6 --- /dev/null +++ b/dev/advection_2d_linearTest/Makefile @@ -0,0 +1,159 @@ + +# Makefile for Clawpack code in this directory. +# This version only sets the local files and frequently changed +# options, and then includes the standard makefile pointed to by CLAWMAKE. +CLAWMAKE = $(CLAW)/clawutil/src/Makefile.common + +# See the above file for details and a list of make options, or type +# make .help +# at the unix prompt. + + +# Adjust these variables if desired: +# ---------------------------------- + +CLAW_PKG = amrclaw # Clawpack package to use +EXE = xamr # Executable to create +SETRUN_FILE = setrun.py # File containing function to make data +OUTDIR = _output # Directory for output +SETPLOT_FILE = setplot.py # File containing function to set plots +PLOTDIR = _plots # Directory for plots + +OVERWRITE ?= True # False ==> make a copy of OUTDIR first +RESTART ?= False # Should = clawdata.restart in setrun + +# Environment variable FC should be set to fortran compiler, e.g. gfortran + +# Compiler flags can be specified here or set as an environment variable +FFLAGS ?= + + +# --------------------------------- +# List of sources for this program: +# --------------------------------- + +AMRLIB = $(CLAW)/amrclaw/src/2d + +MODULES = \ + $(AMRLIB)/amr_module.f90 \ + $(AMRLIB)/regions_module.f90 \ + $(AMRLIB)/gauges_module.f90 + +SOURCES = \ + qinit.f \ + qtrue.f \ + setprob.f \ + bc2amr.f \ + flag2refine2.f90 \ + $(CLAW)/riemann/src/rpn2_advection.f \ + $(CLAW)/riemann/src/rpt2_advection.f \ + $(AMRLIB)/amr2.f90 \ + $(AMRLIB)/setaux.f90 \ + $(AMRLIB)/b4step2.f90 \ + $(AMRLIB)/qad.f \ + $(AMRLIB)/src2.f90 \ + $(AMRLIB)/src1d.f90 \ + $(AMRLIB)/advanc.f \ + $(AMRLIB)/bound.f90 \ + $(AMRLIB)/stepgrid.f \ + $(AMRLIB)/stepgrid_dimSplit.f \ + $(AMRLIB)/auxcoarsen.f \ + $(AMRLIB)/fixcapaq.f \ + $(AMRLIB)/estdt.f \ + $(AMRLIB)/init_iflags.f \ + $(AMRLIB)/igetsp.f \ + $(AMRLIB)/reclam.f \ + $(AMRLIB)/birect.f \ + $(AMRLIB)/cleanup.f \ + $(AMRLIB)/colate2.f \ + $(AMRLIB)/bufnst2.f \ + $(AMRLIB)/spest2.f \ + $(AMRLIB)/flagregions2.f90 \ + $(AMRLIB)/errest.f \ + $(AMRLIB)/errf1.f \ + $(AMRLIB)/gfixup.f \ + $(AMRLIB)/filval.f90 \ + $(AMRLIB)/filpatch2.f90 \ + $(AMRLIB)/prefilp.f90 \ + $(AMRLIB)/flglvl2.f \ + $(AMRLIB)/flagger.f \ + $(AMRLIB)/prepregstep.f \ + $(AMRLIB)/prepbigstep.f \ + $(AMRLIB)/fluxad.f \ + $(AMRLIB)/fluxsv.f \ + $(AMRLIB)/ginit.f \ + $(AMRLIB)/grdfit2.f \ + $(AMRLIB)/intfil.f \ + $(AMRLIB)/moment.f \ + $(AMRLIB)/nestck2.f \ + $(AMRLIB)/prepf.f \ + $(AMRLIB)/prepc.f \ + $(AMRLIB)/projec2.f \ + $(AMRLIB)/signs.f \ + $(AMRLIB)/findcut.f \ + $(AMRLIB)/smartbis.f \ + $(AMRLIB)/putnod.f \ + $(AMRLIB)/putsp.f \ + $(AMRLIB)/regrid.f \ + $(AMRLIB)/setgrd.f \ + $(AMRLIB)/setuse.f \ + $(AMRLIB)/stst1.f \ + $(AMRLIB)/tick.f \ + $(AMRLIB)/trimbd.f90 \ + $(AMRLIB)/update.f \ + $(AMRLIB)/nodget.f \ + $(AMRLIB)/upbnd.f \ + $(AMRLIB)/basic.f \ + $(AMRLIB)/outval.f \ + $(AMRLIB)/copysol.f \ + $(AMRLIB)/outvar.f \ + $(AMRLIB)/outmsh.f \ + $(AMRLIB)/outtre.f \ + $(AMRLIB)/domain.f \ + $(AMRLIB)/setflags.f \ + $(AMRLIB)/shiftset2.f \ + $(AMRLIB)/conck.f \ + $(AMRLIB)/domshrink.f \ + $(AMRLIB)/domprep.f \ + $(AMRLIB)/domup.f \ + $(AMRLIB)/domcopy.f \ + $(AMRLIB)/coarsen.f \ + $(AMRLIB)/intcopy.f \ + $(AMRLIB)/preintcopy.f \ + $(AMRLIB)/icall.f \ + $(AMRLIB)/preicall.f \ + $(AMRLIB)/step2.f90 \ + $(AMRLIB)/step2x.f90 \ + $(AMRLIB)/step2y.f90 \ + $(AMRLIB)/flux2.f \ + $(AMRLIB)/flux2_dimSplit.f \ + $(AMRLIB)/inlinelimiter.f \ + $(AMRLIB)/cstore.f \ + $(AMRLIB)/saveqc.f \ + $(AMRLIB)/valout.f \ + $(AMRLIB)/check.f \ + $(AMRLIB)/restrt.f \ + $(AMRLIB)/quick_sort1.f \ + $(AMRLIB)/opendatafile.f \ + $(AMRLIB)/addflags.f \ + $(AMRLIB)/domgrid.f \ + $(AMRLIB)/drivesort.f \ + $(AMRLIB)/flagcheck.f \ + $(AMRLIB)/setdomflags.f \ + $(AMRLIB)/setIndices.f \ + $(AMRLIB)/setPhysBndry.f \ + $(AMRLIB)/setPhysBndryFlags.f \ + $(AMRLIB)/coarseGridFlagSet.f \ + $(AMRLIB)/griddomcopy.f \ + $(AMRLIB)/griddomshrink.f \ + $(AMRLIB)/griddomup.f \ + $(AMRLIB)/baseCheck.f \ + $(AMRLIB)/init_alloc.f90 \ + $(AMRLIB)/restrt_alloc.f90 \ + $(AMRLIB)/resize_alloc.f90 \ +# $(AMRLIB)/resize_alloc_static.f90 \ + +#------------------------------------------------------------------- +# Include Makefile containing standard definitions and make options: +include $(CLAWMAKE) + diff --git a/dev/advection_2d_linearTest/bc2amr.f b/dev/advection_2d_linearTest/bc2amr.f new file mode 100644 index 000000000..22d3c8094 --- /dev/null +++ b/dev/advection_2d_linearTest/bc2amr.f @@ -0,0 +1,351 @@ +c +c Modified to give inflow boundary conditions at the left edge +c +c ------------------------------------------------------------------ +c + subroutine bc2amr(val,aux,nrow,ncol,meqn,naux, + 1 hx, hy, level, time, + 2 xleft, xright, ybot, ytop, + 3 xlower, ylower,xupper,yupper, + 4 xperiodic, yperiodic,spheredom) + +c +c +c :::::::::: bc2amr ::::::::::::::::::::::::::::::::::::::::::::::; +c +c Take a grid patch with mesh widths hx,hy, of dimensions nrow by +c ncol, and set the values of any piece of +c of the patch which extends outside the physical domain +c using the boundary conditions. +c +c ------------------------------------------------ +c # Standard boundary condition choices for amr2ez in clawpack +c +c # At each boundary k = 1 (left), 2 (right), 3 (bottom), 4 (top): +c # mthbc(k) = 0 for user-supplied BC's (must be inserted!) +c # = 1 for zero-order extrapolation +c # = 2 for periodic boundary conditions +c # = 3 for solid walls, assuming this can be implemented +c # by reflecting the data about the boundary and then +c # negating the 2'nd (for k=1,2) or 3'rd (for k=3,4) +c # component of q. +c # = 5 sphere bcs (left half maps to right half of same +c # side, and vice versa), as if domain folded in half +c ------------------------------------------------ +c +c The corners of the grid patch are at +c (xleft,ybot) -- lower left corner +c (xright,ytop) -- upper right corner +c +c The physical domain itself is a rectangle bounded by +c (xlower,ylower) -- lower left corner +c (xupper,yupper) -- upper right corner +c +c the picture is the following: +c +c _____________________ (xupper,yupper) +c | | +c _________ (xright,ytop) | +c | | | | +c | | | | +c | | | | +c |___|____| | +c (xleft,ybot) | | +c | | +c |_____________________| +c (xlower,ylower) +c +c +c Any cells that lie outside the physical domain are ghost cells whose +c values should be set in this routine. This is tested for by comparing +c xleft with xlower to see if values need to be set at the left, as in +c the figure above, and similarly at the other boundaries. +c +c Patches are guaranteed to have at least 1 row of cells filled +c with interior values so it is possible to extrapolate. +c Fix trimbd if you want more than 1 row pre-set. +c +c Make sure the order the boundaries are specified is correct +c so that diagonal corner cells are also properly taken care of. +c +c Periodic boundaries are set before calling this routine, so if the +c domain is periodic in one direction only you +c can safely extrapolate in the other direction. +c +c Don't overwrite ghost cells in periodic directions! +c +c ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::; + + use amr_module, only: mthbc + + implicit double precision (a-h,o-z) + + common /cparam/ ubar,vbar + + dimension val(meqn,nrow,ncol), aux(naux,nrow,ncol) + logical xperiodic, yperiodic, spheredom + + hxmarg = hx*.01 + hymarg = hy*.01 + + if (xperiodic .and. (yperiodic .or. spheredom)) go to 499 +c +c +c------------------------------------------------------- +c # left boundary: +c------------------------------------------------------- + if (xleft .ge. xlower-hxmarg) then +c # not a physical boundary -- no cells at this edge lies +c # outside the physical bndry. +c # values are set elsewhere in amr code. + go to 199 + endif +c +c # number of grid cells from this patch lying outside physical domain: + nxl = (xlower+hxmarg-xleft)/hx + if (nxl > 2) then + write(6,*) '*** unexpected value nxl = ',nxl + stop + endif +c + go to (100,110,120,130) mthbc(1)+1 +c + 100 continue +! no longer need next test, using eval of exact solution instead +! if (ubar < 0.1d0*vbar) then +! write(6,*) 'inflow BCs at left boundary assume ubar >= vbar/10' +! stop +! endif + xcell1 = xleft + (0.5d0)*hx + xcell2 = xleft + (1.5d0)*hx + do j = 1,ncol + ycell = ybot + (j-0.5d0)*hy + if (nxl >= 1) then + ! first ghost cell: + !tau1 = hx / (2.d0*ubar) + !t1 = time + tau1 ! time at which to evaluate BC + !y1 = ycell + vbar*tau1 ! y at which to evaluate BC + !val(1,nxl,j) = qtrue(0.d0,y1,t1) + val(1,nxl,j) = qtrue(xcell2,ycell,time) + endif + if (nxl == 2) then + ! second ghost cell: + !tau2 = 3.d0 * hx / (2.d0*ubar) + !t2 = time + tau2 ! time at which to evaluate BC + !y2 = ycell + vbar*tau2 ! y at which to evaluate BC + !val(1,1,j) = qtrue(0.d0,y2,t2) + val(1,1,j) = qtrue(xcell1,ycell,time) + endif + enddo + + go to 199 +c + 110 continue +c # zero-order extrapolation: + do 115 j = 1,ncol + do 115 i=1,nxl + do 115 m=1,meqn + val(m,i,j) = val(m,nxl+1,j) + 115 continue + go to 199 + + 120 continue +c # periodic: handled elsewhere in amr + go to 199 + + 130 continue +c # solid wall (assumes 2'nd component is velocity or momentum in x): + do 135 j = 1,ncol + do 135 i=1,nxl + do 135 m=1,meqn + val(m,i,j) = val(m,2*nxl+1-i,j) + 135 continue +c # negate the normal velocity: + do 136 j = 1,ncol + do 136 i=1,nxl + val(2,i,j) = -val(2,i,j) + 136 continue + go to 199 + + 199 continue +c +c------------------------------------------------------- +c # right boundary: +c------------------------------------------------------- + if (xright .le. xupper+hxmarg) then +c # not a physical boundary -- no cells at this edge lies +c # outside the physical bndry. +c # values are set elsewhere in amr code. + go to 299 + endif +c +c # number of grid cells lying outside physical domain: + nxr = (xright - xupper + hxmarg)/hx + ibeg = max0(nrow-nxr+1, 1) +c + go to (200,210,220,230) mthbc(2)+1 +c + 200 continue +c # user-specified boundary conditions go here in place of error output + write(6,*) + & '*** ERROR *** mthbc(2)=0 and no BCs specified in bc2amr' + stop + go to 299 + + 210 continue +c # zero-order extrapolation: + do 215 j = 1,ncol + do 215 i=ibeg,nrow + do 215 m=1,meqn + !val(m,i,j) = val(m,ibeg-1,j) + val(m,i,j) = 2.d0*val(m,ibeg-1,j) - val(m,ibeg-2,j) + 215 continue + go to 299 + + 220 continue +c # periodic: handled elsewhere in amr + go to 299 + + 230 continue +c # solid wall (assumes 2'nd component is velocity or momentum in x): + do 235 j = 1,ncol + do 235 i=ibeg,nrow + do 235 m=1,meqn + val(m,i,j) = val(m,2*ibeg-1-i,j) + 235 continue +c # negate the normal velocity: + do 236 j = 1,ncol + do 236 i=ibeg,nrow + val(2,i,j) = -val(2,i,j) + 236 continue + go to 299 + + 299 continue +c +c------------------------------------------------------- +c # bottom boundary: +c------------------------------------------------------- + if (ybot .ge. ylower-hymarg) then +c # not a physical boundary -- no cells at this edge lies +c # outside the physical bndry. +c # values are set elsewhere in amr code. + go to 399 + endif +c +c # number of grid cells lying outside physical domain: + nyb = (ylower+hymarg-ybot)/hy +c + go to (300,310,320,330) mthbc(3)+1 +c + 300 continue +! if (vbar < 0.1d0*ubar) then +! write(6,*) 'inflow BCs at bottom assume vbar >= ubar/10' +! stop +! endif + ycell1 = ybot + .5d0*hy + ycell2 = ybot + 1.5d0*hy + do i = 1,nrow + xcell = xleft + (i-0.5d0)*hx + if (nyb >= 1) then + ! first ghost cell: + !tau1 = hy / (2.d0*vbar) + !t1 = time + tau1 ! time at which to evaluate BC + !x1 = xcell + vbar*tau1 ! x at which to evaluate BC + !val(1,i,nyb) = qtrue(x1,0.d0,t1) + val(1,i,nyb) = qtrue(xcell,ycell2,time) + endif + if (nyb == 2) then + ! second ghost cell: + !tau2 = 3.d0 * hy / (2.d0*vbar) + !t2 = time + tau2 ! time at which to evaluate BC + !x2 = xcell + ubar*tau2 ! x at which to evaluate BC + !val(1,i,1) = qtrue(x2,0.d0,t2) + val(1,i,1) = qtrue(xcell,ycell1,time) + endif + enddo + go to 399 +c + 310 continue +c # zero-order extrapolation: + do 315 j=1,nyb + do 315 i=1,nrow + do 315 m=1,meqn + val(m,i,j) = val(m,i,nyb+1) + 315 continue + go to 399 + + 320 continue +c # periodic: handled elsewhere in amr + go to 399 + + 330 continue +c # solid wall (assumes 3'rd component is velocity or momentum in y): + do 335 j=1,nyb + do 335 i=1,nrow + do 335 m=1,meqn + val(m,i,j) = val(m,i,2*nyb+1-j) + 335 continue +c # negate the normal velocity: + do 336 j=1,nyb + do 336 i=1,nrow + val(3,i,j) = -val(3,i,j) + 336 continue + go to 399 + + 399 continue +c +c------------------------------------------------------- +c # top boundary: +c------------------------------------------------------- + if (ytop .le. yupper+hymarg) then +c # not a physical boundary -- no cells at this edge lies +c # outside the physical bndry. +c # values are set elsewhere in amr code. + go to 499 + endif +c +c # number of grid cells lying outside physical domain: + nyt = (ytop - yupper + hymarg)/hy + jbeg = max0(ncol-nyt+1, 1) +c + go to (400,410,420,430) mthbc(4)+1 +c + 400 continue +c # user-specified boundary conditions go here in place of error output + write(6,*) + & '*** ERROR *** mthbc(4)=0 and no BCs specified in bc2amr' + stop + go to 499 + + 410 continue +c # zero-order extrapolation: + do 415 j=jbeg,ncol + do 415 i=1,nrow + do 415 m=1,meqn + !val(m,i,j) = val(m,i,jbeg-1) + val(m,i,j) = 2.d0*val(m,i,jbeg-1) - val(m,i,jbeg-2) + 415 continue + go to 499 + + 420 continue +c # periodic: handled elsewhere in amr + go to 499 + + 430 continue +c # solid wall (assumes 3'rd component is velocity or momentum in y): + do 435 j=jbeg,ncol + do 435 i=1,nrow + do 435 m=1,meqn + val(m,i,j) = val(m,i,2*jbeg-1-j) + 435 continue +c # negate the normal velocity: + do 436 j=jbeg,ncol + do 436 i=1,nrow + val(3,i,j) = -val(3,i,j) + 436 continue + go to 499 + + 499 continue + + return + end diff --git a/dev/advection_2d_linearTest/flag2refine2.f90 b/dev/advection_2d_linearTest/flag2refine2.f90 new file mode 100644 index 000000000..5f14caf37 --- /dev/null +++ b/dev/advection_2d_linearTest/flag2refine2.f90 @@ -0,0 +1,107 @@ +! ::::::::::::::::::::: flag2refine :::::::::::::::::::::::::::::::::: +! +! User routine to control flagging of points for refinement. +! +! Default version computes spatial difference dq in each direction and +! for each component of q and flags any point where this is greater than +! the tolerance tolsp. This is consistent with what the routine errsp did in +! earlier versions of amrclaw (4.2 and before). +! +! This routine can be copied to an application directory and modified to +! implement some other desired refinement criterion. +! +! Points may also be flagged for refining based on a Richardson estimate +! of the error, obtained by comparing solutions on the current grid and a +! coarsened grid. Points are flagged if the estimated error is larger than +! the parameter tol in amr2ez.data, provided flag_richardson is .true., +! otherwise the coarsening and Richardson estimation is not performed! +! Points are flagged via Richardson in a separate routine. +! +! Once points are flagged via this routine and/or Richardson, the subroutine +! flagregions is applied to check each point against the min_level and +! max_level of refinement specified in any "region" set by the user. +! So flags set here might be over-ruled by region constraints. +! +! q = grid values including ghost cells (bndry vals at specified +! time have already been set, so can use ghost cell values too) +! +! aux = aux array on this grid patch +! +! amrflags = array to be flagged with either the value +! DONTFLAG (no refinement needed) or +! DOFLAG (refinement desired) +! +! tolsp = tolerance specified by user in input file amr2ez.data, used in default +! version of this routine as a tolerance for spatial differences. +! +! :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: +subroutine flag2refine2(mx,my,mbc,mbuff,meqn,maux,xlower,ylower,dx,dy,t,level, & + tolsp,q,aux,amrflags,DONTFLAG,DOFLAG) + + use regions_module + use amr_module, only: lfine + + implicit none + + ! Subroutine arguments + integer, intent(in) :: mx,my,mbc,meqn,maux,level,mbuff + real(kind=8), intent(in) :: xlower,ylower,dx,dy,t,tolsp + + real(kind=8), intent(in) :: q(meqn,1-mbc:mx+mbc,1-mbc:my+mbc) + real(kind=8), intent(in) :: aux(maux,1-mbc:mx+mbc,1-mbc:my+mbc) + + ! Flagging + real(kind=8),intent(inout) :: amrflags(1-mbuff:mx+mbuff,1-mbuff:my+mbuff) + real(kind=8), intent(in) :: DONTFLAG + real(kind=8), intent(in) :: DOFLAG + + logical :: allowflag + external allowflag + + ! Locals + integer :: i,j,m + real(kind=8) :: x_c,y_c,x_low,y_low,x_hi,y_hi + real(kind=8) :: dqi(meqn), dqj(meqn), dq(meqn) + + ! Initialize flags + amrflags = DONTFLAG + + ! Loop over interior points on this grid + ! (i,j) grid cell is [x_low,x_hi] x [y_low,y_hi], cell center at (x_c,y_c) + ! This information is not needed for the default flagging based on + ! undivided differences, but might be needed in a user's version. + ! Note that if you want to refine only in certain space-time regions, + ! it may be easiest to use the "regions" feature. The flags set here or + ! in the Richardson error estimator are potentially modified by the + ! min_level and max_level specified in any regions. + + y_loop: do j=1,my + y_low = ylower + (j - 1) * dy + y_c = ylower + (j - 0.5d0) * dy + y_hi = ylower + j * dy + + x_loop: do i = 1,mx + x_low = xlower + (i - 1) * dx + x_c = xlower + (i - 0.5d0) * dx + x_hi = xlower + i * dx + + ! ----------------------------------------------------------------- + dq = 0.d0 + dqi = abs(q(:,i+1,j) - q(:,i-1,j)) + dqj = abs(q(:,i,j+1) - q(:,i,j-1)) + dq = max(dq,dqi,dqj) + + ! default checks all components of undivided difference: + do m=1,meqn + if (dq(m) > tolsp) then + amrflags(i,j) = DOFLAG + cycle x_loop + endif + enddo + + enddo x_loop + enddo y_loop + amrflags(mx/2,my/2) = DOFLAG + !if (level .eq. lfine .and. t .gt. 1.) amrflags(mx/2,my/2) = DOFLAG + +end subroutine flag2refine2 diff --git a/dev/advection_2d_linearTest/qinit.f b/dev/advection_2d_linearTest/qinit.f new file mode 100644 index 000000000..026a53e8d --- /dev/null +++ b/dev/advection_2d_linearTest/qinit.f @@ -0,0 +1,25 @@ +c ===================================================== + subroutine qinit(meqn,mbc,mx,my,xlower,ylower, + & dx,dy,q,maux,aux) +c ===================================================== +c +c # Set initial conditions for q. +c # Sample scalar equation with data that is piecewise constant with +c # q = 1.0 if 0.1 < x < 0.6 and 0.1 < y < 0.6 +c # 0.1 otherwise +c + implicit double precision (a-h,o-z) + dimension q(meqn, 1-mbc:mx+mbc, 1-mbc:my+mbc) + dimension aux(maux, 1-mbc:mx+mbc, 1-mbc:my+mbc) + +c # set concentration profile +c --------------------------- + do 20 j=1,my + yj = ylower + (j-0.5d0)*dy + do 20 i=1,mx + xi = xlower + (i-0.5d0)*dx + q(1,i,j) = qtrue(xi,yj,0.d0) + 20 continue + + return + end diff --git a/dev/advection_2d_linearTest/qtrue.f b/dev/advection_2d_linearTest/qtrue.f new file mode 100644 index 000000000..3c6331e0e --- /dev/null +++ b/dev/advection_2d_linearTest/qtrue.f @@ -0,0 +1,30 @@ + + +c ================================== + double precision function qtrue(x,y,t) +c ================================== + implicit double precision (a-h,o-z) + common /cparam/ ubar,vbar + + x0 = x - ubar*t + y0 = y - vbar*t + +c # evaluate desired initial data at (x0,y0): + +!-- r = dsqrt((x0+0.2d0)**2 + (y0-0.4d0)**2) +!-- if (r <= 0.3d0) then +!-- qtrue = 1.d0 +!-- else +!-- qtrue = 0.d0 +!-- endif +!-- +!-- r = dsqrt((x0-0.3d0)**2 + (y0-0.1d0)**2) +!-- qtrue = qtrue - dexp(-15.d0*r**2) + + +c linear test function to see if ghost cell interpolation exact + !qtrue = 2.d0*x0 + 1.0d0*y0 + qtrue = 1.d0*x0 + 1.0d0*y0 + + return + end diff --git a/dev/advection_2d_linearTest/setplot.py b/dev/advection_2d_linearTest/setplot.py new file mode 100644 index 000000000..a8adf864f --- /dev/null +++ b/dev/advection_2d_linearTest/setplot.py @@ -0,0 +1,169 @@ + +""" +Set up the plot figures, axes, and items to be done for each frame. + +This module is imported by the plotting routines and then the +function setplot is called to set the plot parameters. +""" + +from clawpack.clawutil.data import ClawData +probdata = ClawData() +probdata.read('setprob.data', force=True) +ubar = probdata.u +vbar = probdata.v + +#-------------------------- +def setplot(plotdata): +#-------------------------- + + """ + Specify what is to be plotted at each frame. + Input: plotdata, an instance of visclaw.data.ClawPlotData. + Output: a modified version of plotdata. + + """ + + + from clawpack.visclaw import colormaps + + plotdata.clearfigures() # clear any old figures,axes,items data + + + # Figure for pcolor plot + plotfigure = plotdata.new_plotfigure(name='pcolor', figno=0) + + # Set up for axes in this figure: + plotaxes = plotfigure.new_plotaxes() + plotaxes.xlimits = [0,50] + plotaxes.ylimits = [0,50] + plotaxes.title = 'Solution' + plotaxes.scaled = True + + # Set up for item on these axes: + plotitem = plotaxes.new_plotitem(plot_type='2d_pcolor') + plotitem.plot_var = 0 + plotitem.pcolor_cmap = colormaps.red_yellow_blue + plotitem.pcolor_cmin = -150. + plotitem.pcolor_cmax = 150. + plotitem.add_colorbar = True + + plotitem.amr_celledges_show = [0,0] + plotitem.amr_patchedges_show = [0,1] + + + # Figure for pcolor plot + plotfigure = plotdata.new_plotfigure(name='error pcolor', figno=10) + + def error(current_data): + x = current_data.x + y = current_data.y + t = current_data.t + q = current_data.q + qtrue = (x - ubar*t) + (y - vbar*t) + e = q[0,:,:] - qtrue + #import pdb; pdb.set_trace() + return e + + # Set up for axes in this figure: + plotaxes = plotfigure.new_plotaxes() + plotaxes.xlimits = [0,50] + plotaxes.ylimits = [0,50] + plotaxes.title = 'Error' + plotaxes.scaled = True + + # Set up for item on these axes: + plotitem = plotaxes.new_plotitem(plot_type='2d_pcolor') + plotitem.plot_var = error + plotitem.pcolor_cmap = colormaps.red_yellow_blue + plotitem.pcolor_cmin = -0.000001 + plotitem.pcolor_cmax = 0.000001 + plotitem.add_colorbar = True + + plotitem.amr_celledges_show = [0,0] + plotitem.amr_patchedges_show = [0,1] + + + # Figure for contour plot + plotfigure = plotdata.new_plotfigure(name='contour', figno=1) + + # Set up for axes in this figure: + plotaxes = plotfigure.new_plotaxes() + plotaxes.xlimits = [0,50] + plotaxes.ylimits = [0,50] + plotaxes.title = 'Solution' + plotaxes.scaled = True + plotaxes.afteraxes = addgauges + + # Set up for item on these axes: + plotitem = plotaxes.new_plotitem(plot_type='2d_contour') + plotitem.plot_var = 0 + plotitem.contour_nlevels = 40 + plotitem.contour_min = -150.99 + plotitem.contour_max = 150.99 + plotitem.amr_contour_colors = ['r','g','b'] # color on each level + plotitem.amr_patch_bgcolor = ['#ffeeee', '#eeeeff', '#eeffee'] + plotitem.celledges_show = 0 + plotitem.patchedges_show = 1 + + + # Figure for grid cells + plotfigure = plotdata.new_plotfigure(name='cells', figno=2) + + # Set up for axes in this figure: + plotaxes = plotfigure.new_plotaxes() + plotaxes.xlimits = [0,50] + plotaxes.ylimits = [0,50] + plotaxes.title = 'Grid patches' + plotaxes.scaled = True + + # Set up for item on these axes: + plotitem = plotaxes.new_plotitem(plot_type='2d_patch') + plotitem.amr_patch_bgcolor = ['#ffeeee', '#eeeeff', '#eeffee'] + plotitem.amr_celledges_show = [0,0,0] + plotitem.amr_patchedges_show = [1,1] + + #----------------------------------------- + # Figures for gauges + #----------------------------------------- + plotfigure = plotdata.new_plotfigure(name='q', figno=300, \ + type='each_gauge') + plotfigure.clf_each_gauge = True + + # Set up for axes in this figure: + plotaxes = plotfigure.new_plotaxes() + plotaxes.xlimits = 'auto' + plotaxes.ylimits = 'auto' + plotaxes.title = 'q' + + # Plot q as blue curve: + plotitem = plotaxes.new_plotitem(plot_type='1d_plot') + plotitem.plot_var = 0 + plotitem.plotstyle = 'b-' + + + # Parameters used only when creating html and/or latex hardcopy + # e.g., via visclaw.frametools.printframes: + + plotdata.printfigs = True # print figures + plotdata.print_format = 'png' # file format + plotdata.print_framenos = 'all' # list of frames to print + plotdata.print_fignos = 'all' # list of figures to print + plotdata.html = True # create html files of plots? + plotdata.html_homelink = '../README.html' # pointer for top of index + plotdata.html_movie = 'JSAnimation' # new style, or "4.x" for old style + plotdata.latex = True # create latex file of plots? + plotdata.latex_figsperline = 2 # layout of plots + plotdata.latex_framesperline = 1 # layout of plots + plotdata.latex_makepdf = False # also run pdflatex? + + return plotdata + + + +# To plot gauge locations on pcolor or contour plot, use this as +# an afteraxis function: + +def addgauges(current_data): + from clawpack.visclaw import gaugetools + gaugetools.plot_gauge_locations(current_data.plotdata, \ + gaugenos='all', format_string='ko', add_labels=True) diff --git a/dev/advection_2d_linearTest/setprob.f b/dev/advection_2d_linearTest/setprob.f new file mode 100644 index 000000000..392f69c29 --- /dev/null +++ b/dev/advection_2d_linearTest/setprob.f @@ -0,0 +1,22 @@ + subroutine setprob + implicit double precision (a-h,o-z) + character*25 fname + common /cparam/ ubar,vbar +c +c # Set the velocity for scalar advection +c # These values are passed to the Riemann solvers rpn2.f and rpt2.f +c # in a common block +c + +c + iunit = 7 + fname = 'setprob.data' +c # open the unit with new routine from Clawpack 4.4 to skip over +c # comment lines starting with #: + call opendatafile(iunit, fname) + + read(7,*) ubar + read(7,*) vbar + + return + end diff --git a/dev/advection_2d_linearTest/setrun.py b/dev/advection_2d_linearTest/setrun.py new file mode 100644 index 000000000..22dcc6f3d --- /dev/null +++ b/dev/advection_2d_linearTest/setrun.py @@ -0,0 +1,348 @@ +""" +Module to set up run time parameters for Clawpack. + +The values set in the function setrun are then written out to data files +that will be read in by the Fortran code. + +""" + +import os +import numpy as np + +#------------------------------ +def setrun(claw_pkg='amrclaw'): +#------------------------------ + + """ + Define the parameters used for running Clawpack. + + INPUT: + claw_pkg expected to be "amrclaw" for this setrun. + + OUTPUT: + rundata - object of class ClawRunData + + """ + + from clawpack.clawutil import data + + + assert claw_pkg.lower() == 'amrclaw', "Expected claw_pkg = 'amrclaw'" + + num_dim = 2 + rundata = data.ClawRunData(claw_pkg, num_dim) + + #------------------------------------------------------------------ + # Problem-specific parameters to be written to setprob.data: + #------------------------------------------------------------------ + + probdata = rundata.new_UserData(name='probdata',fname='setprob.data') + probdata.add_param('u', .5, 'ubar advection velocity') + probdata.add_param('v', 1. , 'vbar advection velocity') + + #------------------------------------------------------------------ + # Standard Clawpack parameters to be written to claw.data: + # (or to amr2ez.data for AMR) + #------------------------------------------------------------------ + + clawdata = rundata.clawdata # initialized when rundata instantiated + + + # Set single grid parameters first. + # See below for AMR parameters. + + + # --------------- + # Spatial domain: + # --------------- + + # Number of space dimensions: + clawdata.num_dim = num_dim + + # Lower and upper edge of computational domain: + clawdata.lower[0] = 0. # xlower + clawdata.upper[0] = 50. # xupper + clawdata.lower[1] = 0. # ylower + clawdata.upper[1] = 50. # yupper + + # Number of grid cells: + clawdata.num_cells[0] = 50 # mx + clawdata.num_cells[1] = 50 # my + + + # --------------- + # Size of system: + # --------------- + + # Number of equations in the system: + clawdata.num_eqn = 1 + + # Number of auxiliary variables in the aux array (initialized in setaux) + clawdata.num_aux = 0 + + # Index of aux array corresponding to capacity function, if there is one: + clawdata.capa_index = 0 + + + # ------------- + # Initial time: + # ------------- + + clawdata.t0 = 0.0 + + + # Restart from checkpoint file of a previous run? + # Note: If restarting, you must also change the Makefile to set: + # RESTART = True + # If restarting, t0 above should be from original run, and the + # restart_file 'fort.chkNNNNN' specified below should be in + # the OUTDIR indicated in Makefile. + + clawdata.restart = False # True to restart from prior results + clawdata.restart_file = 'fort.chk00006' # File to use for restart data + + + # ------------- + # Output times: + #-------------- + + # Specify at what times the results should be written to fort.q files. + # Note that the time integration stops after the final output time. + + clawdata.output_style = 3 + + if clawdata.output_style==1: + # Output ntimes frames at equally spaced times up to tfinal: + # Can specify num_output_times = 0 for no output + clawdata.num_output_times = 10 + clawdata.tfinal = 2.0 + clawdata.output_t0 = True # output at initial (or restart) time? + + elif clawdata.output_style == 2: + # Specify a list or numpy array of output times: + # Include t0 if you want output at the initial time. + clawdata.output_times = [0., 0.1] + + elif clawdata.output_style == 3: + # Output every step_interval timesteps over total_steps timesteps: + clawdata.output_step_interval = 1 + clawdata.total_steps = 100 + clawdata.output_t0 = True # output at initial (or restart) time? + + + clawdata.output_format = 'ascii' # 'ascii', 'binary', 'netcdf' + + clawdata.output_q_components = 'all' # could be list such as [True,True] + clawdata.output_aux_components = 'none' # could be list + clawdata.output_aux_onlyonce = True # output aux arrays only at t0 + + + # --------------------------------------------------- + # Verbosity of messages to screen during integration: + # --------------------------------------------------- + + # The current t, dt, and cfl will be printed every time step + # at AMR levels <= verbosity. Set verbosity = 0 for no printing. + # (E.g. verbosity == 2 means print only on levels 1 and 2.) + clawdata.verbosity = 0 + + + + # -------------- + # Time stepping: + # -------------- + + # if dt_variable==True: variable time steps used based on cfl_desired, + # if dt_variable==False: fixed time steps dt = dt_initial always used. + clawdata.dt_variable = True + + # Initial time step for variable dt. + # (If dt_variable==0 then dt=dt_initial for all steps) + clawdata.dt_initial = 0.016 + + # Max time step to be allowed if variable dt used: + clawdata.dt_max = 1e+99 + + # Desired Courant number if variable dt used + clawdata.cfl_desired = 0.9 + # max Courant number to allow without retaking step with a smaller dt: + clawdata.cfl_max = 1.0 + + # Maximum number of time steps to allow between output times: + clawdata.steps_max = 100000 + + + # ------------------ + # Method to be used: + # ------------------ + + # Order of accuracy: 1 => Godunov, 2 => Lax-Wendroff plus limiters + clawdata.order = 2 + + # Use dimensional splitting? + clawdata.dimensional_split = 'unsplit' + + # For unsplit method, transverse_waves can be + # 0 or 'none' ==> donor cell (only normal solver used) + # 1 or 'increment' ==> corner transport of waves + # 2 or 'all' ==> corner transport of 2nd order corrections too + clawdata.transverse_waves = 'all' + + + # Number of waves in the Riemann solution: + clawdata.num_waves = 1 + + # List of limiters to use for each wave family: + # Required: len(limiter) == num_waves + # Some options: + # 0 or 'none' ==> no limiter (Lax-Wendroff) + # 1 or 'minmod' ==> minmod + # 2 or 'superbee' ==> superbee + # 3 or 'vanleer' ==> van Leer + # 4 or 'mc' ==> MC limiter + clawdata.limiter = ['none'] + + clawdata.use_fwaves = False # True ==> use f-wave version of algorithms + + # Source terms splitting: + # src_split == 0 or 'none' ==> no source term (src routine never called) + # src_split == 1 or 'godunov' ==> Godunov (1st order) splitting used, + # src_split == 2 or 'strang' ==> Strang (2nd order) splitting used, not recommended. + clawdata.source_split = 'none' + + + # -------------------- + # Boundary conditions: + # -------------------- + + # Number of ghost cells (usually 2) + clawdata.num_ghost = 2 + + # Choice of BCs at xlower and xupper: + # 0 or 'user' => user specified (must modify bcNamr.f to use this option) + # 1 or 'extrap' => extrapolation (non-reflecting outflow) + # 2 or 'periodic' => periodic (must specify this at both boundaries) + # 3 or 'wall' => solid wall for systems where q(2) is normal velocity + + clawdata.bc_lower[0] = 'user' # at xlower + clawdata.bc_upper[0] = 'extrap' # at xupper + + clawdata.bc_lower[1] = 'user' # at ylower + clawdata.bc_upper[1] = 'extrap' # at yupper + + + # --------------- + # Gauges: + # --------------- + rundata.gaugedata.gauges = [] + # for gauges append lines of the form [gaugeno, x, y, t1, t2] + rundata.gaugedata.gauges.append([1, 0.6, 0.6, 0., 10.]) + rundata.gaugedata.gauges.append([2, 0.7, 0.3, 0., 10.]) + + + # -------------- + # Checkpointing: + # -------------- + + # Specify when checkpoint files should be created that can be + # used to restart a computation. + + clawdata.checkpt_style = 0 + + if clawdata.checkpt_style == 0: + # Do not checkpoint at all + pass + + elif clawdata.checkpt_style == 1: + # Checkpoint only at tfinal. + pass + + elif clawdata.checkpt_style == 2: + # Specify a list of checkpoint times. + clawdata.checkpt_times = [0.1,0.15] + + elif clawdata.checkpt_style == 3: + # Checkpoint every checkpt_interval timesteps (on Level 1) + # and at the final time. + clawdata.checkpt_interval = 5 + + # --------------- + # AMR parameters: + # --------------- + amrdata = rundata.amrdata + + + # max number of refinement levels: + amrdata.amr_levels_max = 3 + + # List of refinement ratios at each level (length at least amr_level_max-1) + amrdata.refinement_ratios_x = [6,4,2,8] + amrdata.refinement_ratios_y = [6,4,2,8] + amrdata.refinement_ratios_t = [6,4,2,8] + + + # Specify type of each aux variable in amrdata.auxtype. + # This must be a list of length num_aux, each element of which is one of: + # 'center', 'capacity', 'xleft', or 'yleft' (see documentation). + amrdata.aux_type = [] + + + # Flag for refinement based on Richardson error estimater: + amrdata.flag_richardson = False # use Richardson? + amrdata.flag_richardson_tol = 0.1 # Richardson tolerance + + # Flag for refinement using routine flag2refine: + amrdata.flag2refine = True # use this? + amrdata.flag2refine_tol = 100.55 # tolerance used in this routine + # User can modify flag2refine to change the criterion for flagging. + # Default: check max-norm of difference between q in a cell and + # each of its neighbors. + + # steps to take on each level L between regriddings of level L+1: + amrdata.regrid_interval = 40 + + # width of buffer zone around flagged points: + # (typically the same as regrid_interval so waves don't escape): + amrdata.regrid_buffer_width = 7 + + # clustering alg. cutoff for (# flagged pts) / (total # of cells refined) + # (closer to 1.0 => more small grids may be needed to cover flagged cells) + amrdata.clustering_cutoff = 0.7 + + # print info about each regridding up to this level: + amrdata.verbosity_regrid = 3 + + + # --------------- + # Regions: + # --------------- + regions = rundata.regiondata.regions + # to specify regions of refinement append lines of the form + # [minlevel,maxlevel,t1,t2,x1,x2,y1,y2] + + + # ----- For developers ----- + # Toggle debugging print statements: + amrdata.dprint = False # print domain flags + amrdata.eprint = False # print err est flags + amrdata.edebug = False # even more err est flags + amrdata.gprint = False # grid bisection/clustering + amrdata.nprint = False # proper nesting output + amrdata.pprint = False # proj. of tagged points + amrdata.rprint = True # print regridding summary + amrdata.sprint = False # space/memory output + amrdata.tprint = True # time step reporting each level + amrdata.uprint = False # update/upbnd reporting + + return rundata + + # end of function setrun + # ---------------------- + + +if __name__ == '__main__': + # Set up run-time parameters and write all data files. + import sys + rundata = setrun(*sys.argv[1:]) + rundata.write() + From c6a20b3d1ca82aef2a5bd3ccb2d13709de5086ed Mon Sep 17 00:00:00 2001 From: Marsha Berger Date: Sat, 22 Aug 2015 15:38:31 -0700 Subject: [PATCH 04/53] new test dir - should get zero error for 3d linear solution --- dev/advection_3d_linearTest/Makefile | 142 ++++++ dev/advection_3d_linearTest/bc3amr.f | 453 +++++++++++++++++++ dev/advection_3d_linearTest/flag2refine2.f90 | 113 +++++ dev/advection_3d_linearTest/qinit.f | 26 ++ dev/advection_3d_linearTest/qtrue.f | 20 + dev/advection_3d_linearTest/setaux.f | 36 ++ dev/advection_3d_linearTest/setplot.py | 174 +++++++ dev/advection_3d_linearTest/setprob.f | 23 + dev/advection_3d_linearTest/setrun.py | 368 +++++++++++++++ dev/advection_3d_linearTest/valout.f | 260 +++++++++++ 10 files changed, 1615 insertions(+) create mode 100644 dev/advection_3d_linearTest/Makefile create mode 100644 dev/advection_3d_linearTest/bc3amr.f create mode 100644 dev/advection_3d_linearTest/flag2refine2.f90 create mode 100644 dev/advection_3d_linearTest/qinit.f create mode 100644 dev/advection_3d_linearTest/qtrue.f create mode 100644 dev/advection_3d_linearTest/setaux.f create mode 100644 dev/advection_3d_linearTest/setplot.py create mode 100644 dev/advection_3d_linearTest/setprob.f create mode 100644 dev/advection_3d_linearTest/setrun.py create mode 100644 dev/advection_3d_linearTest/valout.f diff --git a/dev/advection_3d_linearTest/Makefile b/dev/advection_3d_linearTest/Makefile new file mode 100644 index 000000000..0b3015732 --- /dev/null +++ b/dev/advection_3d_linearTest/Makefile @@ -0,0 +1,142 @@ + +# Makefile for Clawpack code in this directory. +# This version only sets the local files and frequently changed +# options, and then includes the standard makefile pointed to by CLAWMAKE. +CLAWMAKE = $(CLAW)/clawutil/src/Makefile.common + +# See the above file for details and a list of make options, or type +# make .help +# at the unix prompt. + + +# Adjust these variables if desired: +# ---------------------------------- + +CLAW_PKG = amrclaw # Clawpack package to use +EXE = xamr # Executable to create +SETRUN_FILE = setrun.py # File containing function to make data +OUTDIR = _output # Directory for output +SETPLOT_FILE = setplot.py # File containing function to set plots +PLOTDIR = _plots # Directory for plots + +OVERWRITE ?= True # False ==> make a copy of OUTDIR first +RESTART ?= False # Should = clawdata.restart in setrun + +# Environment variable FC should be set to fortran compiler, e.g. gfortran + +# Compiler flags can be specified here or set as an environment variable +FFLAGS ?= + +# --------------------------------- +# List of sources for this program: +# --------------------------------- + +AMRLIB = $(CLAW)/amrclaw/src/3d + +MODULES = \ + $(AMRLIB)/amr_module.f90 \ + $(AMRLIB)/gauges_module.f90 + +SOURCES = \ + qinit.f \ + qtrue.f \ + setprob.f \ + bc3amr.f \ + setaux.f \ + valout.f \ + flag2refine2.f90 \ + $(CLAW)/riemann/src/rpn3_vc_advection.f90 \ + $(CLAW)/riemann/src/rpt3_vc_advection.f90 \ + $(CLAW)/riemann/src/rptt3_vc_advection.f90 \ + $(AMRLIB)/amr3.f90 \ + $(AMRLIB)/opendatafile.f \ + $(AMRLIB)/b4step3.f90 \ + $(AMRLIB)/qad.f \ + $(AMRLIB)/src3.f90 \ + $(AMRLIB)/src1d.f90 \ + $(AMRLIB)/advanc.f \ + $(AMRLIB)/bound.f \ + $(AMRLIB)/stepgrid.f \ + $(AMRLIB)/stepgrid_dimSplit.f \ + $(AMRLIB)/auxcoarsen.f \ + $(AMRLIB)/fixcapaq.f \ + $(AMRLIB)/estdt.f \ + $(AMRLIB)/igetsp.f \ + $(AMRLIB)/reclam.f \ + $(AMRLIB)/birect.f \ + $(AMRLIB)/cleanup.f \ + $(AMRLIB)/colate.f \ + $(AMRLIB)/errest.f \ + $(AMRLIB)/allowflag.f \ + $(AMRLIB)/bufnst.f \ + $(AMRLIB)/spest.f \ + $(AMRLIB)/errf1.f \ + $(AMRLIB)/gfixup.f \ + $(AMRLIB)/filval.f \ + $(AMRLIB)/filpatch.f \ + $(AMRLIB)/prefilp.f \ + $(AMRLIB)/flglvl.f \ + $(AMRLIB)/fluxad.f \ + $(AMRLIB)/fluxsv.f \ + $(AMRLIB)/ginit.f \ + $(AMRLIB)/grdfit.f \ + $(AMRLIB)/intfil.f \ + $(AMRLIB)/moment.f \ + $(AMRLIB)/nestck.f \ + $(AMRLIB)/prepc.f \ + $(AMRLIB)/prepf.f \ + $(AMRLIB)/projec.f \ + $(AMRLIB)/signs.f \ + $(AMRLIB)/findcut.f \ + $(AMRLIB)/smartbis.f \ + $(AMRLIB)/putnod.f \ + $(AMRLIB)/putsp.f \ + $(AMRLIB)/regrid.f \ + $(AMRLIB)/setgrd.f \ + $(AMRLIB)/setuse.f \ + $(AMRLIB)/stst1.f \ + $(AMRLIB)/tick.f \ + $(AMRLIB)/trimbd.f \ + $(AMRLIB)/update.f \ + $(AMRLIB)/nodget.f \ + $(AMRLIB)/upbnd.f \ + $(AMRLIB)/basic.f \ + $(AMRLIB)/outval.f \ + $(AMRLIB)/copysol.f \ + $(AMRLIB)/outvar.f \ + $(AMRLIB)/outmsh.f \ + $(AMRLIB)/outtre.f \ + $(AMRLIB)/domain.f \ + $(AMRLIB)/setflags.f \ + $(AMRLIB)/shiftset.f \ + $(AMRLIB)/conck.f \ + $(AMRLIB)/domshrink.f \ + $(AMRLIB)/domprep.f \ + $(AMRLIB)/domup.f \ + $(AMRLIB)/domcopy.f \ + $(AMRLIB)/coarsen.f \ + $(AMRLIB)/intcopy.f \ + $(AMRLIB)/preintcopy.f \ + $(AMRLIB)/icall.f \ + $(AMRLIB)/preicall.f \ + $(AMRLIB)/step3.f \ + $(AMRLIB)/step3x.f \ + $(AMRLIB)/step3y.f \ + $(AMRLIB)/step3z.f \ + $(AMRLIB)/flux3.f \ + $(AMRLIB)/flux3_dimSplit.f \ + $(AMRLIB)/limiter.f \ + $(AMRLIB)/philim.f \ + $(AMRLIB)/cstore.f \ + $(AMRLIB)/saveqc.f \ + $(AMRLIB)/check.f \ + $(AMRLIB)/restrt.f \ + $(AMRLIB)/quick_sort1.f \ + $(AMRLIB)/init_alloc.f90 \ + $(AMRLIB)/resize_alloc.f90 \ + $(AMRLIB)/restrt_alloc.f90 + +#------------------------------------------------------------------- +# Include Makefile containing standard definitions and make options: +include $(CLAWMAKE) + diff --git a/dev/advection_3d_linearTest/bc3amr.f b/dev/advection_3d_linearTest/bc3amr.f new file mode 100644 index 000000000..9747f4453 --- /dev/null +++ b/dev/advection_3d_linearTest/bc3amr.f @@ -0,0 +1,453 @@ +c +c ------------------------------------------------------------------ +c + subroutine bc3amr(val,aux,nrow,ncol,nfil,meqn,naux, + 1 hx, hy, hz, level, time, + 2 xleft, xright, yfront, yrear, + 3 zbot, ztop, + 4 xlower,ylower,zlower, + 5 xupper,yupper,zupper, + 6 xperiodic, yperiodic,zperiodic) + + +c +c +c :::::::::: BC3AMR ::::::::::::::::::::::::::::::::::::::::::::::; +c +c Take a grid patch with mesh widths hx,hy,hz, of dimensions nrow by +c ncol by nfil, and set the values of any piece of +c of the patch which extends outside the physical domain +c using the boundary conditions. +c ------------------------------------------------ +c # Standard boundary condition choices for amr3ez in clawpack +c +c # At each boundary k = 1 (left), 2 (right), 3 (front), 4 (rear), +c 5 (bottom) 6 (top) +c # mthbc(k) = 0 for user-supplied BC's (must be inserted!) +c # = 1 for zero-order extrapolation +c # = 2 for periodic boundary coniditions +c # = 3 for solid walls, assuming this can be implemented +c # by reflecting the data about the boundary and then +c # negating the 2'nd (for k=1,2) or 3'rd (for k=3,4) +c # or 4'th (for k = 5,6) component of q. +c ------------------------------------------------ +c +c The corners of the grid patch are at +c (xleft,yfront,zbot) -- lower front left corner +c (xright,yrear,ztop) -- upper rear right corner +c +c The physical domain itself is a rectangular parallelopiped bounded by +c (xlower,ylower,zlower) -- lower front left corner +c (xupper,yupper,zupper) -- upper rear right corner +c +c the picture is the following: +c +c __________________________(xupper,yupper,zupper) +c / /| +c / / | +c / / | +c /_________________________/ | +c | | | +c | | | +c ___|_____(xright,yrear,ztop) | | +c /___|____/| | | +c | | || | | +c | | || | | +c | | || | | +c |___|____|/ | | +c (xleft,yfront,zbot) | | / +c | | / +c |_________________________|/ +c (xlower,ylower,zlower) +c +c Any cells that lie outside the physical domain are ghost cells whose +c values should be set in this routine. This is tested for by comparing +c xleft with xlower to see if values need to be set at the left, as in +c the figure above, and similarly at the other boundaries. +c +c Patches are guaranteed to have at least 1 row of cells filled +c with interior values so it is possible to extrapolate. +c Fix trimbd if you want more than 1 row pre-set. +c +c Make sure the order the boundaries are specified is correct +c so that diagonal corner cells are also properly taken care of. +c +c Periodic boundaries are set before calling this routine, so if the +c domain is periodic in one direction only you +c can safely extrapolate in the other direction. +c +c Don't overwrite ghost cells in periodic directions! +c +c ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::; + + use amr_module, only: mthbc + implicit none + + real*8 val(meqn,nrow,ncol,nfil), aux(naux,nrow,ncol,nfil) + real*8 ubar, vbar, wbar + logical xperiodic, yperiodic, zperiodic + integer nrow,ncol,nfil,meqn,naux,level + real*8 hx,hy,hz,time,qtrue + real*8 hxmarg,hymarg,hzmarg + real*8 xleft,xright,yfront,yrear,zbot,ztop + real*8 xcell,ycell,zcell + real*8 xlower,ylower,zlower + real*8 xupper,yupper,zupper + integer i,j,k,m,nxl,nxr,ibeg,nyf,nyr,jbeg,nzb,nzt,kbeg + + common /cparam/ ubar,vbar,wbar + + hxmarg = hx*.01d0 + hymarg = hy*.01d0 + hzmarg = hz*.01d0 + + if (xperiodic .and. yperiodic .and. zperiodic) go to 699 +c +c +c------------------------------------------------------- +c # xlower boundary: +c------------------------------------------------------- + if (xleft .ge. xlower-hxmarg) then +c # not a physical boundary -- ghost cells lie within another +c # grid and values are set elsewhere in amr code. + go to 199 + endif +c +c # number of ghost cells lying outside physical domain: + nxl = (xlower+hxmarg-xleft)/hx +c + if (nxl > 2) then + write(6,*) '*** unexpected value nxl = ',nxl + stop + endif +c + go to (100,110,120,130) mthbc(1)+1 +c + 100 continue + do k = 1,nfil + zcell = zbot + (k-0.5d0)*hz + do j = 1,ncol + ycell = yfront + (j-0.5d0)*hy + do i=1,nxl + xcell = xleft + (i-0.5d0)*hx + val(1,i,j,k) = qtrue(xcell,ycell,zcell,time) + enddo + enddo + enddo + go to 199 +c + 110 continue +c # zero-order extrapolation: + do 115 k = 1,nfil + do 115 j = 1,ncol + do 115 i=1,nxl + do 115 m=1,meqn + val(m,i,j,k) = val(m,nxl+1,j,k) + 115 continue + go to 199 + + 120 continue +c # periodic: handled elsewhere in amr + go to 199 + + 130 continue +c # solid wall (assumes 2'nd component is velocity or momentum in x): + do 135 k = 1,nfil + do 135 j = 1,ncol + do 135 i=1,nxl + do 135 m=1,meqn + val(m,i,j,k) = val(m,2*nxl+1-i,j,k) + 135 continue +c # negate the normal velocity: + do 136 i=1,nxl + do 136 j = 1,ncol + do 136 k = 1,nfil + val(2,i,j,k) = -val(2,i,j,k) + 136 continue + go to 199 + + 199 continue +c +c------------------------------------------------------- +c # xupper boundary: +c------------------------------------------------------- + if (xright .le. xupper+hxmarg) then +c # not a physical boundary -- ghost cells lie within another +c # grid and values are set elsewhere in amr code. + go to 299 + endif +c +c # number of ghost cells lying outside physical domain: + nxr = (xright - xupper + hxmarg)/hx + ibeg = max0(nrow-nxr+1, 1) +c + go to (200,210,220,230) mthbc(2)+1 +c + 200 continue +c # user-specified boundary conditions go here in place of error output + write(6,*) + & '*** ERROR *** mthbc(2)=0 and no BCs specified in bc2amr' + stop + go to 299 + + 210 continue +c # linear-order extrapolation: + do 215 k = 1,nfil + do 215 j = 1,ncol + do 215 i=ibeg,nrow + do 215 m=1,meqn + val(m,i,j,k) = 2.d0*val(m,ibeg-1,j,k) - val(m,ibeg-2,j,k) + 215 continue + go to 299 + + 220 continue +c # periodic: handled elsewhere in amr + go to 299 + + 230 continue +c # solid wall (assumes 2'nd component is velocity or momentum in x): + do 235 k = 1,nfil + do 235 j = 1,ncol + do 235 i=ibeg,nrow + do 235 m=1,meqn + val(m,i,j,k) = val(m,2*ibeg-1-i,j,k) + 235 continue +c # negate the normal velocity: + do 236 k = 1,nfil + do 236 j = 1,ncol + do 236 i=ibeg,nrow + val(2,i,j,k) = -val(2,i,j,k) + 236 continue + go to 299 + + 299 continue +c +c------------------------------------------------------- +c # ylower boundary: +c------------------------------------------------------- + if (yfront .ge. ylower-hymarg) then +c # not a physical boundary -- ghost cells lie within another +c # grid and values are set elsewhere in amr code. + go to 399 + endif +c +c # number of ghost cells lying outside physical domain: + nyf = (ylower+hymarg-yfront)/hy +c + go to (300,310,320,330) mthbc(3)+1 +c + 300 continue + do k = 1,nfil + zcell = zbot + (k-0.5d0)*hz + do i = 1,nrow + xcell = xleft + (i-0.5d0)*hx + do j=1,nyf + ycell = yfront + (j-.5d0)*hy + val(1,i,j,k) = qtrue(xcell,ycell,zcell,time) + enddo + enddo + enddo + go to 399 +c + 310 continue +c # zero-order extrapolation: + do 315 k = 1,nfil + do 315 j=1,nyf + do 315 i=1,nrow + do 315 m=1,meqn + val(m,i,j,k) = val(m,i,nyf+1,k) + 315 continue + go to 399 + + 320 continue +c # periodic: handled elsewhere in amr + go to 399 + + 330 continue +c # solid wall (assumes 3'rd component is velocity or momentum in y): + do 335 k = 1,nfil + do 335 j=1,nyf + do 335 i=1,nrow + do 335 m=1,meqn + val(m,i,j,k) = val(m,i,2*nyf+1-j,k) + 335 continue +c # negate the normal velocity: + do 336 k = 1,nfil + do 336 j=1,nyf + do 336 i=1,nrow + val(3,i,j,k) = -val(3,i,j,k) + 336 continue + go to 399 + + 399 continue +c +c------------------------------------------------------- +c # yupper boundary: +c------------------------------------------------------- + if (yrear .le. yupper+hymarg) then +c # not a physical boundary -- ghost cells lie within another +c # grid and values are set elsewhere in amr code. + go to 499 + endif +c +c # number of ghost cells lying outside physical domain: + nyr = (yrear - yupper + hymarg)/hy + jbeg = max0(ncol-nyr+1, 1) +c + go to (400,410,420,430) mthbc(4)+1 +c + 400 continue +c # user-specified boundary conditions go here in place of error output + write(6,*) + & '*** ERROR *** mthbc(4)=0 and no BCs specified in bc3amr' + stop + go to 499 + + 410 continue +c # linear extrapolation: + do 415 k = 1,nfil + do 415 j=jbeg,ncol + do 415 i=1,nrow + do 415 m=1,meqn + val(m,i,j,k) = 2.d0*val(m,i,jbeg-1,k) - val(m,i,jbeg-2,k) + 415 continue + go to 499 + + 420 continue +c # periodic: handled elsewhere in amr + go to 499 + + 430 continue +c # solid wall (assumes 3'rd component is velocity or momentum in y): + do 435 m=1,meqn + do 435 j=jbeg,ncol + do 435 i=1,nrow + do 435 k = 1,nfil + val(m,i,j,k) = val(m,i,2*jbeg-1-j,k) + 435 continue +c # negate the normal velocity: + do 436 j=jbeg,ncol + do 436 i=1,nrow + do 436 k = 1,nfil + val(3,i,j,k) = -val(3,i,j,k) + 436 continue + go to 499 + + 499 continue + +c +c------------------------------------------------------- +c # zlower boundary: +c------------------------------------------------------- + if (zbot .ge. zlower-hzmarg) then +c # not a physical boundary -- ghost cells lie within another +c # grid and values are set elsewhere in amr code. + go to 599 + endif +c +c # number of ghost cells lying outside physical domain: + nzb = (zlower+hzmarg-zbot)/hz +c + go to (500,510,520,530) mthbc(5)+1 +c + 500 continue + do i = 1,nrow + xcell = xleft + (i-0.5d0)*hx + do j = 1,ncol + ycell = yfront + (j-0.5d0)*hy + do k=1,nzb + zcell = zbot + (k-0.5d0)*hz + val(1,i,j,k) = qtrue(xcell,ycell,zcell,time) + enddo + enddo + enddo + go to 599 +c + 510 continue +c # zero-order extrapolation: + do 515 k=1,nzb + do 515 j = 1,ncol + do 515 i=1,nrow + do 515 m=1,meqn + val(m,i,j,k) = val(m,i,j,nzb+1) + 515 continue + go to 599 + + 520 continue +c # periodic: handled elsewhere in amr + go to 599 + + 530 continue +c # solid wall (assumes 3'rd component is velocity or momentum in y): + do 535 k=1,nzb + do 535 j = 1,ncol + do 535 i=1,nrow + do 535 m=1,meqn + val(m,i,j,k) = val(m,i,j,2*nzb+1-k) + 535 continue +c # negate the normal velocity: + do 536 k = 1,nzb + do 536 j = 1,ncol + do 536 i = 1,nrow + val(4,i,j,k) = -val(4,i,j,k) + 536 continue + go to 599 + + 599 continue +c +c------------------------------------------------------- +c # zupper boundary: +c------------------------------------------------------- + if (ztop .le. zupper+hzmarg) then +c # not a physical boundary -- ghost cells lie within another +c # grid and values are set elsewhere in amr code. + go to 699 + endif +c +c # number of ghost cells lying outside physical domain: + nzt = (ztop - zupper + hzmarg)/hz + kbeg = max0(nfil-nzt+1, 1) +c + go to (600,610,620,630) mthbc(6)+1 +c + 600 continue +c # user-specified boundary conditions go here in place of error output + write(6,*) + & '*** ERROR *** mthbc(6)=0 and no BCs specified in bc3amr' + stop + go to 699 + + 610 continue +c # linear extrapolation: + do 615 k = kbeg,nfil + do 615 j = 1,ncol + do 615 i = 1,nrow + do 615 m = 1,meqn + val(m,i,j,k) = 2.d0*val(m,i,j,kbeg-1) - val(m,i,j,kbeg-2) + 615 continue + go to 699 + + 620 continue +c # periodic: handled elsewhere in amr + go to 699 + + 630 continue +c # solid wall (assumes 3'rd component is velocity or momentum in y): + do 635 k = kbeg,nfil + do 635 j = 1,ncol + do 635 i = 1,nrow + do 635 m=1,meqn + val(m,i,j,k) = val(m,i,j,2*kbeg-1-k) + 635 continue +c # negate the normal velocity: + do 636 k=kbeg,nfil + do 636 j = 1,ncol + do 636 i=1,nrow + val(4,i,j,k) = -val(4,i,j,k) + 636 continue + go to 699 + + 699 continue + + return + end diff --git a/dev/advection_3d_linearTest/flag2refine2.f90 b/dev/advection_3d_linearTest/flag2refine2.f90 new file mode 100644 index 000000000..beead4bc8 --- /dev/null +++ b/dev/advection_3d_linearTest/flag2refine2.f90 @@ -0,0 +1,113 @@ +! ::::::::::::::::::::: flag2refine :::::::::::::::::::::::::::::::::: +! +! User routine to control flagging of points for refinement. +! +! Default version computes spatial difference dq in each direction and +! for each component of q and flags any point where this is greater than +! the tolerance tolsp. This is consistent with what the routine errsp did in +! earlier versions of amrclaw (4.2 and before). +! +! This routine can be copied to an application directory and modified to +! implement some other desired refinement criterion. +! +! Points may also be flagged for refining based on a Richardson estimate +! of the error, obtained by comparing solutions on the current grid and a +! coarsened grid. Points are flagged if the estimated error is larger than +! the parameter tol in amr2ez.data, provided flag_richardson is .true., +! otherwise the coarsening and Richardson estimation is not performed! +! Points are flagged via Richardson in a separate routine. +! +! Once points are flagged via this routine and/or Richardson, the subroutine +! flagregions is applied to check each point against the min_level and +! max_level of refinement specified in any "region" set by the user. +! So flags set here might be over-ruled by region constraints. +! +! q = grid values including ghost cells (bndry vals at specified +! time have already been set, so can use ghost cell values too) +! +! aux = aux array on this grid patch +! +! amrflags = array to be flagged with either the value +! DONTFLAG (no refinement needed) or +! DOFLAG (refinement desired) +! +! tolsp = tolerance specified by user in input file amr2ez.data, used in default +! version of this routine as a tolerance for spatial differences. +! +! :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: +subroutine flag2refine(mx,my,mz,mbc,meqn,maux,xlower,ylower,zlower, & + dx,dy,dz,t,level,tolsp,q,aux,amrflags,DONTFLAG,DOFLAG) + + use amr_module, only: lfine + + implicit none + + ! Subroutine arguments + integer, intent(in) :: mx,my,mz,mbc,meqn,maux,level + real(kind=8), intent(in) :: xlower,ylower,zlower,dx,dy,dz,t,tolsp + + real(kind=8), intent(in) :: q(meqn,1-mbc:mx+mbc,1-mbc:my+mbc,1-mbc:mz+mbc) + real(kind=8), intent(in) :: aux(maux,1-mbc:mx+mbc,1-mbc:my+mbc,1-mbc:mz+mbc) + + ! Flagging + real(kind=8),intent(inout) :: amrflags(1-mbc:mx+mbc,1-mbc:my+mbc,1-mbc:mz+mbc) + real(kind=8), intent(in) :: DONTFLAG + real(kind=8), intent(in) :: DOFLAG + + logical :: allowflag + external allowflag + + ! Locals + integer :: i,j,k,m + real(kind=8) :: x_c,y_c,z_c,x_low,y_low,z_low,x_hi,y_hi,z_hi + real(kind=8) :: dqi(meqn), dqj(meqn), dqk(meqn),dq(meqn) + + ! Initialize flags + amrflags = DONTFLAG + + ! Loop over interior points on this grid + ! (i,j) grid cell is [x_low,x_hi] x [y_low,y_hi], cell center at (x_c,y_c) + ! This information is not needed for the default flagging based on + ! undivided differences, but might be needed in a user's version. + ! Note that if you want to refine only in certain space-time regions, + ! it may be easiest to use the "regions" feature. The flags set here or + ! in the Richardson error estimator are potentially modified by the + ! min_level and max_level specified in any regions. + + z_loop: do k=1,mz + z_low = zlower + (k - 1) * dz + z_c = zlower + (k - 0.5d0) * dz + z_hi = zlower + k * dz + + y_loop: do j=1,my + y_low = ylower + (j - 1) * dy + y_c = ylower + (j - 0.5d0) * dy + y_hi = ylower + j * dy + + x_loop: do i = 1,mx + x_low = xlower + (i - 1) * dx + x_c = xlower + (i - 0.5d0) * dx + x_hi = xlower + i * dx + + ! ----------------------------------------------------------------- + dq = 0.d0 + dqi = abs(q(:,i+1,j,k) - q(:,i-1,j,k)) + dqj = abs(q(:,i,j+1,k) - q(:,i,j-1,k)) + dqk = abs(q(:,i,j,k+1) - q(:,i,j,k-1)) + dq = max(dq,dqi,dqj,dqk) + + ! default checks all components of undivided difference: + do m=1,meqn + if (dq(m) > tolsp) then + amrflags(i,j,k) = DOFLAG + cycle x_loop + endif + enddo + + enddo x_loop + enddo y_loop + enddo z_loop + amrflags(mx/2,my/2,mz/2) = DOFLAG + !if (level .eq. lfine .and. t .gt. 1.) amrflags(mx/2,my/2,mz/2) = DOFLAG + +end subroutine flag2refine diff --git a/dev/advection_3d_linearTest/qinit.f b/dev/advection_3d_linearTest/qinit.f new file mode 100644 index 000000000..a5dd7deca --- /dev/null +++ b/dev/advection_3d_linearTest/qinit.f @@ -0,0 +1,26 @@ +c ===================================================== + subroutine qinit(meqn,mbc,mx,my,mz,xlower,ylower,zlower, + & dx,dy,dz,q,maux,aux) +c ===================================================== +c +c # Set initial conditions for q. +c # Sample scalar equation with linear data to test +c # that second order method stays exact. +c + implicit double precision (a-h,o-z) + dimension q(meqn, 1-mbc:mx+mbc, 1-mbc:my+mbc,1-mbc:mz+mbc) + dimension aux(maux, 1-mbc:mx+mbc, 1-mbc:my+mbc,1-mbc:mz+mbc) + +c # set concentration profile +c --------------------------- + do 20 k=1,mz + zk = zlower + (k-0.5d0)*dz + do 20 j=1,my + yj = ylower + (j-0.5d0)*dy + do 20 i=1,mx + xi = xlower + (i-0.5d0)*dx + q(1,i,j,k) = qtrue(xi,yj,zk,0.d0) + 20 continue + + return + end diff --git a/dev/advection_3d_linearTest/qtrue.f b/dev/advection_3d_linearTest/qtrue.f new file mode 100644 index 000000000..671dcb0c6 --- /dev/null +++ b/dev/advection_3d_linearTest/qtrue.f @@ -0,0 +1,20 @@ + + +c ================================== + double precision function qtrue(x,y,z,t) +c ================================== + implicit double precision (a-h,o-z) + common /cparam/ ubar,vbar,wbar + + x0 = x - ubar*t + y0 = y - vbar*t + z0 = z - wbar*t + +c # evaluate desired initial data at (x0,y0,z0): + + +c linear test function to see if ghost cell interpolation exact + qtrue = 1.d0*x0 + 1.0d0*y0 + 1.d0*z0 + + return + end diff --git a/dev/advection_3d_linearTest/setaux.f b/dev/advection_3d_linearTest/setaux.f new file mode 100644 index 000000000..9e78a60b9 --- /dev/null +++ b/dev/advection_3d_linearTest/setaux.f @@ -0,0 +1,36 @@ +c ================================================================== + subroutine setaux(mbc,mx,my,mz,xlower,ylower,zlower, + & dx,dy,dz,maux,aux) +c ================================================================== +c +c # set auxiliary arrays +c +c # advection +c # aux(i,j,k,1) is u velocity on left face of cell +c # aux(i,j,k,2) is v velocity on bottom face of cell +c # aux(i,j,k,3) is w velocity on back face of cell +c +c + implicit none + + integer mx, my, mz, mbc, maux, i,j,k + double precision xlower, ylower, zlower, dx, dy, dz + double precision aux(maux,1-mbc:mx+mbc,1-mbc:my+mbc,1-mbc:mz+mbc) + double precision ubar, vbar, wbar + + common /cparam/ ubar,vbar,wbar + +c + + do k = 1-mbc,mz+mbc + do j = 1-mbc,my+mbc + do i = 1-mbc,mx+mbc + aux(1,i,j,k) = ubar + aux(2,i,j,k) = vbar + aux(3,i,j,k) = wbar + enddo + enddo + enddo + + return + end diff --git a/dev/advection_3d_linearTest/setplot.py b/dev/advection_3d_linearTest/setplot.py new file mode 100644 index 000000000..bbe591f01 --- /dev/null +++ b/dev/advection_3d_linearTest/setplot.py @@ -0,0 +1,174 @@ + +""" +Set up the plot figures, axes, and items to be done for each frame. + +This module is imported by the plotting routines and then the +function setplot is called to set the plot parameters. +""" + +from clawpack.clawutil.data import ClawData +probdata = ClawData() +probdata.read('setprob.data', force=True) +ubar = probdata.u +vbar = probdata.v +vbar = probdata.w + +#-------------------------- +def setplot(plotdata): +#-------------------------- + + """ + Specify what is to be plotted at each frame. + Input: plotdata, an instance of visclaw.data.ClawPlotData. + Output: a modified version of plotdata. + + """ + + + from clawpack.visclaw import colormaps + + plotdata.clearfigures() # clear any old figures,axes,items data + + + # Figure for pcolor plot + plotfigure = plotdata.new_plotfigure(name='pcolor', figno=0) + + # Set up for axes in this figure: + plotaxes = plotfigure.new_plotaxes() + plotaxes.xlimits = [0,50] + plotaxes.ylimits = [0,50] + plotaxes.title = 'Solution' + plotaxes.scaled = True + + # Set up for item on these axes: + plotitem = plotaxes.new_plotitem(plot_type='2d_pcolor') + plotitem.plot_var = 0 + plotitem.pcolor_cmap = colormaps.red_yellow_blue + plotitem.pcolor_cmin = -150. + plotitem.pcolor_cmax = 150. + plotitem.add_colorbar = True + + plotitem.amr_celledges_show = [0,0] + plotitem.amr_patchedges_show = [0,1] + + + # Figure for pcolor plot + plotfigure = plotdata.new_plotfigure(name='error pcolor', figno=10) + + def error(current_data): + x = current_data.x + y = current_data.y + z = current_data.z + t = current_data.t + q = current_data.q + qtrue = (x-ubar*t) + (y-vbar*t) + (z-wbar*t) + e = q[0,:,:,:] - qtrue + #import pdb; pdb.set_trace() + return e + + # Set up for axes in this figure: + plotaxes = plotfigure.new_plotaxes() + plotaxes.xlimits = [0,50] + plotaxes.ylimits = [0,50] + plotaxes.zlimits = [0,50] + plotaxes.title = 'Error' + plotaxes.scaled = True + + # Set up for item on these axes: + plotitem = plotaxes.new_plotitem(plot_type='2d_pcolor') + plotitem.plot_var = error + plotitem.pcolor_cmap = colormaps.red_yellow_blue + plotitem.pcolor_cmin = -0.000001 + plotitem.pcolor_cmax = 0.000001 + plotitem.add_colorbar = True + + plotitem.amr_celledges_show = [0,0] + plotitem.amr_patchedges_show = [0,1] + + + # Figure for contour plot + plotfigure = plotdata.new_plotfigure(name='contour', figno=1) + + # Set up for axes in this figure: + plotaxes = plotfigure.new_plotaxes() + plotaxes.xlimits = [0,50] + plotaxes.ylimits = [0,50] + plotaxes.zlimits = [0,50] + plotaxes.title = 'Solution' + plotaxes.scaled = True + #plotaxes.afteraxes = addgauges + + # Set up for item on these axes: + plotitem = plotaxes.new_plotitem(plot_type='2d_contour') + plotitem.plot_var = 0 + plotitem.contour_nlevels = 40 + plotitem.contour_min = -150.99 + plotitem.contour_max = 150.99 + plotitem.amr_contour_colors = ['r','g','b'] # color on each level + plotitem.amr_patch_bgcolor = ['#ffeeee', '#eeeeff', '#eeffee'] + plotitem.celledges_show = 0 + plotitem.patchedges_show = 1 + + + # Figure for grid cells + plotfigure = plotdata.new_plotfigure(name='cells', figno=2) + + # Set up for axes in this figure: + plotaxes = plotfigure.new_plotaxes() + plotaxes.xlimits = [0,50] + plotaxes.ylimits = [0,50] + plotaxes.title = 'Grid patches' + plotaxes.scaled = True + + # Set up for item on these axes: + plotitem = plotaxes.new_plotitem(plot_type='2d_patch') + plotitem.amr_patch_bgcolor = ['#ffeeee', '#eeeeff', '#eeffee'] + plotitem.amr_celledges_show = [0,0,0] + plotitem.amr_patchedges_show = [1,1] + + #----------------------------------------- + # Figures for gauges + #----------------------------------------- + plotfigure = plotdata.new_plotfigure(name='q', figno=300, \ + type='each_gauge') + plotfigure.clf_each_gauge = True + + # Set up for axes in this figure: + plotaxes = plotfigure.new_plotaxes() + plotaxes.xlimits = 'auto' + plotaxes.ylimits = 'auto' + plotaxes.zlimits = 'auto' + plotaxes.title = 'q' + + # Plot q as blue curve: + plotitem = plotaxes.new_plotitem(plot_type='1d_plot') + plotitem.plot_var = 0 + plotitem.plotstyle = 'b-' + + + # Parameters used only when creating html and/or latex hardcopy + # e.g., via visclaw.frametools.printframes: + + plotdata.printfigs = True # print figures + plotdata.print_format = 'png' # file format + plotdata.print_framenos = 'all' # list of frames to print + plotdata.print_fignos = 'all' # list of figures to print + plotdata.html = True # create html files of plots? + plotdata.html_homelink = '../README.html' # pointer for top of index + plotdata.html_movie = 'JSAnimation' # new style, or "4.x" for old style + plotdata.latex = True # create latex file of plots? + plotdata.latex_figsperline = 2 # layout of plots + plotdata.latex_framesperline = 1 # layout of plots + plotdata.latex_makepdf = False # also run pdflatex? + + return plotdata + + + +# To plot gauge locations on pcolor or contour plot, use this as +# an afteraxis function: + +def addgauges(current_data): + from clawpack.visclaw import gaugetools + gaugetools.plot_gauge_locations(current_data.plotdata, \ + gaugenos='all', format_string='ko', add_labels=True) diff --git a/dev/advection_3d_linearTest/setprob.f b/dev/advection_3d_linearTest/setprob.f new file mode 100644 index 000000000..bb15e9118 --- /dev/null +++ b/dev/advection_3d_linearTest/setprob.f @@ -0,0 +1,23 @@ + subroutine setprob + implicit double precision (a-h,o-z) + character*25 fname + common /cparam/ ubar,vbar,wbar +c +c # Set the velocity for scalar advection +c # These values are passed to the Riemann solvers rpn2.f and rpt2.f +c # in a common block +c + +c + iunit = 7 + fname = 'setprob.data' +c # open the unit with new routine from Clawpack 4.4 to skip over +c # comment lines starting with #: + call opendatafile(iunit, fname) + + read(7,*) ubar + read(7,*) vbar + read(7,*) wbar + + return + end diff --git a/dev/advection_3d_linearTest/setrun.py b/dev/advection_3d_linearTest/setrun.py new file mode 100644 index 000000000..cd6f65d44 --- /dev/null +++ b/dev/advection_3d_linearTest/setrun.py @@ -0,0 +1,368 @@ +""" +Module to set up run time parameters for Clawpack. + +The values set in the function setrun are then written out to data files +that will be read in by the Fortran code. + +""" + +import os +import numpy as np + +# Set which test to do by setting inflow_side to 'x', 'y', or 'z' +inflow_side = 'x' + +# Note: The velocity is set to (u,v,w) = (1,1,1) so the left edge in +# each direction is the inflow boundary. +# With +# clawdata.transverse_waves = 22 +# this method should be stable up to Courant number 1, but does not +# currently seem to be stable with cfl_desired = 0.9. +# This requires more investigation. + + + +#------------------------------ +def setrun(claw_pkg='amrclaw'): +#------------------------------ + + """ + Define the parameters used for running Clawpack. + + INPUT: + claw_pkg expected to be "amrclaw" for this setrun. + + OUTPUT: + rundata - object of class ClawRunData + + """ + + from clawpack.clawutil import data + + + assert claw_pkg.lower() == 'amrclaw', "Expected claw_pkg = 'amrclaw'" + + num_dim = 3 + rundata = data.ClawRunData(claw_pkg, num_dim) + + #------------------------------------------------------------------ + # Problem-specific parameters to be written to setprob.data: + #------------------------------------------------------------------ + + + probdata = rundata.new_UserData(name='probdata',fname='setprob.data') + probdata.add_param('u', 1.0, 'ubar advection velocity') + probdata.add_param('v', 1.0, 'vbar advection velocity') + probdata.add_param('w', 1.0, 'wbar advection velocity') + + #------------------------------------------------------------------ + # Standard Clawpack parameters to be written to claw.data: + # (or to amrclaw.data for AMR) + #------------------------------------------------------------------ + + clawdata = rundata.clawdata # initialized when rundata instantiated + + + # Set single grid parameters first. + # See below for AMR parameters. + + + # --------------- + # Spatial domain: + # --------------- + + # Number of space dimensions: + clawdata.num_dim = num_dim + + # Lower and upper edge of computational domain: + clawdata.lower[0] = 0.000000e+00 # xlower + clawdata.upper[0] = 1.000000e+00 # xupper + clawdata.lower[1] = 0.000000e+00 # ylower + clawdata.upper[1] = 1.000000e+00 # yupper + clawdata.lower[2] = 0.000000e+00 # zlower + clawdata.upper[2] = 1.000000e+00 # zupper + + # Number of grid cells: + clawdata.num_cells[0] = 40 # mx + clawdata.num_cells[1] = 40 # my + clawdata.num_cells[2] = 40 # mz + + + # --------------- + # Size of system: + # --------------- + + # Number of equations in the system: + clawdata.num_eqn = 1 + + # Number of auxiliary variables in the aux array (initialized in setaux) + clawdata.num_aux = 3 + + # Index of aux array corresponding to capacity function, if there is one: + clawdata.capa_index = 0 + + + # ------------- + # Initial time: + # ------------- + + clawdata.t0 = 0.0 + + + # Restart from checkpoint file of a previous run? + # Note: If restarting, you must also change the Makefile to set: + # RESTART = True + # If restarting, t0 above should be from original run, and the + # restart_file 'fort.chkNNNNN' specified below should be in + # the OUTDIR indicated in Makefile. + + clawdata.restart = False # True to restart from prior results + clawdata.restart_file = 'fort.chk00006' # File to use for restart data + + + # ------------- + # Output times: + #-------------- + + # Specify at what times the results should be written to fort.q files. + # Note that the time integration stops after the final output time. + + clawdata.output_style = 3 + + if clawdata.output_style==1: + # Output ntimes frames at equally spaced times up to tfinal: + # Can specify num_output_times = 0 for no output + clawdata.num_output_times = 10 + clawdata.tfinal = 2.0 + clawdata.output_t0 = True # output at initial (or restart) time? + + elif clawdata.output_style == 2: + # Specify a list or numpy array of output times: + # Include t0 if you want output at the initial time. + clawdata.output_times = [0., 0.1] + + elif clawdata.output_style == 3: + # Output every step_interval timesteps over total_steps timesteps: + clawdata.output_step_interval = 1 + clawdata.total_steps = 10 + clawdata.output_t0 = True # output at initial (or restart) time? + + + clawdata.output_format = 'ascii' # 'ascii', 'binary', 'netcdf' + + clawdata.output_q_components = 'all' # could be list such as [True,True] + clawdata.output_aux_components = 'none' # could be list + clawdata.output_aux_onlyonce = True # output aux arrays only at t0 + + + # --------------------------------------------------- + # Verbosity of messages to screen during integration: + # --------------------------------------------------- + + # The current t, dt, and cfl will be printed every time step + # at AMR levels <= verbosity. Set verbosity = 0 for no printing. + # (E.g. verbosity == 2 means print only on levels 1 and 2.) + clawdata.verbosity = 0 + + + + # -------------- + # Time stepping: + # -------------- + + # if dt_variable==True: variable time steps used based on cfl_desired, + # if dt_variable==False: fixed time steps dt = dt_initial always used. + clawdata.dt_variable = True + + # Initial time step for variable dt. + # (If dt_variable==0 then dt=dt_initial for all steps) + clawdata.dt_initial = 0.01 + + # Max time step to be allowed if variable dt used: + clawdata.dt_max = 1e+99 + + # Desired Courant number if variable dt used + clawdata.cfl_desired = 0.9 + # max Courant number to allow without retaking step with a smaller dt: + clawdata.cfl_max = 1.0 + + # Maximum number of time steps to allow between output times: + clawdata.steps_max = 100000 + + + # ------------------ + # Method to be used: + # ------------------ + + # Order of accuracy: 1 => Godunov, 2 => Lax-Wendroff plus limiters + clawdata.order = 2 + + # Use dimensional splitting? + clawdata.dimensional_split = 'unsplit' + + # For unsplit method, transverse_waves can be + # 0 or 'none' ==> donor cell (only normal solver used) + # 1 or 'increment' ==> corner transport of waves + # 2 or 'all' ==> corner transport of 2nd order corrections too + clawdata.transverse_waves = 2 + + + # Number of waves in the Riemann solution: + clawdata.num_waves = 1 + + # List of limiters to use for each wave family: + # Required: len(limiter) == num_waves + # Some options: + # 0 or 'none' ==> no limiter (Lax-Wendroff) + # 1 or 'minmod' ==> minmod + # 2 or 'superbee' ==> superbee + # 3 or 'vanleer' ==> van Leer + # 4 or 'mc' ==> MC limiter + clawdata.limiter = ['none'] + + clawdata.use_fwaves = False # True ==> use f-wave version of algorithms + + # Source terms splitting: + # src_split == 0 or 'none' ==> no source term (src routine never called) + # src_split == 1 or 'godunov' ==> Godunov (1st order) splitting used, + # src_split == 2 or 'strang' ==> Strang (2nd order) splitting used, not recommended. + clawdata.source_split = 'none' + + + # -------------------- + # Boundary conditions: + # -------------------- + + # Number of ghost cells (usually 2) + clawdata.num_ghost = 2 + + # Choice of BCs at xlower and xupper: + # 0 or 'user' => user specified (must modify bcNamr.f to use this option) + # 1 or 'extrap' => extrapolation (non-reflecting outflow) + # 2 or 'periodic' => periodic (must specify this at both boundaries) + # 3 or 'wall' => solid wall for systems where q(2) is normal velocity + + # Note: one of these will be changed below based on inflow_side... + + clawdata.bc_lower[0] = 'user' # at xlower + clawdata.bc_upper[0] = 'extrap' # at xupper + + clawdata.bc_lower[1] = 'user' # at ylower + clawdata.bc_upper[1] = 'extrap' # at yupper + + clawdata.bc_lower[2] = 'user' # at zlower + clawdata.bc_upper[2] = 'extrap' # at zupper + + # --------------- + # Gauges: + # --------------- + rundata.gaugedata.gauges = [] + # for gauges append lines of the form [gaugeno, x, y, t1, t2] + + + # -------------- + # Checkpointing: + # -------------- + + # Specify when checkpoint files should be created that can be + # used to restart a computation. + + clawdata.checkpt_style = 0 + + if clawdata.checkpt_style == 0: + # Do not checkpoint at all + pass + + elif clawdata.checkpt_style == 1: + # Checkpoint only at tfinal. + pass + + elif clawdata.checkpt_style == 2: + # Specify a list of checkpoint times. + clawdata.checkpt_times = [0.1,0.15] + + elif clawdata.checkpt_style == 3: + # Checkpoint every checkpt_interval timesteps (on Level 1) + # and at the final time. + clawdata.checkpt_interval = 5 + + # --------------- + # AMR parameters: + # --------------- + amrdata = rundata.amrdata + + # max number of refinement levels: + amrdata.amr_levels_max = 2 + + # List of refinement ratios at each level (length at least amr_level_max-1) + amrdata.refinement_ratios_x = [4, 2, 2, 2] + amrdata.refinement_ratios_y = [4, 2, 2, 2] + amrdata.refinement_ratios_z = [4, 2, 2, 2] + amrdata.refinement_ratios_t = [4, 2, 2, 2] + + + # Specify type of each aux variable in amrdata.auxtype. + # This must be a list of length num_aux, each element of which is one of: + # 'center', 'capacity', 'xleft', or 'yleft' (see documentation). + amrdata.aux_type = ['xleft', 'yleft', 'zleft'] + + + # Flag for refinement based on Richardson error estimater: + amrdata.flag_richardson = False # use Richardson? + amrdata.flag_richardson_tol = 1.000000e+00 # Richardson tolerance + + # Flag for refinement using routine flag2refine: + amrdata.flag2refine = True # use this? + amrdata.flag2refine_tol = 100.55 # tolerance used in this routine + # User can modify flag2refine to change the criterion for flagging. + # Default: check max-norm of difference between q in a cell and + # each of its neighbors. + + # steps to take on each level L between regriddings of level L+1: + amrdata.regrid_interval = 400000 + + # width of buffer zone around flagged points: + # (typically the same as regrid_interval so waves don't escape): + amrdata.regrid_buffer_width = 7 + + # clustering alg. cutoff for (# flagged pts) / (total # of cells refined) + # (closer to 1.0 => more small grids may be needed to cover flagged cells) + amrdata.clustering_cutoff = 0.7 + + # print info about each regridding up to this level: + amrdata.verbosity_regrid = 3 + + + # --------------- + # Regions: + # --------------- + regions = rundata.regiondata.regions + # to specify regions of refinement append lines of the form + # [minlevel,maxlevel,t1,t2,x1,x2,y1,y2] + + + # ----- For developers ----- + # Toggle debugging print statements: + amrdata.dprint = False # print domain flags + amrdata.eprint = False # print err est flags + amrdata.edebug = False # even more err est flags + amrdata.gprint = False # grid bisection/clustering + amrdata.nprint = False # proper nesting output + amrdata.pprint = False # proj. of tagged points + amrdata.rprint = True # print regridding summary + amrdata.sprint = False # space/memory output + amrdata.tprint = True # time step reporting each level + amrdata.uprint = False # update/upbnd reporting + + + return rundata + + # end of function setrun + # ---------------------- + + +if __name__ == '__main__': + # Set up run-time parameters and write all data files. + import sys + rundata = setrun(*sys.argv[1:]) + rundata.write() diff --git a/dev/advection_3d_linearTest/valout.f b/dev/advection_3d_linearTest/valout.f new file mode 100644 index 000000000..37530dc46 --- /dev/null +++ b/dev/advection_3d_linearTest/valout.f @@ -0,0 +1,260 @@ +c +c ----------------------------------------------------- +c + subroutine valout (lst, lend, time, nvar, naux) +c +c # Output the results for a general system of conservation laws +c # in 3 dimensions +c +c # Write the results to the file fort.q