Changeset 8973

Timestamp:

Sep 12, 2013, 3:18:38 PM (12 years ago)

Author:

steve

Message:

set value of quantity using asc data

Location:

trunk/anuga_core

Files:

• source/anuga/abstract_2d_finite_volumes/quantity.py (modified) (4 diffs)
• source/anuga/abstract_2d_finite_volumes/test_quantity.py (modified) (4 diffs)
• user_manual/anuga_user_manual.pdf (modified) (previous)
• user_manual/anuga_user_manual.tex (modified) (2 diffs)

trunk/anuga_core/source/anuga/abstract_2d_finite_volumes/quantity.py

@@ -1139 +1139 @@
 
             
-        if location == 'centroids':
-            points = self.domain.centroid_coordinates
-        else:
-            points = self.domain.vertex_coordinates
+
             
     
@@ -1238 +1235 @@
         datafile.close()
         
-        #print Z.shape, Z
-    
-        x = num.linspace(xllcorner, xllcorner+cellsize*(nrows-1), nrows)
-        y = num.linspace(yllcorner, yllcorner+cellsize*(ncols-1), ncols)
+        #print Z.shape
+        #print Z
+        
+        # For raster data we need to a flip and transpose
+        Z = numpy.flipud(Z)
+
+        # Transpose z to have y coordinates along the first axis and x coordinates
+        # along the second axis
+        Z = Z.transpose()
+    
+        x = num.linspace(xllcorner, xllcorner+cellsize*(ncols-1), ncols)
+        y = num.linspace(yllcorner, yllcorner+cellsize*(nrows-1), nrows)
+        
+        
+        if location == 'centroids':
+            points = self.domain.centroid_coordinates
+        else:
+            points = self.domain.vertex_coordinates
+            
+        from anuga.geospatial_data.geospatial_data import Geospatial_data,  ensure_absolute
+        points = ensure_absolute(points, geo_reference=self.domain.geo_reference)
+            
+        print numpy.max(points[:,0])
+        print numpy.min(points[:,0])
+        print numpy.max(points[:,1])
+        print numpy.min(points[:,1])
+        
+        print numpy.max(x)
+        print numpy.min(x)
+        print numpy.max(y)
+        print numpy.min(y)
+        
         
         #print x.shape, x
@@ -1275 +1300 @@
                 #assert values.shape[0] == N, msg
 
-                self.vertex_values[:] = values
+                #print values.shape
+                #print self.vertex_values.shape 
+                self.vertex_values[:] = values.reshape((-1,3))
             else:
                 msg = 'Number of values must match number of indices'
@@ -1281 +1308 @@
 
                 # Brute force
-                self.vertex_values[indices] = values
+                self.vertex_values[indices] = values.reshape((-1,3))
                 
     def get_extremum_index(self, mode=None, indices=None):

trunk/anuga_core/source/anuga/abstract_2d_finite_volumes/test_quantity.py

@@ -1388 +1388 @@
 
     def test_set_values_from_asc_vertices(self):
-
-        quantity= Quantity(self.mesh_onslow)
-
-
         """ Format of asc file 
         ncols         11
@@ -1400 +1396 @@
         NODATA_value  -9999
         """
-        
-        # UTM round Onslow 
-        
+
+        #x0 = 314036.58727982
+        #y0 = 6224951.2960092
+        x0 = 0.0
+        y0 = 0.0
+
+        a = [0.0, 0.0]
+        b = [0.0, 2.0]
+        c = [2.0, 0.0]
+        d = [0.0, 4.0]
+        e = [2.0, 2.0]
+        f = [4.0, 0.0]
+
+        points = [a, b, c, d, e, f]
+
+        #bac, bce, ecf, dbe
+        elements = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
+
+        mesh4 = Generic_Domain(points, elements)
+                     #  geo_reference = Geo_reference(56, x0, y0))
+        mesh4.check_integrity()
+        quantity = Quantity(mesh4)
 
 
         ncols = 11  # Nx
         nrows = 12  # Ny
-        xllcorner = 240000
-        yllcorner = 7620000
-        cellsize  = 6000
+        xllcorner = x0
+        yllcorner = y0
+        cellsize  = 1.0
         NODATA_value =  -9999
         
@@ -1457 +1472 @@
 
         # check order of vertices
-        answer =  [ 23298000. , 23358000. , 23118000. ]
-        answer = [ 23298000. ,  23118000. , 23358000 ]
-        print quantity.vertex_values
+
+        answer = [[  6.,   0.,   2.],
+                  [  6.,   2.,   8.],
+                  [  8.,  2.,   4.],
+                  [ 12.,   6.,   8.]]
+        #print quantity.vertex_values
         assert num.allclose(quantity.vertex_values, answer)
         
@@ -1478 +1496 @@
         #print quantity.centroid_values      
         
-        answer =  [ 23258000. ]
+
+        answer = [ 2.66666667,  5.33333333,  4.66666667,  8.66666667]
+        
         assert num.allclose(quantity.centroid_values, answer)
         #Cleanup

trunk/anuga_core/user_manual/anuga_user_manual.tex

@@ -2732 +2732 @@
 
 The \code{Culvert_flow} class takes as input:
-\begin{itemize}
-  \item \code{domain}: a reference to the domain being evolved
-  \item \code{culvert_description_filename}:
-  \item \code{culvert_routine}:
-  \item \code{end_point0}: Coordinates of one opening
-  \item \code{end_point1}: Coordinates of other opening
-  \item \code{enquiry_point0}:
-  \item \code{enquiry_point1}:
-  \item \code{type}: (default is 'box')
-  \item \code{width}:
-  \item \code{height}:
-  \item \code{length}:
-  \item \code{number_of_barrels}: Number of identical pipes in the culvert (default is 1)
-  \item \code{trigger_depth}: (default is 0.01)
-  \item \code{manning}: Mannings Roughness for Culvert
-  \item \code{sum_loss}:
-  \item \code{use_velocity_head}:
-  \item \code{use_momentum_jet}:
-  \item \code{label}: Short text naming the culvert
-  \item \code{description}: Text describing the culvert
-  \item \code{update_interval}:
-  \item \code{log_file}:
-  \item \code{discharge_hydrograph}:
-\end{itemize}
+
+ \code{domain}: a reference to the domain being evolved
+
+\code{culvert_description_filename}:
+
+\code{culvert_routine}:
+
+  \code{end_point0}: Coordinates of one opening
+
+  \code{end_point1}: Coordinates of other opening
+
+  \code{enquiry_point0}:
+
+ \code{enquiry_point1}:
+ 
+ \code{type}: (default is 'box')
+ 
+ \code{width}:
+
+  \code{height}:
+ 
+ \code{length}:
+
+  \code{number_of_barrels}: Number of identical pipes in the culvert (default is 1)
+ 
+ \code{trigger_depth}: (default is 0.01)
+
+  \code{manning}: Mannings Roughness for Culvert
+ 
+ \code{sum_loss}:
+ 
+ \code{use_velocity_head}:
+
+  \code{use_momentum_jet}:
+ 
+   \code{label}: Short text naming the culvert
+
+  \code{description}: Text describing the culvert
+
+    \code{update_interval}:
+ 
+ \code{log_file}:
+
+  \code{discharge_hydrograph}:
+
 
 The user can specify different culvert routines. However, \anuga currently provides only one, namely the
@@ -2792 +2813 @@
 
 
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\section{Fractional Step Operators}\index{fractional step operators}
+
+An alternative way to effect the evolution is via fractional step operators.  The idea is that at each inner timestep we can use 
+a fractional step operator to change the centroid values of our conserved quantities, the stage and the xmomentum and the ymomentum. The elevation can also be changed but more care needs to be applied as the \code{elevation} quantity needs to remain continuous. The moral, play with \code{elevation} at your peril. 
+
+We are tending to move away from forcing terms and towards fractpional step operators as they are easier to parallelize. 
+
+
+Currently predefined fractional step operators:
+
+
+\begin{classdesc}{Set_elevation_operator}{domain,
+                 elevation=None,
+                 indices=None,
+                 polygon=None,
+                 center=None,
+                 radius=None,
+                 description = None,
+                 label = None,
+                 logging = False,
+                 verbose = False}
+Module: \module{anuga.operators.set_elevation_operator}
+
+Set the elevation in a region (careful to maintain continuitiy of elevation)
+
+
+\code{domain} is the domain being evolved.
+
+\code{elevation} a function of $t$ or $x,y$ or $x,y,t$ written to take numpy arrays $x$ and $y$. $t$ is a scalar.
+
+\code{indices} is a list of triangle indices where the function \code{elevation} will be applied.
+
+\code{polygon} is a polygon. The function \code{elevation} will be applied to triangles with in the polygon.
+
+\code{center} is the center of a circle where the function \code{elevation} will be applied. (\code{radius} needs to be set).
+
+\code{radius} is the radius of a circle where the function \code{elevation} will be applied. (\code{center} needs to be set).
+
+\code{description} is a description of this operator.
+
+\code{label} is the lable used when reporting on the operator.
+
+\code{logging} is the boolean flag to set logging to a file. 
+
+\code{verbose} is the boolean flag to setup extended output from the operator. 
+
+Example:
+
+\begin{verbatim}
+from anuga.operators.set_elevation_operator import Set_elevation_operator
+
+op0 =Set_elevation_operator(domain, elevation = lambda t : 0.01*t)
+\end{verbatim}
+would setup an operator which raises the elevation over the whole domain 1 cm per second.
+
+
+While 
+\begin{verbatim}
+from anuga.operators.set_elevation_operator import Set_elevation_operator
+P = [[x0, y0], [x1, y0], [x1, y1], [x0, y1]] 
+
+op0 =Set_elevation_operator(domain, elevation = lambda t : 0.01*t, polygon=P)
+\end{verbatim}
+would setup an operator which raises the elevation over the polygon \code{P} 1 cm per second.
+
+\end{classdesc}
+
+
+
+\subsection{User developed fractional step operator}
+
+Suppose we want to add water to our domain at a steady rate. we can create a class based on the \code{Operator} class which has a \code{__call__} function to do the required update to the centroid values of the stage variable.
+
+Here is some code to do this:
+\begin{verbatim}
+from anuga.operators.base_operator import Operator
+
+class My_operator(Operator):
+    """
+    Add water to the domain at a given rate (m/sec)
+    """
+
+    def __init__(self,
+                 domain,
+                 rate=0.0,
+                 description = None,
+                 label = None,
+                 logging = False,
+                 verbose = False):
+          """Initialize the user defined operator
+          """
+
+        Operator.__init__(self, domain, description, label, logging, verbose)
+       
+        self.rate = rate
+
+
+    def __call__(self):
+        """
+        Apply rate to all triangles 
+        """
+
+        timestep = self.domain.get_timestep()
+
+        self.stage_c[:] = self.stage_c[:]  + self.rate*timestep
+\end{verbatim}
+
+And then in your script you would asosciated \code{My_operator} with the \code{domain} with a call 
+\begin{verbatim}
+My_operator(domain,rate=1.0)
+\end{verbatim}
+
+
+
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 \section{Evolution}\index{evolution}
 \label{sec:evolution}