Changeset 614

Timestamp:

Nov 23, 2004, 5:55:46 PM (20 years ago)

Author:

ole

Message:

Played with wind stress + c-extension

Location:

inundation/ga/storm_surge/pyvolution

Files:

• shallow_water.py (modified) (3 diffs)
• shallow_water_ext.c (modified) (4 diffs)
• wind_example_scalar.py (modified) (1 diff)
• wind_example_variable.py (modified) (5 diffs)

inundation/ga/storm_surge/pyvolution/shallow_water.py

@@ -1033 +1033 @@
                             s_vec, phi_vec, const):
     """Python version of assigning wind field to update vectors.
-    A c version also exists for speed
+    A c version also exists (for speed)
     """
     from math import pi, cos, sin, sqrt
@@ -1062 +1062 @@
     
     FIXME: This class may be incorporated in the generic Wind_stress class
-    SUperseded
+    Superseded
     """
 
@@ -1408 +1408 @@
     #Replace python version with c implementations
     
-    from shallow_water_ext import rotate
+    from shallow_water_ext import rotate, assign_windfield_values
     compute_fluxes = compute_fluxes_c
     gravity = gravity_c

inundation/ga/storm_surge/pyvolution/shallow_water_ext.c

@@ -20 +20 @@
 #include "util_ext.h"
 
+const double pi = 3.14159265358979;
 
 // Computational function for rotation
@@ -285 +286 @@
   return 0;
 }
+
+
+
+int _assign_wind_field_values(int N,
+                              double* xmom_update,
+                              double* ymom_update, 
+                              double* s_vec,
+                              double* phi_vec,
+                              double cw) {
+
+  //Assign windfield values to momentum updates
+
+  int k;
+  double S, s, phi, u, v;
+  
+  for (k=0; k<N; k++) {
+    
+    s = s_vec[k];
+    phi = phi_vec[k];
+
+    //Convert to radians
+    phi = phi*pi/180;
+
+    //Compute velocity vector (u, v)
+    u = s*cos(phi);
+    v = s*sin(phi);
+        
+    //Compute wind stress
+    S = cw * sqrt(u*u + v*v);
+    xmom_update[k] += S*u;
+    ymom_update[k] += S*v;        
+  }
+  return 0;  
+}            
 
 
@@ -673 +708 @@
 
 
+
+
+
+
+PyObject *assign_windfield_values(PyObject *self, PyObject *args) {
+  //
+  //      assign_windfield_values(xmom_update, ymom_update,
+  //                              s_vec, phi_vec, self.const)
+
+  
+
+  PyArrayObject   //(one element per triangle)
+  *s_vec,         //Speeds 
+  *phi_vec,       //Bearings 
+  *xmom_update,   //Momentum updates
+  *ymom_update; 
+
+    
+  int N; 
+  double cw;
+  
+  // Convert Python arguments to C  
+  if (!PyArg_ParseTuple(args, "OOOOd",  
+                        &xmom_update,
+                        &ymom_update,                   
+                        &s_vec, &phi_vec, 
+                        &cw))
+    return NULL;
+
+  N = xmom_update -> dimensions[0];
+    
+  _assign_wind_field_values(N,
+           (double*) xmom_update -> data,
+           (double*) ymom_update -> data, 
+           (double*) s_vec -> data, 
+           (double*) phi_vec -> data,
+           cw);
+ 
+  return Py_BuildValue("");  
+}
+
+
+
+
 //////////////////////////////////////////      
 // Method table for python module
@@ -688 +767 @@
    METH_VARARGS, "Print out"},          
   {"protect", protect, METH_VARARGS | METH_KEYWORDS, "Print out"},   
+  {"assign_windfield_values", assign_windfield_values, 
+   METH_VARARGS | METH_KEYWORDS, "Print out"},     
   //{"distribute_to_vertices_and_edges", 
   // distribute_to_vertices_and_edges, METH_VARARGS},    

inundation/ga/storm_surge/pyvolution/wind_example_scalar.py

@@ -23 +23 @@
 domain.set_quantity('elevation', 0.0)
 domain.set_quantity('level', 1.0)
-domain.set_quantity('friction', 0.1)
+domain.set_quantity('friction', 0.03)
 
 #Constant (quite extreme :-) windfield: 9000 m/s, bearing 120 degrees
-domain.forcing_terms.append( Wind_stress(s=9000, phi=135) )
+domain.forcing_terms.append( Wind_stress(s=9000, phi=120) )
 
 ######################

inundation/ga/storm_surge/pyvolution/wind_example_variable.py

@@ -9 +9 @@
 #
 from mesh_factory import rectangular
-from shallow_water import Domain, Reflective_boundary, Wind_stress
+from shallow_water import Domain, Dirichlet_boundary, Wind_stress
 
 #Create basic mesh
-N = 80
-length = 100
+N = 120
+length = 200
 points, vertices, boundary = rectangular(N, N, length, length)
 
@@ -25 +25 @@
 
 #Set initial conditions
+h = 1.0
 domain.set_quantity('elevation', 0.0)
-domain.set_quantity('level', 2.0)
-domain.set_quantity('friction', 0.07)
+domain.set_quantity('level', h)
+domain.set_quantity('friction', 0.01)
 
 
@@ -36 +37 @@
     """
 
-    from math import sqrt, exp, cos, pi
-
-    N = len(x)
-    s = 0*x  #New array
+    from math import pi
+    from Numeric import sqrt, exp, cos
 
     c = (length/2, length/2)
-    
-    for k in range(N):
-        r = sqrt((x[k] - c[0])**2 + (y[k] - c[1])**2)
+    r = sqrt((x - c[0])**2 + (y - c[1])**2)/length
+    factor = exp( -(r-0.15)**2 )
 
-        r /= max(length,length)
-
-        factor = exp( -(r-0.15)**2 )
-        
-        s[k] = 4000 * factor * (cos(t*2*pi/150) + 2)
-
-    return s
+    #return 9000 * factor
+    return 4000 * factor * (cos(t*2*pi/150) + 2)
 
 
@@ -58 +51 @@
     """Rotating field
     """
-    from math import sqrt, atan, pi
-
-    N = len(x)
-    a = 0*x  #New array
+    
+    from math import pi
+    from Numeric import sqrt, exp, cos, arctan2, choose, less
 
     c = (length/2, length/2)
-    for k in range(N):
-        r = sqrt((x[k] - c[0])**2 + (y[k] - c[1])**2)
+    xx = (x - c[0])/length
+    yy = (y - c[1])/length
+    angle = arctan2(yy,xx)
 
-        xx = (x[k] - c[0])/length
-        yy = (y[k] - c[1])/length
+    #Take normal direction (but reverse direction every 50 seconds)
+    #if sin(t*2*pi/100) < 0:
+    #    sgn = -1
+    #else:
+    #    sgn = 1
+    #angle += sgn*pi/2
+    angle -= pi/2    
 
-        angle = atan(yy/xx)
-
-        if xx < 0:
-            angle+=pi #atan in ]-pi/2; pi/2[
-
-        #Take normal direction
-        angle -= pi/2
-    
-        #Ensure positive radians    
-        if angle < 0:
-            angle += 2*pi
-
-        a[k] = angle/pi*180
-        
-    return a
+    #Convert to degrees and return
+    return angle/pi*180
 
     
@@ -109 +94 @@
 ######################
 # Boundary conditions
-Br = Reflective_boundary(domain)
-domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
+#Br = Reflective_boundary(domain)
+Bd = Dirichlet_boundary([h, 0, 0])
+domain.set_boundary({'left': Bd, 'right': Bd, 'top': Bd, 'bottom': Bd})
 
 ######################