Changeset 3703


Ignore:
Timestamp:
Oct 5, 2006, 5:50:11 PM (17 years ago)
Author:
steve
Message:
 
Location:
anuga_core/source/anuga
Files:
5 edited

Legend:

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

    r3689 r3703  
    751751
    752752            #Yield results
    753             if finaltime is not None and self.time >= finaltime:
    754 
    755                 if self.time > finaltime:
    756                     #FIXME (Ole, 30 April 2006): Do we need this check?
    757                     print 'WARNING (domain.py): time overshot finaltime. Contact Ole.Nielsen@ga.gov.au'
    758                     self.time = finaltime
     753            if finaltime is not None and self.time >= finaltime:
     754
     755                if self.time > finaltime:
     756                              #FIXME (Ole, 30 April 2006): Do we need this check?
     757                      print 'WARNING (domain.py): time overshot finaltime. Contact Ole.Nielsen@ga.gov.au'
     758                      self.time = finaltime
    759759
    760760                # Yield final time and stop
     
    763763
    764764
    765             if self.yieldtime >= yieldstep:
     765            if self.yieldtime >= yieldstep:
    766766                # Yield (intermediate) time and allow inspection of domain
    767767
  • anuga_core/source/anuga/shallow_water/shallow_water_domain.py

    r3689 r3703  
    144144        self.beta_h      = beta_h
    145145
    146 
     146        self.flux_function = flux_function_central
     147        #self.flux_function = flux_function_kinetic
     148       
    147149        self.forcing_terms.append(manning_friction)
    148150        self.forcing_terms.append(gravity)
     
    495497####################################
    496498# Flux computation
    497 def flux_function(normal, ql, qr, zl, zr):
     499def flux_function_central(normal, ql, qr, zl, zr):
    498500    """Compute fluxes between volumes for the shallow water wave equation
    499501    cast in terms of w = h+z using the 'central scheme' as described in
     
    577579
    578580    return edgeflux, max_speed
     581
     582def erfcc(x):
     583
     584    from math import fabs, exp
     585
     586    z=fabs(x)
     587    t=1.0/(1.0+0.5*z)
     588    result=t*exp(-z*z-1.26551223+t*(1.00002368+t*(.37409196+
     589         t*(.09678418+t*(-.18628806+t*(.27886807+t*(-1.13520398+
     590         t*(1.48851587+t*(-.82215223+t*.17087277)))))))))
     591    if x < 0.0:
     592        result = 2.0-result
     593
     594    return result
     595
     596def flux_function_kinetic(normal, ql, qr, zl, zr):
     597    """Compute fluxes between volumes for the shallow water wave equation
     598    cast in terms of w = h+z using the 'central scheme' as described in
     599
     600    Zhang et. al., Advances in Water Resources, 26(6), 2003, 635-647.
     601
     602
     603    Conserved quantities w, uh, vh are stored as elements 0, 1 and 2
     604    in the numerical vectors ql an qr.
     605
     606    Bed elevations zl and zr.
     607    """
     608
     609    from anuga.config import g, epsilon
     610    from math import sqrt
     611    from Numeric import array
     612
     613    #Align momentums with x-axis
     614    q_left  = rotate(ql, normal, direction = 1)
     615    q_right = rotate(qr, normal, direction = 1)
     616
     617    z = (zl+zr)/2 #Take average of field values
     618
     619    w_left  = q_left[0]   #w=h+z
     620    h_left  = w_left-z
     621    uh_left = q_left[1]
     622
     623    if h_left < epsilon:
     624        u_left = 0.0  #Could have been negative
     625        h_left = 0.0
     626    else:
     627        u_left  = uh_left/h_left
     628
     629
     630    w_right  = q_right[0]  #w=h+z
     631    h_right  = w_right-z
     632    uh_right = q_right[1]
     633
     634
     635    if h_right < epsilon:
     636        u_right = 0.0  #Could have been negative
     637        h_right = 0.0
     638    else:
     639        u_right  = uh_right/h_right
     640
     641    vh_left  = q_left[2]
     642    vh_right = q_right[2]
     643
     644    soundspeed_left  = sqrt(g*h_left)
     645    soundspeed_right = sqrt(g*h_right)
     646
     647    #Maximal wave speed
     648    s_max = max(u_left+soundspeed_left, u_right+soundspeed_right, 0)
     649
     650    #Minimal wave speed
     651    s_min = min(u_left-soundspeed_left, u_right-soundspeed_right, 0)
     652
     653    #Flux computation
     654
     655    F_left  = 0.0
     656    F_right = 0.0
     657    from math import sqrt, pi, exp
     658    if h_left > 0.0:
     659        F_left = u_left/sqrt(g*h_left)
     660    if h_right > 0.0:
     661        F_right = u_right/sqrt(g*h_right)
     662
     663    edgeflux = array([0.0, 0.0, 0.0])
     664
     665    edgeflux[0] = h_left*u_left/2.0*erfcc(-F_left) +  \
     666          h_left*sqrt(g*h_left)/2.0/sqrt(pi)*exp(-(F_left**2)) + \
     667          h_right*u_right/2.0*erfcc(F_right) -  \
     668          h_right*sqrt(g*h_right)/2.0/sqrt(pi)*exp(-(F_right**2))
     669
     670    edgeflux[1] = (h_left*u_left**2 + g/2.0*h_left**2)/2.0*erfcc(-F_left) + \
     671          u_left*h_left*sqrt(g*h_left)/2.0/sqrt(pi)*exp(-(F_left**2)) + \
     672          (h_right*u_right**2 + g/2.0*h_right**2)/2.0*erfcc(F_right) -  \
     673          u_right*h_right*sqrt(g*h_right)/2.0/sqrt(pi)*exp(-(F_right**2))
     674
     675    edgeflux[2] = vh_left*u_left/2.0*erfcc(-F_left) + \
     676          vh_left*sqrt(g*h_left)/2.0/sqrt(pi)*exp(-(F_left**2)) + \
     677          vh_right*u_right/2.0*erfcc(F_right) - \
     678          vh_right*sqrt(g*h_right)/2.0/sqrt(pi)*exp(-(F_right**2))
     679
     680
     681    edgeflux = rotate(edgeflux, normal, direction=-1)
     682    max_speed = max(abs(s_max), abs(s_min))
     683
     684    return edgeflux, max_speed
     685
    579686
    580687
     
    648755
    649756            #Flux computation using provided function
    650             edgeflux, max_speed = flux_function(normal, ql, qr, zl, zr)
     757            edgeflux, max_speed = domain.flux_function(normal, ql, qr, zl, zr)
    651758            flux -= edgeflux * domain.edgelengths[k,i]
    652759
     
    712819
    713820    timestep = float(sys.maxint)
    714     from shallow_water_ext import compute_fluxes
    715     domain.timestep = compute_fluxes(timestep, domain.epsilon, domain.g,
     821    from shallow_water_ext import compute_fluxes_ext_central as compute_fluxes_ext
     822    domain.timestep = compute_fluxes_ext(timestep, domain.epsilon, domain.g,
    716823                                     domain.neighbours,
    717824                                     domain.neighbour_edges,
  • anuga_core/source/anuga/shallow_water/shallow_water_ext.c

    r3696 r3703  
    126126
    127127// Computational function for flux computation (using stage w=z+h)
    128 int flux_function(double *q_left, double *q_right,
     128int flux_function_central(double *q_left, double *q_right,
    129129                  double z_left, double z_right,
    130130                  double n1, double n2,
     
    232232  return 0;
    233233}
     234
     235double erfcc(double x){
     236    double z,t,result;
     237
     238    z=fabs(x);
     239    t=1.0/(1.0+0.5*z);
     240    result=t*exp(-z*z-1.26551223+t*(1.00002368+t*(.37409196+
     241         t*(.09678418+t*(-.18628806+t*(.27886807+t*(-1.13520398+
     242         t*(1.48851587+t*(-.82215223+t*.17087277)))))))));
     243    if (x < 0.0) result = 2.0-result;
     244
     245    return result;
     246    }
     247
     248
     249
     250// Computational function for flux computation (using stage w=z+h)
     251int flux_function_kinetic(double *q_left, double *q_right,
     252                  double z_left, double z_right,
     253                  double n1, double n2,
     254                  double epsilon, double g,
     255                  double *edgeflux, double *max_speed) {
     256
     257  /*Compute fluxes between volumes for the shallow water wave equation
     258    cast in terms of the 'stage', w = h+z using
     259    the 'central scheme' as described in
     260
     261    Zhang et. al., Advances in Water Resources, 26(6), 2003, 635-647.
     262  */
     263
     264  int i;
     265
     266  double w_left, h_left, uh_left, vh_left, u_left, F_left;
     267  double w_right, h_right, uh_right, vh_right, u_right, F_right;
     268  double s_min, s_max, soundspeed_left, soundspeed_right;
     269  double z;
     270  double q_left_copy[3], q_right_copy[3];
     271
     272
     273  //Copy conserved quantities to protect from modification
     274  for (i=0; i<3; i++) {
     275    q_left_copy[i] = q_left[i];
     276    q_right_copy[i] = q_right[i];
     277  }
     278
     279  //Align x- and y-momentum with x-axis
     280  _rotate(q_left_copy, n1, n2);
     281  _rotate(q_right_copy, n1, n2);
     282
     283  z = (z_left+z_right)/2; //Take average of field values
     284
     285  //Compute speeds in x-direction
     286  w_left = q_left_copy[0];              // h+z
     287  h_left = w_left-z;
     288  uh_left = q_left_copy[1];
     289
     290  if (h_left < epsilon) {
     291    h_left = 0.0;  //Could have been negative
     292    u_left = 0.0;
     293  } else {
     294    u_left = uh_left/h_left;
     295  }
     296
     297  w_right = q_right_copy[0];
     298  h_right = w_right-z;
     299  uh_right = q_right_copy[1];
     300
     301  if (h_right < epsilon) {
     302    h_right = 0.0; //Could have been negative
     303    u_right = 0.0;
     304  } else {
     305    u_right = uh_right/h_right;
     306  }
     307
     308  //Momentum in y-direction
     309  vh_left  = q_left_copy[2];
     310  vh_right = q_right_copy[2];
     311
     312
     313  //Maximal and minimal wave speeds
     314  soundspeed_left  = sqrt(g*h_left);
     315  soundspeed_right = sqrt(g*h_right);
     316
     317  s_max = max(u_left+soundspeed_left, u_right+soundspeed_right);
     318  if (s_max < 0.0) s_max = 0.0;
     319
     320  s_min = min(u_left-soundspeed_left, u_right-soundspeed_right);
     321  if (s_min > 0.0) s_min = 0.0;
     322
     323
     324  F_left  = 0.0;
     325  F_right = 0.0;
     326  if (h_left > 0.0) F_left = u_left/sqrt(g*h_left);
     327  if (h_right > 0.0) F_right = u_right/sqrt(g*h_right);
     328
     329  for (i=0; i<3; i++) edgeflux[i] = 0.0;
     330  *max_speed = 0.0;
     331
     332  edgeflux[0] = h_left*u_left/2.0*erfcc(-F_left) +  \
     333          h_left*sqrt(g*h_left)/2.0/sqrt(pi)*exp(-(F_left*F_left)) + \
     334          h_right*u_right/2.0*erfcc(F_right) -  \
     335          h_right*sqrt(g*h_right)/2.0/sqrt(pi)*exp(-(F_right*F_right));
     336
     337  edgeflux[1] = (h_left*u_left*u_left + g/2.0*h_left*h_left)/2.0*erfcc(-F_left) + \
     338          u_left*h_left*sqrt(g*h_left)/2.0/sqrt(pi)*exp(-(F_left*F_left)) + \
     339          (h_right*u_right*u_right + g/2.0*h_right*h_right)/2.0*erfcc(F_right) -  \
     340          u_right*h_right*sqrt(g*h_right)/2.0/sqrt(pi)*exp(-(F_right*F_right));
     341
     342  edgeflux[2] = vh_left*u_left/2.0*erfcc(-F_left) + \
     343          vh_left*sqrt(g*h_left)/2.0/sqrt(pi)*exp(-(F_left*F_left)) + \
     344          vh_right*u_right/2.0*erfcc(F_right) - \
     345          vh_right*sqrt(g*h_right)/2.0/sqrt(pi)*exp(-(F_right*F_right));
     346
     347  //Maximal wavespeed
     348  *max_speed = max(fabs(s_max), fabs(s_min));
     349
     350  //Rotate back
     351  _rotate(edgeflux, n1, -n2);
     352
     353  return 0;
     354}
     355
     356
     357
    234358
    235359void _manning_friction(double g, double eps, int N,
     
    9961120}
    9971121
    998 PyObject *compute_fluxes(PyObject *self, PyObject *args) {
     1122PyObject *compute_fluxes_ext_central(PyObject *self, PyObject *args) {
    9991123  /*Compute all fluxes and the timestep suitable for all volumes
    10001124    in domain.
    10011125
    1002     Compute total flux for each conserved quantity using "flux_function"
     1126    Compute total flux for each conserved quantity using "flux_function_central"
    10031127
    10041128    Fluxes across each edge are scaled by edgelengths and summed up
     
    11341258      normal[1] = ((double *) normals -> data)[ki2+1];
    11351259      //Edge flux computation
    1136       flux_function(ql, qr, zl, zr,
     1260      flux_function_central(ql, qr, zl, zr,
    11371261                    normal[0], normal[1],
    11381262                    epsilon, g,
     
    11741298}
    11751299
     1300
     1301PyObject *compute_fluxes_ext_kinetic(PyObject *self, PyObject *args) {
     1302  /*Compute all fluxes and the timestep suitable for all volumes
     1303    in domain.
     1304
     1305    Compute total flux for each conserved quantity using "flux_function_central"
     1306
     1307    Fluxes across each edge are scaled by edgelengths and summed up
     1308    Resulting flux is then scaled by area and stored in
     1309    explicit_update for each of the three conserved quantities
     1310    stage, xmomentum and ymomentum
     1311
     1312    The maximal allowable speed computed by the flux_function for each volume
     1313    is converted to a timestep that must not be exceeded. The minimum of
     1314    those is computed as the next overall timestep.
     1315
     1316    Python call:
     1317    domain.timestep = compute_fluxes(timestep,
     1318                                     domain.epsilon,
     1319                                     domain.g,
     1320                                     domain.neighbours,
     1321                                     domain.neighbour_edges,
     1322                                     domain.normals,
     1323                                     domain.edgelengths,
     1324                                     domain.radii,
     1325                                     domain.areas,
     1326                                     Stage.edge_values,
     1327                                     Xmom.edge_values,
     1328                                     Ymom.edge_values,
     1329                                     Bed.edge_values,
     1330                                     Stage.boundary_values,
     1331                                     Xmom.boundary_values,
     1332                                     Ymom.boundary_values,
     1333                                     Stage.explicit_update,
     1334                                     Xmom.explicit_update,
     1335                                     Ymom.explicit_update,
     1336                                     already_computed_flux)
     1337
     1338
     1339    Post conditions:
     1340      domain.explicit_update is reset to computed flux values
     1341      domain.timestep is set to the largest step satisfying all volumes.
     1342
     1343
     1344  */
     1345
     1346
     1347  PyArrayObject *neighbours, *neighbour_edges,
     1348    *normals, *edgelengths, *radii, *areas,
     1349    *tri_full_flag,
     1350    *stage_edge_values,
     1351    *xmom_edge_values,
     1352    *ymom_edge_values,
     1353    *bed_edge_values,
     1354    *stage_boundary_values,
     1355    *xmom_boundary_values,
     1356    *ymom_boundary_values,
     1357    *stage_explicit_update,
     1358    *xmom_explicit_update,
     1359    *ymom_explicit_update,
     1360    *already_computed_flux;//tracks whether the flux across an edge has already been computed
     1361
     1362
     1363  //Local variables
     1364  double timestep, max_speed, epsilon, g;
     1365  double normal[2], ql[3], qr[3], zl, zr;
     1366  double edgeflux[3]; //Work arrays for summing up fluxes
     1367
     1368  int number_of_elements, k, i, m, n;
     1369  int ki, nm=0, ki2; //Index shorthands
     1370  static long call=1;
     1371
     1372
     1373  // Convert Python arguments to C
     1374  if (!PyArg_ParseTuple(args, "dddOOOOOOOOOOOOOOOOOO",
     1375                        &timestep,
     1376                        &epsilon,
     1377                        &g,
     1378                        &neighbours,
     1379                        &neighbour_edges,
     1380                        &normals,
     1381                        &edgelengths, &radii, &areas,
     1382                        &tri_full_flag,
     1383                        &stage_edge_values,
     1384                        &xmom_edge_values,
     1385                        &ymom_edge_values,
     1386                        &bed_edge_values,
     1387                        &stage_boundary_values,
     1388                        &xmom_boundary_values,
     1389                        &ymom_boundary_values,
     1390                        &stage_explicit_update,
     1391                        &xmom_explicit_update,
     1392                        &ymom_explicit_update,
     1393                        &already_computed_flux)) {
     1394    PyErr_SetString(PyExc_RuntimeError, "Input arguments failed");
     1395    return NULL;
     1396  }
     1397  number_of_elements = stage_edge_values -> dimensions[0];
     1398  call++;//a static local variable to which already_computed_flux is compared
     1399  //set explicit_update to zero for all conserved_quantities.
     1400  //This assumes compute_fluxes called before forcing terms
     1401  for (k=0; k<number_of_elements; k++) {
     1402    ((double *) stage_explicit_update -> data)[k]=0.0;
     1403    ((double *) xmom_explicit_update -> data)[k]=0.0;
     1404    ((double *) ymom_explicit_update -> data)[k]=0.0;
     1405  }
     1406  //Loop through neighbours and compute edge flux for each
     1407  for (k=0; k<number_of_elements; k++) {
     1408    for (i=0; i<3; i++) {
     1409      ki = k*3+i;
     1410      if (((long *) already_computed_flux->data)[ki]==call)//we've already computed the flux across this edge
     1411        continue;
     1412      ql[0] = ((double *) stage_edge_values -> data)[ki];
     1413      ql[1] = ((double *) xmom_edge_values -> data)[ki];
     1414      ql[2] = ((double *) ymom_edge_values -> data)[ki];
     1415      zl =    ((double *) bed_edge_values -> data)[ki];
     1416
     1417      //Quantities at neighbour on nearest face
     1418      n = ((long *) neighbours -> data)[ki];
     1419      if (n < 0) {
     1420        m = -n-1; //Convert negative flag to index
     1421        qr[0] = ((double *) stage_boundary_values -> data)[m];
     1422        qr[1] = ((double *) xmom_boundary_values -> data)[m];
     1423        qr[2] = ((double *) ymom_boundary_values -> data)[m];
     1424        zr = zl; //Extend bed elevation to boundary
     1425      } else {
     1426        m = ((long *) neighbour_edges -> data)[ki];
     1427        nm = n*3+m;
     1428        qr[0] = ((double *) stage_edge_values -> data)[nm];
     1429        qr[1] = ((double *) xmom_edge_values -> data)[nm];
     1430        qr[2] = ((double *) ymom_edge_values -> data)[nm];
     1431        zr =    ((double *) bed_edge_values -> data)[nm];
     1432      }
     1433      // Outward pointing normal vector
     1434      // normal = domain.normals[k, 2*i:2*i+2]
     1435      ki2 = 2*ki; //k*6 + i*2
     1436      normal[0] = ((double *) normals -> data)[ki2];
     1437      normal[1] = ((double *) normals -> data)[ki2+1];
     1438      //Edge flux computation
     1439      flux_function_kinetic(ql, qr, zl, zr,
     1440                    normal[0], normal[1],
     1441                    epsilon, g,
     1442                    edgeflux, &max_speed);
     1443      //update triangle k
     1444      ((long *) already_computed_flux->data)[ki]=call;
     1445      ((double *) stage_explicit_update -> data)[k] -= edgeflux[0]*((double *) edgelengths -> data)[ki];
     1446      ((double *) xmom_explicit_update -> data)[k] -= edgeflux[1]*((double *) edgelengths -> data)[ki];
     1447      ((double *) ymom_explicit_update -> data)[k] -= edgeflux[2]*((double *) edgelengths -> data)[ki];
     1448      //update the neighbour n
     1449      if (n>=0){
     1450        ((long *) already_computed_flux->data)[nm]=call;
     1451        ((double *) stage_explicit_update -> data)[n] += edgeflux[0]*((double *) edgelengths -> data)[nm];
     1452        ((double *) xmom_explicit_update -> data)[n] += edgeflux[1]*((double *) edgelengths -> data)[nm];
     1453        ((double *) ymom_explicit_update -> data)[n] += edgeflux[2]*((double *) edgelengths -> data)[nm];
     1454      }
     1455      ///for (j=0; j<3; j++) {
     1456        ///flux[j] -= edgeflux[j]*((double *) edgelengths -> data)[ki];
     1457        ///}
     1458        //Update timestep
     1459        //timestep = min(timestep, domain.radii[k]/max_speed)
     1460        //FIXME: SR Add parameter for CFL condition
     1461    if ( ((long *) tri_full_flag -> data)[k] == 1) {
     1462            if (max_speed > epsilon) {
     1463                timestep = min(timestep, ((double *) radii -> data)[k]/max_speed);
     1464                //maxspeed in flux_function is calculated as max(|u+a|,|u-a|)
     1465                if (n>=0)
     1466                    timestep = min(timestep, ((double *) radii -> data)[n]/max_speed);
     1467            }
     1468    }
     1469    } // end for i
     1470    //Normalise by area and store for when all conserved
     1471    //quantities get updated
     1472    ((double *) stage_explicit_update -> data)[k] /= ((double *) areas -> data)[k];
     1473    ((double *) xmom_explicit_update -> data)[k] /= ((double *) areas -> data)[k];
     1474    ((double *) ymom_explicit_update -> data)[k] /= ((double *) areas -> data)[k];
     1475  } //end for k
     1476  return Py_BuildValue("d", timestep);
     1477}
     1478
    11761479PyObject *protect(PyObject *self, PyObject *args) {
    11771480  //
     
    14491752  {"rotate", (PyCFunction)rotate, METH_VARARGS | METH_KEYWORDS, "Print out"},
    14501753  {"extrapolate_second_order_sw", extrapolate_second_order_sw, METH_VARARGS, "Print out"},
    1451   {"compute_fluxes", compute_fluxes, METH_VARARGS, "Print out"},
     1754  {"compute_fluxes_ext_central", compute_fluxes_ext_central, METH_VARARGS, "Print out"},
     1755  {"compute_fluxes_ext_kinetic", compute_fluxes_ext_kinetic, METH_VARARGS, "Print out"},
    14521756  {"gravity", gravity, METH_VARARGS, "Print out"},
    14531757  {"manning_friction", manning_friction, METH_VARARGS, "Print out"},
  • anuga_core/source/anuga/shallow_water/shallow_water_kinetic_ext.c

    r3698 r3703  
    11031103      normal[1] = ((double *) normals -> data)[ki2+1];
    11041104      //Edge flux computation
    1105       flux_function_kinetic(ql, qr, zl, zr,
     1105      flux_function_central(ql, qr, zl, zr,
    11061106                    normal[0], normal[1],
    11071107                    epsilon, g,
  • anuga_core/source/anuga/shallow_water/test_shallow_water_domain.py

    r3689 r3703  
    99
    1010from shallow_water_domain import *
    11 
     11from shallow_water_domain import flux_function_central as flux_function
    1212
    1313#Variable windfield implemented using functions
Note: See TracChangeset for help on using the changeset viewer.