Merge pull request #43 from kbarnhart/dclaw4to5_misc

Final changes to make dclaw4 and dclaw5 equivalent.

Author: geoflows

geoclaw/2d/lib/stepgrid_geo.f

@@ -46,7 +46,7 @@ subroutine stepgrid(q,fm,fp,gm,gp,mitot,mjtot,mbc,dt,dtnew,dx,dy,
       dimension work(mwork)
 
       logical    debug,  dump
-      data       debug/.false./,  dump/.false./
+      data       debug/.true./,  dump/.true./
 c
 c     # set tcom = time.  This is in the common block comxyt that could
 c     # be included in the Riemann solver, for example, if t is explicitly
@@ -60,9 +60,8 @@ subroutine stepgrid(q,fm,fp,gm,gp,mitot,mjtot,mbc,dt,dtnew,dx,dy,
          write(outunit,*)" at start of stepgrid: dumping grid ",mptr
          do i = 1, mitot
          do j = 1, mjtot
-            write(outunit,545) i,j,(q(i,j,ivar),ivar=1,nvar),
-     .            aux(i,j,1),aux(i,j,8),aux(i,j,9)
- 545        format(2i4,4e15.7,/,8x,3e15.7,/,8x,3e15.7)
+            write(outunit,545) i,j,(q(i,j,ivar),ivar=1,nvar)
+ 545        format(2i4,7e15.7)
           end do
           end do
       endif
@@ -375,23 +374,11 @@ subroutine stepgrid(q,fm,fp,gm,gp,mitot,mjtot,mbc,dt,dtnew,dx,dy,
             write(outunit,*)" at end of stepgrid: dumping grid ",mptr
             do i = 1, mitot
             do j = 1, mjtot
-                  write(outunit,545) i,j,(q(i,j,ivar),ivar=1,nvar),
-     .            aux(i,j,1),aux(i,j,8),aux(i,j,9)
+                  write(outunit,545) i,j,(q(i,j,ivar),ivar=1,nvar)
             end do
             end do
       endif
 
-c
-!--      if (dump) then
-!--         write(*,*)" at end of stepgrid: dumping grid ",mptr
-!--         do i = 1, mitot
-!--         do j = 1, mjtot
-!--c            write(outunit,545) i,j,(q(i,j,ivar),ivar=1,nvar)
-!--            write(*,545) i,j,(q(i,j,ivar),ivar=1,nvar)
-!--         end do
-!--         end do
-!--      endif
-c
       return
       end
 

geoclaw/2d/lib/valout_geo.f

@@ -111,12 +111,16 @@ subroutine valout (lst, lend, time, nvar, naux)
             endif
             surface = alloc(iadd(i,j,1)) + alloc(iaddaux(i,j,1))
      &              - alloc(iadd(i,j,7))
-            write(matunit1,109) i, j, (alloc(iadd(i,j,ivar)), ivar=1,1),
+!            write(matunit1,110) i, j, (alloc(iadd(i,j,ivar)), ivar=1,1),
+!     &         surface
+             write(matunit1,109) (alloc(iadd(i,j,ivar)), ivar=1,nvar),
      &         surface
          enddo
          write(matunit1,*) ' '
       enddo
-  109       format(2i5, 9e26.16)
+  109       format(9e26.16)
+  110       format(2i5, 9e26.16)
+
 
 
             mptr = node(levelptr, mptr)
@@ -197,8 +201,11 @@ subroutine valout (lst, lend, time, nvar, naux)
                      alloc(iaddaux(i,j,ivar)) = 0.d0
                   endif
                enddo
-               write(matunit3,109) i, j, (alloc(iaddaux(i,j,ivar)),
+!               write(matunit3,110) i, j, (alloc(iaddaux(i,j,ivar)),
+!     &                              ivar=1,naux)
+               write(matunit3,109) (alloc(iaddaux(i,j,ivar)),
      &                              ivar=1,naux)
+
             enddo
             write(matunit3,*) ' '
          enddo

geoclaw/2d/lib_dig/digclaw_mod.f90

@@ -673,6 +673,8 @@ subroutine calc_taudir(maxmx,maxmy,meqn,mbc,mx,my,xlower,ylower,dx,dy,q,maux,aux
             endif
 
             pm = 0.5d0 !does not effect tau. only need tau in different cells
+            ! DIG: pm does impact rho, so it should be calculated here, consistently. should it be q(6,i,j)/h
+
             call admissibleq(h,hu,hv,hm,p,u,v,m,theta)
             call admissibleq(hL,huL,hvL,hmL,pL,uL,vL,mL,theta)
             call admissibleq(hB,huB,hvB,hmL,pB,uB,vB,mB,theta)

geoclaw/2d/lib_dig/riemannsolvers_geo.f

@@ -67,18 +67,18 @@ subroutine riemann_dig2_aug_sswave_ez(ixy,meqn,mwaves,hL,hR,
 
       pmL = chiL
       pmR = chiR
-      pm = 0.5*(chiL + chiR)
+      pm = 0.5d0*(chiL + chiR)
       pm = min(1.0d0,pm)
       pm = max(0.0d0,pm)
-      if (alpha_seg==1.0d0) then
+      if (dabs(alpha_seg-1.d0)<1.d-6) then
          seg = 0.0d0
       else
          seg = 1.0d0
       endif
       !pmtanh01 = seg*(0.5*(tanh(20.0*(pm-0.80))+1.0))
       !pmtanh01 = seg*(0.5*(tanh(40.0*(pm-0.90))+1.0))
       call calc_pmtanh(pm,seg,pmtanh01)
-      rho_fp = (1.0-pmtanh01)*rho_f
+      rho_fp = (1.0d0-pmtanh01)*rho_f
 
 
       if (hL.ge.drytol.and.hR.ge.drytol) then
@@ -114,9 +114,9 @@ subroutine riemann_dig2_aug_sswave_ez(ixy,meqn,mwaves,hL,hR,
       rho = 0.5d0*(rhoL + rhoR)
       tau = 0.5d0*(tauL + tauR)
       theta = 0.5d0*(thetaL + thetaR)
-      gamma = 0.25*(rho_fp + 3.0*rho)/rho
-      gammaL = 0.25*(rho_fp + 3.0*rhoL)/rhoL
-      gammaR = 0.25*(rho_fp + 3.0*rhoR)/rhoR
+      gamma = 0.25d0*(rho_fp + 3.0d0*rho)/rho
+      gammaL = 0.25d0*(rho_fp + 3.0d0*rhoL)/rhoL
+      gammaR = 0.25d0*(rho_fp + 3.0d0*rhoR)/rhoR
 
       eps = kappa + (1.d0-kappa)*gamma
 

geoclaw/2d/lib_dig/src2_geo.f90

@@ -60,6 +60,7 @@ subroutine src2(maxmx,maxmy,meqn,mbc,mx,my,xlower,ylower,dx,dy,q,maux,aux,t,dt)
 
       do i=1-mbc+1,mx+mbc-1
          do j=1-mbc+1,my+mbc-1
+
             theta = 0.d0
             dtheta = 0.d0
             if (bed_normal==1) then
@@ -214,9 +215,9 @@ subroutine src2(maxmx,maxmy,meqn,mbc,mx,my,xlower,ylower,dx,dy,q,maux,aux,t,dt)
                   !write(*,*) '------------'
                   !write(*,*) 'vu',t1bot
                   beta = 1.0d0-m!tanh(10.d0*m) !tan(1.5d0*p/(rho*gmod*h))/14.0d0
-                  gamma= rho*beta2*(vnorm**2)*(beta*gmod*coeff**2)/(tanh(h+1.d0-2.d0)**(1.0d0/3.0d0))
+                  !! DIG not used. segfault.gamma= rho*beta2*(vnorm**2)*(beta*gmod*coeff**2)/(tanh(h+1.d0-2.d0)**(1.0d0/3.0d0))
                   !write(*,*) 'gamma', gamma
-                  t1bot = t1bot + gamma
+                  !! DIG not used. segfault.t1bot = t1bot + gamma
                   t1bot = t1bot + tau!+p*tan(phi)
                   !write(*,*) 'tau',tau
                   t2top = min(t1bot,(1.d0-beta*entrainment_rate)*(tau))
@@ -323,7 +324,7 @@ subroutine src2(maxmx,maxmy,meqn,mbc,mx,my,xlower,ylower,dx,dy,q,maux,aux,t,dt)
                      q(i,j,2)=0.d0
                      q(i,j,3)=0.d0
                   else
-                     beta = 1.0d0-m 
+                     beta = 1.0d0-m
                      gamma= beta*dsqrt(hu**2 + hv**2)*(gmod*coeff**2)/(h**(7.d0/3.d0))
                      dgamma=1.d0 + dt*gamma
                      q(i,j,2)= q(i,j,2)/dgamma

python/pyclaw/data.py

@@ -1594,7 +1594,7 @@ def __init__(self, ndim):
         self.add_attribute("chi_init_val", 0.0, _DIG_ATTRS["chi_init_val"])
         self.add_attribute("kappita_diff", 1.0, _DIG_ATTRS["kappita_diff"])
         self.add_attribute("outaux", False, _DIG_ATTRS["outaux"])
-        self.add_attribute("curvature", 1, _DIG_ATTRS["curvature"])
+        self.add_attribute("curvature", 0, _DIG_ATTRS["curvature"])
         self.add_attribute("momlevel", 1, _DIG_ATTRS["momlevel"])
 
         # self.add_attribute('m_crit2', 0.62, 'critical state value of m (#) for different size species')