Changeset 3678


Ignore:
Timestamp:
Oct 1, 2006, 6:55:46 PM (18 years ago)
Author:
steve
Message:

Added more limiting to cells near dry cells, use beta_*_dry

Files:
17 edited

Legend:

Unmodified
Added
Removed
  • anuga_core/source/anuga/abstract_2d_finite_volumes/domain.py

    r3631 r3678  
    166166        self.set_default_order(1)
    167167        #self.default_order = 1
    168         #self.order = self.default_order
     168        #self._order_ = self.default_order
    169169
    170170        self.smallsteps = 0
     
    222222
    223223        self.default_order = n
    224         self.order = self.default_order
     224        self._order_ = self.default_order
    225225
    226226
     
    695695            yieldstep = float(yieldstep)
    696696
    697         self.order = self.default_order
     697        self._order_ = self.default_order
    698698
    699699
     
    759759            self.yieldtime += self.timestep
    760760            self.number_of_steps += 1
    761             if self.order == 1:
     761            if self._order_ == 1:
    762762                self.number_of_first_order_steps += 1
    763763
     
    846846                self.smallsteps = 0 #Reset
    847847
    848                 if self.order == 1:
     848                if self._order_ == 1:
    849849                    msg = 'WARNING: Too small timestep %.16f reached '\
    850850                          %timestep
     
    857857                else:
    858858                    #Try to overcome situation by switching to 1 order
    859                     self.order = 1
     859                    self._order_ = 1
    860860
    861861        else:
    862862            self.smallsteps = 0
    863             if self.order == 1 and self.default_order == 2:
    864                 self.order = 2
     863            if self._order_ == 1 and self.default_order == 2:
     864                self._order_ = 2
    865865
    866866
     
    933933        for name in self.conserved_quantities:
    934934            Q = self.quantities[name]
    935             if self.order == 1:
     935            if self._order_ == 1:
    936936                Q.extrapolate_first_order()
    937             elif self.order == 2:
     937            elif self._order_ == 2:
    938938                Q.extrapolate_second_order()
    939939                Q.limit()
  • anuga_core/source/anuga/abstract_2d_finite_volumes/mesh_factory.py

    r3560 r3678  
    334334    """
    335335
    336     from anuga.abstract_2d_finite_volumes.util import anglediff
     336    from anuga.utilities.numerical_tools import anglediff
    337337    from math import pi
    338338    import os.path
     
    449449
    450450    from math import pi
    451     from anuga.abstract_2d_finite_volumes.util import anglediff
     451    from anuga.utilities.numerical_tools import anglediff
    452452
    453453
  • anuga_core/source/anuga/abstract_2d_finite_volumes/quantity.py

    r3560 r3678  
    12341234            #Get index of the one neighbour
    12351235            for k0 in surrogate_neighbours[k,:]:
    1236                 if k0 != k: break
     1236                if k0 != k: break
    12371237            assert k0 != k
    12381238
     
    12851285
    12861286    for k in range(N):
    1287         qmax[k] = qmin[k] = qc[k]
     1287        qmax[k] = qc[k]
     1288        qmin[k] = qc[k]
    12881289        for i in range(3):
    12891290            n = quantity.domain.neighbours[k,i]
     
    12931294                qmin[k] = min(qmin[k], qn)
    12941295                qmax[k] = max(qmax[k], qn)
     1296        qmax[k] = min(qmax[k], 2.0*qc[k])
     1297        qmin[k] = max(qmin[k], 0.5*qc[k])
    12951298
    12961299
     
    13261329    #Replace python version with c implementations
    13271330
    1328     from quantity_ext import limit, compute_gradients,\
     1331    from quantity_ext import compute_gradients, limit,\
    13291332    extrapolate_second_order, interpolate_from_vertices_to_edges, update
  • anuga_core/source/anuga/abstract_2d_finite_volumes/quantity_ext.c

    r2676 r3678  
    2626                        long* surrogate_neighbours,
    2727                        double* a,
    28                         double* b) {
    29 
    30   int i, k, k0, k1, k2, index3;
    31   double x0, x1, x2, y0, y1, y2, q0, q1, q2; //, det;
    32 
    33 
    34   for (k=0; k<N; k++) {
    35     index3 = 3*k;
    36 
    37     if (number_of_boundaries[k] < 2) {
    38       //Two or three true neighbours
    39 
    40       //Get indices of neighbours (or self when used as surrogate)
    41       //k0, k1, k2 = surrogate_neighbours[k,:]
    42 
    43       k0 = surrogate_neighbours[index3 + 0];
    44       k1 = surrogate_neighbours[index3 + 1];
    45       k2 = surrogate_neighbours[index3 + 2];
    46 
    47 
    48       if (k0 == k1 || k1 == k2) return -1;
    49 
    50       //Get data
    51       q0 = centroid_values[k0];
    52       q1 = centroid_values[k1];
    53       q2 = centroid_values[k2];
    54 
    55       x0 = centroids[k0*2]; y0 = centroids[k0*2+1];
    56       x1 = centroids[k1*2]; y1 = centroids[k1*2+1];
    57       x2 = centroids[k2*2]; y2 = centroids[k2*2+1];
    58 
    59       //Gradient
    60       _gradient(x0, y0, x1, y1, x2, y2, q0, q1, q2, &a[k], &b[k]);
    61 
    62     } else if (number_of_boundaries[k] == 2) {
    63       //One true neighbour
    64 
    65       //#Get index of the one neighbour
    66       i=0; k0 = k;
    67       while (i<3 && k0==k) {
    68         k0 = surrogate_neighbours[index3 + i];
    69         i++;
    70       }
    71       if (k0 == k) return -1;
    72 
    73       k1 = k; //self
    74 
    75       //Get data
    76       q0 = centroid_values[k0];
    77       q1 = centroid_values[k1];
    78 
    79       x0 = centroids[k0*2]; y0 = centroids[k0*2+1];
    80       x1 = centroids[k1*2]; y1 = centroids[k1*2+1];
    81 
    82       //Two point gradient
    83       _gradient2(x0, y0, x1, y1, q0, q1, &a[k], &b[k]);
    84 
    85 
    86       //Old (wrong code)
    87       //det = x0*y1 - x1*y0;
    88       //if (det != 0.0) {
    89       //        a[k] = (y1*q0 - y0*q1)/det;
    90       //        b[k] = (x0*q1 - x1*q0)/det;
    91       //}
    92     }
    93     //    else:
    94     //        #No true neighbours -
    95     //        #Fall back to first order scheme
    96     //        pass
    97   }
    98   return 0;
     28                        double* b){
     29
     30        int i, k, k0, k1, k2, index3;
     31        double x0, x1, x2, y0, y1, y2, q0, q1, q2; //, det;
     32
     33
     34        for (k=0; k<N; k++) {
     35                index3 = 3*k;
     36
     37                if (number_of_boundaries[k] < 2) {
     38                        //Two or three true neighbours
     39
     40                        //Get indices of neighbours (or self when used as surrogate)
     41                        //k0, k1, k2 = surrogate_neighbours[k,:]
     42
     43                        k0 = surrogate_neighbours[index3 + 0];
     44                        k1 = surrogate_neighbours[index3 + 1];
     45                        k2 = surrogate_neighbours[index3 + 2];
     46
     47
     48                        if (k0 == k1 || k1 == k2) return -1;
     49
     50                        //Get data
     51                        q0 = centroid_values[k0];
     52                        q1 = centroid_values[k1];
     53                        q2 = centroid_values[k2];
     54
     55                        x0 = centroids[k0*2]; y0 = centroids[k0*2+1];
     56                        x1 = centroids[k1*2]; y1 = centroids[k1*2+1];
     57                        x2 = centroids[k2*2]; y2 = centroids[k2*2+1];
     58
     59                        //Gradient
     60                        _gradient(x0, y0, x1, y1, x2, y2, q0, q1, q2, &a[k], &b[k]);
     61
     62                } else if (number_of_boundaries[k] == 2) {
     63                        //One true neighbour
     64
     65                        //Get index of the one neighbour
     66                        i=0; k0 = k;
     67                        while (i<3 && k0==k) {
     68                                k0 = surrogate_neighbours[index3 + i];
     69                                i++;
     70                        }
     71                        if (k0 == k) return -1;
     72
     73                        k1 = k; //self
     74
     75                        //Get data
     76                        q0 = centroid_values[k0];
     77                        q1 = centroid_values[k1];
     78
     79                        x0 = centroids[k0*2]; y0 = centroids[k0*2+1];
     80                        x1 = centroids[k1*2]; y1 = centroids[k1*2+1];
     81
     82                        //Two point gradient
     83                        _gradient2(x0, y0, x1, y1, q0, q1, &a[k], &b[k]);
     84
     85
     86                        //Old (wrong code)
     87                        //det = x0*y1 - x1*y0;
     88                        //if (det != 0.0) {
     89                        //      a[k] = (y1*q0 - y0*q1)/det;
     90                        //      b[k] = (x0*q1 - x1*q0)/det;
     91                        //}
     92                }
     93                //    else:
     94                //        #No true neighbours -
     95                //        #Fall back to first order scheme
     96                //        pass
     97        }
     98        return 0;
    9999}
    100100
     
    108108                 double* b) {
    109109
    110   int k, k2, k3, k6;
    111   double x, y, x0, y0, x1, y1, x2, y2;
    112 
    113   for (k=0; k<N; k++) {
    114     k6 = 6*k;
    115     k3 = 3*k;
    116     k2 = 2*k;
    117 
    118     //Centroid coordinates
    119     x = centroids[k2]; y = centroids[k2+1];
    120 
    121     //vertex coordinates
    122     //x0, y0, x1, y1, x2, y2 = X[k,:]
    123     x0 = vertex_coordinates[k6 + 0];
    124     y0 = vertex_coordinates[k6 + 1];
    125     x1 = vertex_coordinates[k6 + 2];
    126     y1 = vertex_coordinates[k6 + 3];
    127     x2 = vertex_coordinates[k6 + 4];
    128     y2 = vertex_coordinates[k6 + 5];
    129 
    130     //Extrapolate
    131     vertex_values[k3+0] = centroid_values[k] + a[k]*(x0-x) + b[k]*(y0-y);
    132     vertex_values[k3+1] = centroid_values[k] + a[k]*(x1-x) + b[k]*(y1-y);
    133     vertex_values[k3+2] = centroid_values[k] + a[k]*(x2-x) + b[k]*(y2-y);
    134 
    135   }
    136   return 0;
     110        int k, k2, k3, k6;
     111        double x, y, x0, y0, x1, y1, x2, y2;
     112
     113        for (k=0; k<N; k++){
     114                k6 = 6*k;
     115                k3 = 3*k;
     116                k2 = 2*k;
     117
     118                //Centroid coordinates
     119                x = centroids[k2]; y = centroids[k2+1];
     120
     121                //vertex coordinates
     122                //x0, y0, x1, y1, x2, y2 = X[k,:]
     123                x0 = vertex_coordinates[k6 + 0];
     124                y0 = vertex_coordinates[k6 + 1];
     125                x1 = vertex_coordinates[k6 + 2];
     126                y1 = vertex_coordinates[k6 + 3];
     127                x2 = vertex_coordinates[k6 + 4];
     128                y2 = vertex_coordinates[k6 + 5];
     129
     130                //Extrapolate
     131                vertex_values[k3+0] = centroid_values[k] + a[k]*(x0-x) + b[k]*(y0-y);
     132                vertex_values[k3+1] = centroid_values[k] + a[k]*(x1-x) + b[k]*(y1-y);
     133                vertex_values[k3+2] = centroid_values[k] + a[k]*(x2-x) + b[k]*(y2-y);
     134
     135        }
     136        return 0;
    137137}
    138138
     
    144144                 double* edge_values) {
    145145
    146   int k, k3;
    147   double q0, q1, q2;
    148 
    149 
    150   for (k=0; k<N; k++) {
    151     k3 = 3*k;
    152 
    153     q0 = vertex_values[k3 + 0];
    154     q1 = vertex_values[k3 + 1];
    155     q2 = vertex_values[k3 + 2];
    156 
    157     //printf("%f, %f, %f\n", q0, q1, q2);
    158     edge_values[k3 + 0] = 0.5*(q1+q2);
    159     edge_values[k3 + 1] = 0.5*(q0+q2);
    160     edge_values[k3 + 2] = 0.5*(q0+q1);
    161   }
    162   return 0;
     146        int k, k3;
     147        double q0, q1, q2;
     148
     149
     150        for (k=0; k<N; k++) {
     151                k3 = 3*k;
     152
     153                q0 = vertex_values[k3 + 0];
     154                q1 = vertex_values[k3 + 1];
     155                q2 = vertex_values[k3 + 2];
     156
     157                //printf("%f, %f, %f\n", q0, q1, q2);
     158                edge_values[k3 + 0] = 0.5*(q1+q2);
     159                edge_values[k3 + 1] = 0.5*(q0+q2);
     160                edge_values[k3 + 2] = 0.5*(q0+q1);
     161        }
     162        return 0;
    163163}
    164164
     
    168168            double* explicit_update,
    169169            double* semi_implicit_update) {
    170   //Update centroid values based on values stored in
    171   //explicit_update and semi_implicit_update as well as given timestep
    172 
    173 
    174   int k;
    175   double denominator, x;
    176 
    177 
    178   //Divide semi_implicit update by conserved quantity
    179   for (k=0; k<N; k++) {
    180     x = centroid_values[k];
    181     if (x == 0.0) {
    182       semi_implicit_update[k] = 0.0;
    183     } else {
    184       semi_implicit_update[k] /= x;
    185     }
    186   }
    187 
    188 
    189   //Semi implicit updates
    190   for (k=0; k<N; k++) {
    191     denominator = 1.0 - timestep*semi_implicit_update[k];
    192     if (denominator == 0.0) {
    193       return -1;
    194     } else {
    195       //Update conserved_quantities from semi implicit updates
    196       centroid_values[k] /= denominator;
    197     }
    198   }
    199 
    200 /*  for (k=0; k<N; k++) {*/
    201 /*    centroid_values[k] = exp(timestep*semi_implicit_update[k])*centroid_values[k];*/
    202 /*  }*/
    203 
    204 
    205   //Explicit updates
    206   for (k=0; k<N; k++) {
    207     centroid_values[k] += timestep*explicit_update[k];
    208   }
    209 
    210 
    211   //MH080605 set semi_impliit_update[k] to 0.0 here, rather than in update_conserved_quantities.py
    212   for (k=0;k<N;k++){
    213     semi_implicit_update[k]=0.0;
    214   }
    215 
    216   return 0;
     170        //Update centroid values based on values stored in
     171        //explicit_update and semi_implicit_update as well as given timestep
     172
     173
     174        int k;
     175        double denominator, x;
     176
     177
     178        //Divide semi_implicit update by conserved quantity
     179        for (k=0; k<N; k++) {
     180                x = centroid_values[k];
     181                if (x == 0.0) {
     182                        semi_implicit_update[k] = 0.0;
     183                } else {
     184                        semi_implicit_update[k] /= x;
     185                }
     186        }
     187
     188
     189        //Semi implicit updates
     190        for (k=0; k<N; k++) {
     191                denominator = 1.0 - timestep*semi_implicit_update[k];
     192                if (denominator == 0.0) {
     193                        return -1;
     194                } else {
     195                        //Update conserved_quantities from semi implicit updates
     196                        centroid_values[k] /= denominator;
     197                }
     198        }
     199
     200        /*  for (k=0; k<N; k++) {*/
     201        /*    centroid_values[k] = exp(timestep*semi_implicit_update[k])*centroid_values[k];*/
     202        /*  }*/
     203
     204
     205        //Explicit updates
     206        for (k=0; k<N; k++) {
     207                centroid_values[k] += timestep*explicit_update[k];
     208        }
     209
     210
     211        //MH080605 set semi_impliit_update[k] to 0.0 here, rather than in update_conserved_quantities.py
     212        for (k=0;k<N;k++){
     213                semi_implicit_update[k]=0.0;
     214        }
     215
     216        return 0;
    217217}
    218218
     
    222222PyObject *update(PyObject *self, PyObject *args) {
    223223
    224   PyObject *quantity;
    225   PyArrayObject *centroid_values, *explicit_update, *semi_implicit_update;
    226 
    227   double timestep;
    228   int N, err;
    229 
    230 
    231   // Convert Python arguments to C
    232   if (!PyArg_ParseTuple(args, "Od", &quantity, &timestep))
    233     return NULL;
    234 
    235   centroid_values = get_consecutive_array(quantity, "centroid_values");
    236   explicit_update = get_consecutive_array(quantity, "explicit_update");
    237   semi_implicit_update = get_consecutive_array(quantity, "semi_implicit_update");
    238 
    239   N = centroid_values -> dimensions[0];
    240 
    241   err = _update(N, timestep,
    242                 (double*) centroid_values -> data,
    243                 (double*) explicit_update -> data,
    244                 (double*) semi_implicit_update -> data);
    245 
    246 
    247   if (err != 0) {
    248     PyErr_SetString(PyExc_RuntimeError,
    249                     "Zero division in semi implicit update - call Stephen :)");
    250     return NULL;
    251   }
    252 
    253   //Release and return
    254   Py_DECREF(centroid_values);
    255   Py_DECREF(explicit_update);
    256   Py_DECREF(semi_implicit_update);
    257 
    258   return Py_BuildValue("");
     224        PyObject *quantity;
     225        PyArrayObject *centroid_values, *explicit_update, *semi_implicit_update;
     226
     227        double timestep;
     228        int N, err;
     229
     230
     231        // Convert Python arguments to C
     232        if (!PyArg_ParseTuple(args, "Od", &quantity, &timestep))
     233                return NULL;
     234
     235        centroid_values = get_consecutive_array(quantity, "centroid_values");
     236        explicit_update = get_consecutive_array(quantity, "explicit_update");
     237        semi_implicit_update = get_consecutive_array(quantity, "semi_implicit_update");
     238
     239        N = centroid_values -> dimensions[0];
     240
     241        err = _update(N, timestep,
     242                        (double*) centroid_values -> data,
     243                        (double*) explicit_update -> data,
     244                        (double*) semi_implicit_update -> data);
     245
     246
     247        if (err != 0) {
     248                PyErr_SetString(PyExc_RuntimeError,
     249                        "Zero division in semi implicit update - call Stephen :)");
     250                return NULL;
     251        }
     252
     253        //Release and return
     254        Py_DECREF(centroid_values);
     255        Py_DECREF(explicit_update);
     256        Py_DECREF(semi_implicit_update);
     257
     258        return Py_BuildValue("");
    259259}
    260260
     
    262262PyObject *interpolate_from_vertices_to_edges(PyObject *self, PyObject *args) {
    263263
    264   PyObject *quantity;
    265   PyArrayObject *vertex_values, *edge_values;
    266 
    267   int N, err;
    268 
    269   // Convert Python arguments to C
    270   if (!PyArg_ParseTuple(args, "O", &quantity))
    271     return NULL;
    272 
    273   vertex_values = get_consecutive_array(quantity, "vertex_values");
    274   edge_values = get_consecutive_array(quantity, "edge_values");
    275 
    276   N = vertex_values -> dimensions[0];
    277 
    278   err = _interpolate(N,
     264        PyObject *quantity;
     265        PyArrayObject *vertex_values, *edge_values;
     266
     267        int N, err;
     268
     269        // Convert Python arguments to C
     270        if (!PyArg_ParseTuple(args, "O", &quantity))
     271                return NULL;
     272
     273        vertex_values = get_consecutive_array(quantity, "vertex_values");
     274        edge_values = get_consecutive_array(quantity, "edge_values");
     275
     276        N = vertex_values -> dimensions[0];
     277
     278        err = _interpolate(N,
    279279                     (double*) vertex_values -> data,
    280280                     (double*) edge_values -> data);
    281281
    282   if (err != 0) {
    283     PyErr_SetString(PyExc_RuntimeError, "Interpolate could not be computed");
    284     return NULL;
    285   }
    286 
    287   //Release and return
    288   Py_DECREF(vertex_values);
    289   Py_DECREF(edge_values);
    290 
    291   return Py_BuildValue("");
     282        if (err != 0) {
     283                PyErr_SetString(PyExc_RuntimeError, "Interpolate could not be computed");
     284                return NULL;
     285        }
     286
     287        //Release and return
     288        Py_DECREF(vertex_values);
     289        Py_DECREF(edge_values);
     290
     291        return Py_BuildValue("");
    292292}
    293293
     
    295295PyObject *compute_gradients(PyObject *self, PyObject *args) {
    296296
    297   PyObject *quantity, *domain, *R;
    298   PyArrayObject
    299     *centroids,            //Coordinates at centroids
    300     *centroid_values,      //Values at centroids
    301     *number_of_boundaries, //Number of boundaries for each triangle
    302     *surrogate_neighbours, //True neighbours or - if one missing - self
    303     *a, *b;                //Return values
    304 
    305   int dimensions[1], N, err;
    306 
    307   // Convert Python arguments to C
    308   if (!PyArg_ParseTuple(args, "O", &quantity))
    309     return NULL;
    310 
    311   domain = PyObject_GetAttrString(quantity, "domain");
    312   if (!domain)
    313     return NULL;
    314 
    315   //Get pertinent variables
    316 
    317   centroids = get_consecutive_array(domain, "centroid_coordinates");
    318   centroid_values = get_consecutive_array(quantity, "centroid_values");
    319   surrogate_neighbours = get_consecutive_array(domain, "surrogate_neighbours");
    320   number_of_boundaries = get_consecutive_array(domain, "number_of_boundaries");
    321 
    322   N = centroid_values -> dimensions[0];
    323 
    324   //Release
    325   Py_DECREF(domain);
    326 
    327   //Allocate space for return vectors a and b (don't DECREF)
    328   dimensions[0] = N;
    329   a = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE);
    330   b = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE);
    331 
    332 
    333 
    334   err = _compute_gradients(N,
    335                            (double*) centroids -> data,
    336                            (double*) centroid_values -> data,
    337                            (long*) number_of_boundaries -> data,
    338                            (long*) surrogate_neighbours -> data,
    339                            (double*) a -> data,
    340                            (double*) b -> data);
    341 
    342   if (err != 0) {
    343     PyErr_SetString(PyExc_RuntimeError, "Gradient could not be computed");
    344     return NULL;
    345   }
    346 
    347   //Release
    348   Py_DECREF(centroids);
    349   Py_DECREF(centroid_values);
    350   Py_DECREF(number_of_boundaries);
    351   Py_DECREF(surrogate_neighbours);
    352 
    353   //Build result, release and return
    354   R = Py_BuildValue("OO", PyArray_Return(a), PyArray_Return(b));
    355   Py_DECREF(a);
    356   Py_DECREF(b);
    357   return R;
     297        PyObject *quantity, *domain, *R;
     298        PyArrayObject
     299                *centroids,            //Coordinates at centroids
     300                *centroid_values,      //Values at centroids
     301                *number_of_boundaries, //Number of boundaries for each triangle
     302                *surrogate_neighbours, //True neighbours or - if one missing - self
     303                *a, *b;                //Return values
     304
     305        int dimensions[1], N, err;
     306
     307        // Convert Python arguments to C
     308        if (!PyArg_ParseTuple(args, "O", &quantity))
     309                return NULL;
     310
     311        domain = PyObject_GetAttrString(quantity, "domain");
     312        if (!domain)
     313                return NULL;
     314
     315        //Get pertinent variables
     316
     317        centroids = get_consecutive_array(domain, "centroid_coordinates");
     318        centroid_values = get_consecutive_array(quantity, "centroid_values");
     319        surrogate_neighbours = get_consecutive_array(domain, "surrogate_neighbours");
     320        number_of_boundaries = get_consecutive_array(domain, "number_of_boundaries");
     321
     322        N = centroid_values -> dimensions[0];
     323
     324        //Release
     325        Py_DECREF(domain);
     326
     327        //Allocate space for return vectors a and b (don't DECREF)
     328        dimensions[0] = N;
     329        a = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE);
     330        b = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE);
     331
     332
     333
     334        err = _compute_gradients(N,
     335                        (double*) centroids -> data,
     336                        (double*) centroid_values -> data,
     337                        (long*) number_of_boundaries -> data,
     338                        (long*) surrogate_neighbours -> data,
     339                        (double*) a -> data,
     340                        (double*) b -> data);
     341
     342        if (err != 0) {
     343                PyErr_SetString(PyExc_RuntimeError, "Gradient could not be computed");
     344                return NULL;
     345        }
     346
     347        //Release
     348        Py_DECREF(centroids);
     349        Py_DECREF(centroid_values);
     350        Py_DECREF(number_of_boundaries);
     351        Py_DECREF(surrogate_neighbours);
     352
     353        //Build result, release and return
     354        R = Py_BuildValue("OO", PyArray_Return(a), PyArray_Return(b));
     355        Py_DECREF(a);
     356        Py_DECREF(b);
     357        return R;
    358358}
    359359
     
    362362PyObject *extrapolate_second_order(PyObject *self, PyObject *args) {
    363363
    364   PyObject *quantity, *domain;
    365   PyArrayObject
    366     *centroids,            //Coordinates at centroids
    367     *centroid_values,      //Values at centroids
    368     *vertex_coordinates,   //Coordinates at vertices
    369     *vertex_values,        //Values at vertices
    370     *number_of_boundaries, //Number of boundaries for each triangle
    371     *surrogate_neighbours, //True neighbours or - if one missing - self
    372     *a, *b;                //Gradients
    373 
    374   //int N, err;
    375   int dimensions[1], N, err;
    376   //double *a, *b;  //Gradients
    377 
    378   // Convert Python arguments to C
    379   if (!PyArg_ParseTuple(args, "O", &quantity))
    380     return NULL;
    381 
    382   domain = PyObject_GetAttrString(quantity, "domain");
    383   if (!domain)
    384     return NULL;
    385 
    386   //Get pertinent variables
    387   centroids = get_consecutive_array(domain, "centroid_coordinates");
    388   centroid_values = get_consecutive_array(quantity, "centroid_values");
    389   surrogate_neighbours = get_consecutive_array(domain, "surrogate_neighbours");
    390   number_of_boundaries = get_consecutive_array(domain, "number_of_boundaries");
    391   vertex_coordinates = get_consecutive_array(domain, "vertex_coordinates");
    392   vertex_values = get_consecutive_array(quantity, "vertex_values");
    393 
    394   N = centroid_values -> dimensions[0];
    395 
    396   //Release
    397   Py_DECREF(domain);
    398 
    399   //Allocate space for return vectors a and b (don't DECREF)
    400   dimensions[0] = N;
    401   a = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE);
    402   b = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE);
    403 
    404   //FIXME: Odd that I couldn't use normal arrays
    405   //Allocate space for return vectors a and b (don't DECREF)
    406   //a = (double*) malloc(N * sizeof(double));
    407   //if (!a) return NULL;
    408   //b = (double*) malloc(N * sizeof(double));
    409   //if (!b) return NULL;
    410 
    411 
    412   err = _compute_gradients(N,
    413                            (double*) centroids -> data,
    414                            (double*) centroid_values -> data,
    415                            (long*) number_of_boundaries -> data,
    416                            (long*) surrogate_neighbours -> data,
    417                            (double*) a -> data,
    418                            (double*) b -> data);
    419 
    420   if (err != 0) {
    421     PyErr_SetString(PyExc_RuntimeError, "Gradient could not be computed");
    422     return NULL;
    423   }
    424 
    425   err = _extrapolate(N,
    426                      (double*) centroids -> data,
    427                      (double*) centroid_values -> data,
    428                      (double*) vertex_coordinates -> data,
    429                      (double*) vertex_values -> data,
    430                      (double*) a -> data,
    431                      (double*) b -> data);
    432 
    433 
    434   if (err != 0) {
    435     PyErr_SetString(PyExc_RuntimeError,
    436                     "Internal function _extrapolate failed");
    437     return NULL;
    438   }
    439 
    440 
    441 
    442   //Release
    443   Py_DECREF(centroids);
    444   Py_DECREF(centroid_values);
    445   Py_DECREF(number_of_boundaries);
    446   Py_DECREF(surrogate_neighbours);
    447   Py_DECREF(vertex_coordinates);
    448   Py_DECREF(vertex_values);
    449   Py_DECREF(a);
    450   Py_DECREF(b);
    451 
    452   return Py_BuildValue("");
     364        PyObject *quantity, *domain;
     365        PyArrayObject
     366                *centroids,            //Coordinates at centroids
     367                *centroid_values,      //Values at centroids
     368                *vertex_coordinates,   //Coordinates at vertices
     369                *vertex_values,        //Values at vertices
     370                *number_of_boundaries, //Number of boundaries for each triangle
     371                *surrogate_neighbours, //True neighbours or - if one missing - self
     372                *a, *b;                //Gradients
     373
     374        //int N, err;
     375        int dimensions[1], N, err;
     376        //double *a, *b;  //Gradients
     377
     378        // Convert Python arguments to C
     379        if (!PyArg_ParseTuple(args, "O", &quantity))
     380                return NULL;
     381
     382        domain = PyObject_GetAttrString(quantity, "domain");
     383        if (!domain)
     384                return NULL;
     385
     386        //Get pertinent variables
     387        centroids = get_consecutive_array(domain, "centroid_coordinates");
     388        centroid_values = get_consecutive_array(quantity, "centroid_values");
     389        surrogate_neighbours = get_consecutive_array(domain, "surrogate_neighbours");
     390        number_of_boundaries = get_consecutive_array(domain, "number_of_boundaries");
     391        vertex_coordinates = get_consecutive_array(domain, "vertex_coordinates");
     392        vertex_values = get_consecutive_array(quantity, "vertex_values");
     393
     394        N = centroid_values -> dimensions[0];
     395
     396        //Release
     397        Py_DECREF(domain);
     398
     399        //Allocate space for return vectors a and b (don't DECREF)
     400        dimensions[0] = N;
     401        a = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE);
     402        b = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE);
     403
     404        //FIXME: Odd that I couldn't use normal arrays
     405        //Allocate space for return vectors a and b (don't DECREF)
     406        //a = (double*) malloc(N * sizeof(double));
     407        //if (!a) return NULL;
     408        //b = (double*) malloc(N * sizeof(double));
     409        //if (!b) return NULL;
     410
     411
     412        err = _compute_gradients(N,
     413                        (double*) centroids -> data,
     414                        (double*) centroid_values -> data,
     415                        (long*) number_of_boundaries -> data,
     416                        (long*) surrogate_neighbours -> data,
     417                        (double*) a -> data,
     418                        (double*) b -> data);
     419
     420        if (err != 0) {
     421                PyErr_SetString(PyExc_RuntimeError, "Gradient could not be computed");
     422                return NULL;
     423        }
     424
     425        err = _extrapolate(N,
     426                        (double*) centroids -> data,
     427                        (double*) centroid_values -> data,
     428                        (double*) vertex_coordinates -> data,
     429                        (double*) vertex_values -> data,
     430                        (double*) a -> data,
     431                        (double*) b -> data);
     432
     433
     434        if (err != 0) {
     435                PyErr_SetString(PyExc_RuntimeError,
     436                        "Internal function _extrapolate failed");
     437                return NULL;
     438        }
     439
     440
     441
     442        //Release
     443        Py_DECREF(centroids);
     444        Py_DECREF(centroid_values);
     445        Py_DECREF(number_of_boundaries);
     446        Py_DECREF(surrogate_neighbours);
     447        Py_DECREF(vertex_coordinates);
     448        Py_DECREF(vertex_values);
     449        Py_DECREF(a);
     450        Py_DECREF(b);
     451
     452        return Py_BuildValue("");
    453453}
    454454
     
    457457PyObject *limit(PyObject *self, PyObject *args) {
    458458
    459   PyObject *quantity, *domain, *Tmp;
    460   PyArrayObject
    461     *qv, //Conserved quantities at vertices
    462     *qc, //Conserved quantities at centroids
    463     *neighbours;
    464 
    465   int k, i, n, N, k3;
    466   double beta_w; //Safety factor
    467   double *qmin, *qmax, qn;
    468 
    469   // Convert Python arguments to C
    470   if (!PyArg_ParseTuple(args, "O", &quantity))
    471     return NULL;
    472 
    473   domain = PyObject_GetAttrString(quantity, "domain");
    474   if (!domain)
    475     return NULL;
    476 
    477   //neighbours = (PyArrayObject*) PyObject_GetAttrString(domain, "neighbours");
    478   neighbours = get_consecutive_array(domain, "neighbours");
    479 
    480   //Get safety factor beta_w
    481   Tmp = PyObject_GetAttrString(domain, "beta_w");
    482   if (!Tmp)
    483     return NULL;
    484 
    485   beta_w = PyFloat_AsDouble(Tmp);
    486 
    487   Py_DECREF(Tmp);
    488   Py_DECREF(domain);
    489 
    490 
    491   qc = get_consecutive_array(quantity, "centroid_values");
    492   qv = get_consecutive_array(quantity, "vertex_values");
    493 
    494 
    495   N = qc -> dimensions[0];
    496 
    497   //Find min and max of this and neighbour's centroid values
    498   qmin = malloc(N * sizeof(double));
    499   qmax = malloc(N * sizeof(double));
    500   for (k=0; k<N; k++) {
    501     k3=k*3;
    502 
    503     qmin[k] = ((double*) qc -> data)[k];
    504     qmax[k] = qmin[k];
    505 
    506     for (i=0; i<3; i++) {
    507       n = ((long*) neighbours -> data)[k3+i];
    508       if (n >= 0) {
    509         qn = ((double*) qc -> data)[n]; //Neighbour's centroid value
    510 
    511         qmin[k] = min(qmin[k], qn);
    512         qmax[k] = max(qmax[k], qn);
    513       }
    514     }
    515   }
    516 
    517   // Call underlying routine
    518   _limit(N, beta_w, (double*) qc -> data, (double*) qv -> data, qmin, qmax);
    519 
    520   free(qmin);
    521   free(qmax);
    522   return Py_BuildValue("");
     459        PyObject *quantity, *domain, *Tmp;
     460        PyArrayObject
     461                *qv, //Conserved quantities at vertices
     462                *qc, //Conserved quantities at centroids
     463                *neighbours;
     464
     465        int k, i, n, N, k3;
     466        double beta_w; //Safety factor
     467        double *qmin, *qmax, qn;
     468
     469        // Convert Python arguments to C
     470        if (!PyArg_ParseTuple(args, "O", &quantity))
     471                return NULL;
     472
     473        domain = PyObject_GetAttrString(quantity, "domain");
     474        if (!domain)
     475                return NULL;
     476
     477        //neighbours = (PyArrayObject*) PyObject_GetAttrString(domain, "neighbours");
     478        neighbours = get_consecutive_array(domain, "neighbours");
     479
     480        //Get safety factor beta_w
     481        Tmp = PyObject_GetAttrString(domain, "beta_w");
     482        if (!Tmp)
     483                return NULL;
     484
     485        beta_w = PyFloat_AsDouble(Tmp);
     486
     487        Py_DECREF(Tmp);
     488        Py_DECREF(domain);
     489
     490
     491        qc = get_consecutive_array(quantity, "centroid_values");
     492        qv = get_consecutive_array(quantity, "vertex_values");
     493
     494
     495        N = qc -> dimensions[0];
     496
     497        //Find min and max of this and neighbour's centroid values
     498        qmin = malloc(N * sizeof(double));
     499        qmax = malloc(N * sizeof(double));
     500        for (k=0; k<N; k++) {
     501                k3=k*3;
     502
     503                qmin[k] = ((double*) qc -> data)[k];
     504                qmax[k] = qmin[k];
     505
     506                for (i=0; i<3; i++) {
     507                        n = ((long*) neighbours -> data)[k3+i];
     508                        if (n >= 0) {
     509                                qn = ((double*) qc -> data)[n]; //Neighbour's centroid value
     510
     511                                qmin[k] = min(qmin[k], qn);
     512                                qmax[k] = max(qmax[k], qn);
     513                        }
     514                        //qmin[k] = max(qmin[k],0.5*((double*) qc -> data)[k]);
     515                        //qmax[k] = min(qmax[k],2.0*((double*) qc -> data)[k]);
     516                }
     517        }
     518
     519        // Call underlying routine
     520        _limit(N, beta_w, (double*) qc -> data, (double*) qv -> data, qmin, qmax);
     521
     522        free(qmin);
     523        free(qmax);
     524        return Py_BuildValue("");
    523525}
    524526
     
    527529// Method table for python module
    528530static struct PyMethodDef MethodTable[] = {
    529   {"limit", limit, METH_VARARGS, "Print out"},
    530   {"update", update, METH_VARARGS, "Print out"},
    531   {"compute_gradients", compute_gradients, METH_VARARGS, "Print out"},
    532   {"extrapolate_second_order", extrapolate_second_order,
    533    METH_VARARGS, "Print out"},
    534   {"interpolate_from_vertices_to_edges",
    535    interpolate_from_vertices_to_edges,
    536    METH_VARARGS, "Print out"},
    537   {NULL, NULL, 0, NULL}   // sentinel
     531        {"limit", limit, METH_VARARGS, "Print out"},
     532        {"update", update, METH_VARARGS, "Print out"},
     533        {"compute_gradients", compute_gradients, METH_VARARGS, "Print out"},
     534        {"extrapolate_second_order", extrapolate_second_order,
     535                METH_VARARGS, "Print out"},
     536        {"interpolate_from_vertices_to_edges",
     537                interpolate_from_vertices_to_edges,
     538                METH_VARARGS, "Print out"},
     539        {NULL, NULL, 0, NULL}   // sentinel
    538540};
    539541
  • anuga_core/source/anuga/abstract_2d_finite_volumes/test_quantity.py

    r3514 r3678  
    941941        quantity.limit()
    942942
    943 
     943        # limited value for beta_w = 0.9
    944944        assert allclose(quantity.vertex_values[1,:], [2.2, 4.9, 4.9])
     945        # limited values for beta_w = 0.5
     946        #assert allclose(quantity.vertex_values[1,:], [3.0, 4.5, 4.5])
    945947
    946948
  • anuga_core/source/anuga/compile_all.py

    r3560 r3678  
    1313
    1414os.chdir('..')
     15os.chdir('shallow_water')
     16execfile('..' + os.sep + 'utilities' + os.sep + 'compile.py')
     17
     18os.chdir('..')
    1519os.chdir('mesh_engine')
    1620execfile('compile.py')
  • anuga_core/source/anuga/config.py

    r3642 r3678  
    3939max_smallsteps = 50  #Max number of degenerate steps allowed b4 trying first order
    4040
    41 manning = 0.3  #Manning's friction coefficient
     41manning = 0.03  #Manning's friction coefficient
    4242#g = 9.80665       #Gravity
    4343g = 9.8
     
    6767#
    6868#
    69 #There are separate betas for the w-limiter and the h-limiter
    70 #
    71 #
    72 #
     69#There are separate betas for the w, uh, vh and h limiters
    7370#
    7471#Good values are:
    75 #beta_w = 0.9
    76 #beta_h = 0.2
     72beta_w      = 0.9
     73beta_w_dry  = 0.9
     74beta_uh     = 0.9
     75beta_uh_dry = 0.9
     76beta_vh     = 0.9
     77beta_vh_dry = 0.9
     78beta_h      = 0.2
     79
     80# I think these are better SR but they conflict with the unit tests!
     81# beta_w      = 1.0
     82# beta_w_dry  = 0.2
     83# beta_uh     = 1.0
     84# beta_uh_dry = 0.2
     85# beta_vh     = 1.0
     86# beta_vh_dry = 0.2
     87# beta_h      = 0.2
    7788
    7889
    79 
    80 beta_w = 0.9
    81 beta_h = 0.2
    8290CFL = 1.0  #FIXME (ole): Is this in use yet??
    8391           #(Steve) yes, change domain.CFL to
  • anuga_core/source/anuga/examples/netherlands.py

    r3659 r3678  
    1515     Transmissive_boundary, Constant_height, Constant_stage
    1616
    17 from anuga.mesh_factory import rectangular_cross
     17from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
    1818from Numeric import array
    19 #from vtk_realtime_visualiser import Visualiser
     19from anuga.visualiser.vtk_realtime_visualiser import Visualiser
    2020
    2121class Weir:
     
    8383N = 130 #size = 33800
    8484N = 600 #Size = 720000
    85 N = 50
     85N = 100
    8686
    8787
     
    9696
    9797domain.check_integrity()
     98
     99#Setup order and all the beta's for the limiters (these should become defaults
    98100domain.default_order = 2
    99 #domain.beta_h=0
     101domain.beta_w      = 1.0
     102domain.beta_w_dry  = 0.2
     103domain.beta_uh     = 1.0
     104domain.beta_uh_dry = 0.2
     105domain.beta_vh     = 1.0
     106domain.beta_vh_dry = 0.2
    100107
    101108#Output params
  • anuga_core/source/anuga/examples/sww_file_visualiser_example.py

    r3670 r3678  
    1010# The argument to OfflineVisualiser is the path to a sww file
    1111o = OfflineVisualiser("../../swollen_viewer/tests/cylinders.sww")
     12#o = OfflineVisualiser("../../../../anuga_validation/analytical solutions/circular_second_order.sww")
    1213
    1314# Specify the height-based-quantities to render.
  • anuga_core/source/anuga/shallow_water/__init__.py

    r3659 r3678  
    88# Make selected classes available directly
    99from shallow_water_domain import Domain,\
    10      Transmissive_boundary, Reflective_boundary,\
    11      Dirichlet_boundary, Time_boundary, File_boundary,\
    12      Transmissive_Momentum_Set_Stage_boundary,\
     10    Transmissive_boundary, Reflective_boundary,\
     11    Dirichlet_boundary, Time_boundary, File_boundary,\
     12    Transmissive_Momentum_Set_Stage_boundary,\
     13    Dirichlet_Discharge_boundary,\
    1314    Constant_stage, Constant_height
    1415
    1516
    16 
  • anuga_core/source/anuga/shallow_water/data_manager.py

    r3664 r3678  
    724724
    725725
    726 #### NBED national exposure database
     726#### NEED national exposure database
    727727   
    728728LAT_TITLE = 'LATITUDE'
  • anuga_core/source/anuga/shallow_water/shallow_water_domain.py

    r3642 r3678  
    127127
    128128
    129         from anuga.config import minimum_allowed_height, maximum_allowed_speed, g
     129        from anuga.config import *
    130130        self.minimum_allowed_height = minimum_allowed_height
    131131        self.maximum_allowed_speed = maximum_allowed_speed
    132132        self.g = g
     133        self.beta_w      = beta_w
     134        self.beta_w_dry  = beta_w_dry
     135        self.beta_uh     = beta_uh
     136        self.beta_uh_dry = beta_uh_dry
     137        self.beta_vh     = beta_vh
     138        self.beta_vh_dry = beta_vh_dry
     139        self.beta_h      = beta_h
    133140
    134141
     
    150157        self.minimum_storable_height = minimum_storable_height
    151158        self.quantities_to_be_stored = ['stage','xmomentum','ymomentum']
     159               
    152160
    153161
     
    301309
    302310        msg = 'Parameter beta_h must be in the interval [0, 1['
    303         assert 0 <= self.beta_h < 1.0, msg
     311        assert 0 <= self.beta_h <= 1.0, msg
    304312        msg = 'Parameter beta_w must be in the interval [0, 1['
    305         assert 0 <= self.beta_w < 1.0, msg
     313        assert 0 <= self.beta_w <= 1.0, msg
    306314
    307315
     
    612620    Xmom = domain.quantities['xmomentum']
    613621    Ymom = domain.quantities['ymomentum']
     622    Elevation = domain.quantities['elevation']
    614623    from shallow_water_ext import extrapolate_second_order_sw
    615     extrapolate_second_order_sw(domain,domain.surrogate_neighbours,
     624    extrapolate_second_order_sw(domain,
     625                                domain.surrogate_neighbours,
    616626                                domain.number_of_boundaries,
    617627                                domain.centroid_coordinates,
     
    619629                                Xmom.centroid_values,
    620630                                Ymom.centroid_values,
     631                                Elevation.centroid_values,
    621632                                domain.vertex_coordinates,
    622633                                Stage.vertex_values,
    623634                                Xmom.vertex_values,
    624                                 Ymom.vertex_values)
     635                                Ymom.vertex_values,
     636                                Elevation.vertex_values)
    625637
    626638def compute_fluxes_c(domain):
     
    699711        #all of the conserved quantitie
    700712
    701         if (domain.order == 1):
     713        if (domain._order_ == 1):
    702714            for name in domain.conserved_quantities:
    703715                Q = domain.quantities[name]
    704716                Q.extrapolate_first_order()
    705         elif domain.order == 2:
     717        elif domain._order_ == 2:
    706718            domain.extrapolate_second_order_sw()
    707719        else:
     
    711723        for name in domain.conserved_quantities:
    712724            Q = domain.quantities[name]
    713             if domain.order == 1:
     725            if domain._order_ == 1:
    714726                Q.extrapolate_first_order()
    715             elif domain.order == 2:
     727            elif domain._order_ == 2:
    716728                Q.extrapolate_second_order()
    717729                Q.limit()
  • anuga_core/source/anuga/shallow_water/shallow_water_ext.c

    r3563 r3678  
    622622    *xmom_centroid_values,
    623623    *ymom_centroid_values,
     624        *elevation_centroid_values,
    624625    *vertex_coordinates,
    625626    *stage_vertex_values,
    626627    *xmom_vertex_values,
    627     *ymom_vertex_values;
    628   PyObject *domain, *Tmp;
     628    *ymom_vertex_values,
     629        *elevation_vertex_values;
     630  PyObject *domain, *Tmp_w, *Tmp_w_dry, *Tmp_uh, *Tmp_uh_dry, *Tmp_vh, *Tmp_vh_dry;
    629631  //Local variables
    630632  double a, b;//gradient vector, not stored but used to calculate vertex values from centroids
     
    632634  double x,y,x0,y0,x1,y1,x2,y2,xv0,yv0,xv1,yv1,xv2,yv2;//vertices of the auxiliary triangle
    633635  double dx1,dx2,dy1,dy2,dxv0,dxv1,dxv2,dyv0,dyv1,dyv2,dq0,dq1,dq2,area2;
    634   double dqv[3], qmin, qmax, beta_w;//provisional jumps from centroids to v'tices and safety factor re limiting
     636  double dqv[3], qmin, qmax, hmin;
     637  double hc, h0, h1, h2;
     638  double beta_w, beta_w_dry, beta_uh, beta_uh_dry, beta_vh, beta_vh_dry, beta_tmp;
     639  //provisional jumps from centroids to v'tices and safety factor re limiting
    635640  //by which these jumps are limited
    636641  // Convert Python arguments to C
    637   if (!PyArg_ParseTuple(args, "OOOOOOOOOOO",
     642  if (!PyArg_ParseTuple(args, "OOOOOOOOOOOOO",
    638643                        &domain,
    639644                        &surrogate_neighbours,
     
    643648                        &xmom_centroid_values,
    644649                        &ymom_centroid_values,
     650                        &elevation_centroid_values,
    645651                        &vertex_coordinates,
    646652                        &stage_vertex_values,
    647653                        &xmom_vertex_values,
    648                         &ymom_vertex_values)) {
     654                        &ymom_vertex_values,
     655                        &elevation_vertex_values)) {
    649656    PyErr_SetString(PyExc_RuntimeError, "Input arguments failed");
    650657    return NULL;
     
    652659
    653660  //get the safety factor beta_w, set in the config.py file. This is used in the limiting process
    654   Tmp = PyObject_GetAttrString(domain, "beta_w");
    655   if (!Tmp)
    656     return NULL;
    657   beta_w = PyFloat_AsDouble(Tmp);
    658   Py_DECREF(Tmp);
     661  Tmp_w = PyObject_GetAttrString(domain, "beta_w");
     662  if (!Tmp_w)
     663    return NULL;
     664  beta_w = PyFloat_AsDouble(Tmp_w);
     665  Py_DECREF(Tmp_w);
     666 
     667  Tmp_w_dry = PyObject_GetAttrString(domain, "beta_w_dry");
     668  if (!Tmp_w_dry)
     669    return NULL;
     670  beta_w_dry = PyFloat_AsDouble(Tmp_w_dry);
     671  Py_DECREF(Tmp_w_dry);
     672 
     673  Tmp_uh = PyObject_GetAttrString(domain, "beta_uh");
     674  if (!Tmp_uh)
     675    return NULL;
     676  beta_uh = PyFloat_AsDouble(Tmp_uh);
     677  Py_DECREF(Tmp_uh);
     678 
     679  Tmp_uh_dry = PyObject_GetAttrString(domain, "beta_uh_dry");
     680  if (!Tmp_uh_dry)
     681    return NULL;
     682  beta_uh_dry = PyFloat_AsDouble(Tmp_uh_dry);
     683  Py_DECREF(Tmp_uh_dry);
     684
     685  Tmp_vh = PyObject_GetAttrString(domain, "beta_vh");
     686  if (!Tmp_vh)
     687    return NULL;
     688  beta_vh = PyFloat_AsDouble(Tmp_vh);
     689  Py_DECREF(Tmp_vh);
     690 
     691  Tmp_vh_dry = PyObject_GetAttrString(domain, "beta_vh_dry");
     692  if (!Tmp_vh_dry)
     693    return NULL;
     694  beta_vh_dry = PyFloat_AsDouble(Tmp_vh_dry);
     695  Py_DECREF(Tmp_vh_dry);
     696 
    659697  number_of_elements = stage_centroid_values -> dimensions[0];
    660698  for (k=0; k<number_of_elements; k++) {
     
    710748      //If the mesh is 'weird' near the boundary, the trianlge might be flat or clockwise:
    711749      if (area2<=0)
    712         return NULL;
    713 
     750                return NULL;
     751
     752          //### Calculate heights of neighbouring cells
     753          hc = ((double *)stage_centroid_values->data)[k]  - ((double *)elevation_centroid_values->data)[k];
     754          h0 = ((double *)stage_centroid_values->data)[k0] - ((double *)elevation_centroid_values->data)[k0];
     755          h1 = ((double *)stage_centroid_values->data)[k1] - ((double *)elevation_centroid_values->data)[k1];
     756          h2 = ((double *)stage_centroid_values->data)[k2] - ((double *)elevation_centroid_values->data)[k2];
     757          hmin = min(hc,min(h0,min(h1,h2)));
    714758      //### stage ###
    715759      //calculate the difference between vertex 0 of the auxiliary triangle and the FV triangle centroid
     
    730774      //and compute jumps from the centroid to the min and max
    731775      find_qmin_and_qmax(dq0,dq1,dq2,&qmin,&qmax);
    732       limit_gradient(dqv,qmin,qmax,beta_w);//the gradient will be limited
     776          // Playing with dry wet interface
     777          hmin = qmin;
     778          beta_tmp = beta_w;
     779          if (hmin<0.001)
     780                beta_tmp = beta_w_dry;
     781      limit_gradient(dqv,qmin,qmax,beta_tmp);//the gradient will be limited
    733782      for (i=0;i<3;i++)
    734783        ((double *)stage_vertex_values->data)[k3+i]=((double *)stage_centroid_values->data)[k]+dqv[i];
     
    752801      //and compute jumps from the centroid to the min and max
    753802      find_qmin_and_qmax(dq0,dq1,dq2,&qmin,&qmax);
    754       limit_gradient(dqv,qmin,qmax,beta_w);//the gradient will be limited
     803          beta_tmp = beta_uh;
     804          if (hmin<0.001)
     805                beta_tmp = beta_uh_dry;
     806      limit_gradient(dqv,qmin,qmax,beta_tmp);//the gradient will be limited
    755807      for (i=0;i<3;i++)
    756808        ((double *)xmom_vertex_values->data)[k3+i]=((double *)xmom_centroid_values->data)[k]+dqv[i];
     
    774826      //and compute jumps from the centroid to the min and max
    775827      find_qmin_and_qmax(dq0,dq1,dq2,&qmin,&qmax);
    776       limit_gradient(dqv,qmin,qmax,beta_w);//the gradient will be limited
     828          beta_tmp = beta_vh;
     829          if (hmin<0.001)
     830                beta_tmp = beta_vh_dry;
     831      limit_gradient(dqv,qmin,qmax,beta_tmp);//the gradient will be limited
    777832      for (i=0;i<3;i++)
    778833        ((double *)ymom_vertex_values->data)[k3+i]=((double *)ymom_centroid_values->data)[k]+dqv[i];
  • anuga_core/source/anuga/shallow_water/test_shallow_water_domain.py

    r3632 r3678  
    11881188
    11891189
    1190         domain.order = 1
     1190        domain._order_ = 1
    11911191        domain.distribute_to_vertices_and_edges()
    11921192
     
    12341234        assert not alltrue(alltrue(greater_equal(L,E-epsilon)))
    12351235
    1236         domain.order = 1
     1236        domain._order_ = 1
    12371237        domain.distribute_to_vertices_and_edges()
    12381238
     
    12691269
    12701270        #First order
    1271         domain.order = 1
     1271        domain._order_ = 1
    12721272        domain.distribute_to_vertices_and_edges()
    12731273        assert allclose(L[1], val1)
    12741274
    12751275        #Second order
    1276         domain.order = 2
     1276        domain._order_ = 2
    12771277        domain.distribute_to_vertices_and_edges()
    12781278        assert allclose(L[1], [2.2, 4.9, 4.9])
     
    13071307        assert allclose(b[1], 0.0)
    13081308
    1309         domain.order = 1
     1309        domain._order_ = 1
    13101310        domain.distribute_to_vertices_and_edges()
    13111311        assert allclose(L[1], 1.77777778)
    13121312
    1313         domain.order = 2
     1313        domain._order_ = 2
    13141314        domain.distribute_to_vertices_and_edges()
    13151315        assert allclose(L[1], [0.57777777, 2.37777778, 2.37777778])
     
    13451345        assert allclose(b[1], 3.33333333)
    13461346
    1347         domain.order = 1
     1347        domain._order_ = 1
    13481348        domain.distribute_to_vertices_and_edges()
    13491349        assert allclose(L[1], 4.9382716)
    13501350
    1351         domain.order = 2
     1351        domain._order_ = 2
    13521352        domain.distribute_to_vertices_and_edges()
    13531353        assert allclose(L[1], [1.07160494, 6.46058131, 7.28262855])
     
    13911391
    13921392        #print E
    1393         domain.order = 1
     1393        domain._order_ = 1
    13941394        domain.distribute_to_vertices_and_edges()
    13951395        ##assert allclose(L[1], [0.19999999, 20.05, 20.05])
    13961396        assert allclose(L[1], [0.1, 20.1, 20.1])
    13971397
    1398         domain.order = 2
     1398        domain._order_ = 2
    13991399        domain.distribute_to_vertices_and_edges()
    14001400        assert allclose(L[1], [0.1, 20.1, 20.1])
     
    14361436
    14371437        #print E
    1438         domain.order = 1
     1438        domain._order_ = 1
    14391439        domain.distribute_to_vertices_and_edges()
    14401440        ##assert allclose(L[1], [4.19999999, 16.07142857, 20.02857143])
    14411441        assert allclose(L[1], [4.1, 16.1, 20.1])
    14421442
    1443         domain.order = 2
     1443        domain._order_ = 2
    14441444        domain.distribute_to_vertices_and_edges()
    14451445        assert allclose(L[1], [4.1, 16.1, 20.1])
     
    14881488
    14891489        #print E
    1490         domain.order = 2
     1490        domain._order_ = 2
    14911491        domain.beta_h = 0.0 #Use first order in h-limiter
    14921492        domain.distribute_to_vertices_and_edges()
     
    22322232        domain.smooth = False
    22332233        domain.visualise = False
    2234         domain.default_order=domain.order=2
     2234        domain.default_order=domain._order_=2
    22352235
    22362236        # Boundary conditions
  • anuga_core/source/anuga/test_all.py

    r3573 r3678  
    7676        for file in exclude_files:
    7777            print 'WARNING: File '+ file + ' to be excluded from testing'
    78             try:   
    79                 files.remove(file)
    80             except ValueError, e:
    81                 msg = 'File "%s" was not found in test suite.\n' %file
    82                 msg += 'Original error is "%s"\n' %e
    83                 msg += 'Perhaps it should be removed from exclude list?'
    84                 raise Exception, msg
     78        try:   
     79            files.remove(file)
     80        except ValueError, e:
     81            msg = 'File "%s" was not found in test suite.\n' %file
     82            msg += 'Original error is "%s"\n' %e
     83            msg += 'Perhaps it should be removed from exclude list?'
     84            raise Exception, msg
    8585
    8686    filenameToModuleName = lambda f: os.path.splitext(f)[0]
     
    101101    runner = unittest.TextTestRunner() #verbosity=2
    102102    runner.run(suite)
     103   
  • anuga_core/source/anuga_parallel/run_parallel_sw_rectangle.py

    r3579 r3678  
    3535processor_name = pypar.Get_processor_name()
    3636
    37 M = 22
     37M = 50
    3838N = M*numprocs
    3939
     
    6969
    7070
    71 domain.default_order = 1
     71
     72
    7273
    7374#Boundaries
     
    8586    """
    8687
    87     def __init__(self, x0=0.25, x1=0.75, y0=0.0, y1=1.0, h=5.0):
     88    def __init__(self, x0=0.25, x1=0.75, y0=0.0, y1=1.0, h=5.0, h0=0.0):
    8889        self.x0 = x0
    8990        self.x1 = x1
     
    9192        self.y1 = y1
    9293        self.h  = h
     94        self.h0 = h0
    9395
    9496    def __call__(self, x, y):
    95         return self.h*((x>self.x0)&(x<self.x1)&(y>self.y0)&(y<self.y1))
     97        return self.h0 + self.h*((x>self.x0)&(x<self.x1)&(y>self.y0)&(y<self.y1))
    9698
    97 domain.set_quantity('stage', Set_Stage(0.2,0.4,0.25, 0.75, 1.0))
     99domain.set_quantity('stage', Set_Stage(0.2, 0.4, 0.25, 0.75, 1.0, 0.00))
    98100
    99101if myid == 0:
     
    106108rect = [0.0, 0.0, 1.0, 1.0]
    107109domain.initialise_visualiser()
     110
     111domain.default_order = 2
     112domain.beta_w      = 1.0
     113domain.beta_w_dry  = 0.2
     114domain.beta_uh     = 1.0
     115domain.beta_uh_dry = 0.2
     116domain.beta_vh     = 1.0
     117domain.beta_vh_dry = 0.2
    108118
    109119yieldstep = 0.005
  • anuga_validation/analytical solutions/Analytical_solution_circular_hydraulic_jump.py

    r3514 r3678  
    88
    99#-------------------------------
    10 # Set up path and module imports
    11 import sys
    12 from os import sep
    13 sys.path.append('..'+sep+'pyvolution')
     10# Setup modules
    1411
    15 from shallow_water import Domain, Dirichlet_Discharge_boundary
    16 from shallow_water import Transmissive_Momentum_Set_Stage_boundary, Dirichlet_boundary
     12from anuga.shallow_water import Domain, Dirichlet_Discharge_boundary
     13from anuga.shallow_water import Transmissive_Momentum_Set_Stage_boundary, Dirichlet_boundary
    1714from math import pi, sqrt
    18 from mesh_factory import strang_mesh
     15from anuga.abstract_2d_finite_volumes.mesh_factory import strang_mesh
    1916
    2017
     
    6663#------------------------------------------
    6764# Reduction operation for get_vertex_values
    68 from anuga.pyvolution.util import mean
    69 domain.reduction = mean
     65from anuga.utilities.numerical_tools import mean
     66#domain.reduction = mean
    7067#domain.reduction = min  #Looks better near steep slopes
    7168
     
    116113#------------------
    117114# Order of accuracy
    118 domain.default_order = 1
    119 domain.CFL = 0.75
    120 #domain.beta_w = 0.5
    121 #domain.beta_h = 0.2
    122 domain.smooth = True
     115domain.default_order = 2
     116domain.beta_w      = 1.0
     117domain.beta_w_dry  = 0.2
     118domain.beta_uh     = 1.0
     119domain.beta_uh_dry = 0.2
     120domain.beta_vh     = 1.0
     121domain.beta_vh_dry = 0.2
     122domain.CFL = 0.5
     123
     124#domain.smooth = True
    123125
    124126
     
    135137#
    136138#
     139
     140domain.initialise_visualiser()
     141#from anuga.visualiser.vtk_realtime_visualiser import Visualiser
     142
    137143#from realtime_visualisation_new import Visualiser
     144#vis = Visualiser(domain,title="stage")
     145#vis.setup['elevation'] = True
     146#vis.updating['stage'] = True
     147#vis.qcolor['stage'] = (0.0,0.0,0.8)
     148#vis.coloring['stage']= True
    138149##vxmom = Visualiser(domain,title='xmomentum',scale_z=10.0)
    139150##vymom = Visualiser(domain,title='ymomentum',scale_z=10.0)
     
    148159
    149160t0 = time.time()
    150 for t in domain.evolve(yieldstep = .1, finaltime = 3.0):
     161for t in domain.evolve(yieldstep = .02, finaltime = 3.0):
    151162    domain.write_time()
    152 
     163    #vis.update()
    153164    exp = '(xmomentum**2 + ymomentum**2)**0.5'
    154165    radial_momentum = domain.create_quantity_from_expression(exp)
     
    167178
    168179#    f.write('time = %25.15e wall clock time %g \n' % (domain.time, time.time()))
    169     f.write('%10.3f %25.15e %25.15e %25.15e %25.15e \n' % (domain.time, inner_stage, inner_radial_mom, outer_stage, outer_radial_mom))
     180#    f.write('%10.3f %25.15e %25.15e %25.15e %25.15e \n' % (domain.time, inner_stage, inner_radial_mom, outer_stage, outer_radial_mom))
    170181
    171182    f.write('time = %25.15e wall clock time %g \n' % (domain.time, time.time()))
Note: See TracChangeset for help on using the changeset viewer.