Changeset 8135

Timestamp:

Mar 10, 2011, 3:31:46 PM (14 years ago)

Author:

steve

Message:

Changing over kv to use quantities

Location:

trunk/anuga_core/source/anuga/operators

Files:

• kinematic_viscosity.py (modified) (7 diffs)
• kinematic_viscosity_ext.c (modified) (9 diffs)
• test_kinematic_viscosity.py (modified) (2 diffs)

trunk/anuga_core/source/anuga/operators/kinematic_viscosity.py

@@ -1 +1 @@
 
 from anuga import Domain
+from anuga import Quantity
 from anuga.utilities.sparse import Sparse, Sparse_CSR
 from anuga.utilities.cg_solve import conjugate_gradient
 import anuga.abstract_2d_finite_volumes.neighbour_mesh as neighbour_mesh
-from anuga.abstract_2d_finite_volumes.generic_boundary_conditions import Dirichlet_boundary
+from anuga import Dirichlet_boundary
 import numpy as num
 import kinematic_viscosity_ext
@@ -11 +12 @@
 class Kinematic_Viscosity:
     """
-    Class for setting up structures and matrices for kinematic viscosity
+    Class for setting up structures and matrices for kinematic viscosity differential
+    operator using centroid values.
+
+    div ( diffusivity grad )
+
+    where diffusvity is scalar quantity (defaults to quantity with values = 1)
+    boundary values of f are used to setup entries associated with cells with boundaries
+
+    There are procedures to apply this operator, ie
+
+    (1) Calculate div( diffusivity grad u )
+    using boundary values stored in u
+
+    (2) Calculate ( u + dt div( diffusivity grad u )
+    using boundary values stored in u
+
+    (3) Solve div( diffusivity grad u ) = f
+    for quantity f and using boundary values stored in u
+
+    (4) Solve ( u + dt div( diffusivity grad u ) = f
+    for quantity f using boundary values stored in u
+
     """
 
-    def __init__(self, domain, triangle_areas=True, verbose=False):
+    def __init__(self, domain, diffusivity = None, triangle_areas=True, verbose=False):
         if verbose: log.critical('Kinematic Viscosity: Beginning Initialisation')
         #Expose the domain attributes
@@ -21 +43 @@
         self.mesh = domain.mesh
         self.boundary = domain.boundary
+
+        # Setup diffusivity quantity
+        if diffusivity is None:
+            diffusivity = Quantity(domain)
+            diffusivity.set_values(1.0)
+        #check that diffusivity is a quantity associated with domain
+        assert diffusivity.domain == domain
+
+        self.diffusivity = diffusivity
+        self.diffusivity_bdry_data = self.diffusivity.boundary_values
+        self.diffusivity_data = self.diffusivity.centroid_values
+
+
         self.n = len(self.domain)
         self.dt = 1e-6 #Need to set to domain.timestep
         self.boundary_len = len(domain.boundary)
         self.tot_len = self.n + self.boundary_len
+
+        # FIXME SR: maybe this should already exist in mesh!
         self.boundary_enum = self.enumerate_boundary()
+        
         self.verbose = verbose
 
         #Geometric Information
         if verbose: log.critical('Kinematic Viscosity: Building geometric structure')
+
         self.geo_structure_indices = num.zeros((self.n, 3), num.int)
         self.geo_structure_values = num.zeros((self.n, 3), num.float)
+
+        # Only needs to built once, doesn't change
         kinematic_viscosity_ext.build_geo_structure(self, self.n, self.tot_len)
+
+        # Setup type of scaling
         self.apply_triangle_areas = triangle_areas
-        self.triangle_areas = Sparse(self.n, self.n)
+
+        # FIXME SR: should this really be a matrix?
+        temp  = Sparse(self.n, self.n)
         for i in range(self.n):
-            self.triangle_areas[i, i] = 1.0 / self.mesh.areas[i]
+            temp[i, i] = 1.0 / self.mesh.areas[i]
             
-        self.triangle_areas = Sparse_CSR(self.triangle_areas)
-
-        #Application Information
+        self.triangle_areas = Sparse_CSR(temp)
+
+        # FIXME SR: More to do with solving equation
         self.qty_considered = 1 #1 or 2 (uh or vh respectively)
+
+        #FIXME SR: Do we need Sparse or should we just go with Sparse_CSR
         self.operator_matrix = Sparse(self.n, self.tot_len)
 
@@ -49 +96 @@
         self.operator_colind = num.zeros((4 * self.n, ), num.int)
         #Sparse_CSR.rowptr (4 entries in every row, we know this already) = [0,4,8,...,4*n]
-        self.operator_rowptr = 4 * num.indices((self.n + 1, ))[0, :]
-        self.stage_heights = num.zeros((self.n, 1), num.float)
-        self.stage_heights_scaling = Sparse(self.n, self.n)
-        self.boundary_vector = num.zeros((self.n, ), num.float)
+        self.operator_rowptr = 4 * num.arange(self.n + 1)
+
+        # Build matrix self.elliptic_operator_matrix [A B]
+        self.build_elliptic_operator()
+
+        # Build self.boundary_term
+        self.build_operator_boundary_term()
 
         self.parabolic_solve = False #Are we doing a parabolic solve at the moment?
@@ -61 +111 @@
     def enumerate_boundary(self):
         #Enumerate by boundary index
+        # FIXME SR: Probably should be part of mesh
         enumeration = {}
         for i in range(self.n):
@@ -70 +121 @@
 
     def set_qty_considered(self, qty):
+        # FIXME SR: Probably should just be set by quantity to which operation is applied
         if qty == 1 or qty == 'u':
             self.qty_considered = 1
@@ -78 +130 @@
             assert 0 == 1, msg
 
-    def apply_stage_heights(self, h):
-        msg = "(KV_Operator.apply_stage_heights) h vector has incorrect length"
-        assert h.size == self.n, msg
-
-        self.operator_matrix = Sparse(self.n, self.tot_len)
-
-        #Evaluate the boundary stage heights - use generic_domain.update_boundary? (check enumeration)
-        boundary_heights = num.zeros((self.boundary_len, 1), num.float)
-        for key in self.boundary_enum.keys():
-            boundary_heights[self.boundary_enum[key]] = self.boundary[key].evaluate()[0]
-
-        kinematic_viscosity_ext.build_operator_matrix(self, self.n, self.tot_len, h, boundary_heights)
-        self.operator_matrix = \
-            Sparse_CSR(None, self.operator_data, self.operator_colind, self.operator_rowptr, self.n, self.tot_len)
-
-        self.stage_heights = h
-        #Set up the scaling matrix
-        data = h
-        num.putmask(data, data != 0, 1 / data) #take the reciprocal of each entry unless it is zero
-        self.stage_heights_scaling = \
-            Sparse_CSR(None, data, num.arange(self.n), num.append(num.arange(self.n), self.n), self.n, self.n)
-
-        self.build_boundary_vector()
-
-        return self.operator_matrix
-
-    def build_boundary_vector(self):
-        #Require stage heights applied
-        X = num.zeros((self.tot_len, 2), num.float)
-        for key in self.boundary_enum.keys():
-            quantities = self.boundary[key].evaluate()
-            h_i = quantities[0]
-            if h_i == 0.0:
-                X[self.n + self.boundary_enum[key], :] = num.array([0.0, 0.0])
-            else:
-                X[self.n + self.boundary_enum[key], :] = quantities[1:] / h_i
-        self.boundary_vector = self.operator_matrix * X
+    def build_elliptic_operator(self):
+        """
+        Builds matrix representing
+
+        div ( diffusivity grad )
+
+        which has the form [ A B ]
+        """
+
+        #Arrays self.operator_data, self.operator_colind, self.operator_rowptr
+        # are setup via this call
+        kinematic_viscosity_ext.build_operator_matrix(self,self.diffusivity_data, self.diffusivity_bdry_data)
+
+        self.elliptic_operator_matrix = Sparse_CSR(None, \
+                self.operator_data, self.operator_colind, self.operator_rowptr, \
+                self.n, self.tot_len)
+
+#        #Set up the scaling matrix
+#        data = h
+#        num.putmask(data, data != 0, 1 / data) #take the reciprocal of each entry unless it is zero
+#        self.stage_heights_scaling = \
+#            Sparse_CSR(None, self.diffusivity_data, num.arange(self.n), num.arange(self.n +1), self.n, self.n)
+
+    def update_elliptic_operator(self):
+        """
+        Updates the data values of matrix representing
+
+        div ( diffusivity grad )
+
+        (as when diffusivitiy has changed)
+        """
+
+        #Array self.operator_data is changed by this call, which should flow
+        # through to the Sparse_CSR matrix.
+
+        kinematic_viscosity_ext.update_operator_matrix(self, self.n, self.tot_len, \
+          self.diffusivity_data, self.diffusivity_bdry_data)
+
+
+
+
+    def build_operator_boundary_term(self):
+        """
+        Operator has form [A B] and U = [ u ; b]
+
+        This procedure calculates B b which can be calculated as
+
+        [A B] [ 0 ; b]
+        """
+
+        n = self.n
+        tot_len = self.tot_len
+
+        self.boundary_term = num.zeros((n, ), num.float)
+        
+        X = num.zeros((tot_len,), num.float)
+
+        X[n:] = self.diffusivity_bdry_data
+        self.boundary_term[:] = self.operator_matrix * X
+
         #Tidy up
         if self.apply_triangle_areas:
-            self.boundary_vector = self.triangle_areas * self.boundary_vector
-
-        return self.boundary_vector
+            self.boundary_term[:] = self.triangle_areas * self.boundary_term
+
 
     def elliptic_multiply(self, V, qty_considered=None, include_boundary=True):

trunk/anuga_core/source/anuga/operators/kinematic_viscosity_ext.c

@@ -85 +85 @@
 }
 
-int build_operator_matrix(int n, 
-                          int tot_len, 
-                          long *geo_indices, 
-                          double *geo_values, 
-                          double *h, 
-                          double *boundary_heights, 
-                          double *data, 
+
+int build_operator_matrix(int n,
+                          int tot_len,
+                          long *geo_indices,
+                          double *geo_values,
+                          double *cell_data,
+                          double *bdry_data,
+                          double *data,
                           long *colind) {
         int i, k, edge, j[4], *sorted_j, this_index;
@@ -102 +103 @@
                         j[edge] = geo_indices[3*i+edge];
                         if (j[edge]<n) { //interior
-                                h_j = h[j[edge]];
+                                h_j = cell_data[j[edge]];
                         } else { //boundary
-                                h_j = boundary_heights[j[edge]-n];
-                        }
-                        v[edge] = -0.5*(h[i] + h_j)*geo_values[3*i+edge]; //the negative of the individual interaction
-                        v_i += 0.5*(h[i] + h_j)*geo_values[3*i+edge]; //sum the three interactions
+                                h_j = bdry_data[j[edge]-n];
+                        }
+                        v[edge] = -0.5*(cell_data[i] + h_j)*geo_values[3*i+edge]; //the negative of the individual interaction
+                        v_i += 0.5*(cell_data[i] + h_j)*geo_values[3*i+edge]; //sum the three interactions
                 }
                 //Organise the set of 4 values/indices into the data and colind arrays
@@ -131 +132 @@
 }
 
+int update_operator_matrix(int n,
+                          int tot_len,
+                          long *geo_indices,
+                          double *geo_values,
+                          double *cell_data,
+                          double *bdry_data,
+                          double *data,
+                          long *colind) {
+        int i, k, edge, j[4], *sorted_j, this_index;
+        double h_j, v[3], v_i; //v[k] = value of the interaction of edge k in a given triangle, v_i = (i,i) entry
+        for (i=0; i<n; i++) {
+                v_i = 0.0;
+                j[3] = i;
+                //Get the values of each interaction, and the column index at which they occur
+                for (edge=0; edge<3; edge++) {
+                        j[edge] = geo_indices[3*i+edge];
+                        if (j[edge]<n) { //interior
+                                h_j = cell_data[j[edge]];
+                        } else { //boundary
+                                h_j = bdry_data[j[edge]-n];
+                        }
+                        v[edge] = -0.5*(cell_data[i] + h_j)*geo_values[3*i+edge]; //the negative of the individual interaction
+                        v_i += 0.5*(cell_data[i] + h_j)*geo_values[3*i+edge]; //sum the three interactions
+                }
+                //Organise the set of 4 values/indices into the data and colind arrays
+                sorted_j = quicksort((int *)j,4);
+                for (k=0; k<4; k++) { //loop through the nonzero indices
+                        this_index = sorted_j[k];
+                        if (this_index == i) {
+                                data[4*i+k] = v_i;
+                                colind[4*i+k] = i;
+                        } else if (this_index == j[0]) {
+                                data[4*i+k] = v[0];
+                                colind[4*i+k] = j[0];
+                        } else if (this_index == j[1]) {
+                                data[4*i+k] = v[1];
+                                colind[4*i+k] = j[1];
+                        } else { //this_index == j[2]
+                                data[4*i+k] = v[2];
+                                colind[4*i+k] = j[2];
+                        }
+                }
+        }
+        return 0;
+}
+
 /**
 * Python wrapper methods
@@ -147 +194 @@
     
     //Convert Python arguments to C
-    if (!PyArg_ParseTuple(args, "Oii", &kv_operator, &n, &tot_len)) {
+    if (!PyArg_ParseTuple(args, "Oii", &kv_operator)) {
         PyErr_SetString(PyExc_RuntimeError, "build_geo_structure could not parse input");
         return NULL;
@@ -158 +205 @@
 
     //Extract parameters
+    n       = get_python_integer(kv_operator,"n");
+    tot_len = get_python_integer(kv_operator,"tot_len");
+
     centroid_coordinates = get_consecutive_array(mesh,"centroid_coordinates");
     neighbours = get_consecutive_array(mesh,"neighbours");
@@ -195 +245 @@
         int n, tot_len, err;
         PyArrayObject
-                *h,
-                *boundary_heights,
+                *cell_data,
+                *bdry_data,
                 *geo_indices,
                 *geo_values,
@@ -203 +253 @@
         
         //Convert Python arguments to C
-    if (!PyArg_ParseTuple(args, "OiiOO", &kv_operator, &n, &tot_len, &h, &boundary_heights)) {
-        PyErr_SetString(PyExc_RuntimeError, "get_stage_height_interactions could not parse input");
-        return NULL;
-    }
-        
+    if (!PyArg_ParseTuple(args, "OOO", &kv_operator, &cell_data, &bdry_data)) {
+        PyErr_SetString(PyExc_RuntimeError, "build_operator_matrix could not parse input");
+        return NULL;
+    }
+
+
+    n       = get_python_integer(kv_operator,"n");
+    tot_len = get_python_integer(kv_operator,"tot_len");
+
+
         geo_indices = get_consecutive_array(kv_operator,"geo_structure_indices");
         geo_values = get_consecutive_array(kv_operator,"geo_structure_values");
@@ -216 +271 @@
                                                                 (long *)geo_indices -> data, 
                                                                 (double *)geo_values -> data, 
-                                                                (double *)h -> data, 
-                                                                (double *)boundary_heights -> data, 
+                                                                (double *)cell_data -> data,
+                                                                (double *)bdry_data -> data,
                                                                 (double *)_data -> data, 
                                                                 (long *)colind -> data);
@@ -233 +288 @@
 }
 
+
+static PyObject *py_update_operator_matrix(PyObject *self, PyObject *args) {
+        PyObject *kv_operator;
+        int n, tot_len, err;
+        PyArrayObject
+                *cell_data,
+                *bdry_data,
+                *geo_indices,
+                *geo_values,
+                *_data,
+                *colind;
+
+        //Convert Python arguments to C
+    if (!PyArg_ParseTuple(args, "OOO", &kv_operator, &cell_data, &bdry_data)) {
+        PyErr_SetString(PyExc_RuntimeError, "update_operator_matrix could not parse input");
+        return NULL;
+    }
+
+
+    n       = get_python_integer(kv_operator,"n");
+    tot_len = get_python_integer(kv_operator,"tot_len");
+
+
+        geo_indices = get_consecutive_array(kv_operator,"geo_structure_indices");
+        geo_values = get_consecutive_array(kv_operator,"geo_structure_values");
+        _data = get_consecutive_array(kv_operator,"operator_data");
+        colind = get_consecutive_array(kv_operator,"operator_colind");
+
+        err = update_operator_matrix(n,tot_len,
+                                                                (long *)geo_indices -> data,
+                                                                (double *)geo_values -> data,
+                                                                (double *)cell_data -> data,
+                                                                (double *)bdry_data -> data,
+                                                                (double *)_data -> data,
+                                                                (long *)colind -> data);
+    if (err != 0) {
+        PyErr_SetString(PyExc_RuntimeError, "Could not get stage height interactions");
+        return NULL;
+    }
+
+        Py_DECREF(geo_indices);
+        Py_DECREF(geo_values);
+        Py_DECREF(_data);
+        Py_DECREF(colind);
+
+    return Py_BuildValue("");
+
+}
+
 static struct PyMethodDef MethodTable[] = {
         {"build_geo_structure",py_build_geo_structure,METH_VARARGS,"Print out"},
                 {"build_operator_matrix",py_build_operator_matrix,METH_VARARGS,"Print out"},
+        {"update_operator_matrix",py_update_operator_matrix,METH_VARARGS,"Print out"},
         {NULL,NULL,0,NULL} // sentinel
 };

trunk/anuga_core/source/anuga/operators/test_kinematic_viscosity.py

@@ -29 +29 @@
         domain.set_boundary({'edge0': D0, 'edge1': D1, 'edge2': D2})
 
+        domain.set_quantity('stage', lambda x,y : x+2*y )
+        domain.set_quantity('elevation', lambda x,y : 3*x+5*y )
+
+
+        #print domain.quantities['stage'].vertex_values
+
+        #print domain.quantities['stage'].edge_values
+
         domain.update_boundary()
 
 
-        print domain.quantities['elevation'].boundary_values
+        #print domain.quantities['stage'].boundary_values
         
         return Kinematic_Viscosity(domain)
@@ -88 +96 @@
         assert num.allclose(values, num.array([[-3.0,-6.0/sqrt(5),-6.0/sqrt(5)],[-6.0/sqrt(5),-3.0,-6.0/sqrt(5)]]))
 
-    def test_apply_heights(self):
+    def test_elliptic_matrix(self):
         operator1 = self.operator1()
         operator2 = self.operator2()
 
-        h = num.array([1.0])
-        A = operator1.apply_stage_heights(h)
+        domain1 = operator1.domain
+        domain2 = operator2.domain
+
+
+        A = operator1.elliptic_operator_matrix
+
+        print A
+        print num.array([-6.0-12.0/sqrt(5), 6.0,  6.0/sqrt(5), 6.0/sqrt(5)])
+        
         assert num.allclose(A.todense(), num.array([-6.0-12.0/sqrt(5), 6.0,  6.0/sqrt(5), 6.0/sqrt(5)]))
 
-        h = num.array([2.0])
-        A = operator1.apply_stage_heights(h)
+
+        A = operator1.elliptic_operator_matrix
         assert num.allclose(A.todense(), 1.5*num.array([-6.0-12.0/sqrt(5), 6.0, 6.0/sqrt(5), 6.0/sqrt(5)]))