diff --git a/dev/advection_2d_linearTest/Makefile b/dev/advection_2d_linearTest/Makefile new file mode 100644 index 000000000..cb48321cf --- /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)/domain.f \ + $(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 \ + ./valout.f \ + $(AMRLIB)/flagregions2.f90 \ + $(AMRLIB)/errest.f \ + $(AMRLIB)/errf1.f \ + $(AMRLIB)/gfixup.f \ + $(AMRLIB)/filval.f90 \ + $(AMRLIB)/filpatch.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)/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)/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..2c5353c46 --- /dev/null +++ b/dev/advection_2d_linearTest/bc2amr.f @@ -0,0 +1,348 @@ +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 implementedsrt +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 + do i=1,nxl + xcell = xleft +(i-0.5d0)*hx + do j = 1,ncol + ycell = ybot + (j-0.5d0)*hy + val(1,i,j) = qtrue(xcell,ycell,time) +! 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 +! 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 +! endif + enddo + 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 # first-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 + do i = 1,nrow + xcell = xleft + (i-0.5d0)*hx + do j = 1, nyb + ycell = ybot + (j-0.5d0)*hy + val(1,i,nyb) = qtrue(xcell,ycell,time) + + ! 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 +! 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 +! endif + enddo + 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 # first-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..de0c9f418 --- /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..36a8471a5 --- /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..c139d8e67 --- /dev/null +++ b/dev/advection_2d_linearTest/setplot.py @@ -0,0 +1,177 @@ + +""" +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 + +#import os +#import clawpack.clawutil.data as clawutil +#clawdata = clawutil.ClawInputData(2) +#clawdata.read(os.path.join(plotdata.outdir,'claw.data')) +#xlimits = [clawdata.lower[0],clawdata.upper[0]] +#ylimits = [clawdata.lower[1],clawdata.upper[1]] +xlimits = [0, 26] +ylimits = [0, 28] +#-------------------------- +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 = xlimits + plotaxes.ylimits = ylimits + 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 = 2*(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 = xlimits + plotaxes.ylimits = ylimits + 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 = xlimits + plotaxes.ylimits = ylimits + 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 = xlimits + plotaxes.ylimits = ylimits + 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..8a6f60664 --- /dev/null +++ b/dev/advection_2d_linearTest/setrun.py @@ -0,0 +1,352 @@ +""" +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] = 26. # xupper + clawdata.lower[1] = 0. # ylower + clawdata.upper[1] = 28. # yupper + + # Number of grid cells: + clawdata.num_cells[0] = 52 # mx + clawdata.num_cells[1] = 56 # 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.chk00016' # 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 = 1 + 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 = 1 + + # --------------- + # 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 = [4,4,2,8] + amrdata.refinement_ratios_y = [4,4,2,8] + amrdata.refinement_ratios_t = [4,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 = 4 + + # width of buffer zone around flagged points: + # (typically the same as regrid_interval so waves don't escape): + # donna amrdata.regrid_buffer_width = 7 + amrdata.regrid_buffer_width = 0 + + # 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] + regions.append([3,3,-1.0,1,.149625,.359375,0,1]) + regions.append([2,3,-1.0,1,0.0,.5625,0,1]) + + + + # ----- 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 = False # print regridding summary + amrdata.sprint = False # space/memory output + amrdata.tprint = False # 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_2d_linearTest/valout.f b/dev/advection_2d_linearTest/valout.f new file mode 100644 index 000000000..03170ed2a --- /dev/null +++ b/dev/advection_2d_linearTest/valout.f @@ -0,0 +1,277 @@ +c +c ----------------------------------------------------- +c + subroutine valout (lst, lend, time, nvar, naux) +c + use amr_module + implicit double precision (a-h,o-z) + character*10 fname1, fname2, fname3, fname4, fname5 + +c # Output the results for a general system of conservation laws +c # in 2 dimensions +c +c # Write the results to the file fort.q