Reverse ordering of corners of a triangular fault

.travis.yml

@@ -7,7 +7,7 @@ env:
 
 before_install:
   - sudo apt-get update
-  - sudo apt-get install gfortran liblapack-dev
+  - sudo apt-get install gfortran liblapack-pic liblapack-dev
   # Print NumPy version that is already installed by Travis CI:
   - python -c "import numpy; print(numpy.__version__)"
   # - pip install matplotlib

src/2d/shallow/Makefile.geoclaw

@@ -35,7 +35,6 @@ COMMON_SOURCES += \
   $(AMRLIB)/quick_sort1.f \
   $(AMRLIB)/estdt.f \
   $(AMRLIB)/check4nans.f90 \
-  $(AMRLIB)/spest2.f  \
   $(AMRLIB)/init_iflags.f \
   $(AMRLIB)/igetsp.f \
   $(AMRLIB)/reclam.f \

src/2d/shallow/flag2refine2.f90

@@ -21,9 +21,9 @@
 !
 ! ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
 subroutine flag2refine2(mx,my,mbc,mbuff,meqn,maux,xlower,ylower,dx,dy,t,level, &
-                       tolsp,q,aux,amrflags,DONTFLAG,DOFLAG)
+                       tolsp,q,aux,amrflags)
 
-    use amr_module, only: mxnest, t0
+    use amr_module, only: mxnest, t0, DOFLAG, UNSET
     use geoclaw_module, only:dry_tolerance, sea_level
     use geoclaw_module, only: spherical_distance, coordinate_system
 
@@ -53,8 +53,6 @@ subroutine flag2refine2(mx,my,mbc,mbuff,meqn,maux,xlower,ylower,dx,dy,t,level, &
 
     ! Flagging
     real(kind=8), intent(in out) :: 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
@@ -67,8 +65,9 @@ subroutine flag2refine2(mx,my,mbc,mbuff,meqn,maux,xlower,ylower,dx,dy,t,level, &
     ! Storm specific variables
     real(kind=8) :: R_eye(2), wind_speed
 
-    ! Initialize flags
-    amrflags = DONTFLAG
+    ! Don't initialize flags, since they were already
+    ! flagged by flagregions2
+    ! 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)
@@ -97,7 +96,7 @@ subroutine flag2refine2(mx,my,mbc,mbuff,meqn,maux,xlower,ylower,dx,dy,t,level, &
                         ds = sqrt((x_c - R_eye(1))**2 + (y_c - R_eye(2))**2)
                     end if
 
-                    if ( ds < R_refine(m) .and. level <= m ) then
+                    if ( ds < R_refine(m) .and. level <= m .and. amrflags(i,j) == UNSET) then
                         amrflags(i,j) = DOFLAG
                         cycle x_loop
                     endif
@@ -107,7 +106,8 @@ subroutine flag2refine2(mx,my,mbc,mbuff,meqn,maux,xlower,ylower,dx,dy,t,level, &
                 if (wind_forcing) then
                     wind_speed = sqrt(aux(wind_index,i,j)**2 + aux(wind_index+1,i,j)**2)
                     do m=1,size(wind_refine,1)
-                        if ((wind_speed > wind_refine(m)) .and. (level <= m)) then
+                        if ((wind_speed > wind_refine(m)) .and. (level <= m) &
+                             .and. amrflags(i,j) == UNSET) then
                             amrflags(i,j) = DOFLAG
                             cycle x_loop
                         endif
@@ -120,20 +120,8 @@ subroutine flag2refine2(mx,my,mbc,mbuff,meqn,maux,xlower,ylower,dx,dy,t,level, &
             do m=1,mtopofiles
                 if (level < minleveltopo(m) .and. t >= tlowtopo(m) .and. t <= thitopo(m)) then
                     if (  x_hi > xlowtopo(m) .and. x_low < xhitopo(m) .and. &
-                          y_hi > ylowtopo(m) .and. y_low < yhitopo(m) ) then
-
-                        amrflags(i,j) = DOFLAG
-                        cycle x_loop
-                    endif
-                endif
-            enddo
-
-            ! Check to see if refinement is forced in any other region:
-            do m=1,num_regions
-                if (level < regions(m)%min_level .and. &
-                    t >= regions(m)%t_low .and. t <= regions(m)%t_hi) then
-                    if (x_hi > regions(m)%x_low .and. x_low < regions(m)%x_hi .and. &
-                        y_hi > regions(m)%y_low .and. y_low < regions(m)%y_hi ) then
+                          y_hi > ylowtopo(m) .and. y_low < yhitopo(m) &
+                           .and. amrflags(i,j) == UNSET) then
 
                         amrflags(i,j) = DOFLAG
                         cycle x_loop
@@ -147,7 +135,8 @@ subroutine flag2refine2(mx,my,mbc,mbuff,meqn,maux,xlower,ylower,dx,dy,t,level, &
                 if (level < minleveldtopo(m).and. &
                     t <= tfdtopo(m) .and. & !t.ge.t0dtopo(m).and.
                     x_hi > xlowdtopo(m) .and. x_low < xhidtopo(m).and. &
-                    y_hi > ylowdtopo(m) .and. y_low < yhidtopo(m)) then
+                    y_hi > ylowdtopo(m) .and. y_low < yhidtopo(m) &
+                    .and. amrflags(i,j) == UNSET) then
 
                     amrflags(i,j) = DOFLAG
                     cycle x_loop
@@ -158,7 +147,7 @@ subroutine flag2refine2(mx,my,mbc,mbuff,meqn,maux,xlower,ylower,dx,dy,t,level, &
             ! specified and need to force refinement:
             ! This assumes that t0 = 0.d0, should really be t0 but we do
             ! not have access to that parameter in this routine
-            if (qinit_type > 0 .and. t == t0) then
+            if (qinit_type > 0 .and. t == t0 .and. amrflags(i,j) == UNSET) then
                 if (level < min_level_qinit .and. &
                     x_hi > x_low_qinit .and. x_low < x_hi_qinit .and. &
                     y_hi > y_low_qinit .and. y_low < y_hi_qinit) then
@@ -169,9 +158,12 @@ subroutine flag2refine2(mx,my,mbc,mbuff,meqn,maux,xlower,ylower,dx,dy,t,level, &
             endif
 
             ! -----------------------------------------------------------------
-            ! Refinement not forced, so check if it is allowed and if so,
-            ! check if there is a reason to flag this point:
-            if (allowflag(x_c,y_c,t,level)) then
+            ! Refinement not forced, so check if it is allowed 
+            ! and if the flag is still UNSET. If so,
+            ! check if there is a reason to flag this point.
+            ! If flag == DONTFLAG then refinement is forbidden by a region,
+            ! if flag == DOFLAG checking is not needed
+            if (allowflag(x_c,y_c,t,level) .and. amrflags(i,j) == UNSET) then
 
                 if (q(1,i,j) > dry_tolerance) then
                     eta = q(1,i,j) + aux(1,i,j)

src/2d/shallow/multilayer/flag2refine2.f90

@@ -21,9 +21,9 @@
 !
 ! ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
 subroutine flag2refine2(mx,my,mbc,mbuff,meqn,maux,xlower,ylower,dx,dy,t,level, &
-                       tolsp,q,aux,amrflags,DONTFLAG,DOFLAG)
+                       tolsp,q,aux,amrflags)
 
-    use amr_module, only: mxnest, t0
+    use amr_module, only: mxnest, t0, DONTFLAG, DOFLAG, UNSET
     use geoclaw_module, only: sea_level, rho
     use geoclaw_module, only: spherical_distance, coordinate_system
 
@@ -55,8 +55,6 @@ subroutine flag2refine2(mx,my,mbc,mbuff,meqn,maux,xlower,ylower,dx,dy,t,level, &
 
     ! Flagging
     real(kind=8), intent(in out) :: 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
@@ -69,8 +67,9 @@ subroutine flag2refine2(mx,my,mbc,mbuff,meqn,maux,xlower,ylower,dx,dy,t,level, &
     ! Storm specific variables
     real(kind=8) :: R_eye(2), wind_speed
 
-    ! Initialize flags
-    amrflags = DONTFLAG
+    ! Don't initialize flags, since they were already
+    ! flagged by flagregions2
+    ! 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)
@@ -99,7 +98,7 @@ subroutine flag2refine2(mx,my,mbc,mbuff,meqn,maux,xlower,ylower,dx,dy,t,level, &
                         ds = sqrt((x_c - R_eye(1))**2 + (y_c - R_eye(2))**2)
                     end if
 
-                    if ( ds < R_refine(m) .and. level <= m ) then
+                    if ( ds < R_refine(m) .and. level <= m .and. amrflags(i,j) == UNSET) then
                         amrflags(i,j) = DOFLAG
                         cycle x_loop
                     endif
@@ -109,7 +108,8 @@ subroutine flag2refine2(mx,my,mbc,mbuff,meqn,maux,xlower,ylower,dx,dy,t,level, &
                 if (wind_forcing) then
                     wind_speed = sqrt(aux(wind_index,i,j)**2 + aux(wind_index+1,i,j)**2)
                     do m=1,size(wind_refine,1)
-                        if ((wind_speed > wind_refine(m)) .and. (level <= m)) then
+                        if ((wind_speed > wind_refine(m)) .and. (level <= m) &
+                              .and. amrflags(i,j) == UNSET) then
                             amrflags(i,j) = DOFLAG
                             cycle x_loop
                         endif
@@ -122,20 +122,8 @@ subroutine flag2refine2(mx,my,mbc,mbuff,meqn,maux,xlower,ylower,dx,dy,t,level, &
             do m=1,mtopofiles
                 if (level < minleveltopo(m) .and. t >= tlowtopo(m) .and. t <= thitopo(m)) then
                     if (  x_hi > xlowtopo(m) .and. x_low < xhitopo(m) .and. &
-                          y_hi > ylowtopo(m) .and. y_low < yhitopo(m) ) then
-
-                        amrflags(i,j) = DOFLAG
-                        cycle x_loop
-                    endif
-                endif
-            enddo
-
-            ! Check to see if refinement is forced in any other region:
-            do m=1,num_regions
-                if (level < regions(m)%min_level .and. &
-                    t >= regions(m)%t_low .and. t <= regions(m)%t_hi) then
-                    if (x_hi > regions(m)%x_low .and. x_low < regions(m)%x_hi .and. &
-                        y_hi > regions(m)%y_low .and. y_low < regions(m)%y_hi ) then
+                          y_hi > ylowtopo(m) .and. y_low < yhitopo(m) &
+                           .and. amrflags(i,j) == UNSET) then
 
                         amrflags(i,j) = DOFLAG
                         cycle x_loop
@@ -149,7 +137,8 @@ subroutine flag2refine2(mx,my,mbc,mbuff,meqn,maux,xlower,ylower,dx,dy,t,level, &
                 if (level < minleveldtopo(m).and. &
                     t <= tfdtopo(m) .and. & !t.ge.t0dtopo(m).and.
                     x_hi > xlowdtopo(m) .and. x_low < xhidtopo(m).and. &
-                    y_hi > ylowdtopo(m) .and. y_low < yhidtopo(m)) then
+                    y_hi > ylowdtopo(m) .and. y_low < yhidtopo(m) &
+                     .and. amrflags(i,j) == UNSET) then
 
                     amrflags(i,j) = DOFLAG
                     cycle x_loop
@@ -160,7 +149,7 @@ subroutine flag2refine2(mx,my,mbc,mbuff,meqn,maux,xlower,ylower,dx,dy,t,level, &
             ! specified and need to force refinement:
             ! This assumes that t0 = 0.d0, should really be t0 but we do
             ! not have access to that parameter in this routine
-            if (qinit_type > 0 .and. t == t0) then
+            if (qinit_type > 0 .and. t == t0 .and. amrflags(i,j) == UNSET) then
                 if (level < min_level_qinit .and. &
                     x_hi > x_low_qinit .and. x_low < x_hi_qinit .and. &
                     y_hi > y_low_qinit .and. y_low < y_hi_qinit) then
@@ -173,7 +162,7 @@ subroutine flag2refine2(mx,my,mbc,mbuff,meqn,maux,xlower,ylower,dx,dy,t,level, &
             ! -----------------------------------------------------------------
             ! Refinement not forced, so check if it is allowed and if so,
             ! check if there is a reason to flag this point:
-            if (allowflag(x_c,y_c,t,level)) then
+            if (allowflag(x_c,y_c,t,level) .and. amrflags(i,j) == UNSET) then
                 b = aux(1,i,j)
                 do layer = num_layers, 1, -1
                     if (q(3*layer-2,i,j) / rho(layer) > dry_tolerance(layer)) then

src/python/geoclaw/dtopotools.py

@@ -1316,6 +1316,7 @@ def __init__(self):
         self._corners = None
         self._gauss_pts = None
         self.n_gauss_pts = 4
+        self._fix_orientation = False
 
 
     def convert_to_standard_units(self, input_units, verbose=False):
@@ -1523,17 +1524,17 @@ def calculate_geometry(self):
                                             self._centers[0][1] 
                                                                  - up_strike[1])
 
-    def calculate_geometry_triangles(self,rake=90.):
+    def calculate_geometry_triangles(self):
         r"""
         Calculate geometry for triangular subfaults
 
         - Uses *corners* to calculate *centers*, *longitude*, *latitude*,
-          *depth*, *strike*, *dip*, *rake*, *length*, *width*.
+          *depth*, *strike*, *dip*, *length*, *width*.
 
         - sets *coordinate_specification* as "triangular"
 
         - Note that calculate_geometry() computes 
-          long/lat/strik/dip/rake/length/width to calculate centers/corners
+          long/lat/strike/dip/length/width to calculate centers/corners
 
         """
 
@@ -1552,9 +1553,8 @@ def calculate_geometry_triangles(self,rake=90.):
 
             x0 = numpy.array(self.corners)
             x = x0.copy()
-            x = numpy.array(x0)
 
-            x[:,2] = - numpy.abs(x[:,2]) # set depth to be always negative(lazy)
+            x[:,2] = - numpy.abs(x[:,2]) # set depth to be always negative
 
             if 0:
                 # old coordinate transform
@@ -1581,6 +1581,7 @@ def calculate_geometry_triangles(self,rake=90.):
             normal = cross(v1,v2)
             if normal[2] < 0:
                 normal = -normal
+                self._fix_orientation = True
             strikev = cross(normal,e3)   # vector in strike direction
 
             a = normal[0]
@@ -1589,19 +1590,16 @@ def calculate_geometry_triangles(self,rake=90.):
 
             #Compute strike
             strike_deg = rad2deg(numpy.arctan(-b/a))
-            print('+++ initial strike_deg = %g' % strike_deg)
 
             #Compute dip
             beta = deg2rad(strike_deg + 90)
             m = numpy.array([sin(beta),cos(beta),0]) #Points in dip direction
             n = numpy.array([a,b,c]) #Normal to the plane
-            #dip_deg = abs(rad2deg(numpy.arcsin(m.dot(n)/(norm(m)*norm(n)))))
 
             if abs(c) < 1e-8:
                 dip_deg = 90.   # vertical fault
             else:
                 dip_deg = rad2deg(numpy.arcsin(m.dot(n)/(norm(m)*norm(n))))
-            print('+++ initial dip_deg = %g' % dip_deg)
 
             # dip should be between 0 and 90. If negative, reverse strike:
             if dip_deg < 0:
@@ -1616,7 +1614,6 @@ def calculate_geometry_triangles(self,rake=90.):
 
             self.strike = strike_deg
             self.dip = dip_deg
-            self.rake = rake     # set default rake to 90 degrees
 
             # find the center line
             xx = numpy.zeros((3,3))
@@ -1649,18 +1646,22 @@ def calculate_geometry_triangles(self,rake=90.):
             raise ValueError("Invalid coordinate specification %s." \
                                                 % self.coordinate_specification)
 
-    def set_corners(self,corners,rake=90.,projection_zone=None):
+    def set_corners(self,corners,projection_zone=None):
         r"""
-            set three corners for a triangular fault.
-            Input *corners* should be iterable of length 3.
+        Set three corners for a triangular fault.
+        :Inputs: 
+
+        - *corners* should be iterable of length 3.
+        - *rake* should be between -180. and 180.
+
         """
 
         if len(corners) == 3:
             self._corners =\
                 [corners[0],corners[1],corners[2]]
             self._projection_zone = projection_zone
             self.coordinate_specification = 'triangular'
-            self.calculate_geometry_triangles(rake=rake)
+            self.calculate_geometry_triangles()
         else:
             raise ValueError("Expected input of length 3")
 
@@ -1873,9 +1874,13 @@ def okada(self, x, y):
                     Odepth = abs(O2_list[k][2])
 
                 if reverse_list[k]:
-                    sgn = (-1.)**(floor(j/3)+1)
+                    sgn = (-1.)**(floor(j/3))
                 else:
-                    sgn = (-1.)**floor(j/3)
+                    sgn = (-1.)**floor(j/3+1)
+
+                # fix orientation 
+                if self._fix_orientation:
+                    sgn *= -1.
 
                 Y1,Y2,Y3,Z1,Z2,Z3,Yb1,Yb2,Yb3,Zb1,Zb2,Zb3 = \
                 self._get_halfspace_coords(X1,X2,X3,alpha,beta,Olong,Olat,Odepth)
@@ -1906,7 +1911,7 @@ def okada(self, x, y):
             dZ = -v31*burgersv[0] - v32*burgersv[1] + v33*burgersv[2]
 
             dtopo = DTopography()
-            dtopo.X = X1    # DR: X1, X2 varname confusing?
+            dtopo.X = X1    
             dtopo.Y = X2
             dtopo.dX = numpy.array(dX, ndmin=3)
             dtopo.dY = numpy.array(dY, ndmin=3)
@@ -1940,7 +1945,7 @@ def _get_leg_angles(self):
         y[:,1] = LAT2METER * x[:,1]
         y[:,2] = - numpy.abs(x[:,2])    # force sign
 
-        v_list = [y[1,:] - y[0,:], y[2,:] - y[1,:], y[0,:] - y[2,:]]
+        v_list = [y[0,:] - y[1,:], y[1,:] - y[2,:], y[2,:] - y[0,:]]
 
         e3 = numpy.array([0.,0.,-1.])
 
@@ -1961,10 +1966,8 @@ def _get_leg_angles(self):
                 l = j
                 reverse_list[j] = True
 
-            O1 = x[k,:].tolist()     # set origin for the vector v
-            O2 = x[l,:].tolist()     # set dest.  for the vector v
-            O1_list.append(O1)
-            O2_list.append(O2)
+            O1_list.append(x[k,:].copy())  # set origin for the vector v
+            O2_list.append(x[l,:].copy())  # set dest.  for the vector v
 
             alpha = numpy.arctan2(vn[0],vn[1])
             alpha_list.append(alpha)