Changeset 5206
 Timestamp:
 Apr 11, 2008, 3:34:48 PM (15 years ago)
 Location:
 anuga_core/source/anuga/shallow_water
 Files:

 3 added
 2 edited
Legend:
 Unmodified
 Added
 Removed

anuga_core/source/anuga/shallow_water/test_tsunami_okada.py
r5204 r5206 18 18 from anuga.shallow_water import Domain 19 19 from anuga.abstract_2d_finite_volumes.quantity import Quantity 20 21 """ 22 Pick the test you want to do p; T= 0 test a point source,20 from anuga.utilities.system_tools import get_pathname_from_package 21 """ 22 Pick the test you want to do; T= 0 test a point source, 23 23 T= 1 test single rectangular source, T= 2 test multiple 24 24 rectangular sources 25 25 """ 26 26 #get path where this test is run 27 path= get_pathname_from_package('anuga.shallow_water') 27 28 #choose what test to proceed 28 T= 229 30 home = getenv('INUNDATIONHOME') + sep +'sandpits'+sep+'hdamlami'+sep+'test_okada'+sep 29 T=1 30 31 31 32 32 33 if T==0: 33 34 # fotran output file 34 filename = home+'fullokada_SP.txt' 35 36 filename = path+sep+'fullokada_SP.txt' 35 37 # initial condition of earthquake for multiple source 36 38 x0 = 7000.0 … … 47 49 elif T==1: 48 50 # fotran output file 49 filename = home+'fullokada_SS.txt'51 filename = path+sep+'fullokada_SS.txt' 50 52 # initial condition of earthquake for multiple source 51 53 x0 = 7000.0 … … 64 66 65 67 # fotran output file 66 filename = home+'fullokada_MS.txt'68 filename = path+sep+'fullokada_MS.txt' 67 69 # initial condition of earthquake for multiple source 68 70 x0 = [7000.0,10000.0] … … 115 117 116 118 # call okada 117 Ts= Okada_func(ns ,NSMAX,length=length, width=width, dip=dip, \119 Ts= Okada_func(ns=ns, NSMAX=NSMAX,length=length, width=width, dip=dip, \ 118 120 x0=x0, y0=y0, strike=strike, depth=depth, \ 119 121 slip=slip, rake=rake,zrec=zrec) … … 148 150 from anuga.shallow_water import Domain 149 151 from anuga.abstract_2d_finite_volumes.quantity import Quantity 150 152 from anuga.utilities.system_tools import get_pathname_from_package 151 153 """ 152 154 Pick the test you want to dop; T= 0 test a point source, … … 154 156 rectangular sources 155 157 """ 156 158 159 #get path where this test is run 160 path= get_pathname_from_package('anuga.shallow_water') 161 157 162 #choose what test to proceed 158 T=2 159 160 home = getenv('INUNDATIONHOME') + sep +'sandpits'+sep+'hdamlami'+sep+'test_okada'+sep 163 T=1 161 164 162 165 if T==0: 163 # fotran output file 164 filename = home+'fullokada_SP.txt' 166 # fotran output file 167 168 filename = path+sep+'fullokada_SP.txt' 165 169 # initial condition of earthquake for multiple source 166 170 x0 = 7000.0 … … 177 181 elif T==1: 178 182 # fotran output file 179 filename = home+'fullokada_SS.txt'183 filename = path+sep+'fullokada_SS.txt' 180 184 # initial condition of earthquake for multiple source 181 185 x0 = 7000.0 … … 194 198 195 199 # fotran output file 196 filename = home+'fullokada_MS.txt'200 filename = path+sep+'fullokada_MS.txt' 197 201 # initial condition of earthquake for multiple source 198 202 x0 = [7000.0,10000.0] … … 237 241 domain.set_name('test') 238 242 domain.set_quantity('elevation',topography) 239 240 earthquake_tsunami(ns,NSMAX,length, width, strike, depth, \ 241 dip, x0, y0, slip=1.0, rake=90.,\ 242 domain=domain, verbose=False) 243 243 zrec = Quantity(domain) 244 zrec.set_values(topography) 245 246 Ts = earthquake_tsunami(ns=ns,NSMAX=NSMAX,length=length, width=width, strike=strike,\ 247 depth=depth,dip=dip, xi=x0, yi=y0, slip=slip, rake=rake,\ 248 zrec=zrec,domain=domain, verbose=False) 249 250 #create a variable to store vertical displacement throughout the domain 251 tsunami = Quantity(domain) 252 tsunami.set_values(Ts) 244 253 k=0.0 245 254 for i in range(0,6): 
anuga_core/source/anuga/shallow_water/tsunami_okada.py
r5204 r5206 27 27 """ 28 28 29 def earthquake_tsunami(ns,NSMAX,length, width, strike, depth, 30 dip, x 0, y0, slip=1.0, rake=90.,\31 domain=None, verbose=False):29 def earthquake_tsunami(ns,NSMAX,length, width, strike, depth,\ 30 dip, xi, yi, slip, rake,\ 31 zrec,domain=None, verbose=False): 32 32 33 33 from math import sin, radians 34 34 from Numeric import zeros, Float 35 x0= zeros(ns,Float) 36 y0= zeros(ns,Float) 37 if ns ==1: 38 x0[0]=xi 39 y0[0]=yi 40 else: 41 x0=xi 42 y10=yi 35 43 if domain is not None: 36 44 xllcorner = domain.geo_reference.get_xllcorner() 37 45 yllcorner = domain.geo_reference.get_yllcorner() 38 x0 = x0  xllcorner # fault origin (relative) 39 y0 = y0  yllcorner 40 zrec=Quantity(domain) 41 zrec.set_values(topography) 46 for i in range(0,ns): 47 x0[i] = x0[i]  xllcorner # fault origin (relative) 48 y0[i] = y0[i]  yllcorner 42 49 #a few temporary print statements 43 50 if verbose is True: … … 62 69 raise Exception, msg 63 70 64 return Okada_func(ns ,NSMAX,length=length, width=width, dip=dip, \71 return Okada_func(ns=ns,NSMAX=NSMAX,length=length, width=width, dip=dip, \ 65 72 x0=x0, y0=y0, strike=strike, depth=depth, \ 66 73 slip=slip, rake=rake, zrec=zrec) … … 105 112 self.ns=ns 106 113 self.zrec=zrec 107 print 'hello'108 114 109 115 def __call__(self, x, y): … … 118 124 from math import sin, cos, radians, exp, cosh 119 125 from Numeric import zeros, Float 120 #from anuga.abstract_2d_finite_volumes.quantity import Quantity121 #from okada import okadatest122 126 #ensure vectors x and y have the same length 123 127 N = len(x) … … 128 132 length = self.length 129 133 width = self.width 130 y 1= self.x0131 x 1= self.y0134 y0 = self.x0 135 x0 = self.y0 132 136 strike = self.strike 133 137 dip = self.dip … … 164 168 widths[0]=width 165 169 dips[0]=dip 166 xs[0]=x 1167 ys[0]=y 1170 xs[0]=x0 171 ys[0]=y0 168 172 else: 169 173 dislocations=dislocation … … 173 177 widths=width 174 178 dips=dip 175 xs=x 1176 ys=y 1179 xs=x0 180 ys=y0 177 181 depths=depth 178 182 #double Gaussian calculation assumes water displacement is oriented … … 199 203 for ist in range(0,ns): 200 204 #st = radians(strikes[ist]) 201 st= strikes[ist]*rad205 st=radians(strikes[ist]) 202 206 csst = cos(st) 203 207 ssst = sin(st) … … 206 210 # 207 211 #di = radians(dips[ist]) 208 di=rad *dips[ist]212 di=radians(dips[ist]) 209 213 csdi =cos(di) 210 214 ssdi=sin(di) 211 215 # 212 216 #ra= radians(rakes[ist]) 213 ra=rad *rakes[ist]217 ra=radians(rakes[ist]) 214 218 csra=cos(ra) 215 219 ssra=sin(ra)
Note: See TracChangeset
for help on using the changeset viewer.