Changeset 5162


Ignore:
Timestamp:
Mar 11, 2008, 8:43:22 PM (16 years ago)
Author:
steve
Message:

Updated some methods for quantity. Looks like we can use old
limiting system with larger values of beta.

Files:
2 added
14 edited

Legend:

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

    r5080 r5162  
    1010from Numeric import allclose, argmax, zeros, Float
    1111from anuga.config import epsilon
     12from anuga.config import beta_euler, beta_rk2
    1213
    1314from anuga.abstract_2d_finite_volumes.neighbour_mesh import Mesh
     
    301302        self.default_order = n
    302303        self._order_ = self.default_order
     304
     305        if self.default_order == 1:
     306            self.set_timestepping_method('euler')
     307            #self.set_all_limiters(beta_euler)
     308
     309        if self.default_order == 2:
     310            self.set_timestepping_method('rk2')
     311            #self.set_all_limiters(beta_rk2)
    303312       
    304313
     
    15501559            elif self._order_ == 2:
    15511560                Q.extrapolate_second_order()
    1552                 Q.limit()
     1561                #Q.limit()
    15531562            else:
    15541563                raise 'Unknown order'
    1555             Q.interpolate_from_vertices_to_edges()
     1564            #Q.interpolate_from_vertices_to_edges()
    15561565
    15571566
  • anuga_core/source/anuga/abstract_2d_finite_volumes/quantity.py

    r4990 r5162  
    6262        #Allocate space for Gradient
    6363        self.x_gradient = zeros(N, Float)
    64         self.y_gradient = zeros(N, Float)       
     64        self.y_gradient = zeros(N, Float)
     65
     66        #Allocate space for Limiter Phi
     67        self.phi = zeros(N, Float)       
    6568
    6669        #Intialise centroid and edge_values
     
    7477        #flux calculations and forcing functions
    7578
    76         N = len(domain) # number_of_triangles
     79        #Allocate space for update fields
    7780        self.explicit_update = zeros(N, Float )
    7881        self.semi_implicit_update = zeros(N, Float )
    79         self.centroid_backup_values = zeros(N, Float)       
     82        self.centroid_backup_values = zeros(N, Float)
     83
     84        self.beta = 1.0
    8085
    8186
     
    14031408        compute_gradients(self)
    14041409        extrapolate_from_gradient(self)
    1405 
     1410       
    14061411    def extrapolate_second_order_and_limit(self):
    14071412        #Call correct module function
    14081413        #(either from this module or C-extension)
    1409         extrapolate_second_order_and_limit(self)       
     1414        extrapolate_second_order_and_limit(self)
     1415
     1416    def bound_vertices_below_by_constant(self, bound):
     1417        #Call correct module function
     1418        #(either from this module or C-extension)
     1419        bound_vertices_below_by_constant(self, bound)
     1420
     1421    def bound_vertices_below_by_quantity(self, quantity):
     1422        #Call correct module function
     1423        #(either from this module or C-extension)
     1424
     1425        #check consistency
     1426        assert self.domain == quantity.domain
     1427        bound_vertices_below_by_quantity(self, quantity)                       
    14101428
    14111429    def backup_centroid_values(self):
     
    14501468         limit_gradient_by_neighbour,\
    14511469         extrapolate_from_gradient,\
     1470         extrapolate_second_order_and_limit,\
     1471         bound_vertices_below_by_constant,\
     1472         bound_vertices_below_by_quantity,\
    14521473         interpolate_from_vertices_to_edges,\
    14531474         interpolate_from_edges_to_vertices,\
  • anuga_core/source/anuga/abstract_2d_finite_volumes/quantity_ext.c

    r4990 r5162  
    145145                                         double* vertex_values,
    146146                                         double* edge_values,
     147                                         double* phi,
    147148                                         double* x_gradient,
    148149                                         double* y_gradient) {
     
    151152        double x, y, x0, y0, x1, y1, x2, y2;
    152153        long n;
    153         double qmin, qmax, qn, qc;
    154         double dq, dqa[3], phi, r;
    155 
     154        double qmin, qmax, qc;
     155        double qn[3];
     156        double dq, dqa[3], r;
    156157
    157158        for (k=0; k<N; k++){
     
    161162
    162163                // Centroid coordinates
    163                 x = centroids[k2]; y = centroids[k2+1];
     164                x = centroids[k2+0];
     165                y = centroids[k2+1];
    164166
    165167                // vertex coordinates
    166168                // x0, y0, x1, y1, x2, y2 = X[k,:]
    167                 x0 = vertex_coordinates[k6 + 0];
     169                x0 = vertex_coordinates[k6 + 0];
    168170                y0 = vertex_coordinates[k6 + 1];
    169171                x1 = vertex_coordinates[k6 + 2];
     
    172174                y2 = vertex_coordinates[k6 + 5];
    173175
    174                 //Centroid Value
    175                 qc = centroid_values[k];
    176 
    177176                // Extrapolate to Vertices
    178177                vertex_values[k3+0] = centroid_values[k] + x_gradient[k]*(x0-x) + y_gradient[k]*(y0-y);
     
    184183                edge_values[k3+1] = 0.5*(vertex_values[k3 + 2]+vertex_values[k3 + 0]);
    185184                edge_values[k3+2] = 0.5*(vertex_values[k3 + 0]+vertex_values[k3 + 1]);
    186 
    187 
    188                 phi = 1.0;
     185        }
     186
     187
     188
     189        for (k=0; k<N; k++){
     190                k6 = 6*k;
     191                k3 = 3*k;
     192                k2 = 2*k;
     193
     194
     195                qc = centroid_values[k];
     196               
     197                qmin = qc;
     198                qmax = qc;
     199               
     200                for (i=0; i<3; i++) {
     201                    n = neighbours[k3+i];
     202                    if (n < 0) {
     203                        qn[i] = qc;
     204                    } else {
     205                        qn[i] = centroid_values[n];
     206                    }
     207
     208                    qmin = min(qmin, qn[i]);
     209                    qmax = max(qmax, qn[i]);
     210                }
     211
     212                //qtmin = min(min(min(qn[0],qn[1]),qn[2]),qc);
     213                //qtmax = max(max(max(qn[0],qn[1]),qn[2]),qc);
     214
     215/*              for (i=0; i<3; i++) { */
     216/*                  n = neighbours[k3+i]; */
     217/*                  if (n < 0) { */
     218/*                      qn[i] = qc; */
     219/*                      qmin[i] = qtmin; */
     220/*                      qmax[i] = qtmax; */
     221/*                  }  */
     222/*              } */
     223                       
     224                phi[k] = 1.0;
    189225
    190226                for (i=0; i<3; i++) {
     
    192228                    dqa[i] = dq;                      //Save dq for use in updating vertex values
    193229
    194                     n = neighbours[k3+i];
    195                     if (n >= 0) {
    196                         qn = centroid_values[n];     //Neighbour's centroid value
    197 
    198                         qmin = min(qc, qn);
    199                         qmax = max(qc, qn);
    200 
    201 
    202                         r = 1.0;
     230                    r = 1.0;
    203231     
    204                         if (dq > 0.0) r = (qmax - qc)/dq;
    205                         if (dq < 0.0) r = (qmin - qc)/dq;     
     232                    if (dq > 0.0) r = (qmax - qc)/dq;
     233                    if (dq < 0.0) r = (qmin - qc)/dq;
    206234                   
    207                         phi = min( min(r*beta, 1.0), phi);   
    208                     }
     235                    phi[k] = min( min(r*beta, 1.0), phi[k]);
     236                   
    209237                }
    210238
    211239
     240
    212241                //Update gradient, edge and vertex values using phi limiter
    213                 x_gradient[k] = x_gradient[k]*phi;
    214                 y_gradient[k] = y_gradient[k]*phi;
    215 
    216                 edge_values[k3+0] = qc + phi*dqa[0];
    217                 edge_values[k3+1] = qc + phi*dqa[1];
    218                 edge_values[k3+2] = qc + phi*dqa[2];
     242                x_gradient[k] = x_gradient[k]*phi[k];
     243                y_gradient[k] = y_gradient[k]*phi[k];
     244
     245                edge_values[k3+0] = qc + phi[k]*dqa[0];
     246                edge_values[k3+1] = qc + phi[k]*dqa[1];
     247                edge_values[k3+2] = qc + phi[k]*dqa[2];
    219248
    220249
     
    222251                vertex_values[k3+1] = edge_values[k3+2] + edge_values[k3+0] - edge_values[k3+1];
    223252                vertex_values[k3+2] = edge_values[k3+0] + edge_values[k3+1] - edge_values[k3+2];
     253               
    224254
    225255        }
     
    290320}
    291321
     322
     323
     324
    292325int _limit_edges_by_all_neighbours(int N, double beta,
    293                      double* centroid_values,
    294                      double* vertex_values,
    295                      double* edge_values,
    296                      long*   neighbours) {
     326                                   double* centroid_values,
     327                                   double* vertex_values,
     328                                   double* edge_values,
     329                                   long*   neighbours,
     330                                   double* x_gradient,
     331                                   double* y_gradient) {
    297332
    298333        int i, k, k2, k3, k6;
     
    334369                }
    335370   
    336                 //Update vertex and edge values using phi limiter
     371                //Update gradient, vertex and edge values using phi limiter
     372                x_gradient[k] = x_gradient[k]*phi;
     373                y_gradient[k] = y_gradient[k]*phi;
     374
    337375                edge_values[k3+0] = qc + phi*dqa[0];
    338376                edge_values[k3+1] = qc + phi*dqa[1];
     
    353391                     double* vertex_values,
    354392                     double* edge_values,
     393                     long*   neighbours) {
     394
     395        int i, k, k2, k3, k6;
     396        long n;
     397        double qmin, qmax, qn, qc;
     398        double dq, dqa[3], phi, r;
     399
     400        for (k=0; k<N; k++){
     401                k6 = 6*k;
     402                k3 = 3*k;
     403                k2 = 2*k;
     404
     405                qc = centroid_values[k];
     406                phi = 1.0;
     407
     408                for (i=0; i<3; i++) {
     409                    dq = edge_values[k3+i] - qc;     //Delta between edge and centroid values
     410                    dqa[i] = dq;                      //Save dqa for use in updating vertex values
     411
     412                    n = neighbours[k3+i];
     413                    qn = qc;
     414                    if (n >= 0)  qn = centroid_values[n]; //Neighbour's centroid value
     415
     416                    qmin = min(qc, qn);
     417                    qmax = max(qc, qn);
     418
     419                    r = 1.0;
     420     
     421                    if (dq > 0.0) r = (qmax - qc)/dq;
     422                    if (dq < 0.0) r = (qmin - qc)/dq;     
     423                   
     424                    phi = min( min(r*beta, 1.0), phi);   
     425                   
     426                }
     427
     428
     429                //Update edge and vertex values using phi limiter
     430                edge_values[k3+0] = qc + phi*dqa[0];
     431                edge_values[k3+1] = qc + phi*dqa[1];
     432                edge_values[k3+2] = qc + phi*dqa[2];
     433               
     434                vertex_values[k3+0] = edge_values[k3+1] + edge_values[k3+2] - edge_values[k3+0];
     435                vertex_values[k3+1] = edge_values[k3+2] + edge_values[k3+0] - edge_values[k3+1];
     436                vertex_values[k3+2] = edge_values[k3+0] + edge_values[k3+1] - edge_values[k3+2];
     437
     438        }
     439
     440        return 0;
     441}
     442
     443
     444int _limit_gradient_by_neighbour(int N, double beta,
     445                     double* centroid_values,
     446                     double* vertex_values,
     447                     double* edge_values,
     448                     double* x_gradient,
     449                     double* y_gradient,
    355450                     long*   neighbours) {
    356451
     
    403498}
    404499
    405 
    406 int _limit_gradient_by_neighbour(int N, double beta,
     500int _bound_vertices_below_by_constant(int N, double bound,
    407501                     double* centroid_values,
    408502                     double* vertex_values,
    409503                     double* edge_values,
    410504                     double* x_gradient,
    411                      double* y_gradient,
    412                      long*   neighbours) {
     505                     double* y_gradient) {
    413506
    414507        int i, k, k2, k3, k6;
    415         long n;
    416         double qmin, qmax, qn, qc;
     508        double qmin, qc;
    417509        double dq, dqa[3], phi, r;
    418510
     
    423515
    424516                qc = centroid_values[k];
     517                qmin = bound;
     518
     519
    425520                phi = 1.0;
    426 
    427                 for (i=0; i<3; i++) {
    428                     dq = edge_values[k3+i] - qc;     //Delta between edge and centroid values
     521                for (i=0; i<3; i++) {   
     522                    r = 1.0;
     523     
     524                    dq = vertex_values[k3+i] - qc;    //Delta between vertex and centroid values
    429525                    dqa[i] = dq;                      //Save dq for use in updating vertex values
    430 
    431                     n = neighbours[k3+i];
    432                     if (n >= 0) {
    433                         qn = centroid_values[n]; //Neighbour's centroid value
    434 
    435                         qmin = min(qc, qn);
    436                         qmax = max(qc, qn);
    437 
    438                         r = 1.0;
    439526     
    440                         if (dq > 0.0) r = (qmax - qc)/dq;
    441                         if (dq < 0.0) r = (qmin - qc)/dq;     
     527                    if (dq < 0.0) r = (qmin - qc)/dq;     
     528 
    442529                   
    443                         phi = min( min(r*beta, 1.0), phi);   
    444                     }
     530                    phi = min( min(r, 1.0), phi);   
    445531                }
    446 
    447 
    448                 //Update edge and vertex values using phi limiter
    449                 edge_values[k3+0] = qc + phi*dqa[0];
    450                 edge_values[k3+1] = qc + phi*dqa[1];
    451                 edge_values[k3+2] = qc + phi*dqa[2];
     532   
     533
     534                //Update gradient, vertex and edge values using phi limiter
     535                x_gradient[k] = x_gradient[k]*phi;
     536                y_gradient[k] = y_gradient[k]*phi;
     537
     538                vertex_values[k3+0] = qc + phi*dqa[0];
     539                vertex_values[k3+1] = qc + phi*dqa[1];
     540                vertex_values[k3+2] = qc + phi*dqa[2];
    452541               
    453                 vertex_values[k3+0] = edge_values[k3+1] + edge_values[k3+2] - edge_values[k3+0];
    454                 vertex_values[k3+1] = edge_values[k3+2] + edge_values[k3+0] - edge_values[k3+1];
    455                 vertex_values[k3+2] = edge_values[k3+0] + edge_values[k3+1] - edge_values[k3+2];
     542                edge_values[k3+0] = 0.5*(vertex_values[k3+1] + vertex_values[k3+2]);
     543                edge_values[k3+1] = 0.5*(vertex_values[k3+2] + vertex_values[k3+0]);
     544                edge_values[k3+2] = 0.5*(vertex_values[k3+0] + vertex_values[k3+1]);
    456545
    457546        }
     
    460549}
    461550
    462 
     551int _bound_vertices_below_by_quantity(int N,
     552                                      double* bound_vertex_values,
     553                                      double* centroid_values,
     554                                      double* vertex_values,
     555                                      double* edge_values,
     556                                      double* x_gradient,
     557                                      double* y_gradient) {
     558
     559        int i, k, k2, k3, k6;
     560        double qmin, qc;
     561        double dq, dqa[3], phi, r;
     562
     563        for (k=0; k<N; k++){
     564                k6 = 6*k;
     565                k3 = 3*k;
     566                k2 = 2*k;
     567
     568                qc = centroid_values[k];
     569
     570                phi = 1.0;
     571                for (i=0; i<3; i++) {   
     572                    r = 1.0;
     573     
     574                    dq = vertex_values[k3+i] - qc;    //Delta between vertex and centroid values
     575                    dqa[i] = dq;                      //Save dq for use in updating vertex values
     576     
     577                    qmin = bound_vertex_values[k3+i];
     578                    if (dq < 0.0) r = (qmin - qc)/dq;     
     579 
     580                   
     581                    phi = min( min(r, 1.0), phi);   
     582                }
     583   
     584
     585                //Update gradient, vertex and edge values using phi limiter
     586                x_gradient[k] = x_gradient[k]*phi;
     587                y_gradient[k] = y_gradient[k]*phi;
     588
     589                vertex_values[k3+0] = qc + phi*dqa[0];
     590                vertex_values[k3+1] = qc + phi*dqa[1];
     591                vertex_values[k3+2] = qc + phi*dqa[2];
     592               
     593                edge_values[k3+0] = 0.5*(vertex_values[k3+1] + vertex_values[k3+2]);
     594                edge_values[k3+1] = 0.5*(vertex_values[k3+2] + vertex_values[k3+0]);
     595                edge_values[k3+2] = 0.5*(vertex_values[k3+0] + vertex_values[k3+1]);
     596
     597        }
     598
     599        return 0;
     600}
    463601
    464602int _interpolate_from_vertices_to_edges(int N,
     
    9061044
    9071045
    908 /* PyObject *compute_gradients(PyObject *self, PyObject *args) { */
    909 /*   //"""Compute gradients of triangle surfaces defined by centroids of */
    910 /*   //neighbouring volumes. */
    911 /*   //If one edge is on the boundary, use own centroid as neighbour centroid. */
    912 /*   //If two or more are on the boundary, fall back to first order scheme. */
    913 /*   //""" */
    914 
    915 
    916 /*      PyObject *quantity, *domain, *R; */
    917 /*      PyArrayObject */
    918 /*              *centroids,            //Coordinates at centroids */
    919 /*              *centroid_values,      //Values at centroids */
    920 /*              *number_of_boundaries, //Number of boundaries for each triangle */
    921 /*              *surrogate_neighbours, //True neighbours or - if one missing - self */
    922 /*              *a, *b;                //Return values */
    923 
    924 /*      int dimensions[1], N, err; */
    925 
    926 /*      // Convert Python arguments to C */
    927 /*      if (!PyArg_ParseTuple(args, "O", &quantity)) { */
    928 /*        PyErr_SetString(PyExc_RuntimeError,  */
    929 /*                        "quantity_ext.c: compute_gradients could not parse input");    */
    930 /*        return NULL; */
    931 /*      } */
    932 
    933 /*      domain = PyObject_GetAttrString(quantity, "domain"); */
    934 /*      if (!domain) { */
    935 /*        PyErr_SetString(PyExc_RuntimeError,  */
    936 /*                        "compute_gradients could not obtain domain object from quantity"); */
    937 /*        return NULL; */
    938 /*      } */
    939 
    940 /*      // Get pertinent variables */
    941 
    942 /*      centroids = get_consecutive_array(domain, "centroid_coordinates"); */
    943 /*      centroid_values = get_consecutive_array(quantity, "centroid_values"); */
    944 /*      surrogate_neighbours = get_consecutive_array(domain, "surrogate_neighbours"); */
    945 /*      number_of_boundaries = get_consecutive_array(domain, "number_of_boundaries"); */
    946 
    947 /*      N = centroid_values -> dimensions[0]; */
    948 
    949 /*      // Release */
    950 /*      Py_DECREF(domain); */
    951 
    952 /*      // Allocate space for return vectors a and b (don't DECREF) */
    953 /*      dimensions[0] = N; */
    954 /*      a = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE); */
    955 /*      b = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE); */
    956 
    957 
    958 
    959 /*      err = _compute_gradients(N, */
    960 /*                      (double*) centroids -> data, */
    961 /*                      (double*) centroid_values -> data, */
    962 /*                      (long*) number_of_boundaries -> data, */
    963 /*                      (long*) surrogate_neighbours -> data, */
    964 /*                      (double*) a -> data, */
    965 /*                      (double*) b -> data); */
    966 
    967 /*      if (err != 0) { */
    968 /*        PyErr_SetString(PyExc_RuntimeError, "Gradient could not be computed"); */
    969 /*        return NULL; */
    970 /*      } */
    971 
    972 /*      // Release */
    973 /*      Py_DECREF(centroids); */
    974 /*      Py_DECREF(centroid_values); */
    975 /*      Py_DECREF(number_of_boundaries); */
    976 /*      Py_DECREF(surrogate_neighbours); */
    977 
    978 /*      // Build result, release and return */
    979 /*      R = Py_BuildValue("OO", PyArray_Return(a), PyArray_Return(b)); */
    980 /*      Py_DECREF(a); */
    981 /*      Py_DECREF(b); */
    982 /*      return R; */
    983 /* } */
    984 
    985 
    986 
    9871046PyObject *extrapolate_from_gradient(PyObject *self, PyObject *args) {
    9881047
     
    10871146        *quantity_vertex_values,        //Values at vertices
    10881147        *quantity_edge_values,          //Values at edges
     1148        *quantity_phi,                  //limiter phi values
    10891149        *quantity_x_gradient,           //x gradient
    10901150        *quantity_y_gradient;           //y gradient
     
    10971157   
    10981158  // Convert Python arguments to C
    1099   if (!PyArg_ParseTuple(args, "OOd",
    1100                         &domain,
    1101                         &quantity,
    1102                         &beta)) {
     1159  if (!PyArg_ParseTuple(args, "O",&quantity)) {
    11031160      PyErr_SetString(PyExc_RuntimeError,
    11041161                      "quantity_ext.c: extrapolate_second_order_and_limit could not parse input");     
     1162      return NULL;
     1163  }
     1164
     1165  domain = PyObject_GetAttrString(quantity, "domain");
     1166  if (!domain) {
     1167      PyErr_SetString(PyExc_RuntimeError,
     1168                      "quantity_ext.c: extrapolate_second_order_and_limit could not obtain domain object from quantity");       
    11051169      return NULL;
    11061170  }
     
    11171181  quantity_vertex_values        = get_consecutive_array(quantity, "vertex_values");
    11181182  quantity_edge_values          = get_consecutive_array(quantity, "edge_values");
     1183  quantity_phi                  = get_consecutive_array(quantity, "phi");
    11191184  quantity_x_gradient           = get_consecutive_array(quantity, "x_gradient");
    11201185  quantity_y_gradient           = get_consecutive_array(quantity, "y_gradient");
    11211186
     1187  beta = get_python_double(quantity,"beta");
     1188
    11221189  ntri = quantity_centroid_values -> dimensions[0];
    1123 
    1124 
    11251190
    11261191  err = _compute_gradients(ntri,
     
    11391204
    11401205
    1141   err = _extrapolate_and_limit_from_gradient(ntri, beta,
    1142                                              (double*) domain_centroids -> data,
    1143                                              (long*)   domain_neighbours -> data,
    1144                                              (double*) quantity_centroid_values -> data,
    1145                                              (double*) domain_vertex_coordinates -> data,
    1146                                              (double*) quantity_vertex_values -> data,
    1147                                              (double*) quantity_edge_values -> data,
    1148                                              (double*) quantity_x_gradient -> data,
    1149                                              (double*) quantity_y_gradient -> data);
     1206  err = _extrapolate_from_gradient(ntri,
     1207                                   (double*) domain_centroids -> data,
     1208                                   (double*) quantity_centroid_values -> data,
     1209                                   (double*) domain_vertex_coordinates -> data,
     1210                                   (double*) quantity_vertex_values -> data,
     1211                                   (double*) quantity_edge_values -> data,
     1212                                   (double*) quantity_x_gradient -> data,
     1213                                   (double*) quantity_y_gradient -> data);
    11501214
    11511215  if (err != 0) {
    11521216      PyErr_SetString(PyExc_RuntimeError,
    1153                           "quantity_ext.c: Internal function _extrapolate_and_limit_from_gradient failed");
     1217                          "quantity_ext.c: Internal function _extrapolate_from_gradient failed");
     1218      return NULL;
     1219  }
     1220
     1221
     1222  err = _limit_edges_by_all_neighbours(ntri, beta,
     1223                                       (double*) quantity_centroid_values -> data,
     1224                                       (double*) quantity_vertex_values -> data,
     1225                                       (double*) quantity_edge_values -> data,
     1226                                       (long*)   domain_neighbours -> data,
     1227                                       (double*) quantity_x_gradient -> data,
     1228                                       (double*) quantity_y_gradient -> data);
     1229
     1230  if (err != 0) {
     1231      PyErr_SetString(PyExc_RuntimeError,
     1232                          "quantity_ext.c: Internal function _limit_edges_by_all_neighbours failed");
    11541233      return NULL;
    11551234  }
     
    11651244  Py_DECREF(quantity_vertex_values);
    11661245  Py_DECREF(quantity_edge_values);
     1246  Py_DECREF(quantity_phi);
    11671247  Py_DECREF(quantity_x_gradient);
    11681248  Py_DECREF(quantity_y_gradient);
     
    14451525  //
    14461526
    1447         PyObject *quantity, *domain, *Tmp;
     1527        PyObject *quantity, *domain;
    14481528        PyArrayObject
    14491529            *vertex_values,   //Conserved quantities at vertices
    14501530            *centroid_values, //Conserved quantities at centroids
    14511531            *edge_values,     //Conserved quantities at edges
     1532            *x_gradient,
     1533            *y_gradient,
    14521534            *neighbours;
    14531535
     
    14631545        }
    14641546
    1465         domain = PyObject_GetAttrString(quantity, "domain");
    1466         if (!domain) {
    1467           PyErr_SetString(PyExc_RuntimeError,
    1468                           "quantity_ext.c: limit_edges_by_all_neighbours could not obtain domain object from quantity");                       
    1469          
    1470           return NULL;
    1471         }
    1472 
    1473         // Get safety factor beta_w
    1474         Tmp = PyObject_GetAttrString(domain, "beta_w");
    1475         if (!Tmp) {
    1476           PyErr_SetString(PyExc_RuntimeError,
    1477                           "quantity_ext.c: limit_edges_by_all_neighbours could not obtain beta_w object from domain");                 
    1478          
    1479           return NULL;
    1480         }       
     1547        domain = get_python_object(quantity, "domain");
    14811548
    14821549
     
    14861553        vertex_values    = get_consecutive_array(quantity, "vertex_values");
    14871554        edge_values      = get_consecutive_array(quantity, "edge_values");
    1488         beta_w           = PyFloat_AsDouble(Tmp);
    1489 
     1555        x_gradient       = get_consecutive_array(quantity, "x_gradient");
     1556        y_gradient       = get_consecutive_array(quantity, "y_gradient");
     1557        beta_w           = get_python_double(domain,"beta_w");
     1558
     1559        Py_DECREF(domain);
    14901560
    14911561        N = centroid_values -> dimensions[0];
     
    14951565                                             (double*) vertex_values -> data,
    14961566                                             (double*) edge_values -> data,
    1497                                              (long*)   neighbours -> data);
     1567                                             (long*)   neighbours -> data,
     1568                                             (double*) x_gradient -> data,
     1569                                             (double*) y_gradient -> data);
    14981570
    14991571        if (err != 0) {
    15001572          PyErr_SetString(PyExc_RuntimeError,
    1501                           "Internal function _limit_by_vertex failed");
     1573                          "quantity_ect.c: limit_edges_by_all_neighbours internal function _limit_edges_by_all_neighbours failed");
    15021574          return NULL;
    15031575        }       
     
    15091581        Py_DECREF(vertex_values);
    15101582        Py_DECREF(edge_values);
    1511         Py_DECREF(Tmp);
     1583
     1584
     1585
     1586        return Py_BuildValue("");
     1587}
     1588
     1589PyObject *bound_vertices_below_by_constant(PyObject *self, PyObject *args) {
     1590  //Bound a quantity below by a contant (useful for ensuring positivity
     1591  //precondition:
     1592  //  vertex values are already calulated, gradient consistent
     1593  //postcondition:
     1594  //  gradient, vertex and edge values are updated
     1595  //
     1596
     1597        PyObject *quantity, *domain;
     1598        PyArrayObject
     1599            *vertex_values,   //Conserved quantities at vertices
     1600            *centroid_values, //Conserved quantities at centroids
     1601            *edge_values,     //Conserved quantities at edges
     1602            *x_gradient,
     1603            *y_gradient;
     1604
     1605        double bound; //Safety factor
     1606        int N, err;
     1607
     1608
     1609        // Convert Python arguments to C
     1610        if (!PyArg_ParseTuple(args, "Od", &quantity, &bound)) {
     1611          PyErr_SetString(PyExc_RuntimeError,
     1612                          "quantity_ext.c: bound_vertices_below_by_constant could not parse input");
     1613          return NULL;
     1614        }
     1615
     1616        domain = get_python_object(quantity, "domain");
     1617
     1618        // Get pertinent variables
     1619        centroid_values  = get_consecutive_array(quantity, "centroid_values");
     1620        vertex_values    = get_consecutive_array(quantity, "vertex_values");
     1621        edge_values      = get_consecutive_array(quantity, "edge_values");
     1622        x_gradient       = get_consecutive_array(quantity, "x_gradient");
     1623        y_gradient       = get_consecutive_array(quantity, "y_gradient");
     1624
     1625
     1626
     1627        Py_DECREF(domain);
     1628
     1629        N = centroid_values -> dimensions[0];
     1630
     1631        err = _bound_vertices_below_by_constant(N, bound,
     1632                                             (double*) centroid_values -> data,
     1633                                             (double*) vertex_values -> data,
     1634                                             (double*) edge_values -> data,
     1635                                             (double*) x_gradient -> data,
     1636                                             (double*) y_gradient -> data);
     1637
     1638        if (err != 0) {
     1639          PyErr_SetString(PyExc_RuntimeError,
     1640                          "quantity_ect.c: bound_vertices_below_by_constant internal function _bound_vertices_below_by_constant failed");
     1641          return NULL;
     1642        }       
     1643
     1644
     1645        // Release
     1646        Py_DECREF(centroid_values);
     1647        Py_DECREF(vertex_values);
     1648        Py_DECREF(edge_values);
     1649        Py_DECREF(x_gradient);
     1650        Py_DECREF(y_gradient);
     1651
     1652
     1653
     1654        return Py_BuildValue("");
     1655}
     1656
     1657
     1658PyObject *bound_vertices_below_by_quantity(PyObject *self, PyObject *args) {
     1659  //Bound a quantity below by a contant (useful for ensuring positivity
     1660  //precondition:
     1661  //  vertex values are already calulated, gradient consistent
     1662  //postcondition:
     1663  //  gradient, vertex and edge values are updated
     1664  //
     1665
     1666        PyObject *quantity, *bounding_quantity, *domain;
     1667        PyArrayObject
     1668            *vertex_values,   //Conserved quantities at vertices
     1669            *centroid_values, //Conserved quantities at centroids
     1670            *edge_values,     //Conserved quantities at edges
     1671            *x_gradient,
     1672            *y_gradient,
     1673            *bound_vertex_values;
     1674
     1675        int N, err;
     1676
     1677
     1678        // Convert Python arguments to C
     1679        if (!PyArg_ParseTuple(args, "OO", &quantity, &bounding_quantity)) {
     1680          PyErr_SetString(PyExc_RuntimeError,
     1681                          "quantity_ext.c: bound_vertices_below_by_quantity could not parse input");
     1682          return NULL;
     1683        }
     1684
     1685        domain = get_python_object(quantity, "domain");
     1686
     1687        // Get pertinent variables
     1688        centroid_values     = get_consecutive_array(quantity, "centroid_values");
     1689        vertex_values       = get_consecutive_array(quantity, "vertex_values");
     1690        edge_values         = get_consecutive_array(quantity, "edge_values");
     1691        x_gradient          = get_consecutive_array(quantity, "x_gradient");
     1692        y_gradient          = get_consecutive_array(quantity, "y_gradient");
     1693        bound_vertex_values = get_consecutive_array(bounding_quantity, "vertex_values");
     1694
     1695
     1696
     1697        Py_DECREF(domain);
     1698
     1699        N = centroid_values -> dimensions[0];
     1700
     1701        err = _bound_vertices_below_by_quantity(N,
     1702                                                (double*) bound_vertex_values -> data,
     1703                                                (double*) centroid_values -> data,
     1704                                                (double*) vertex_values -> data,
     1705                                                (double*) edge_values -> data,
     1706                                                (double*) x_gradient -> data,
     1707                                                (double*) y_gradient -> data);
     1708
     1709        if (err != 0) {
     1710          PyErr_SetString(PyExc_RuntimeError,
     1711                          "quantity_ect.c: bound_vertices_below_by_quantity internal function _bound_vertices_below_by_quantity failed");
     1712          return NULL;
     1713        }       
     1714
     1715
     1716        // Release
     1717        Py_DECREF(centroid_values);
     1718        Py_DECREF(vertex_values);
     1719        Py_DECREF(edge_values);
     1720        Py_DECREF(x_gradient);
     1721        Py_DECREF(y_gradient);
     1722        Py_DECREF(bound_vertex_values);
     1723
    15121724
    15131725
     
    15891801
    15901802        // Release
     1803        Py_DECREF(domain);
    15911804        Py_DECREF(neighbours);
    15921805        Py_DECREF(centroid_values);
     
    17001913        {"limit_edges_by_neighbour", limit_edges_by_neighbour, METH_VARARGS, "Print out"},
    17011914        {"limit_gradient_by_neighbour", limit_gradient_by_neighbour, METH_VARARGS, "Print out"},
     1915        {"bound_vertices_below_by_constant", bound_vertices_below_by_constant, METH_VARARGS, "Print out"},
     1916        {"bound_vertices_below_by_quantity", bound_vertices_below_by_quantity, METH_VARARGS, "Print out"},
    17021917        {"update", update, METH_VARARGS, "Print out"},
    17031918        {"backup_centroid_values", backup_centroid_values, METH_VARARGS, "Print out"},
  • anuga_core/source/anuga/advection/advection_ext.c

    r4978 r5162  
    22//
    33// To compile (Python2.X):
    4 // use python ../utilities/compile.py
    5 //
    6 // See the module advection.py
     4// use python ../utilities/compile.py to
     5// compile all C files within a directory
    76//
    87//
     
    1918
    2019
     20//-------------------------------------------
     21// Low level routines (called from wrappers)
     22//------------------------------------------
    2123
    2224double _compute_fluxes(
     
    3840
    3941 
    40         //Loop
     42        //Local Variables
    4143
    4244        double qr,ql;
     
    5153        int k3;
    5254       
     55        //Loop through triangles
     56
    5357        timestep = max_timestep;
    5458
     
    5963            for (i=0; i<3; i++){
    6064                k_i = k3+i;
    61                 //Quantities inside volume facing neighbour i
     65                //Quantities inside triangle facing neighbour i
    6266                ql = quantity_edge[k_i];
    6367
     
    144148
    145149  PyObject *domain, *quantity;
     150 
    146151  PyArrayObject
    147152    * quantity_update,
  • anuga_core/source/anuga/config.py

    r4837 r5162  
    100100#timestepping_method = 'rk2'   # 2nd Order TVD scheme
    101101
     102# rk2 is a little more stable than euler, so rk2 timestepping
     103# can deal with a larger beta when slope limiting the reconstructed
     104# solution. The large beta is needed if solving problems sensitive
     105# to numerical diffusion, like a small forced wave in an ocean
     106beta_euler = 1.0
     107beta_rk2   = 1.6
     108
     109
     110
    102111# Option to search for signatures where isolated triangles are
    103112# responsible for a small global timestep.
  • anuga_core/source/anuga/shallow_water/shallow_water_domain.py

    r5089 r5162  
    103103from anuga.config import alpha_balance
    104104from anuga.config import optimise_dry_cells
     105from anuga.config import optimised_gradient_limiter
    105106
    106107#---------------------
     
    176177        self.minimum_storable_height = minimum_storable_height
    177178        self.quantities_to_be_stored = ['stage','xmomentum','ymomentum']
    178        
     179
     180        # Limiters
     181        self.use_old_limiter = True
     182
     183        self.optimised_gradient_limiter = optimised_gradient_limiter
    179184               
    180185
     
    187192        self.beta_w      = beta
    188193        self.beta_w_dry  = beta
     194        self.quantities['stage'].beta = beta
     195       
    189196        self.beta_uh     = beta
    190197        self.beta_uh_dry = beta
     198        self.quantities['xmomentum'].beta = beta
     199       
    191200        self.beta_vh     = beta
    192201        self.beta_vh_dry = beta
     202        self.quantities['ymomentum'].beta = beta
     203       
    193204        self.beta_h      = beta
    194205       
     
    384395        # Call correct module function
    385396        # (either from this module or C-extension)
    386         distribute_to_vertices_and_edges(self)
     397        if self.use_old_limiter:
     398            distribute_to_vertices_and_edges(self)
     399        else:
     400            for name in self.conserved_quantities:
     401                Q = self.quantities[name]
     402                if self._order_ == 1:
     403                    Q.extrapolate_first_order()
     404                elif self._order_ == 2:
     405                    Q.extrapolate_second_order_and_limit()
     406                    if name == 'stage':
     407                        Q.bound_vertices_below_by_quantity(self.quantities['elevation'])
     408                else:
     409                    raise 'Unknown order'
     410 
    387411
    388412
     
    399423        # self.check_integrity()
    400424
    401         msg = 'Parameter beta_h must be in the interval [0, 1['
    402         assert 0 <= self.beta_h <= 1.0, msg
    403         msg = 'Parameter beta_w must be in the interval [0, 1['
    404         assert 0 <= self.beta_w <= 1.0, msg
     425        msg = 'Parameter beta_h must be in the interval [0, 2['
     426        assert 0 <= self.beta_h <= 2.0, msg
     427        msg = 'Parameter beta_w must be in the interval [0, 2['
     428        assert 0 <= self.beta_w <= 2.0, msg
    405429
    406430
     
    737761    """
    738762
    739     from anuga.config import optimised_gradient_limiter
     763   
    740764
    741765    # Remove very thin layers of water
     
    743767
    744768    # Extrapolate all conserved quantities
    745     if optimised_gradient_limiter:
     769    if domain.optimised_gradient_limiter:
    746770        # MH090605 if second order,
    747771        # perform the extrapolation and limiting on
     
    764788                Q.extrapolate_first_order()
    765789            elif domain._order_ == 2:
    766 
    767                 # Experiment
    768                 #if name == 'stage':
    769                 #    #print name, 'second'
    770                 #    Q.extrapolate_second_order()
    771                 #    Q.limit()
    772                 #else:
    773                 #    #print name, 'first'               
    774                 #    Q.extrapolate_first_order()
    775                 #    #Q.extrapolate_second_order()
    776                 #    #Q.limit()               
    777                
    778790                Q.extrapolate_second_order()
    779791                Q.limit()
  • anuga_core/source/anuga/shallow_water/shallow_water_ext.c

    r4978 r5162  
    10881088  double dx1,dx2,dy1,dy2,dxv0,dxv1,dxv2,dyv0,dyv1,dyv2,dq0,dq1,dq2,area2;
    10891089  double dqv[3], qmin, qmax, hmin, hmax;
    1090   double hc, h0, h1, h2, beta_tmp;
     1090  double hc, h0, h1, h2, beta_tmp, hfactor;
    10911091
    10921092       
     
    11851185      h1 = stage_centroid_values[k1] - elevation_centroid_values[k1];
    11861186      h2 = stage_centroid_values[k2] - elevation_centroid_values[k2];
    1187       hmin = min(h0,min(h1,h2));
    1188      
     1187      hmin = min(min(h0,min(h1,h2)),hc);
     1188      //hfactor = hc/(hc + 1.0);
     1189
     1190      hfactor = 0.0;
     1191      if (hmin > 0.1 ) {
     1192          hfactor = (hmin-0.1)/(hmin+0.4);
     1193      }
    11891194     
    11901195      if (optimise_dry_cells) {     
     
    12301235     
    12311236      // Playing with dry wet interface
    1232       hmin = qmin;
    1233       beta_tmp = beta_w;
    1234       if (hmin<minimum_allowed_height)
    1235         beta_tmp = beta_w_dry;
     1237      //hmin = qmin;
     1238      //beta_tmp = beta_w_dry;
     1239      //if (hmin>minimum_allowed_height)
     1240      beta_tmp = beta_w_dry + (beta_w - beta_w_dry) * hfactor;
    12361241       
     1242      //printf("min_alled_height = %f\n",minimum_allowed_height);
     1243      //printf("hmin = %f\n",hmin);
     1244      //printf("beta_w = %f\n",beta_w);
     1245      //printf("beta_tmp = %f\n",beta_tmp);
    12371246      // Limit the gradient
    12381247      limit_gradient(dqv,qmin,qmax,beta_tmp);
     
    12711280      // from the centroid to the min and max
    12721281      find_qmin_and_qmax(dq0,dq1,dq2,&qmin,&qmax);
    1273       beta_tmp = beta_uh;
    1274       if (hmin<minimum_allowed_height)
    1275         beta_tmp = beta_uh_dry;
    1276        
     1282      //beta_tmp = beta_uh;
     1283      //if (hmin<minimum_allowed_height)
     1284      //beta_tmp = beta_uh_dry;
     1285      beta_tmp = beta_uh_dry + (beta_uh - beta_uh_dry) * hfactor;
     1286
    12771287      // Limit the gradient
    12781288      limit_gradient(dqv,qmin,qmax,beta_tmp);
     
    13111321      // from the centroid to the min and max
    13121322      find_qmin_and_qmax(dq0,dq1,dq2,&qmin,&qmax);
    1313       beta_tmp = beta_vh;
    1314      
    1315       if (hmin<minimum_allowed_height)
    1316         beta_tmp = beta_vh_dry;
    1317        
     1323     
     1324      //beta_tmp = beta_vh;
     1325      //
     1326      //if (hmin<minimum_allowed_height)
     1327      //beta_tmp = beta_vh_dry;
     1328      beta_tmp = beta_vh_dry + (beta_vh - beta_vh_dry) * hfactor;       
     1329
    13181330      // Limit the gradient
    13191331      limit_gradient(dqv,qmin,qmax,beta_tmp);
     
    15631575 
    15641576
    1565   beta_w      = get_double(domain,"beta_w");
    1566   beta_w_dry  = get_double(domain,"beta_w_dry");
    1567   beta_uh     = get_double(domain,"beta_uh");
    1568   beta_uh_dry = get_double(domain,"beta_uh_dry");
    1569   beta_vh     = get_double(domain,"beta_vh");
    1570   beta_vh_dry = get_double(domain,"beta_vh_dry"); 
    1571 
    1572   minimum_allowed_height = get_double(domain,"minimum_allowed_height");
    1573   epsilon = get_double(domain,"epsilon");
     1577  beta_w                 = get_python_double(domain,"beta_w");
     1578  beta_w_dry             = get_python_double(domain,"beta_w_dry");
     1579  beta_uh                = get_python_double(domain,"beta_uh");
     1580  beta_uh_dry            = get_python_double(domain,"beta_uh_dry");
     1581  beta_vh                = get_python_double(domain,"beta_vh");
     1582  beta_vh_dry            = get_python_double(domain,"beta_vh_dry"); 
     1583
     1584  minimum_allowed_height = get_python_double(domain,"minimum_allowed_height");
     1585  epsilon                = get_python_double(domain,"epsilon");
    15741586
    15751587  number_of_elements = stage_centroid_values -> dimensions[0]; 
  • anuga_core/source/anuga/shallow_water/test_data_manager.py

    r4924 r5162  
    248248       
    249249        range = fid.variables['xmomentum_range'][:]
    250         assert allclose(range,[0,0.4695096]) or\
    251                allclose(range,[0,0.47790655]) # Old slope limiters
     250
     251
     252
     253        assert allclose(range,[0,0.4695096]) or \
     254               allclose(range,[0,0.47790655]) or\
     255               allclose(range,[0,0.46941957])
    252256       
    253257        range = fid.variables['ymomentum_range'][:]
    254258        #assert allclose(range,[0,0.02174380])
     259       
    255260        assert allclose(range,[0,0.02174439]) or\
    256                allclose(range,[0,0.02283983]) # Old slope limiters
     261               allclose(range,[0,0.02283983]) or\
     262               allclose(range,[0,0.0217342]) # Old slope limiters
    257263       
    258264        fid.close()
     
    324330
    325331        extrema = fid.variables['xmomentum.extrema'][:]
    326         assert allclose(extrema,[-0.06062178, 0.47886313])
     332        assert allclose(extrema,[-0.06062178, 0.47873023])
    327333       
    328334        extrema = fid.variables['ymomentum.extrema'][:]
    329         assert allclose(extrema,[0.00, 0.06241221])       
     335        assert allclose(extrema,[0.00, 0.0625786])       
    330336
    331337        time_interval = fid.variables['extrema.time_interval'][:]
  • anuga_core/source/anuga/utilities/util_ext.h

    r4978 r5162  
    261261   
    262262  B = (PyArrayObject*) PyObject_GetAttrString(O, name);
    263   if (!B) return NULL;     
     263  if (!B) {
     264    PyErr_SetString(PyExc_RuntimeError, "util_ext.h: get_consecutive_array could not obtain array");
     265    return NULL;
     266  }     
    264267 
    265268  //Convert to consecutive array
     
    269272  Py_DECREF(B); //FIXME: Is this really needed??
    270273 
    271   if (!A) return NULL;
     274  if (!A) {
     275    PyErr_SetString(PyExc_RuntimeError, "util_ext.h: get_consecutive_array could not obtain array");
     276    return NULL;
     277  }
    272278  return A;
    273279}
    274280
    275 double get_double(PyObject *O, char *name) {
     281double get_python_double(PyObject *O, char *name) {
    276282  PyObject *TObject;
    277283  double tmp;
     
    281287  TObject = PyObject_GetAttrString(O, name);
    282288  if (!TObject) {
    283     PyErr_SetString(PyExc_RuntimeError, "util_ext.h: get_double could not obtain double from object");
     289    PyErr_SetString(PyExc_RuntimeError, "util_ext.h: get_python_double could not obtain double from object");
    284290    return 0.0;
    285291  } 
     
    293299
    294300
    295 
     301PyObject *get_python_object(PyObject *O, char *name) {
     302  PyObject *Oout;
     303
     304  Oout = PyObject_GetAttrString(O, name);
     305  if (!Oout) {
     306    PyErr_SetString(PyExc_RuntimeError, "util_ext.h: get_python_object could not obtain object");
     307    return NULL;
     308  }
     309
     310  return Oout;
     311
     312}
     313
     314
     315
  • anuga_core/source/anuga/utilities/xml_tools.py

    r5022 r5162  
    239239        fid = xml
    240240
    241     try:
    242         dom = parse(fid)
    243     except Exception, e:
    244         # Throw filename into dom exception
    245         msg = 'XML file "%s" could not be parsed: ' %fid.name
    246         msg += str(e)
    247         raise Exception, msg
     241    dom = parse(fid)
     242
     243##     try:
     244##         dom = parse(fid)
     245##     except Exception, e:
     246##         # Throw filename into dom exception
     247##         msg = 'XML file "%s" could not be parsed: ' %fid.name
     248##         msg += str(e)
     249##         raise Exception, msg
    248250
    249251    return dom2object(dom)
  • anuga_validation/convergence_study/animatesww2d_alt.py

    r4838 r5162  
    1 def animatesww2d_alt(sww_filename, movie_filename, range):
    2     """Plot cross section of model output
    3     """
     1## def animatesww2d_alt(sww_filename, movie_filename, range):
     2##     """Plot cross section of model output
     3##     """
    44
    5     # Read SWW file     
    6     from Scientific.IO.NetCDF import NetCDFFile
    7     fid = NetCDFFile(sww_filename, 'r')
     5##     # Read SWW file 
     6##     from Scientific.IO.NetCDF import NetCDFFile
     7##     fid = NetCDFFile(sww_filename, 'r')
    88   
    9     x = fid.variables['x']
    10     y = fid.variables['y']
    11     volumes = fid.variables['volumes']
     9##     x = fid.variables['x']
     10##     y = fid.variables['y']
     11##     volumes = fid.variables['volumes']
    1212   
    13     elevation = fid.variables['elevation']
    14     time = fid.variables['time']
    15     stage = fid.variables['stage']
    16     xmomentum = fid.variables['xmomentum']
    17     ymomentum = fid.variables['ymomentum']
     13##     elevation = fid.variables['elevation']
     14##     time = fid.variables['time']
     15##     stage = fid.variables['stage']
     16##     xmomentum = fid.variables['xmomentum']
     17##     ymomentum = fid.variables['ymomentum']
    1818   
    1919def animatesww2d(swwfile):
     
    3434    # interpolation points are for y = 0
    3535    # number of increments in x
    36     m = 100
     36    m = 1000
    3737    max_x = 100000.
    3838    x = 0
  • anuga_validation/convergence_study/convergence_structured.py

    r4992 r5162  
    1717from anuga.shallow_water import Transmissive_boundary
    1818from anuga.shallow_water import Transmissive_Momentum_Set_Stage_boundary
    19 from anuga.geospatial_data.geospatial_data import *
     19#from anuga.geospatial_data.geospatial_data import *
    2020from math import cos
     21from Numeric import zeros, Float
    2122
    2223#------------------------------------------------------------------------------
    2324# Setup computational domain
    2425#------------------------------------------------------------------------------
    25 dx = 1000.
     26dx = 200.
    2627dy = dx
    2728L = 100000.
     
    3334domain = Domain(points, vertices, boundary)
    3435
    35 domain.set_timestepping_method('euler')
     36domain.set_timestepping_method('rk2')
    3637domain.set_default_order(2)
    3738domain.set_name('myexample9')               
    3839domain.set_datadir('.')                     # Use current directory for output
    3940
    40 domain.beta_w      = 1.0
    41 domain.beta_w_dry  = 0.2
    42 domain.beta_uh     = 1.0
    43 domain.beta_uh_dry = 0.2
    44 domain.beta_vh     = 1.0
    45 domain.beta_vh_dry = 0.2
    46 domain.beta_h      = 1.0
     41domain.set_all_limiters(0.9)
     42
     43print domain.beta_w
     44domain.use_old_limiter = False
     45domain.CFL = 1.0
    4746
    4847#------------------------------------------------------------------------------
     
    5049#------------------------------------------------------------------------------
    5150#domain.set_quantity('elevation', topography) # Use function for elevation
    52 domain.set_quantity('elevation',-100)
     51domain.set_quantity('elevation',0.0)
    5352domain.set_quantity('friction', 0.00)
    54 domain.set_quantity('stage', 0.0)           
     53
     54h0 = 10.0
     55h1 = 1.0
     56
     57def height(x,y):
     58    z = zeros(len(x),Float)
     59    for i in range(len(x)):
     60        if x[i]<=50000.0:
     61            z[i] = h0
     62        else:
     63            z[i] = h1
     64    return z
     65
     66
     67domain.set_quantity('stage', height)
     68#domain.set_quantity('stage', 0.0)           
    5569
    5670#-----------------------------------------------------------------------------
     
    6579Bw = Time_boundary(domain=domain,     # Time dependent boundary 
    6680## Sine wave
    67                    f=lambda t: [(-amplitude*sin((1./300.)*t*2*pi)), 0.0, 0.0])
     81#                   f=lambda t: [(-amplitude*sin((1./300.)*t*2*pi)), 0.0, 0.0])
     82## Single wave
     83                   f=lambda t: [h0, 0.0, 0.0])
    6884## Sawtooth?
    6985#                   f=lambda t: [(-8.0*(sin((1./180.)*t*2*pi))+(1./2.)*sin((2./180.)*t*2*pi)+(1./3.)*sin((3./180.)*t*2*pi)), 0.0, 0.0])
     
    92108#------------------------------------------------------------------------------
    93109
    94 for t in domain.evolve(yieldstep = 20.0, finaltime = 10*40*60.):
     110for t in domain.evolve(yieldstep = 50.0, finaltime = 10*40*60.):
    95111    domain.write_time()
    96112    vis.update()
  • anuga_work/development/analytical_solutions/oblique_shock.py

    r5030 r5162  
    2929
    3030
    31 from shallow_water import Domain, Constant_height
    32 from shallow_water import Transmissive_boundary, Reflective_boundary,\
     31from anuga.shallow_water import Domain
     32from anuga.shallow_water import Transmissive_boundary, Reflective_boundary,\
    3333     Dirichlet_boundary
    34 
     34from anuga.visualiser import RealtimeVisualiser
    3535from math import sqrt, cos, sin, pi
    36 from mesh_factory import oblique
     36from mesh_factory import oblique_cross
    3737
    3838
     
    4444leny = 30.
    4545lenx = 40.
    46 n = 50
    47 m = 60
     46n = 100
     47m = 120
    4848
    49 points, elements, boundary = oblique(m, n, lenx, leny)
     49points, elements, boundary = oblique_cross(m, n, lenx, leny)
    5050domain = Domain(points, elements, boundary)
    5151
     
    5454
    5555# Store output
    56 domain.store=True
     56#domain.store=True
    5757
    5858# Output format
    59 domain.format="sww" #NET.CDF binary format
     59#domain.format="sww" #NET.CDF binary format
    6060                    # "dat" for ASCII
    6161
     
    8787#Initial condition
    8888h = 0.5
    89 domain.set_quantity('stage', Constant_height(x_slope, h) )
     89domain.set_quantity('stage', expression = 'elevation + %d'%h )
    9090
     91#---------------------------------
     92# Setup visualization
     93#---------------------------------
     94vis = RealtimeVisualiser(domain)
     95vis.render_quantity_height("elevation", dynamic=False)
     96vis.render_quantity_height("stage", zScale=10, dynamic=True)
     97vis.colour_height_quantity('stage', (0.75, 0.5, 0.5))
     98vis.start()
    9199
    92100######################
     
    96104for t in domain.evolve(yieldstep = 0.5, finaltime = 50):
    97105    domain.write_time()
     106    vis.update()
    98107
    99108print 'That took %.2f seconds' %(time.time()-t0)
     109
     110
     111vis.evolveFinished()
     112vis.join()
    100113
    101114#FIXME: Compute average water depth on either side of shock and compare
  • anuga_work/development/near_shore_PMD/run_beach.py

    r5154 r5162  
    167167    domain.set_boundary( {'wall': Br, 'wave': Bwp_velocity} )
    168168
     169
    169170    #-------------------------------------------------------------------------
    170171    # Evolve system through time
     
    174175    for t in domain.evolve(yieldstep, finaltime):   
    175176        domain.write_time()
     177
    176178    print 'That took %.2f seconds' %(time.time()-t0)
    177179    print 'finished'
Note: See TracChangeset for help on using the changeset viewer.