Vectorization for Riemann solvers for the single-layer shallow water equations

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:

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
+      !$OMP 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
+         !$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:
+            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
+      !$OMP 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.
+          !$OMP 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:
+      !$OMP 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:
+      !$OMP 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