Changeset 8371


Ignore:
Timestamp:
Mar 27, 2012, 8:46:05 PM (13 years ago)
Author:
steve
Message:

working on

Location:
trunk/anuga_core/source/anuga
Files:
6 edited

Legend:

Unmodified
Added
Removed

  • trunk/anuga_core/source/anuga/abstract_2d_finite_volumes/generic_domain.py

    r8360 r8371  
    243243        from anuga.config import default_order
    244244        from anuga.config import max_timestep, min_timestep
    245 
     245        from anuga.config import g
     246
     247        self.g = g
    246248        self.beta_w = beta_w
    247249        self.epsilon = epsilon

  • trunk/anuga_core/source/anuga/shallow_water/shallow_water_domain.py

    r8366 r8371  
    169169        #-------------------------------
    170170        self.forcing_terms.append(manning_friction_implicit)
    171         self.forcing_terms.append(gravity)
     171        #self.forcing_terms.append(gravity)
     172        self.forcing_terms.append(gravity_new)
    172173
    173174
    … …  
    12961297################################################################################
    12971298
     1299#def gravity(domain):
     1300#    """Apply gravitational pull in the presence of bed slope
     1301#    Wrapper calls underlying C implementation
     1302#    """
     1303#
     1304#    from shallow_water_ext import gravity as gravity_c
     1305#
     1306#    xmom_update = domain.quantities['xmomentum'].explicit_update
     1307#    ymom_update = domain.quantities['ymomentum'].explicit_update
     1308#
     1309#    stage = domain.quantities['stage']
     1310#    elevation = domain.quantities['elevation']
     1311#
     1312#    #FIXME SR Should avoid allocating memory!
     1313#    height = stage.centroid_values - elevation.centroid_values
     1314#
     1315#
     1316#    elevation = elevation.vertex_values
     1317#
     1318#    point = domain.get_vertex_coordinates()
     1319#
     1320#    gravity_c(domain.g, height, elevation, point, xmom_update, ymom_update)
     1321
     1322
     1323
     1324def gravity_new(domain):
     1325    """Apply gravitational pull in the presence of bed slope
     1326    Wrapper calls underlying C implementation
     1327    """
     1328
     1329    from shallow_water_ext import gravity_new as gravity_c
     1330
     1331    gravity_c(domain)
     1332
     1333
    12981334def gravity(domain):
    12991335    """Apply gravitational pull in the presence of bed slope
    … …  
    13031339    from shallow_water_ext import gravity as gravity_c
    13041340
    1305     xmom_update = domain.quantities['xmomentum'].explicit_update
    1306     ymom_update = domain.quantities['ymomentum'].explicit_update
    1307 
    1308     stage = domain.quantities['stage']
    1309     elevation = domain.quantities['elevation']
    1310 
    1311     #FIXME SR Should avoid allocating memory!
    1312     height = stage.centroid_values - elevation.centroid_values
    1313     elevation = elevation.vertex_values
    1314 
    1315     point = domain.get_vertex_coordinates()
    1316 
    1317     gravity_c(domain.g, height, elevation, point, xmom_update, ymom_update)
     1341    gravity_c(domain)
    13181342
    13191343def manning_friction_implicit(domain):

  • trunk/anuga_core/source/anuga/shallow_water/shallow_water_ext.c

    r8367 r8371  
    2222// Shared code snippets
    2323#include "util_ext.h"
     24#include "sw_domain.h"
    2425
    2526
    … …  
    10771078
    10781079
     1080/*
    10791081PyObject *gravity(PyObject *self, PyObject *args) {
    10801082  //
    … …  
    11131115    z2 = ((double*) z -> data)[k3 + 2];
    11141116
     1117    //printf("z0 %g, z1 %g, z2 %g \n",z0,z1,z2);
     1118
    11151119    // Optimise for flat bed
    11161120    // Note (Ole): This didn't produce measurable speed up.
    … …  
    11231127    avg_h = ((double *) h -> data)[k];
    11241128
     1129    //printf("avg_h  %g \n",avg_h);
     1130
    11251131    // Compute bed slope
    11261132    k6 = 6*k;  // base index
    … …  
    11331139    y2 = ((double*) x -> data)[k6 + 5];
    11341140
    1135 
     1141    //printf("x0 %g, y0 %g, x1 %g, y1 %g, x2 %g, y2 %g \n",x0,y0,x1,y1,x2,y2);
     1142   
    11361143    _gradient(x0, y0, x1, y1, x2, y2, z0, z1, z2, &zx, &zy);
    11371144
     1145    //printf("zx %g, zy %g \n",zx,zy);
    11381146    // Update momentum
    11391147    ((double*) xmom -> data)[k] += -g*zx*avg_h;
    … …  
    11431151  return Py_BuildValue("");
    11441152}
    1145 
    1146 
     1153*/
     1154
     1155//-------------------------------------------------------------------------------
     1156// gravity term using new get_python_domain interface
     1157//
     1158// FIXME SR: Probably should do this calculation together with hte compute
     1159// fluxes
     1160//-------------------------------------------------------------------------------
     1161PyObject *gravity(PyObject *self, PyObject *args) {
     1162  //
     1163  //  gravity(domain)
     1164  //
     1165
     1166   
     1167    struct domain D;
     1168    PyObject *domain;
     1169
     1170    int k, N, k3, k6;
     1171    double g, avg_h, zx, zy;
     1172    double x0, y0, x1, y1, x2, y2, z0, z1, z2;
     1173    //double h0,h1,h2;
     1174    //double epsilon;
     1175
     1176    if (!PyArg_ParseTuple(args, "O", &domain)) {
     1177        report_python_error(AT, "could not parse input arguments");
     1178        return NULL;
     1179    }
     1180
     1181    // populate the C domain structure with pointers
     1182    // to the python domain data
     1183    get_python_domain(&D,domain);
     1184
     1185
     1186    g = D.g;
     1187   
     1188    N = D.number_of_elements;
     1189    for (k=0; k<N; k++) {
     1190        k3 = 3*k;  // base index
     1191
     1192        // Get bathymetry
     1193        z0 = D.bed_vertex_values[k3 + 0];
     1194        z1 = D.bed_vertex_values[k3 + 1];
     1195        z2 = D.bed_vertex_values[k3 + 2];
     1196
     1197        //printf("z0 %g, z1 %g, z2 %g \n",z0,z1,z2);
     1198       
     1199        // Get average depth from centroid values
     1200        avg_h = D.stage_centroid_values[k] - D.bed_centroid_values[k];
     1201
     1202        //printf("avg_h  %g \n",avg_h);
     1203        // Compute bed slope
     1204        k6 = 6*k;  // base index
     1205
     1206        x0 = D.vertex_coordinates[k6 + 0];
     1207        y0 = D.vertex_coordinates[k6 + 1];
     1208        x1 = D.vertex_coordinates[k6 + 2];
     1209        y1 = D.vertex_coordinates[k6 + 3];
     1210        x2 = D.vertex_coordinates[k6 + 4];
     1211        y2 = D.vertex_coordinates[k6 + 5];
     1212
     1213        //printf("x0 %g, y0 %g, x1 %g, y1 %g, x2 %g, y2 %g \n",x0,y0,x1,y1,x2,y2);
     1214        _gradient(x0, y0, x1, y1, x2, y2, z0, z1, z2, &zx, &zy);
     1215
     1216        //printf("zx %g, zy %g \n",zx,zy);
     1217       
     1218        // Update momentum
     1219        D.xmom_explicit_update[k] += -g*zx*avg_h;
     1220        D.ymom_explicit_update[k] += -g*zy*avg_h;
     1221    }
     1222
     1223    return Py_BuildValue("");
     1224}
     1225
     1226
     1227//-------------------------------------------------------------------------------
     1228// New well balanced gravity term using new get_python_domain interface
     1229//
     1230// FIXME SR: Probably should do this calculation together with hte compute
     1231// fluxes
     1232//-------------------------------------------------------------------------------
    11471233PyObject *gravity_new(PyObject *self, PyObject *args) {
    11481234  //
    1149   //  gravity_new(g, h, v, x, xmom, ymom)
     1235  //  gravity_new(domain)
    11501236  //
    11511237
    11521238
    1153   PyArrayObject *h, *z, *x, *xmom, *ymom;
    1154   int k, N, k3, k6;
    1155   double g, avg_h, zx, zy;
    1156   double x0, y0, x1, y1, x2, y2, z0, z1, z2;
    1157   //double epsilon;
    1158 
    1159   if (!PyArg_ParseTuple(args, "dOOOOO",
    1160             &g, &h, &z, &x,
    1161             &xmom, &ymom)) {
    1162     //&epsilon)) {
    1163       report_python_error(AT, "could not parse input arguments");
    1164       return NULL;
    1165   }
    1166 
    1167   // check that numpy array objects arrays are C contiguous memory
    1168   CHECK_C_CONTIG(h);
    1169   CHECK_C_CONTIG(z);
    1170   CHECK_C_CONTIG(x);
    1171   CHECK_C_CONTIG(xmom);
    1172   CHECK_C_CONTIG(ymom);
    1173 
    1174   N = h -> dimensions[0];
    1175   for (k=0; k<N; k++) {
    1176     k3 = 3*k;  // base index
    1177 
    1178     // Get bathymetry
    1179     z0 = ((double*) z -> data)[k3 + 0];
    1180     z1 = ((double*) z -> data)[k3 + 1];
    1181     z2 = ((double*) z -> data)[k3 + 2];
    1182 
    1183     // Optimise for flat bed
    1184     // Note (Ole): This didn't produce measurable speed up.
    1185     // Revisit later
    1186     //if (fabs(z0-z1)<epsilon && fabs(z1-z2)<epsilon) {
    1187     //  continue;
    1188     //}
    1189 
    1190     // Get average depth from centroid values
    1191     avg_h = ((double *) h -> data)[k];
    1192 
    1193     // Compute bed slope
    1194     k6 = 6*k;  // base index
    1195 
    1196     x0 = ((double*) x -> data)[k6 + 0];
    1197     y0 = ((double*) x -> data)[k6 + 1];
    1198     x1 = ((double*) x -> data)[k6 + 2];
    1199     y1 = ((double*) x -> data)[k6 + 3];
    1200     x2 = ((double*) x -> data)[k6 + 4];
    1201     y2 = ((double*) x -> data)[k6 + 5];
    1202 
    1203 
    1204     _gradient(x0, y0, x1, y1, x2, y2, z0, z1, z2, &zx, &zy);
    1205 
    1206     // Update momentum
    1207     ((double*) xmom -> data)[k] += -g*zx*avg_h;
    1208     ((double*) ymom -> data)[k] += -g*zy*avg_h;
    1209   }
    1210 
    1211   return Py_BuildValue("");
    1212 }
    1213 
     1239    struct domain D;
     1240    PyObject *domain;
     1241
     1242    int i, k, N, k3, k6;
     1243    double g, avg_h, wx, wy, fact;
     1244    double x0, y0, x1, y1, x2, y2;
     1245    double n[2];
     1246    double h[3];
     1247    double w[3];
     1248    double side[2];
     1249    //double epsilon;
     1250
     1251    if (!PyArg_ParseTuple(args, "O", &domain)) {
     1252        report_python_error(AT, "could not parse input arguments");
     1253        return NULL;
     1254    }
     1255
     1256    // populate the C domain structure with pointers
     1257    // to the python domain data
     1258    get_python_domain(&D,domain);
     1259
     1260
     1261    g = D.g;
     1262
     1263    N = D.number_of_elements;
     1264    for (k=0; k<N; k++) {
     1265        k3 = 3*k;  // base index
     1266
     1267        //------------------------------------
     1268        // Calculate side terms -ghw_x term
     1269        //------------------------------------
     1270
     1271        // Get vertex bed values for gradient calculation
     1272        w[0] = D.stage_vertex_values[k3 + 0];
     1273        w[1] = D.stage_vertex_values[k3 + 1];
     1274        w[2] = D.stage_vertex_values[k3 + 2];
     1275
     1276        // Compute stage slope
     1277        k6 = 6*k;  // base index
     1278
     1279        x0 = D.vertex_coordinates[k6 + 0];
     1280        y0 = D.vertex_coordinates[k6 + 1];
     1281        x1 = D.vertex_coordinates[k6 + 2];
     1282        y1 = D.vertex_coordinates[k6 + 3];
     1283        x2 = D.vertex_coordinates[k6 + 4];
     1284        y2 = D.vertex_coordinates[k6 + 5];
     1285
     1286        //printf("x0 %g, y0 %g, x1 %g, y1 %g, x2 %g, y2 %g \n",x0,y0,x1,y1,x2,y2);
     1287        _gradient(x0, y0, x1, y1, x2, y2, w[0], w[1], w[2], &wx, &wy);
     1288
     1289        avg_h = D.stage_centroid_values[k] - D.bed_centroid_values[k];
     1290
     1291        // Update using -ghw_x term
     1292        D.xmom_explicit_update[k] += -g*wx*avg_h;
     1293        D.ymom_explicit_update[k] += -g*wy*avg_h;
     1294
     1295        //------------------------------------
     1296        // Calculate side terms \sum_i 0.5 g l_i h_i^2 n_i
     1297        //------------------------------------
     1298       
     1299        // Get edge depths
     1300        h[0] = D.stage_edge_values[k3 + 0] - D.bed_edge_values[k3 + 0];
     1301        h[1] = D.stage_edge_values[k3 + 1] - D.bed_edge_values[k3 + 1];
     1302        h[2] = D.stage_edge_values[k3 + 2] - D.bed_edge_values[k3 + 2];
     1303
     1304        // Calculate the side correction term
     1305        side[0] = 0.0;
     1306        side[1] = 0.0;
     1307        for (i=0; i<3; i++) {
     1308            n[0] = D.normals[k6+2*i];
     1309            n[1] = D.normals[k6+2*i+1];
     1310            fact = 0.5*g*D.edgelengths[k3+i]*h[i]*h[i];
     1311            side[0] += fact*n[0];
     1312            side[1] += fact*n[1];
     1313        }
     1314
     1315        // Update momentum with side terms
     1316        D.xmom_explicit_update[k] += side[0];
     1317        D.ymom_explicit_update[k] += side[1];
     1318
     1319    }
     1320
     1321    return Py_BuildValue("");
     1322}
    12141323
    12151324
    … …  
    30233132  {"compute_fluxes_ext_kinetic", compute_fluxes_ext_kinetic, METH_VARARGS, "Print out"},
    30243133  {"gravity", gravity, METH_VARARGS, "Print out"},
    3025   {"gravity_new", gravity, METH_VARARGS, "Print out"},
     3134  {"gravity_new", gravity_new, METH_VARARGS, "Print out"},
    30263135  {"manning_friction_flat", manning_friction_flat, METH_VARARGS, "Print out"},
    30273136  {"manning_friction_sloped", manning_friction_sloped, METH_VARARGS, "Print out"},

  • trunk/anuga_core/source/anuga/shallow_water/sw_domain.h

    r8369 r8371  
    1111// structure
    1212struct domain {
    13     double  timestep;
    1413    long    number_of_elements;
    1514    double  epsilon;
    … …  
    2524    long*   tri_full_flag;
    2625    long*   already_computed_flux;
    27     double* max_speed_array;
     26
     27    double* vertex_coordinates;
    2828
    2929    double* stage_edge_values;
    … …  
    3131    double* ymom_edge_values;
    3232    double* bed_edge_values;
     33
     34    double* stage_centroid_values;
     35    double* xmom_centroid_values;
     36    double* ymom_centroid_values;
     37    double* bed_centroid_values;
     38
     39    double* stage_vertex_values;
     40    double* xmom_vertex_values;
     41    double* ymom_vertex_values;
     42    double* bed_vertex_values;
     43
     44
    3345    double* stage_boundary_values;
    3446    double* xmom_boundary_values;
    3547    double* ymom_boundary_values;
     48    double* bed_boundary_values;
     49
    3650    double* stage_explicit_update;
    3751    double* xmom_explicit_update;
    … …  
    4155
    4256
    43 struct domain* get_python_domain(struct domain *D, PyObject *domain, double timestep) {
     57struct domain* get_python_domain(struct domain *D, PyObject *domain) {
    4458    PyArrayObject
    4559            *neighbours,
    … …  
    5064            *areas,
    5165            *tri_full_flag,
    52             *max_speed_array,
    5366            *already_computed_flux,
     67            *vertex_coordinates;
    5468
    55             *stage_edge_values,
    56             *xmom_edge_values,
    57             *ymom_edge_values,
    58             *bed_edge_values,
    59             *stage_boundary_values,
    60             *xmom_boundary_values,
    61             *ymom_boundary_values,
    62             *stage_explicit_update,
    63             *xmom_explicit_update,
    64             *ymom_explicit_update;
     69    PyObject *quantities;
    6570
    66 
    67 
    68     D->timestep = timestep;
    6971    D->number_of_elements = get_python_integer(domain, "number_of_elements");
    7072    D->epsilon = get_python_double(domain, "epsilon");
    7173    D->H0 = get_python_double(domain, "H0");
    7274    D->g = get_python_double(domain, "g");
    73 
    7475
    7576    neighbours = get_consecutive_array(domain, "neighbours");
    … …  
    9495    D->tri_full_flag = (long *) tri_full_flag->data;
    9596
    96     max_speed_array = get_consecutive_array(domain, "max_speed_array");
    97     D->max_speed_array = (double *) max_speed_array->data;
     97    already_computed_flux = get_consecutive_array(domain, "already_computed_flux");
     98    D->already_computed_flux = (long *) already_computed_flux->data;
    9899
    99     already_computed_flux = get_consecutive_array(domain, "already_computed_flux");
    100     D->already_computed_flux = (double *) already_computed_flux->data;
     100    vertex_coordinates = get_consecutive_array(domain, "vertex_coordinates");
     101    D->vertex_coordinates = (double *) vertex_coordinates->data;
    101102
    102     quantities = get_python_object(domain, "quantities")
     103    quantities = get_python_object(domain, "quantities");
    103104
    104105    D->stage_edge_values     = get_python_array_data_from_dict(quantities, "stage",     "edge_values");
    … …  
    106107    D->ymom_edge_values      = get_python_array_data_from_dict(quantities, "ymomentum", "edge_values");
    107108    D->bed_edge_values       = get_python_array_data_from_dict(quantities, "elevation", "edge_values");
     109
     110    D->stage_centroid_values     = get_python_array_data_from_dict(quantities, "stage",     "centroid_values");
     111    D->xmom_centroid_values      = get_python_array_data_from_dict(quantities, "xmomentum", "centroid_values");
     112    D->ymom_centroid_values      = get_python_array_data_from_dict(quantities, "ymomentum", "centroid_values");
     113    D->bed_centroid_values       = get_python_array_data_from_dict(quantities, "elevation", "centroid_values");
     114
     115    D->stage_vertex_values     = get_python_array_data_from_dict(quantities, "stage",     "vertex_values");
     116    D->xmom_vertex_values      = get_python_array_data_from_dict(quantities, "xmomentum", "vertex_values");
     117    D->ymom_vertex_values      = get_python_array_data_from_dict(quantities, "ymomentum", "vertex_values");
     118    D->bed_vertex_values       = get_python_array_data_from_dict(quantities, "elevation", "vertex_values");
     119
    108120    D->stage_boundary_values = get_python_array_data_from_dict(quantities, "stage",     "boundary_values");
    109     D->stage_boundary_values = get_python_array_data_from_dict(quantities, "xmomentum", "boundary_values");
    110     D->stage_boundary_values = get_python_array_data_from_dict(quantities, "ymomentum", "boundary_values");
     121    D->xmom_boundary_values  = get_python_array_data_from_dict(quantities, "xmomentum", "boundary_values");
     122    D->ymom_boundary_values  = get_python_array_data_from_dict(quantities, "ymomentum", "boundary_values");
     123    D->bed_boundary_values   = get_python_array_data_from_dict(quantities, "elevation", "boundary_values");
     124
    111125    D->stage_explicit_update = get_python_array_data_from_dict(quantities, "stage",     "explicit_update");
    112126    D->xmom_explicit_update  = get_python_array_data_from_dict(quantities, "xmomentum", "explicit_update");
    … …  
    114128
    115129
    116 
    117 
    118   Py_DECREF(neighbours);
    119   Py_DECREF(neighbour_vertices);
    120   Py_DECREF(normals);
    121   Py_DECREF(areas);
    122   Py_DECREF(max_speed_array);
    123 
    124   return D;
     130    return D;
    125131}

  • trunk/anuga_core/source/anuga/shallow_water/test_forcing.py

    r8126 r8371  
    19521952        domain.compute_forcing_terms()
    19531953
     1954
     1955
     1956       
    19541957        assert num.allclose(domain.quantities['stage'].explicit_update, 0)
    1955         assert num.allclose(domain.quantities['xmomentum'].explicit_update,
    1956                             -g*h*3)
     1958        assert num.allclose(domain.quantities['xmomentum'].explicit_update, -g*h*3)
    19571959        assert num.allclose(domain.quantities['ymomentum'].explicit_update, 0)
    19581960

  • trunk/anuga_core/source/anuga/utilities/util_ext.h

    r8370 r8371  
    1515
    1616#include <stdio.h>
     17
     18
     19#ifndef ANUGA_UTIL_EXT_H
     20#define ANUGA_UTIL_EXT_H
     21
     22
     23
     24
    1725#define STRINGIFY(x) #x
    1826#define TOSTRING(x) STRINGIFY(x)
    … …  
    401409                                    return NULL; \
    402410                                }
     411
     412
     413#endif
Note: See TracChangeset for help on using the changeset viewer.