diff --git a/.travis.yml b/.travis.yml index 64ec49223..107eead81 100644 --- a/.travis.yml +++ b/.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 diff --git a/src/2d/shallow/Makefile.geoclaw b/src/2d/shallow/Makefile.geoclaw index 1bd276b13..f9c55e768 100644 --- a/src/2d/shallow/Makefile.geoclaw +++ b/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 \ diff --git a/src/2d/shallow/flag2refine2.f90 b/src/2d/shallow/flag2refine2.f90 index 6e3c20beb..1cb2a88c2 100644 --- a/src/2d/shallow/flag2refine2.f90 +++ b/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) diff --git a/src/2d/shallow/multilayer/flag2refine2.f90 b/src/2d/shallow/multilayer/flag2refine2.f90 index bcef062d5..4e9c49cba 100644 --- a/src/2d/shallow/multilayer/flag2refine2.f90 +++ b/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 diff --git a/src/python/geoclaw/dtopotools.py b/src/python/geoclaw/dtopotools.py index 55c0009bb..7efb3a84a 100644 --- a/src/python/geoclaw/dtopotools.py +++ b/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,10 +1646,14 @@ 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: @@ -1660,7 +1661,7 @@ def set_corners(self,corners,rake=90.,projection_zone=None): [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)