Changeset 5206

Timestamp:

Apr 11, 2008, 3:34:48 PM (15 years ago)

Author:

herve

Message:

modification of test_tsunami_okada.py to make txt file available.

Location:

anuga_core/source/anuga/shallow_water

Files:

• fullokada_MS.txt (added)
• fullokada_SP.txt (added)
• fullokada_SS.txt (added)
• test_tsunami_okada.py (modified) (9 diffs)
• tsunami_okada.py (modified) (9 diffs)

anuga_core/source/anuga/shallow_water/test_tsunami_okada.py

@@ -18 +18 @@
         from anuga.shallow_water import Domain
         from anuga.abstract_2d_finite_volumes.quantity import Quantity
-
-        """
-        Pick the test you want to dop; T= 0 test a point source,
+        from anuga.utilities.system_tools import get_pathname_from_package
+        """
+        Pick the test you want to do; T= 0 test a point source,
         T= 1  test single rectangular source, T= 2 test multiple
         rectangular sources
         """
-
+        #get path where this test is run
+        path= get_pathname_from_package('anuga.shallow_water')
         #choose what test to proceed
-        T=2
-
-        home = getenv('INUNDATIONHOME') + sep +'sandpits'+sep+'hdamlami'+sep+'test_okada'+sep 
+        T=1
+        
+
 
         if T==0:
             # fotran output file
-            filename = home+'fullokada_SP.txt'
+            
+            filename = path+sep+'fullokada_SP.txt'
         # initial condition of earthquake for multiple source
             x0 = 7000.0
@@ -47 +49 @@
         elif T==1:
         # fotran output file
-            filename = home+'fullokada_SS.txt'
+            filename = path+sep+'fullokada_SS.txt'
         # initial condition of earthquake for multiple source
             x0 = 7000.0
@@ -64 +66 @@
 
         # fotran output file
-            filename = home+'fullokada_MS.txt'
+            filename = path+sep+'fullokada_MS.txt'
         # initial condition of earthquake for multiple source
             x0 = [7000.0,10000.0]
@@ -115 +117 @@
         
         # call okada
-        Ts= Okada_func(ns, NSMAX,length=length, width=width, dip=dip, \
+        Ts= Okada_func(ns=ns, NSMAX=NSMAX,length=length, width=width, dip=dip, \
                        x0=x0, y0=y0, strike=strike, depth=depth, \
                        slip=slip, rake=rake,zrec=zrec)
@@ -148 +150 @@
         from anuga.shallow_water import Domain
         from anuga.abstract_2d_finite_volumes.quantity import Quantity
-
+        from anuga.utilities.system_tools import get_pathname_from_package
         """
         Pick the test you want to dop; T= 0 test a point source,
@@ -154 +156 @@
         rectangular sources
         """
-
+        
+        #get path where this test is run
+        path= get_pathname_from_package('anuga.shallow_water')
+        
         #choose what test to proceed
-        T=2
-
-        home = getenv('INUNDATIONHOME') + sep +'sandpits'+sep+'hdamlami'+sep+'test_okada'+sep 
+        T=1
 
         if T==0:
-            # fotran output file
-            filename = home+'fullokada_SP.txt'
+        # fotran output file
+            
+            filename = path+sep+'fullokada_SP.txt'
         # initial condition of earthquake for multiple source
             x0 = 7000.0
@@ -177 +181 @@
         elif T==1:
         # fotran output file
-            filename = home+'fullokada_SS.txt'
+            filename = path+sep+'fullokada_SS.txt'
         # initial condition of earthquake for multiple source
             x0 = 7000.0
@@ -194 +198 @@
 
         # fotran output file
-            filename = home+'fullokada_MS.txt'
+            filename = path+sep+'fullokada_MS.txt'
         # initial condition of earthquake for multiple source
             x0 = [7000.0,10000.0]
@@ -237 +241 @@
         domain.set_name('test')
         domain.set_quantity('elevation',topography)
-        
-        earthquake_tsunami(ns,NSMAX,length, width, strike, depth, \
-                       dip, x0, y0, slip=1.0, rake=90.,\
-                       domain=domain, verbose=False)
-
+        zrec = Quantity(domain)
+        zrec.set_values(topography)
+        
+        Ts = earthquake_tsunami(ns=ns,NSMAX=NSMAX,length=length, width=width, strike=strike,\
+                                depth=depth,dip=dip, xi=x0, yi=y0, slip=slip, rake=rake,\
+                                zrec=zrec,domain=domain, verbose=False)
+        
+        #create a variable to store vertical displacement throughout the domain
+        tsunami = Quantity(domain)
+        tsunami.set_values(Ts)
         k=0.0
         for i in range(0,6):

anuga_core/source/anuga/shallow_water/tsunami_okada.py

@@ -27 +27 @@
 """
 
-def earthquake_tsunami(ns,NSMAX,length, width, strike, depth, \
-                       dip, x0, y0, slip=1.0, rake=90.,\
-                       domain=None, verbose=False):
+def earthquake_tsunami(ns,NSMAX,length, width, strike, depth,\
+                       dip, xi, yi, slip, rake,\
+                       zrec,domain=None, verbose=False):
 
     from math import sin, radians
-    
+    from Numeric import zeros, Float
+    x0= zeros(ns,Float)
+    y0= zeros(ns,Float)
+    if ns ==1:
+        x0[0]=xi
+        y0[0]=yi
+    else:
+        x0=xi
+        y10=yi
     if domain is not None:
         xllcorner = domain.geo_reference.get_xllcorner()
         yllcorner = domain.geo_reference.get_yllcorner()
-        x0 = x0 - xllcorner  # fault origin (relative)
-        y0 = y0 - yllcorner
-        zrec=Quantity(domain)
-        zrec.set_values(topography)
+        for i in range(0,ns):
+            x0[i] = x0[i] - xllcorner  # fault origin (relative)
+            y0[i] = y0[i] - yllcorner
     #a few temporary print statements
     if verbose is True:
@@ -62 +69 @@
         raise Exception, msg
 
-    return Okada_func(ns,NSMAX,length=length, width=width, dip=dip, \
+    return Okada_func(ns=ns,NSMAX=NSMAX,length=length, width=width, dip=dip, \
                       x0=x0, y0=y0, strike=strike, depth=depth, \
                       slip=slip, rake=rake, zrec=zrec)
@@ -105 +112 @@
         self.ns=ns
         self.zrec=zrec
-        print 'hello'
         
     def __call__(self, x, y):
@@ -118 +124 @@
         from math import sin, cos, radians, exp, cosh
         from Numeric import zeros, Float
-        #from anuga.abstract_2d_finite_volumes.quantity import Quantity
-        #from okada import okadatest
         #ensure vectors x and y have the same length
         N = len(x)
@@ -128 +132 @@
         length = self.length
         width = self.width
-        y1 = self.x0
-        x1 = self.y0
+        y0 = self.x0
+        x0 = self.y0
         strike = self.strike
         dip = self.dip
@@ -164 +168 @@
             widths[0]=width
             dips[0]=dip
-            xs[0]=x1
-            ys[0]=y1
+            xs[0]=x0
+            ys[0]=y0
         else:
             dislocations=dislocation
@@ -173 +177 @@
             widths=width
             dips=dip
-            xs=x1
-            ys=y1
+            xs=x0
+            ys=y0
             depths=depth
         #double Gaussian calculation assumes water displacement is oriented
@@ -199 +203 @@
             for ist in range(0,ns):
                 #st = radians(strikes[ist])
-                st=strikes[ist]*rad
+                st=radians(strikes[ist])
                 csst = cos(st)
                 ssst = sin(st)
@@ -206 +210 @@
 #
                 #di = radians(dips[ist])
-                di=rad*dips[ist]
+                di=radians(dips[ist])
                 csdi =cos(di)
                 ssdi=sin(di)
 #
                 #ra= radians(rakes[ist])
-                ra=rad*rakes[ist]
+                ra=radians(rakes[ist])
                 csra=cos(ra)
                 ssra=sin(ra)