From 7cafbe908f7526ae90840032b0c99e7f8c0a0e39 Mon Sep 17 00:00:00 2001 From: Chaulio Ferreira Date: Tue, 10 Apr 2018 15:51:14 +0200 Subject: [PATCH 1/3] Vectorization for Riemann solvers for the single-layer shallow water equations (chile2010 example). To use it, define GEOCLAW_VEC=1 before compiling. More info in the chile2010 Makefile. --- examples/tsunami/chile2010/Makefile | 30 ++- src/2d/shallow/flux2fw_vec.f | 270 ++++++++++++++++++++ src/2d/shallow/inlinelimiter_vec.f | 110 ++++++++ src/2d/shallow/qad_vec.f | 374 ++++++++++++++++++++++++++++ src/2d/shallow/step2_vec.f90 | 235 +++++++++++++++++ 5 files changed, 1015 insertions(+), 4 deletions(-) create mode 100644 src/2d/shallow/flux2fw_vec.f create mode 100644 src/2d/shallow/inlinelimiter_vec.f create mode 100644 src/2d/shallow/qad_vec.f create mode 100644 src/2d/shallow/step2_vec.f90 diff --git a/examples/tsunami/chile2010/Makefile b/examples/tsunami/chile2010/Makefile index 2106f068a..5204bba29 100644 --- a/examples/tsunami/chile2010/Makefile +++ b/examples/tsunami/chile2010/Makefile @@ -7,6 +7,10 @@ CLAWMAKE = $(CLAW)/clawutil/src/Makefile.common # make .help # at the unix prompt. +# Define GEOCLAW_VEC as "1" to use the code with vectorized Riemann solvers. +# The following compiler flags are suggested for vectorization (Intel Compiler): +# -O2 -fopenmp -align array64byte -xHost + # Adjust these variables if desired: # ---------------------------------- @@ -47,10 +51,28 @@ EXCLUDE_SOURCES = \ MODULES = \ -SOURCES = \ - $(CLAW)/riemann/src/rpn2_geoclaw.f \ - $(CLAW)/riemann/src/rpt2_geoclaw.f \ - $(CLAW)/riemann/src/geoclaw_riemann_utils.f \ +ifneq ($(GEOCLAW_VEC),1) + SOURCES = \ + $(CLAW)/riemann/src/rpn2_geoclaw.f \ + $(CLAW)/riemann/src/rpt2_geoclaw.f \ + $(CLAW)/riemann/src/geoclaw_riemann_utils.f +else + EXCLUDE_SOURCES = \ + $(GEOLIB)/qad.f \ + $(GEOLIB)/flux2fw.f \ + $(GEOLIB)/step2.f \ + $(AMRLIB)/inlinelimiter.f + + SOURCES = \ + $(CLAW)/riemann/src/rpn2_geoclaw_vec.f90 \ + $(CLAW)/riemann/src/rpt2_geoclaw_vec.f90 \ + $(CLAW)/riemann/src/geoclaw_riemann_utils_vec.f90 \ + $(GEOLIB)/qad_vec.f \ + $(GEOLIB)/flux2fw_vec.f \ + $(GEOLIB)/step2_vec.f90 \ + $(GEOLIB)/inlinelimiter_vec.f +endif + #------------------------------------------------------------------- # Include Makefile containing standard definitions and make options: diff --git a/src/2d/shallow/flux2fw_vec.f b/src/2d/shallow/flux2fw_vec.f new file mode 100644 index 000000000..6b27e5d58 --- /dev/null +++ b/src/2d/shallow/flux2fw_vec.f @@ -0,0 +1,270 @@ +! This code has been adapted for vectorization. +! For more details on how to use it, see the Makefile of the chile2010 example. + +c +c +c ===================================================== + subroutine flux2(ixy,maxm,meqn,maux,mbc,mx, + & q1d,dtdx1d,aux1,aux2,aux3, + & faddm,faddp,gaddm,gaddp,cfl1d, + & rpn2,rpt2) +!! & fwave,s,amdq,apdq,rpn2,rpt2) +c ===================================================== +c +c # clawpack routine ... modified for AMRCLAW +c +c-------------------------------------------------------------------- +c # flux2fw is a modified version of flux2 to use fwave instead of wave. +c # A modified Riemann solver rp2n must be used in conjunction with this +c # routine, which returns fwave's instead of wave's. +c # See http://amath.washington.edu/~claw/fwave.html +c +c # Limiters are applied to the fwave's, and the only significant +c # modification of this code is in the "do 119" loop, for the +c # second order corrections. +c +c-------------------------------------------------------------------- +c +c +c # Compute the modification to fluxes f and g that are generated by +c # all interfaces along a 1D slice of the 2D grid. +c # ixy = 1 if it is a slice in x +c # 2 if it is a slice in y +c # This value is passed into the Riemann solvers. The flux modifications +c # go into the arrays fadd and gadd. The notation is written assuming +c # we are solving along a 1D slice in the x-direction. +c +c # fadd(i,.) modifies F to the left of cell i +c # gadd(i,.,1) modifies G below cell i +c # gadd(i,.,2) modifies G above cell i +c +c # The method used is specified by method(2:3): +c +c method(2) = 1 if only first order increment waves are to be used. +c = 2 if second order correction terms are to be added, with +c a flux limiter as specified by mthlim. +c +c method(3) = 0 if no transverse propagation is to be applied. +c Increment and perhaps correction waves are propagated +c normal to the interface. +c = 1 if transverse propagation of increment waves +c (but not correction waves, if any) is to be applied. +c = 2 if transverse propagation of correction waves is also +c to be included. +c +c Note that if mcapa>0 then the capa array comes into the second +c order correction terms, and is already included in dtdx1d: +c If ixy = 1 then +c dtdx1d(i) = dt/dx if mcapa= 0 +c = dt/(dx*aux(i,jcom,mcapa)) if mcapa = 1 +c If ixy = 2 then +c dtdx1d(j) = dt/dy if mcapa = 0 +c = dt/(dy*aux(icom,j,mcapa)) if mcapa = 1 +c +c Notation: +c The jump in q (q1d(i,:)-q1d(i-1,:)) is split by rpn2 into +c amdq = the left-going flux difference A^- Delta q +c apdq = the right-going flux difference A^+ Delta q +c Each of these is split by rpt2 into +c bmasdq = the down-going transverse flux difference B^- A^* Delta q +c bpasdq = the up-going transverse flux difference B^+ A^* Delta q +c where A^* represents either A^- or A^+. +c + +c # modifications for GeoClaw + +c--------------------------flux2fw_geo.f-------------------------- +c This version of flux2fw.f is modified slightly to be used with +c step2_geo.f The only modification is for the first-order +c mass fluxes, faddm(i,j,1) and faddp(i,j,1), so that those terms are true +c interface fluxes. +c +c The only change is in loop 40 +c to revert to the original version, set relimit = .false. +c---------------------last modified 1/04/05----------------------------- + + use amr_module, only: use_fwaves, mwaves, method + use amr_module, only: mthlim + use geoclaw_module, only: coordinate_system, earth_radius, deg2rad + + implicit double precision (a-h,o-z) + + external rpn2, rpt2 + dimension q1d(1-mbc:maxm+mbc, meqn) + + dimension bmasdq(1-mbc:maxm+mbc, meqn) + dimension bpasdq(1-mbc:maxm+mbc,meqn) + dimension cqxx(1-mbc:maxm+mbc, meqn) + dimension faddm(1-mbc:maxm+mbc,meqn) + dimension faddp(1-mbc:maxm+mbc,meqn) + dimension gaddm(1-mbc:maxm+mbc,meqn, 2) + dimension gaddp(1-mbc:maxm+mbc,meqn, 2) + dimension dtdx1d(1-mbc:maxm+mbc) + dimension aux1(1-mbc:maxm+mbc, maux) + dimension aux2(1-mbc:maxm+mbc, maux) + dimension aux3(1-mbc:maxm+mbc, maux) +c + dimension s(1-mbc:maxm+mbc, mwaves) + dimension fwave(1-mbc:maxm+mbc, meqn, mwaves) + dimension amdq(1-mbc:maxm+mbc, meqn) + dimension apdq(1-mbc:maxm+mbc, meqn) +c + logical limit, relimit + + relimit = .false. +c + limit = .false. + do 5 mw=1,mwaves + if (mthlim(mw) .gt. 0) limit = .true. + 5 continue +c +c # initialize flux increments: +c ----------------------------- +c + + do 30 jside=1,2 + do 20 m=1,meqn + do 10 i = 1-mbc, mx+mbc + faddm(i,m) = 0.d0 + faddp(i,m) = 0.d0 + gaddm(i,m,jside) = 0.d0 + gaddp(i,m,jside) = 0.d0 + 10 continue + 20 continue + 30 continue + +c +c +c # solve Riemann problem at each interface and compute Godunov updates +c --------------------------------------------------------------------- +c + call rpn2(ixy,maxm,meqn,mwaves,maux,mbc,mx,q1d,q1d, + & aux2,aux2,fwave,s,amdq,apdq) + +c +c # Set fadd for the donor-cell upwind method (Godunov) + if (ixy.eq.1) mu=2 + if (ixy.eq.2) mu=3 + !DIR$ SIMD PRIVATE(dxdc) + do 40 i=1-mbc+1,mx+mbc-1 + if (coordinate_system.eq.2) then + if (ixy.eq.1) then + dxdc=earth_radius*deg2rad + else + dxdc=earth_radius*cos(aux2(i,3))*deg2rad +! if (ixy.eq.2) dxdc=earth_radius*cos(aux2(3,i))*deg2rad !why test again + endif + else + dxdc=1.d0 + endif + + do m=1,meqn + faddp(i,m) = faddp(i,m) - apdq(i,m) + faddm(i,m) = faddm(i,m) + amdq(i,m) + enddo + if (relimit) then + faddp(i,1) = faddp(i,1) + dxdc*q1d(i,mu) + faddm(i,1) = faddp(i,1) + endif + 40 continue +c +c # compute maximum wave speed for checking Courant number: + cfl1d = 0.d0 + do 50 mw=1,mwaves + !DIR$ SIMD + do 50 i=1,mx+1 +c # if s>0 use dtdx1d(i) to compute CFL, +c # if s<0 use dtdx1d(i-1) to compute CFL: + cfl1d = dmax1(cfl1d, dtdx1d(i)*s(i,mw), + & -dtdx1d(i-1)*s(i,mw)) + + 50 continue +c + if (method(2).eq.1) go to 130 +c +c # modify F fluxes for second order q_{xx} correction terms: +c ----------------------------------------------------------- +c +c # apply limiter to fwaves: + if (limit) call limiter(maxm,meqn,mwaves,mbc,mx,fwave,s,mthlim) +c + !DIR$ SIMD + do 120 i = 1, mx+1 +c +c # For correction terms below, need average of dtdx in cell +c # i-1 and i. Compute these and overwrite dtdx1d: +c + dtdx1d(i-1) = 0.5d0 * (dtdx1d(i-1) + dtdx1d(i)) +c + do 120 m=1,meqn + cqxx(i,m) = 0.d0 + do 119 mw=1,mwaves +c +c # second order corrections: + cqxx(i,m) = cqxx(i,m) + dsign(1.d0,s(i,mw)) + & * (1.d0 - dabs(s(i,mw))*dtdx1d(i-1)) * fwave(i,m,mw) +c + 119 continue + faddm(i,m) = faddm(i,m) + 0.5d0 * cqxx(i,m) + faddp(i,m) = faddp(i,m) + 0.5d0 * cqxx(i,m) + 120 continue +c +c + 130 continue +c + if (method(3).eq.0) go to 999 !# no transverse propagation +c + if (method(2).gt.1 .and. method(3).eq.2) then +c # incorporate cqxx into amdq and apdq so that it is split also. + !DIR$ SIMD + do 150 i = 1, mx+1 + do 150 m=1,meqn + amdq(i,m) = amdq(i,m) + cqxx(i,m) + apdq(i,m) = apdq(i,m) - cqxx(i,m) + 150 continue + endif +c +c +c # modify G fluxes for transverse propagation +c -------------------------------------------- +c +c +c # split the left-going flux difference into down-going and up-going: + call rpt2(ixy,1,maxm,meqn,mwaves,maux,mbc,mx, + & q1d,q1d,aux1,aux2,aux3, + & amdq,bmasdq,bpasdq) +c +c # modify flux below and above by B^- A^- Delta q and B^+ A^- Delta q: + !DIR$ SIMD PRIVATE(gupdate) + do 160 i = 1, mx+1 + do 160 m=1,meqn + gupdate = 0.5d0*dtdx1d(i-1) * bmasdq(i,m) + gaddm(i-1,m,1) = gaddm(i-1,m,1) - gupdate + gaddp(i-1,m,1) = gaddp(i-1,m,1) - gupdate +c + gupdate = 0.5d0*dtdx1d(i-1) * bpasdq(i,m) + gaddm(i-1,m,2) = gaddm(i-1,m,2) - gupdate + gaddp(i-1,m,2) = gaddp(i-1,m,2) - gupdate + 160 continue +c +c # split the right-going flux difference into down-going and up-going: + call rpt2(ixy,2,maxm,meqn,mwaves,maux,mbc,mx, + & q1d,q1d,aux1,aux2,aux3, + & apdq,bmasdq,bpasdq) +c +c # modify flux below and above by B^- A^+ Delta q and B^+ A^+ Delta q: + !DIR$ SIMD PRIVATE(gupdate) + do 180 i = 1, mx+1 + do 180 m=1,meqn + gupdate = 0.5d0*dtdx1d(i-1) * bmasdq(i,m) + gaddm(i,m,1) = gaddm(i,m,1) - gupdate + gaddp(i,m,1) = gaddp(i,m,1) - gupdate +c + gupdate = 0.5d0*dtdx1d(i-1) * bpasdq(i,m) + gaddm(i,m,2) = gaddm(i,m,2) - gupdate + gaddp(i,m,2) = gaddp(i,m,2) - gupdate + 180 continue +c + 999 continue + return + end diff --git a/src/2d/shallow/inlinelimiter_vec.f b/src/2d/shallow/inlinelimiter_vec.f new file mode 100644 index 000000000..07941a4ae --- /dev/null +++ b/src/2d/shallow/inlinelimiter_vec.f @@ -0,0 +1,110 @@ +! This code has been adapted for vectorization. +! For more details on how to use it, see the Makefile of the chile2010 example. + +c +c ===================================================== + subroutine limiter(maxm,meqn,mwaves,mbc,mx,wave,s,mthlim) +c ===================================================== +c +c # Apply a limiter to the waves. +c +c # Version of December, 2002. +c # Modified from the original CLAWPACK routine to eliminate calls +c # to philim. Since philim was called for every wave at each cell +c # interface, this was adding substantial overhead in some cases. +c +c # The limiter is computed by comparing the 2-norm of each wave with +c # the projection of the wave from the interface to the left or +c # right onto the current wave. For a linear system this would +c # correspond to comparing the norms of the two waves. For a +c # nonlinear problem the eigenvectors are not colinear and so the +c # projection is needed to provide more limiting in the case where the +c # neighboring wave has large norm but points in a different direction +c # in phase space. +c +c # The specific limiter used in each family is determined by the +c # value of the corresponding element of the array mthlim. +c # Note that a different limiter may be used in each wave family. +c +c # dotl and dotr denote the inner product of wave with the wave to +c # the left or right. The norm of the projections onto the wave are then +c # given by dotl/wnorm2 and dotr/wnorm2, where wnorm2 is the 2-norm +c # of wave. +c + implicit double precision (a-h,o-z) + dimension mthlim(mwaves) + dimension wave(1-mbc:maxm+mbc, meqn, mwaves) + dimension s(1-mbc:maxm+mbc, mwaves) +c +c + do 200 mw=1,mwaves + if (mthlim(mw) .eq. 0) go to 200 + dotr = 0.d0 + do 190 i = 0, mx+1 + wnorm2 = 0.d0 + dotl = dotr + dotr = 0.d0 + do 5 m=1,meqn + wnorm2 = wnorm2 + wave(i,m,mw)**2 + dotr = dotr + wave(i,m,mw)*wave(i+1,m,mw) + 5 continue + if (i.eq.0) go to 190 + if (wnorm2.eq.0.d0) go to 190 +c + if (s(i,mw) .gt. 0.d0) then + r = dotl / wnorm2 + else + r = dotr / wnorm2 + endif +c + go to (10,20,30,40,50) mthlim(mw) +c + 10 continue +c -------- +c # minmod +c -------- + wlimitr = dmax1(0.d0, dmin1(1.d0, r)) + go to 170 +c + 20 continue +c ---------- +c # superbee +c ---------- + wlimitr = dmax1(0.d0, dmin1(1.d0, 2.d0*r), dmin1(2.d0, r)) + go to 170 +c + 30 continue +c ---------- +c # van Leer +c ---------- + wlimitr = (r + dabs(r)) / (1.d0 + dabs(r)) + go to 170 +c + 40 continue +c ------------------------------ +c # monotinized centered +c ------------------------------ + c = (1.d0 + r)/2.d0 + wlimitr = dmax1(0.d0, dmin1(c, 2.d0, 2.d0*r)) + go to 170 +c + 50 continue +c ------------------------------ +c # Beam-Warming +c ------------------------------ + wlimitr = r + go to 170 +c + 170 continue +c +c # apply limiter to waves: +c + do 180 m=1,meqn + wave(i,m,mw) = wlimitr * wave(i,m,mw) + 180 continue + + 190 continue + 200 continue +c + return + end diff --git a/src/2d/shallow/qad_vec.f b/src/2d/shallow/qad_vec.f new file mode 100644 index 000000000..230dff0f1 --- /dev/null +++ b/src/2d/shallow/qad_vec.f @@ -0,0 +1,374 @@ +! This code has been adapted for vectorization. +! For more details on how to use it, see the Makefile of the chile2010 example. + +c +c ------------------------------------------------------------- +c + subroutine qad(valbig,mitot,mjtot,nvar, + . svdflx,qc1d,lenbc,lratiox,lratioy,hx,hy, + . maux,aux,auxc1d,delt,mptr) + + use amr_module + implicit double precision (a-h, o-z) + + + logical qprint + + dimension valbig(nvar,mitot,mjtot) + dimension qc1d(nvar,lenbc) + dimension svdflx(nvar,lenbc) + dimension aux(maux,mitot,mjtot) + dimension auxc1d(maux,lenbc) + +c +c ::::::::::::::::::::::::::: QAD :::::::::::::::::::::::::::::::::: +c solve RP between ghost cell value on fine grid and coarse grid +c value that ghost cell overlaps. The resulting fluctuations +c are added in to coarse grid value, as a conservation fixup. +c Done each fine grid time step. If source terms are present, the +c coarse grid value is advanced by source terms each fine time step too. + +c No change needed in this sub. for spherical mapping: correctly +c mapped vals already in bcs on this fine grid and coarse saved +c vals also properly prepared +c +c Side 1 is the left side of the fine grid patch. Then go around clockwise. +c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: +c +c # local storage +c # note that dimension here are bigger than dimensions used +c # in rp2, but shouldn't matter since wave is not used in qad +c # and for other arrays it is only the last parameter that is wrong +c # ok as long as meqn, mwaves < maxvar + + parameter (max1dp1 = max1d+1) + dimension ql(max1dp1,nvar), qr(max1dp1,nvar) + dimension wave(max1dp1,nvar,mwaves), s(max1dp1,mwaves) + dimension amdq(max1dp1,nvar), apdq(max1dp1,nvar) + dimension auxl(maxaux*max1dp1), auxr(maxaux*max1dp1) +c +c WARNING: auxl,auxr dimensioned at max possible, but used as if +c they were dimensioned as the real maux by max1dp1. Would be better +c of course to dimension by maux by max1dp1 but this wont work if maux=0 +c So need to access using your own indexing into auxl,auxr. + iaddaux(iaux,i) = i + max1dp1*(iaux-1) + + data qprint/.false./ +c +c aux is auxiliary array with user parameters needed in Riemann solvers +c on fine grid corresponding to valbig +c auxc1d is coarse grid stuff from around boundary, same format as qc1d +c auxl, auxr are work arrays needed to pass stuff to rpn2 +c maux is the number of aux variables, which may be zero. + + +c --------------------------------------------------------------- +c # for GeoClaw we skip this routine --- need to fix someday: + +c # Is this routine properly implemented with dry states and +c # interpolation based on surface eta? + +c return + +c --------------------------------------------------------------- + + + if (qprint) write(dbugunit,*)" working on grid ",mptr + tgrid = rnode(timemult, mptr) + nc = mjtot-2*nghost + nr = mitot-2*nghost + level = node(nestlevel, mptr) + index = 0 +c +c-------- +c side 1 +c-------- +c + do 10 j = nghost+1, mjtot-nghost + if (maux.gt.0) then + do 5 ma = 1,maux + if (auxtype(ma).eq."xleft") then +c # Assuming velocity at left-face, this fix +c # preserves conservation in incompressible flow: + auxl(iaddaux(ma,j-nghost+1)) = aux(ma,nghost+1,j) + else +c # Normal case -- we set the aux arrays +c # from the cell corresponding to q + auxl(iaddaux(ma,j-nghost+1)) = aux(ma,nghost,j) + endif + 5 continue + endif + do 10 ivar = 1, nvar + ql(j-nghost+1,ivar) = valbig(ivar,nghost,j) + 10 continue + + lind = 0 + ncrse = (mjtot-2*nghost)/lratioy + do 20 jc = 1, ncrse + index = index + 1 + do 25 l = 1, lratioy + lind = lind + 1 + if (maux.gt.0) then + do 24 ma=1,maux + auxr(iaddaux(ma,lind)) = auxc1d(ma,index) + 24 continue + endif + do 25 ivar = 1, nvar + 25 qr(lind,ivar) = qc1d(ivar,index) + 20 continue + index = ncrse + + if (qprint) then + write(dbugunit,*) 'side 1, ql and qr:' + do i=2,nc + write(dbugunit,4101) i,qr(i-1,1),ql(i,1) + enddo + 4101 format(i3,4e16.6) + if (maux .gt. 0) then + write(dbugunit,*) 'side 1, auxr:' + do i=2,nc + write(dbugunit,4101) i,(auxr(iaddaux(ma,i-1)),ma=1,maux) + enddo + write(dbugunit,*) 'side 1, auxl:' + do i=2,nc + write(dbugunit,4101) i,(auxl(iaddaux(ma,i)),ma=1,maux) + enddo + endif + endif + + call rpn2(1,max1dp1-2*nghost,nvar,mwaves,maux,nghost, + . nc+1-2*nghost,ql,qr,auxl,auxr,wave,s,amdq,apdq) +c +c we have the wave. for side 1 add into sdflxm +c + influx = 0 + do 30 j = 1, nc/lratioy + influx = influx + 1 + jfine = (j-1)*lratioy + do 40 ivar = 1, nvar + do 50 l = 1, lratioy + svdflx(ivar,influx) = svdflx(ivar,influx) + . + amdq(jfine+l+1,ivar) * hy * delt + . + apdq(jfine+l+1,ivar) * hy * delt + 50 continue + 40 continue + 30 continue + +c-------- +c side 2 +c-------- +c + if (mjtot .eq. 2*nghost+1) then +c # a single row of interior cells only happens when using the +c # 2d amrclaw code to do a 1d problem with refinement. +c # (feature added in Version 4.3) +c # skip over sides 2 and 4 in this case + go to 299 + endif + + do 210 i = nghost+1, mitot-nghost + if (maux.gt.0) then + do 205 ma = 1,maux + auxr(iaddaux(ma,i-nghost)) = aux(ma,i,mjtot-nghost+1) + 205 continue + endif + do 210 ivar = 1, nvar + qr(i-nghost,ivar) = valbig(ivar,i,mjtot-nghost+1) + 210 continue + + lind = 0 + ncrse = (mitot-2*nghost)/lratiox + do 220 ic = 1, ncrse + index = index + 1 + do 225 l = 1, lratiox + lind = lind + 1 + if (maux.gt.0) then + do 224 ma=1,maux + if (auxtype(ma).eq."yleft") then +c # Assuming velocity at bottom-face, this fix +c # preserves conservation in incompressible flow: + ifine = (ic-1)*lratiox + nghost + l + auxl(iaddaux(ma,lind+1)) = aux(ma,ifine,mjtot-nghost+1) + else + auxl(iaddaux(ma,lind+1)) = auxc1d(ma,index) + endif + 224 continue + endif + do 225 ivar = 1, nvar + 225 ql(lind+1,ivar) = qc1d(ivar,index) + 220 continue + + if (qprint) then + write(dbugunit,*) 'side 2, ql and qr:' + do i=1,nr + write(dbugunit,4101) i,ql(i+1,1),qr(i,1) + enddo + if (maux .gt. 0) then + write(dbugunit,*) 'side 2, auxr:' + do i = 1, nr + write(dbugunit,4101) i, (auxr(iaddaux(ma,i)),ma=1,maux) + enddo + write(dbugunit,*) 'side 2, auxl:' + do i = 1, nr + write(dbugunit,4101) i, (auxl(iaddaux(ma,i)),ma=1,maux) + enddo + endif + endif + call rpn2(2,max1dp1-2*nghost,nvar,mwaves,maux,nghost, + . nr+1-2*nghost,ql,qr,auxl,auxr,wave,s,amdq,apdq) +c +c we have the wave. for side 2. add into sdflxp +c + do 230 i = 1, nr/lratiox + influx = influx + 1 + ifine = (i-1)*lratiox + do 240 ivar = 1, nvar + do 250 l = 1, lratiox + svdflx(ivar,influx) = svdflx(ivar,influx) + . - amdq(ifine+l+1,ivar) * hx * delt + . - apdq(ifine+l+1,ivar) * hx * delt + 250 continue + 240 continue + 230 continue + + 299 continue + +c-------- +c side 3 +c-------- +c + do 310 j = nghost+1, mjtot-nghost + if (maux.gt.0) then + do 305 ma = 1,maux + auxr(iaddaux(ma,j-nghost)) = aux(ma,mitot-nghost+1,j) + 305 continue + endif + do 310 ivar = 1, nvar + qr(j-nghost,ivar) = valbig(ivar,mitot-nghost+1,j) + 310 continue + + lind = 0 + ncrse = (mjtot-2*nghost)/lratioy + do 320 jc = 1, ncrse + index = index + 1 + do 325 l = 1, lratioy + lind = lind + 1 + if (maux.gt.0) then + do 324 ma=1,maux + if (auxtype(ma).eq."xleft") then +c # Assuming velocity at left-face, this fix +c # preserves conservation in incompressible flow: + jfine = (jc-1)*lratioy + nghost + l + auxl(iaddaux(ma,lind+1)) = aux(ma,mitot-nghost+1,jfine) + else + auxl(iaddaux(ma,lind+1)) = auxc1d(ma,index) + endif + 324 continue + endif + do 325 ivar = 1, nvar + 325 ql(lind+1,ivar) = qc1d(ivar,index) + 320 continue + + if (qprint) then + write(dbugunit,*) 'side 3, ql and qr:' + do i=1,nc + write(dbugunit,4101) i,ql(i,1),qr(i,1) + enddo + endif + call rpn2(1,max1dp1-2*nghost,nvar,mwaves,maux,nghost, + . nc+1-2*nghost,ql,qr,auxl,auxr,wave,s,amdq,apdq) +c +c we have the wave. for side 3 add into sdflxp +C + do 330 j = 1, nc/lratioy + influx = influx + 1 + jfine = (j-1)*lratioy + do 340 ivar = 1, nvar + do 350 l = 1, lratioy + svdflx(ivar,influx) = svdflx(ivar,influx) + . - amdq(jfine+l+1,ivar) * hy * delt + . - apdq(jfine+l+1,ivar) * hy * delt + 350 continue + 340 continue + 330 continue + +c-------- +c side 4 +c-------- +c + if (mjtot .eq. 2*nghost+1) then +c # a single row of interior cells only happens when using the +c # 2d amrclaw code to do a 1d problem with refinement. +c # (feature added in Version 4.3) +c # skip over sides 2 and 4 in this case + go to 499 + endif +c + do 410 i = nghost+1, mitot-nghost + if (maux.gt.0) then + do 405 ma = 1,maux + if (auxtype(ma).eq."yleft") then +c # Assuming velocity at bottom-face, this fix +c # preserves conservation in incompressible flow: + auxl(iaddaux(ma,i-nghost+1)) = aux(ma,i,nghost+1) + else + auxl(iaddaux(ma,i-nghost+1)) = aux(ma,i,nghost) + endif + 405 continue + endif + do 410 ivar = 1, nvar + ql(i-nghost+1,ivar) = valbig(ivar,i,nghost) + 410 continue + + lind = 0 + ncrse = (mitot-2*nghost)/lratiox + do 420 ic = 1, ncrse + index = index + 1 + do 425 l = 1, lratiox + lind = lind + 1 + if (maux.gt.0) then + do 424 ma=1,maux + auxr(iaddaux(ma,lind)) = auxc1d(ma,index) + 424 continue + endif + do 425 ivar = 1, nvar + 425 qr(lind,ivar) = qc1d(ivar,index) + 420 continue + + if (qprint) then + write(dbugunit,*) 'side 4, ql and qr:' + do i=1,nr + write(dbugunit,4101) i, ql(i,1),qr(i,1) + enddo + endif + call rpn2(2,max1dp1-2*nghost,nvar,mwaves,maux,nghost, + . nr+1-2*nghost,ql,qr,auxl,auxr,wave,s,amdq,apdq) +c +c we have the wave. for side 4. add into sdflxm +c + do 430 i = 1, nr/lratiox + influx = influx + 1 + ifine = (i-1)*lratiox + do 440 ivar = 1, nvar + do 450 l = 1, lratiox + svdflx(ivar,influx) = svdflx(ivar,influx) + . + amdq(ifine+l+1,ivar) * hx * delt + . + apdq(ifine+l+1,ivar) * hx * delt + 450 continue + 440 continue + 430 continue + + 499 continue + +c # for source terms: +c # This is commented out for now --- updating the momenumtum due +c # to friction/Coriolis source terms affects the resulting mass flux +c # and leads to loss of conservation when src1d is used. +c # It would be better to incorporate these source terms into the +c # the f-wave formulation? +c if (method(5) .ne. 0) then +c call src1d(nvar,nghost,lenbc,qc1d,maux,auxc1d,tgrid,delt) +c endif + + return + end diff --git a/src/2d/shallow/step2_vec.f90 b/src/2d/shallow/step2_vec.f90 new file mode 100644 index 000000000..5b44e45f7 --- /dev/null +++ b/src/2d/shallow/step2_vec.f90 @@ -0,0 +1,235 @@ +! This code has been adapted for vectorization. +! For more details on how to use it, see the Makefile of the chile2010 example. + +subroutine step2(maxm,meqn,maux,mbc,mx,my, & + qold,aux,dx,dy,dt,cflgrid, & + fm,fp,gm,gp,rpn2,rpt2) +! ========================================================== +! +! # clawpack routine ... modified for AMRCLAW +! +! # Take one time step, updating q. +! # On entry, qold gives +! # initial data for this step +! # and is unchanged in this version. +! +! # fm, fp are fluxes to left and right of single cell edge +! # See the flux2 documentation for more information. +! +! # modified again for GeoClaw +!------------------step2_geo.f----------------------- +! This version of step2 relimits the fluxes in order to +! maintain positivity. +! to do so set relimit=.true. +! The only modification is the 101 loop +!------------last modified 12/30/04-------------------------- +! + + use geoclaw_module, only: dry_tolerance + use amr_module, only: mwaves, mcapa + + implicit none + + external rpn2, rpt2 + + ! Arguments + integer, intent(in) :: maxm,meqn,maux,mbc,mx,my + real(kind=8), intent(in) :: dx,dy,dt + real(kind=8), intent(inout) :: cflgrid + real(kind=8), intent(inout) :: qold(meqn, 1-mbc:mx+mbc, 1-mbc:my+mbc) + real(kind=8), intent(inout) :: aux(maux,1-mbc:mx+mbc, 1-mbc:my+mbc) + real(kind=8), intent(inout) :: fm(meqn, 1-mbc:mx+mbc, 1-mbc:my+mbc) + real(kind=8), intent(inout) :: fp(meqn,1-mbc:mx+mbc, 1-mbc:my+mbc) + real(kind=8), intent(inout) :: gm(meqn,1-mbc:mx+mbc, 1-mbc:my+mbc) + real(kind=8), intent(inout) :: gp(meqn,1-mbc:mx+mbc, 1-mbc:my+mbc) + + ! Local storage for flux accumulation + real(kind=8) :: faddm(1-mbc:maxm+mbc,meqn) + real(kind=8) :: faddp(1-mbc:maxm+mbc,meqn) + real(kind=8) :: gaddm(1-mbc:maxm+mbc,meqn,2) + real(kind=8) :: gaddp(1-mbc:maxm+mbc,meqn,2) + + ! Scratch storage for Sweeps and Riemann problems + real(kind=8) :: q1d(1-mbc:maxm+mbc,meqn) + real(kind=8) :: aux1(1-mbc:maxm+mbc,maux) + real(kind=8) :: aux2(1-mbc:maxm+mbc,maux) + real(kind=8) :: aux3(1-mbc:maxm+mbc,maux) + real(kind=8) :: dtdx1d(1-mbc:maxm+mbc) + real(kind=8) :: dtdy1d(1-mbc:maxm+mbc) + + ! real(kind=8) :: wave(meqn, mwaves, 1-mbc:maxm+mbc) + ! real(kind=8) :: s(mwaves, 1-mbc:maxm + mbc) + ! real(kind=8) :: amdq(meqn,1-mbc:maxm + mbc) + ! real(kind=8) :: apdq(meqn,1-mbc:maxm + mbc) +! real(kind=8) :: cqxx(meqn,1-mbc:maxm + mbc) +! real(kind=8) :: bmadq(meqn,1-mbc:maxm + mbc) +! real(kind=8) :: bpadq(meqn,1-mbc:maxm + mbc) + + ! Looping scalar storage + integer :: i,j,m,thread_num + real(kind=8) :: dtdx,dtdy,cfl1d,p,phi,cm,dtdxij,dtdyij + + ! Common block storage + integer :: icom,jcom + + ! Parameters + ! Relimit fluxes to maintain positivity + logical, parameter :: relimit = .false. + + cflgrid = 0.d0 + dtdx = dt/dx + dtdy = dt/dy + + fm = 0.d0 + fp = 0.d0 + gm = 0.d0 + gp = 0.d0 + + ! ========================================================================== + ! Perform X-Sweeps + do j = 0,my+1 + + ! Copy old q into 1d slice + do i = 1-mbc,mx+mbc + do m = 1,meqn + q1d(i,m) = qold(m,i,j) + enddo + enddo + + ! Set dtdx slice if a capacity array exists + if (mcapa > 0) then + dtdx1d(1-mbc:mx+mbc) = dtdx / aux(mcapa,1-mbc:mx+mbc,j) + else + dtdx1d = dtdx + endif + + ! Copy aux array into slices + if (maux > 0) then + do i = 1-mbc,mx+mbc + do m = 1,maux + aux1(i,m) = aux(m,i,j-1) + aux2(i,m) = aux(m,i,j ) + aux3(i,m) = aux(m,i,j+1) + enddo + enddo + endif + + ! Store value of j along the slice into common block + ! *** WARNING *** This may not working with threading + jcom = j + + ! Compute modifications fadd and gadd to fluxes along this slice: + call flux2(1,maxm,meqn,maux,mbc,mx,q1d,dtdx1d,aux1,aux2,aux3, & + faddm,faddp,gaddm,gaddp,cfl1d,rpn2,rpt2) + + cflgrid = max(cflgrid,cfl1d) + ! write(53,*) 'x-sweep: ',cfl1d,cflgrid + + ! Update fluxes + do i=1,meqn + fm(i,1:mx+1,j) = fm(i,1:mx+1,j) + faddm(1:mx+1,i) + fp(i,1:mx+1,j) = fp(i,1:mx+1,j) + faddp(1:mx+1,i) + gm(i,1:mx+1,j) = gm(i,1:mx+1,j) + gaddm(1:mx+1,i,1) + gp(i,1:mx+1,j) = gp(i,1:mx+1,j) + gaddp(1:mx+1,i,1) + gm(i,1:mx+1,j+1) = gm(i,1:mx+1,j+1) + gaddm(1:mx+1,i,2) + gp(i,1:mx+1,j+1) = gp(i,1:mx+1,j+1) + gaddp(1:mx+1,i,2) + end do + + enddo + + ! ============================================================================ + ! y-sweeps + ! + do i = 0, mx+1 + + ! Copy data along a slice into 1d arrays: + do j = 1-mbc, my+mbc + do m=1, meqn + q1d(j,m) = qold(m,i,j) + enddo + enddo + + ! Set dt/dy ratio in slice + if (mcapa > 0) then + dtdy1d(1-mbc:my+mbc) = dtdy / aux(mcapa,i,1-mbc:my+mbc) + else + dtdy1d = dtdy + endif + + ! Copy aux slices + if (maux .gt. 0) then + do j = 1-mbc, my+mbc + do m=1, maux + aux1(j,m) = aux(m,i-1,j) + aux2(j,m) = aux(m,i ,j) + aux3(j,m) = aux(m,i+1,j) + enddo + enddo + endif + + ! Store the value of i along this slice in the common block + ! *** WARNING *** This may not working with threading + icom = i + + ! Compute modifications fadd and gadd to fluxes along this slice + call flux2(2,maxm,meqn,maux,mbc,my,q1d,dtdy1d,aux1,aux2,aux3, & + faddm,faddp,gaddm,gaddp,cfl1d,rpn2,rpt2) + + cflgrid = max(cflgrid,cfl1d) + ! write(53,*) 'y-sweep: ',cfl1d,cflgrid + + ! Update fluxes + do j=1,meqn + gm(j,i,1:my+1) = gm(j,i,1:my+1) + faddm(1:my+1,j) + gp(j,i,1:my+1) = gp(j,i,1:my+1) + faddp(1:my+1,j) + fm(j,i,1:my+1) = fm(j,i,1:my+1) + gaddm(1:my+1,j,1) + fp(j,i,1:my+1) = fp(j,i,1:my+1) + gaddp(1:my+1,j,1) + fm(j,i+1,1:my+1) = fm(j,i+1,1:my+1) + gaddm(1:my+1,j,2) + fp(j,i+1,1:my+1) = fp(j,i+1,1:my+1) + gaddp(1:my+1,j,2) + end do + + end do + + ! Relimit correction fluxes if they drive a cell negative + if (relimit) then + dtdxij = dtdx + dtdyij = dtdy + do i=1,mx + do j=1,my + if (mcapa > 0) then + dtdxij = dtdx / aux(mcapa,i,j) + dtdyij = dtdy / aux(mcapa,i,j) + endif + p = max(0.d0,dtdxij*fm(1,i+1,j)) + max(0.d0,dtdyij*gm(1,i,j+1)) & + - min(0.d0,dtdxij*fp(1,i,j)) - min(0.d0,dtdyij*gp(1,i,j)) + phi = min(1.d0,abs(qold(1,i,j) / (p+dry_tolerance))) + + if (phi < 1.d0) then + do m=1,meqn + if (fp(1,i,j) < 0.d0) then + cm = fp(m,i,j) - fm(m,i,j) + fm(m,i,j) = phi * fm(m,i,j) + fp(m,i,j) = fm(m,i,j) + cm + endif + if (gp(1,i,j) < 0.d0) then + cm = gp(m,i,j) - gm(m,i,j) + gm(m,i,j) = phi * gm(m,i,j) + gp(m,i,j) = gm(m,i,j) + cm + endif + if (fm(1,i+1,j) > 0.d0) then + cm = fp(m,i+1,j) - fm(m,i+1,j) + fp(m,i+1,j) = phi * fp(m,i+1,j) + fm(m,i+1,j) = fp(m,i+1,j) - cm + endif + if (gm(1,i,j+1) > 0.d0) then + cm = gp(m,i,j+1) - gm(m,i,j+1) + gp(m,i,j+1) = phi * gp(m,i,j+1) + gm(m,i,j+1) = gp(m,i,j+1) - cm + endif + end do + endif + enddo + enddo + endif + +end subroutine step2 From ddf66fffb2342b7f20b245b67bc602535bd0feea Mon Sep 17 00:00:00 2001 From: Chaulio Ferreira Date: Wed, 11 Apr 2018 13:37:31 +0200 Subject: [PATCH 2/3] flux2fw_vec.f: replaced !DIR$ SIMD directives with !$OMP SIMD --- src/2d/shallow/flux2fw_vec.f | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/2d/shallow/flux2fw_vec.f b/src/2d/shallow/flux2fw_vec.f index 6b27e5d58..ab61405d3 100644 --- a/src/2d/shallow/flux2fw_vec.f +++ b/src/2d/shallow/flux2fw_vec.f @@ -145,7 +145,7 @@ subroutine flux2(ixy,maxm,meqn,maux,mbc,mx, c # Set fadd for the donor-cell upwind method (Godunov) if (ixy.eq.1) mu=2 if (ixy.eq.2) mu=3 - !DIR$ SIMD PRIVATE(dxdc) + !$OMP SIMD PRIVATE(dxdc) do 40 i=1-mbc+1,mx+mbc-1 if (coordinate_system.eq.2) then if (ixy.eq.1) then @@ -171,7 +171,7 @@ subroutine flux2(ixy,maxm,meqn,maux,mbc,mx, c # compute maximum wave speed for checking Courant number: cfl1d = 0.d0 do 50 mw=1,mwaves - !DIR$ SIMD + !$OMP SIMD do 50 i=1,mx+1 c # if s>0 use dtdx1d(i) to compute CFL, c # if s<0 use dtdx1d(i-1) to compute CFL: @@ -188,7 +188,7 @@ subroutine flux2(ixy,maxm,meqn,maux,mbc,mx, c # apply limiter to fwaves: if (limit) call limiter(maxm,meqn,mwaves,mbc,mx,fwave,s,mthlim) c - !DIR$ SIMD + !$OMP SIMD do 120 i = 1, mx+1 c c # For correction terms below, need average of dtdx in cell @@ -216,7 +216,7 @@ subroutine flux2(ixy,maxm,meqn,maux,mbc,mx, c if (method(2).gt.1 .and. method(3).eq.2) then c # incorporate cqxx into amdq and apdq so that it is split also. - !DIR$ SIMD + !$OMP SIMD do 150 i = 1, mx+1 do 150 m=1,meqn amdq(i,m) = amdq(i,m) + cqxx(i,m) @@ -235,7 +235,7 @@ subroutine flux2(ixy,maxm,meqn,maux,mbc,mx, & amdq,bmasdq,bpasdq) c c # modify flux below and above by B^- A^- Delta q and B^+ A^- Delta q: - !DIR$ SIMD PRIVATE(gupdate) + !$OMP SIMD PRIVATE(gupdate) do 160 i = 1, mx+1 do 160 m=1,meqn gupdate = 0.5d0*dtdx1d(i-1) * bmasdq(i,m) @@ -253,7 +253,7 @@ subroutine flux2(ixy,maxm,meqn,maux,mbc,mx, & apdq,bmasdq,bpasdq) c c # modify flux below and above by B^- A^+ Delta q and B^+ A^+ Delta q: - !DIR$ SIMD PRIVATE(gupdate) + !$OMP SIMD PRIVATE(gupdate) do 180 i = 1, mx+1 do 180 m=1,meqn gupdate = 0.5d0*dtdx1d(i-1) * bmasdq(i,m) From 4f07e98cdf0d8c568749715dbd5c7861f13eec74 Mon Sep 17 00:00:00 2001 From: Chaulio Ferreira Date: Thu, 17 May 2018 18:09:07 +0200 Subject: [PATCH 3/3] Created test for vectorized Riemann solvers on a separate directory (bowl_slosh_vec) --- tests/bowl_slosh_vec/Makefile | 90 ++++ tests/bowl_slosh_vec/__init__.py | 0 tests/bowl_slosh_vec/make_fgmax_grid.py | 33 ++ tests/bowl_slosh_vec/maketopo.py | 45 ++ tests/bowl_slosh_vec/plot_fgmax.py | 38 ++ tests/bowl_slosh_vec/qinit.f90 | 37 ++ .../regression_data/claw_git_status.txt | 94 ++++ .../regression_data/gauge00001.txt | 57 +++ .../regression_data/regression_data_fgmax.txt | 2 + tests/bowl_slosh_vec/regression_tests.py | 75 ++++ tests/bowl_slosh_vec/setplot.py | 177 ++++++++ tests/bowl_slosh_vec/setrun.py | 404 ++++++++++++++++++ 12 files changed, 1052 insertions(+) create mode 100644 tests/bowl_slosh_vec/Makefile create mode 100644 tests/bowl_slosh_vec/__init__.py create mode 100644 tests/bowl_slosh_vec/make_fgmax_grid.py create mode 100644 tests/bowl_slosh_vec/maketopo.py create mode 100644 tests/bowl_slosh_vec/plot_fgmax.py create mode 100644 tests/bowl_slosh_vec/qinit.f90 create mode 100644 tests/bowl_slosh_vec/regression_data/claw_git_status.txt create mode 100644 tests/bowl_slosh_vec/regression_data/gauge00001.txt create mode 100644 tests/bowl_slosh_vec/regression_data/regression_data_fgmax.txt create mode 100644 tests/bowl_slosh_vec/regression_tests.py create mode 100644 tests/bowl_slosh_vec/setplot.py create mode 100644 tests/bowl_slosh_vec/setrun.py diff --git a/tests/bowl_slosh_vec/Makefile b/tests/bowl_slosh_vec/Makefile new file mode 100644 index 000000000..2fdd92b2b --- /dev/null +++ b/tests/bowl_slosh_vec/Makefile @@ -0,0 +1,90 @@ +# 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 = geoclaw # Clawpack package to use +EXE = xgeoclaw # 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 + +RESTART = False + +# 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 ?= + +FFLAGS += -O2 -fopenmp + +ifeq ($(FC), ifort) + FFLAGS += -align array64byte -xHost +else + $(warning Warning: For vectorization the Intel Compiler is recommended) +endif + +# --------------------------------- +# package sources for this program: +# --------------------------------- + +GEOLIB = $(CLAW)/geoclaw/src/2d/shallow +include $(GEOLIB)/Makefile.geoclaw + +# --------------------------------------- +# package sources specifically to exclude +# (i.e. if a custom replacement source +# under a different name is provided) +# --------------------------------------- + +EXCLUDE_MODULES = \ + +EXCLUDE_SOURCES = \ + $(GEOLIB)/qad.f \ + $(GEOLIB)/flux2fw.f \ + $(GEOLIB)/step2.f \ + $(AMRLIB)/inlinelimiter.f \ + $(GEOLIB)/fgmax_interpolate0.f90 \ + +# ---------------------------------------- +# List of custom sources for this program: +# ---------------------------------------- + +MODULES = \ + +# This is where this test differs mostly from the original bowl-slosh one: the source files + +SOURCES = \ + qinit.f90 \ + $(GEOLIB)/fgmax_interpolate.f90 \ + $(CLAW)/riemann/src/rpn2_geoclaw_vec.f90 \ + $(CLAW)/riemann/src/rpt2_geoclaw_vec.f90 \ + $(CLAW)/riemann/src/geoclaw_riemann_utils_vec.f90 \ + $(GEOLIB)/qad_vec.f \ + $(GEOLIB)/flux2fw_vec.f \ + $(GEOLIB)/step2_vec.f90 \ + $(GEOLIB)/inlinelimiter_vec.f + +#------------------------------------------------------------------- +# Include Makefile containing standard definitions and make options: +include $(CLAWMAKE) + +# Construct the topography data +.PHONY: topo all +topo: + python maketopo.py + +all: + $(MAKE) topo + $(MAKE) .plots + $(MAKE) .htmls + diff --git a/tests/bowl_slosh_vec/__init__.py b/tests/bowl_slosh_vec/__init__.py new file mode 100644 index 000000000..e69de29bb diff --git a/tests/bowl_slosh_vec/make_fgmax_grid.py b/tests/bowl_slosh_vec/make_fgmax_grid.py new file mode 100644 index 000000000..f814aa756 --- /dev/null +++ b/tests/bowl_slosh_vec/make_fgmax_grid.py @@ -0,0 +1,33 @@ + +""" +Create fgmax_grid.txt input file +""" + +from __future__ import absolute_import +from clawpack.geoclaw import fgmax_tools +import os + + +def make_fgmax_grid1(datadir): + fg = fgmax_tools.FGmaxGrid() + fg.point_style = 2 # will specify a 2d grid of points + fg.x1 = -2. + fg.x2 = 2. + fg.y1 = -2. + fg.y2 = 2. + fg.dx = 0.1 + fg.tstart_max = 0. # when to start monitoring max values + fg.tend_max = 1.e10 # when to stop monitoring max values + fg.dt_check = 0.1 # target time (sec) increment between updating + # max values + fg.min_level_check = 2 # which levels to monitor max on + fg.arrival_tol = 1.e-2 # tolerance for flagging arrival + + fg.input_file_name = os.path.join(datadir, 'fgmax1.txt') + fg.write_input_data() + + +if __name__ == "__main__": + make_fgmax_grid1('.') + + diff --git a/tests/bowl_slosh_vec/maketopo.py b/tests/bowl_slosh_vec/maketopo.py new file mode 100644 index 000000000..a3f71dff7 --- /dev/null +++ b/tests/bowl_slosh_vec/maketopo.py @@ -0,0 +1,45 @@ + +""" +Module to create topo and qinit data files for this example. +""" + +from __future__ import absolute_import +from clawpack.geoclaw.topotools import Topography + +from numpy import * + +#from pyclaw.data import Data +#probdata = Data('setprob.data') + +a = 1. +sigma = 0.5 +h0 = 0.1 +grav = 9.81 +omega = sqrt(2.*grav*h0) / a + +def maketopo(): + """ + Output topography file for the entire domain + """ + nxpoints=200 + nypoints=200 + xupper=2.e0 + yupper=2.e0 + xlower = -2.e0 + ylower = -2.e0 + outfile= "bowl.topotype2" + topography = Topography(topo_func=topo) + topography.x = linspace(xlower,xupper,nxpoints) + topography.y = linspace(ylower,yupper,nypoints) + topography.write(outfile, topo_type=2, Z_format="%22.15e") + +def topo(x,y): + """ + Parabolic bowl + """ + z = h0*(x**2 + y**2)/a**2 - h0 + return z + + +if __name__=='__main__': + maketopo() diff --git a/tests/bowl_slosh_vec/plot_fgmax.py b/tests/bowl_slosh_vec/plot_fgmax.py new file mode 100644 index 000000000..aa4431f4f --- /dev/null +++ b/tests/bowl_slosh_vec/plot_fgmax.py @@ -0,0 +1,38 @@ + +""" +Plot fgmax output from GeoClaw run. + +""" + +from __future__ import absolute_import +from pylab import * +from numpy import ma +from clawpack.geoclaw import fgmax_tools + +fg = fgmax_tools.FGmaxGrid() +fg.read_input_data('fgmax1.txt') +fg.read_output() + +figure(1) +clf() +surface = ma.masked_where(fg.h < 0.001, fg.h + fg.B) +contourf(fg.X,fg.Y,surface,10) +cb = colorbar() +cb.set_label('meters') +title('Max surface elevation') + +figure(2) +clf() +#s = ma.masked_where(fg.s<-1e10, fg.s) +contourf(fg.X,fg.Y,fg.s,10) +cb = colorbar() +cb.set_label('meters / sec') +title('Max speed') + +figure(3) +clf() +contourf(fg.X,fg.Y,fg.arrival_time,10) +cb = colorbar() +cb.set_label('seconds') +title('Arrival time') + diff --git a/tests/bowl_slosh_vec/qinit.f90 b/tests/bowl_slosh_vec/qinit.f90 new file mode 100644 index 000000000..3fd1d902e --- /dev/null +++ b/tests/bowl_slosh_vec/qinit.f90 @@ -0,0 +1,37 @@ +! qinit routine for parabolic bowl problem, only single layer +subroutine qinit(meqn,mbc,mx,my,xlower,ylower,dx,dy,q,maux,aux) + + use geoclaw_module, only: grav + + implicit none + + ! Subroutine arguments + integer, intent(in) :: meqn,mbc,mx,my,maux + real(kind=8), intent(in) :: xlower,ylower,dx,dy + real(kind=8), intent(inout) :: q(meqn,1-mbc:mx+mbc,1-mbc:my+mbc) + real(kind=8), intent(inout) :: aux(maux,1-mbc:mx+mbc,1-mbc:my+mbc) + + ! Parameters for problem + real(kind=8), parameter :: a = 1.d0 + real(kind=8), parameter :: sigma = 0.5d0 + real(kind=8), parameter :: h0 = 0.1d0 + + ! Other storage + integer :: i,j + real(kind=8) :: omega,x,y,eta + + omega = sqrt(2.d0 * grav * h0) / a + + do i=1-mbc,mx+mbc + x = xlower + (i - 0.5d0)*dx + do j=1-mbc,my+mbc + y = ylower + (j - 0.5d0) * dx + eta = sigma * h0 / a**2 * (2.d0 * x - sigma) + + q(1,i,j) = max(0.d0,eta - aux(1,i,j)) + q(2,i,j) = 0.d0 + q(3,i,j) = sigma * omega * q(1,i,j) + enddo + enddo + +end subroutine qinit diff --git a/tests/bowl_slosh_vec/regression_data/claw_git_status.txt b/tests/bowl_slosh_vec/regression_data/claw_git_status.txt new file mode 100644 index 000000000..47d020a1b --- /dev/null +++ b/tests/bowl_slosh_vec/regression_data/claw_git_status.txt @@ -0,0 +1,94 @@ +Clawpack Git Status +Diffs can be found in /Users/rjl/git/clawpack/geoclaw/tests/bowl_slosh/regression_data/claw_git_diffs.txt + +Tue, 20 Dec 2016 22:18:41 PST +$CLAW = /Users/rjl/git/clawpack +$FC = gfortran + + +=========== +classic +=========== +/Users/rjl/git/clawpack/classic + +--- last commit --- +8e6d0c2 make plotdata an optional argument in setplot + +--- branch and status --- +## setplot_fixes + + +=========== +amrclaw +=========== +/Users/rjl/git/clawpack/amrclaw + +--- last commit --- +5c81870 set xlimits and ylimits in quadrants example + +--- branch and status --- +## setplot_fixes + + +=========== +clawutil +=========== +/Users/rjl/git/clawpack/clawutil + +--- last commit --- +c3950d1 Fix call to StringIO to match new way it is imported + +--- branch and status --- +## claw_git_status_StringIO + + +=========== +pyclaw +=========== +/Users/rjl/git/clawpack/pyclaw + +--- last commit --- +695af06 Merge pull request #550 from clawpack/python3 + +--- branch and status --- +## master...origin/master + + +=========== +visclaw +=========== +/Users/rjl/git/clawpack/visclaw + +--- last commit --- +6833be6 Merge pull request #190 from clawpack/python3 + +--- branch and status --- +## master...origin/master + + +=========== +riemann +=========== +/Users/rjl/git/clawpack/riemann + +--- last commit --- +b68d73d Merge pull request #111 from rjleveque/geoclaw_riemann_utils_update + +--- branch and status --- +## master...origin/master + + +=========== +geoclaw +=========== +/Users/rjl/git/clawpack/geoclaw + +--- last commit --- +824ca39 fix regression data in tests/multilayer and remove claw_git_diffs.txt: full or recursive garbage + +--- branch and status --- +## regression_data_update_riemann + M examples/tsunami/bowl-radial/setplot.py + M tests/bowl_slosh/regression_data/claw_git_status.txt + M tests/bowl_slosh/regression_data/gauge00001.txt + M tests/bowl_slosh/regression_tests.py diff --git a/tests/bowl_slosh_vec/regression_data/gauge00001.txt b/tests/bowl_slosh_vec/regression_data/gauge00001.txt new file mode 100644 index 000000000..2383062d8 --- /dev/null +++ b/tests/bowl_slosh_vec/regression_data/gauge00001.txt @@ -0,0 +1,57 @@ +# gauge_id= 1 location=( 0.5000000E+00 0.5000000E+00 ) num_var= 4 + # level, time, q[ 1 2 3], eta, aux[] + 02 0.0000000E+00 0.7497749E-01 0.0000000E+00 0.5251101E-01 0.2500000E-01 + 02 0.1000000E-03 0.7498449E-01 -0.7355608E-05 0.5251591E-01 0.2500700E-01 + 02 0.8759850E-02 0.7558763E-01 -0.6468797E-03 0.5293694E-01 0.2561014E-01 + 02 0.1741970E-01 0.7618225E-01 -0.1296497E-02 0.5334473E-01 0.2620476E-01 + 02 0.2607955E-01 0.7676870E-01 -0.1955904E-02 0.5373907E-01 0.2679122E-01 + 02 0.3473940E-01 0.7734737E-01 -0.2624783E-02 0.5411970E-01 0.2736988E-01 + 02 0.4339925E-01 0.7791933E-01 -0.3302832E-02 0.5448613E-01 0.2794184E-01 + 02 0.5205910E-01 0.7848445E-01 -0.3989757E-02 0.5483802E-01 0.2850697E-01 + 02 0.6071895E-01 0.7904271E-01 -0.4685290E-02 0.5517516E-01 0.2906522E-01 + 02 0.6937880E-01 0.7959405E-01 -0.5389169E-02 0.5549738E-01 0.2961656E-01 + 02 0.7803865E-01 0.8013854E-01 -0.6101133E-02 0.5580453E-01 0.3016105E-01 + 02 0.8669850E-01 0.8067617E-01 -0.6820904E-02 0.5609644E-01 0.3069868E-01 + 02 0.9753442E-01 0.8133911E-01 -0.7731392E-02 0.5644092E-01 0.3136163E-01 + 02 0.1083703E+00 0.8199077E-01 -0.8653056E-02 0.5676093E-01 0.3201328E-01 + 02 0.1192063E+00 0.8263056E-01 -0.9585221E-02 0.5705613E-01 0.3265308E-01 + 02 0.1300422E+00 0.8325786E-01 -0.1052719E-01 0.5732622E-01 0.3328037E-01 + 02 0.1387318E+00 0.8375138E-01 -0.1128967E-01 0.5752390E-01 0.3377389E-01 + 02 0.1474215E+00 0.8423610E-01 -0.1205759E-01 0.5770533E-01 0.3425861E-01 + 02 0.1561111E+00 0.8471167E-01 -0.1283054E-01 0.5787048E-01 0.3473418E-01 + 02 0.1648008E+00 0.8517772E-01 -0.1360810E-01 0.5801935E-01 0.3520023E-01 + 02 0.1734904E+00 0.8563385E-01 -0.1438981E-01 0.5815192E-01 0.3565637E-01 + 02 0.1822098E+00 0.8608118E-01 -0.1517786E-01 0.5826857E-01 0.3610370E-01 + 02 0.1909292E+00 0.8651772E-01 -0.1596909E-01 0.5836880E-01 0.3654024E-01 + 02 0.1996485E+00 0.8694305E-01 -0.1676296E-01 0.5845259E-01 0.3696557E-01 + 02 0.2083679E+00 0.8735670E-01 -0.1755891E-01 0.5851981E-01 0.3737921E-01 + 02 0.2170872E+00 0.8775823E-01 -0.1835636E-01 0.5857028E-01 0.3778074E-01 + 02 0.2280465E+00 0.8824585E-01 -0.1935947E-01 0.5861106E-01 0.3826837E-01 + 02 0.2390058E+00 0.8871490E-01 -0.2036338E-01 0.5862758E-01 0.3873741E-01 + 02 0.2499651E+00 0.8916683E-01 -0.2136755E-01 0.5862245E-01 0.3918934E-01 + 02 0.2609243E+00 0.8960144E-01 -0.2237112E-01 0.5859525E-01 0.3962396E-01 + 02 0.2697571E+00 0.8993482E-01 -0.2317859E-01 0.5855007E-01 0.3995734E-01 + 02 0.2785899E+00 0.9024317E-01 -0.2398073E-01 0.5846799E-01 0.4026568E-01 + 02 0.2874227E+00 0.9054041E-01 -0.2477994E-01 0.5837046E-01 0.4056293E-01 + 02 0.2962555E+00 0.9083315E-01 -0.2557727E-01 0.5826789E-01 0.4085567E-01 + 02 0.3050883E+00 0.9112000E-01 -0.2637444E-01 0.5815740E-01 0.4114251E-01 + 02 0.3139929E+00 0.9139976E-01 -0.2717637E-01 0.5803128E-01 0.4142227E-01 + 02 0.3228975E+00 0.9166768E-01 -0.2797485E-01 0.5788578E-01 0.4169020E-01 + 02 0.3318020E+00 0.9192412E-01 -0.2876905E-01 0.5772053E-01 0.4194663E-01 + 02 0.3407066E+00 0.9216946E-01 -0.2955895E-01 0.5753511E-01 0.4219198E-01 + 02 0.3496112E+00 0.9240469E-01 -0.3034469E-01 0.5732996E-01 0.4242720E-01 + 02 0.3608577E+00 0.9268918E-01 -0.3133151E-01 0.5704410E-01 0.4271169E-01 + 02 0.3721041E+00 0.9296115E-01 -0.3231229E-01 0.5672919E-01 0.4298366E-01 + 02 0.3833506E+00 0.9322199E-01 -0.3328731E-01 0.5638607E-01 0.4324451E-01 + 02 0.3945971E+00 0.9347234E-01 -0.3425662E-01 0.5601494E-01 0.4349485E-01 + 02 0.4037091E+00 0.9366758E-01 -0.3503783E-01 0.5569340E-01 0.4369010E-01 + 02 0.4128211E+00 0.9385575E-01 -0.3581491E-01 0.5535343E-01 0.4387826E-01 + 02 0.4219331E+00 0.9403638E-01 -0.3658737E-01 0.5499511E-01 0.4405890E-01 + 02 0.4310451E+00 0.9420891E-01 -0.3735457E-01 0.5461874E-01 0.4423142E-01 + 02 0.4401570E+00 0.9437262E-01 -0.3811573E-01 0.5422468E-01 0.4439514E-01 + 02 0.4493903E+00 0.9452873E-01 -0.3887991E-01 0.5380784E-01 0.4455124E-01 + 02 0.4586237E+00 0.9467415E-01 -0.3963594E-01 0.5337384E-01 0.4469667E-01 + 02 0.4678570E+00 0.9480809E-01 -0.4038279E-01 0.5292318E-01 0.4483060E-01 + 02 0.4770903E+00 0.9492982E-01 -0.4111949E-01 0.5245631E-01 0.4495233E-01 + 02 0.4863236E+00 0.9503879E-01 -0.4184516E-01 0.5197365E-01 0.4506130E-01 + 02 0.4931618E+00 0.9511096E-01 -0.4237519E-01 0.5160611E-01 0.4513348E-01 diff --git a/tests/bowl_slosh_vec/regression_data/regression_data_fgmax.txt b/tests/bowl_slosh_vec/regression_data/regression_data_fgmax.txt new file mode 100644 index 000000000..9a00deb36 --- /dev/null +++ b/tests/bowl_slosh_vec/regression_data/regression_data_fgmax.txt @@ -0,0 +1,2 @@ +2.028603161028199864e+01 +2.669859327530000428e+02 diff --git a/tests/bowl_slosh_vec/regression_tests.py b/tests/bowl_slosh_vec/regression_tests.py new file mode 100644 index 000000000..1ffcee608 --- /dev/null +++ b/tests/bowl_slosh_vec/regression_tests.py @@ -0,0 +1,75 @@ +#!/usr/bin/env python + +r"""Bowl-Slosh regression test for GeoClaw + +To create new regression data use + `python regression_tests.py True` +""" + +from __future__ import absolute_import +import os +import sys +import unittest + +import numpy + +import clawpack.geoclaw.test as test +import clawpack.geoclaw.topotools as topotools + + +class BowlSloshTest(test.GeoClawRegressionTest): + + r"""Bowl-Slosh regression test for GeoClaw""" + + def setUp(self): + + super(BowlSloshTest, self).setUp() + + # Make topography + a = 1. + h0 = 0.1 + topo_func = lambda x,y: h0 * (x**2 + y**2) / a**2 - h0 + + topo = topotools.Topography(topo_func=topo_func) + topo.topo_type = 2 + topo.x = numpy.linspace(-2.0, 2.0, 200) + topo.y = numpy.linspace(-2.0, 2.0, 200) + topo.write(os.path.join(self.temp_path, "bowl.topotype2"), \ + topo_type=2, Z_format="%22.15e") + + from . import make_fgmax_grid + make_fgmax_grid.make_fgmax_grid1(self.temp_path) + + + def runTest(self, save=False, indices=(2, 3)): + r"""Test bowl-slosh example with vectorization + + Note that this stub really only runs the code and performs no tests. + + """ + + # Write out data files + self.load_rundata() + self.write_rundata_objects() + + # Run code + self.run_code() + + # Perform tests + self.check_gauges(save=save, gauge_id=1, indices=(2, 3)) + self.check_fgmax(save=save) + self.success = True + + +if __name__=="__main__": + if len(sys.argv) > 1: + if bool(sys.argv[1]): + # Fake the setup and save out output + test = BowlSloshTest() + try: + test.setUp() + test.runTest(save=True) + finally: + test.tearDown() + sys.exit(0) + unittest.main() diff --git a/tests/bowl_slosh_vec/setplot.py b/tests/bowl_slosh_vec/setplot.py new file mode 100644 index 000000000..2a9f7a80e --- /dev/null +++ b/tests/bowl_slosh_vec/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 __future__ import absolute_import +import numpy +a = 1. +sigma = 0.5 +h0 = 0.1 +grav = 9.81 +omega = numpy.sqrt(2.*grav*h0) / a + +#-------------------------- +def setplot(plotdata): +#-------------------------- + + """ + Specify what is to be plotted at each frame. + Input: plotdata, an instance of pyclaw.plotters.data.ClawPlotData. + Output: a modified version of plotdata. + + """ + + + from clawpack.visclaw import colormaps, geoplot + + plotdata.clearfigures() # clear any old figures,axes,items data + + def set_drytol(current_data): + # The drytol parameter is used in masking land and water and + # affects what color map is used for cells with small water depth h. + # The cell will be plotted as dry if h < drytol. + # The best value to use often depends on the application and can + # be set here (measured in meters): + current_data.user["drytol"] = 1.e-3 + + plotdata.beforeframe = set_drytol + + #----------------------------------------- + # Figure for pcolor plot + #----------------------------------------- + plotfigure = plotdata.new_plotfigure(name='pcolor', figno=0) + + # Set up for axes in this figure: + plotaxes = plotfigure.new_plotaxes('pcolor') + plotaxes.title = 'Surface' + plotaxes.scaled = True + + # Water + plotitem = plotaxes.new_plotitem(plot_type='2d_pcolor') + plotitem.plot_var = geoplot.surface + plotitem.pcolor_cmap = geoplot.tsunami_colormap + plotitem.pcolor_cmin = -0.1 + plotitem.pcolor_cmax = 0.1 + plotitem.add_colorbar = True + plotitem.amr_celledges_show = [0,0,0] + plotitem.patchedges_show = 1 + + # Land + plotitem = plotaxes.new_plotitem(plot_type='2d_pcolor') + plotitem.plot_var = geoplot.land + plotitem.pcolor_cmap = geoplot.land_colors + plotitem.pcolor_cmin = 0.0 + plotitem.pcolor_cmax = 100.0 + plotitem.add_colorbar = False + plotitem.amr_celledges_show = [0,0,0] + plotitem.patchedges_show = 1 + plotaxes.xlimits = [-2,2] + plotaxes.ylimits = [-2,2] + + # Add contour lines of bathymetry: + plotitem = plotaxes.new_plotitem(plot_type='2d_contour') + plotitem.plot_var = geoplot.topo + from numpy import arange, linspace + plotitem.contour_levels = linspace(-.1, 0.5, 20) + plotitem.amr_contour_colors = ['k'] # color on each level + plotitem.kwargs = {'linestyles':'solid'} + plotitem.amr_contour_show = [1] + plotitem.celledges_show = 0 + plotitem.patchedges_show = 0 + plotitem.show = True + + #----------------------------------------- + # Figure for cross section + #----------------------------------------- + plotfigure = plotdata.new_plotfigure(name='cross-section', figno=1) + + # Set up for axes in this figure: + plotaxes = plotfigure.new_plotaxes() + plotaxes.xlimits = [-2,2] + plotaxes.ylimits = [-0.15,0.3] + plotaxes.title = 'Cross section at y=0' + def plot_topo_xsec(current_data): + from pylab import plot, hold, cos,sin,where,legend,nan + t = current_data.t + + hold(True) + x = linspace(-2,2,201) + y = 0. + B = h0*(x**2 + y**2)/a**2 - h0 + eta1 = sigma*h0/a**2 * (2.*x*cos(omega*t) + 2.*y*sin(omega*t) -sigma) + etatrue = where(eta1>B, eta1, nan) + plot(x, etatrue, 'r', label="true solution", linewidth=2) + plot(x, B, 'g', label="bathymetry") + ## plot([0],[-1],'kx',label="Level 1") # shouldn't show up in plots, + ## plot([0],[-1],'bo',label="Level 2") # but will produced desired legend + plot([0],[-1],'bo',label="Computed") ## need to fix plotstyle + legend() + hold(False) + plotaxes.afteraxes = plot_topo_xsec + + plotitem = plotaxes.new_plotitem(plot_type='1d_from_2d_data') + + def xsec(current_data): + # Return x value and surface eta at this point, along y=0 + from pylab import find,ravel + x = current_data.x + y = current_data.y + dy = current_data.dy + q = current_data.q + + ij = find((y <= dy/2.) & (y > -dy/2.)) + x_slice = ravel(x)[ij] + eta_slice = ravel(q[3,:,:])[ij] + return x_slice, eta_slice + + plotitem.map_2d_to_1d = xsec + plotitem.plotstyle = 'kx' ## need to be able to set amr_plotstyle + plotitem.kwargs = {'markersize':3} + plotitem.amr_show = [1] # plot on all levels + + + #----------------------------------------- + # Figure for grids alone + #----------------------------------------- + plotfigure = plotdata.new_plotfigure(name='grids', figno=2) + plotfigure.show = True + + # Set up for axes in this figure: + plotaxes = plotfigure.new_plotaxes() + plotaxes.xlimits = [-2,2] + plotaxes.ylimits = [-2,2] + plotaxes.title = 'grids' + 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 = [1,1,0] + plotitem.amr_patchedges_show = [1] + + + #----------------------------------------- + + # Parameters used only when creating html and/or latex hardcopy + # e.g., via pyclaw.plotters.frametools.printframes: + + plotdata.printfigs = True # print figures + plotdata.print_format = 'png' # file format + plotdata.print_framenos = 'all' # list of frames to print + plotdata.print_gaugenos = [] # list of gauges 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.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 + + diff --git a/tests/bowl_slosh_vec/setrun.py b/tests/bowl_slosh_vec/setrun.py new file mode 100644 index 000000000..067fc1ddc --- /dev/null +++ b/tests/bowl_slosh_vec/setrun.py @@ -0,0 +1,404 @@ +""" +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. + +""" + +from __future__ import absolute_import +from __future__ import print_function +import os +import numpy as np + + +#------------------------------ +def setrun(claw_pkg='geoclaw'): +#------------------------------ + + """ + Define the parameters used for running Clawpack. + + INPUT: + claw_pkg expected to be "geoclaw" for this setrun. + + OUTPUT: + rundata - object of class ClawRunData + + """ + + from clawpack.clawutil import data + + assert claw_pkg.lower() == 'geoclaw', "Expected claw_pkg = 'geoclaw'" + + num_dim = 2 + rundata = data.ClawRunData(claw_pkg, num_dim) + + #------------------------------------------------------------------ + # GeoClaw specific parameters: + #------------------------------------------------------------------ + rundata = setgeo(rundata) + + #------------------------------------------------------------------ + # 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] = -2.0 + clawdata.upper[0] = 2.0 + + clawdata.lower[1] = -2.0 + clawdata.upper[1] = 2.0 + + + + # Number of grid cells: Coarsest grid + clawdata.num_cells[0] = 41 + clawdata.num_cells[1] = 41 + + + # --------------- + # Size of system: + # --------------- + + # Number of equations in the system: + clawdata.num_eqn = 3 + + # Number of auxiliary variables in the aux array (initialized in setaux) + clawdata.num_aux = 1 + + # 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. + # The solution at initial time t0 is always written in addition. + + clawdata.output_style = 1 + + if clawdata.output_style == 1: + # Output nout frames at equally spaced times up to tfinal: + clawdata.num_output_times = 1 + clawdata.tfinal = 0.5 + clawdata.output_t0 = True # output at initial (or restart) time? + + elif clawdata.output_style == 2: + # Specify a list of output times. + clawdata.output_times = [0.5, 1.0] + + elif clawdata.output_style == 3: + # Output every iout timesteps with a total of ntot time steps: + clawdata.output_step_interval = 1 + clawdata.total_steps = 1 + clawdata.output_t0 = True + + + clawdata.output_format = 'ascii' # 'ascii' or '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 = 3 + + + + # -------------- + # Time stepping: + # -------------- + + # if dt_variable==1: variable time steps used based on cfl_desired, + # if dt_variable==0: fixed time steps dt = dt_initial will always be 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.0001 + + # Max time step to be allowed if variable dt used: + clawdata.dt_max = 1e+99 + + # Desired Courant number if variable dt used, and max to allow without + # retaking step with a smaller dt: + clawdata.cfl_desired = 0.75 + clawdata.cfl_max = 1.0 + + # Maximum number of time steps to allow between output times: + clawdata.steps_max = 5000 + + + + + # ------------------ + # Method to be used: + # ------------------ + + # Order of accuracy: 1 => Godunov, 2 => Lax-Wendroff plus limiters + clawdata.order = 2 + + # Use dimensional splitting? (not yet available for AMR) + 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 = 3 + + # 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 'mc' ==> MC limiter + # 4 or 'vanleer' ==> van Leer + clawdata.limiter = ['mc', 'mc', 'mc'] + + clawdata.use_fwaves = True # 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 = 'godunov' + + + # -------------------- + # Boundary conditions: + # -------------------- + + # Number of ghost cells (usually 2) + clawdata.num_ghost = 2 + + # Choice of BCs at xlower and xupper: + # 0 => user specified (must modify bcN.f to use this option) + # 1 => extrapolation (non-reflecting outflow) + # 2 => periodic (must specify this at both boundaries) + # 3 => solid wall for systems where q(2) is normal velocity + + clawdata.bc_lower[0] = 'extrap' + clawdata.bc_upper[0] = 'extrap' + + clawdata.bc_lower[1] = 'extrap' + clawdata.bc_upper[1] = 'extrap' + + # 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 mxnest-1) + amrdata.refinement_ratios_x = [4,4] + amrdata.refinement_ratios_y = [4,4] + amrdata.refinement_ratios_t = [2,6] + + + # Specify type of each aux variable in amrdata.auxtype. + # This must be a list of length maux, each element of which is one of: + # 'center', 'capacity', 'xleft', or 'yleft' (see documentation). + + amrdata.aux_type = ['center'] + + + # Flag using refinement routine flag2refine rather than richardson error + amrdata.flag_richardson = False # use Richardson? + amrdata.flag2refine = True + + # steps to take on each level L between regriddings of level L+1: + amrdata.regrid_interval = 3 + + # width of buffer zone around flagged points: + # (typically the same as regrid_interval so waves don't escape): + amrdata.regrid_buffer_width = 2 + + # 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.700000 + + # print info about each regridding up to this level: + amrdata.verbosity_regrid = 0 + + + # ----- 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 + + # More AMR parameters can be set -- see the defaults in pyclaw/data.py + + # == setregions.data values == + regions = rundata.regiondata.regions + # to specify regions of refinement append lines of the form + # [minlevel,maxlevel,t1,t2,x1,x2,y1,y2] + + # == setgauges.data values == + # for gauges append lines of the form [gaugeno, x, y, t1, t2] + rundata.gaugedata.gauges.append([1,0.5,0.5,0,1e10]) + + return rundata + # end of function setrun + # ---------------------- + + +#------------------- +def setgeo(rundata): +#------------------- + """ + Set GeoClaw specific runtime parameters. + For documentation see .... + """ + + try: + geo_data = rundata.geo_data + except: + print("*** Error, this rundata has no geo_data attribute") + raise AttributeError("Missing geo_data attribute") + + + # == Physics == + geo_data.gravity = 9.81 + geo_data.coordinate_system = 1 + geo_data.earth_radius = 6367.5e3 + + # == Forcing Options + geo_data.coriolis_forcing = False + + # == Algorithm and Initial Conditions == + geo_data.sea_level = -10.0 + geo_data.dry_tolerance = 1.e-3 + geo_data.friction_forcing = False + geo_data.manning_coefficient = 0.0 + geo_data.friction_depth = 1.e6 + + # Refinement data + refinement_data = rundata.refinement_data + refinement_data.wave_tolerance = 1.e-2 + refinement_data.deep_depth = 1e2 + refinement_data.max_level_deep = 3 + refinement_data.variable_dt_refinement_ratios = True + + # == settopo.data values == + topo_data = rundata.topo_data + # for topography, append lines of the form + # [topotype, minlevel, maxlevel, t1, t2, fname] + topo_data.topofiles.append([2, 1, 10, 0., 1.e10, 'bowl.topotype2']) + + # == setdtopo.data values == + dtopo_data = rundata.dtopo_data + # for moving topography, append lines of the form : (<= 1 allowed for now!) + # [topotype, minlevel,maxlevel,fname] + + # == setqinit.data values == + rundata.qinit_data.qinit_type = 0 + rundata.qinit_data.qinitfiles = [] + # for qinit perturbations, append lines of the form: (<= 1 allowed for now!) + # [minlev, maxlev, fname] + + # == setfixedgrids.data values == + fixedgrids = rundata.fixed_grid_data + # for fixed grids append lines of the form + # [t1,t2,noutput,x1,x2,y1,y2,xpoints,ypoints,\ + # ioutarrivaltimes,ioutsurfacemax] + + # == fgmax.data values == + fgmax_files = rundata.fgmax_data.fgmax_files + # for fixed grids append to this list names of any fgmax input files + fgmax_files.append('fgmax1.txt') + rundata.fgmax_data.num_fgmax_val = 2 + + return rundata + # end of function setgeo + # ---------------------- + + + +if __name__ == '__main__': + # Set up run-time parameters and write all data files. + import sys + rundata = setrun(*sys.argv[1:]) + rundata.write() +