Changeset 8385


Ignore:
Timestamp:
Apr 6, 2012, 11:30:59 AM (13 years ago)
Author:
davies
Message:

balanced_dev: Experimenting with flux_based boundary conditions

Location:
trunk/anuga_work/development/gareth
Files:
5 edited

Legend:

Unmodified
Added
Removed

  • trunk/anuga_work/development/gareth/experimental/balanced_dev/__init__.py

    r8353 r8385  
    99from balanced_dev.swb2_domain import Domain as Domain
    1010from balanced_dev.swb2_domain import Domain
     11from balanced_dev.swb2_boundary_conditions import zero_mass_flux_transmissive_momentum_flux_boundary as zero_mass_flux_transmissive_momentum_flux_boundary
    1112#, \
    1213#     Transmissive_boundary,\

  • trunk/anuga_work/development/gareth/experimental/balanced_dev/swb2_boundary_conditions.py

    r8353 r8385  
    22from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\
    33     import Boundary
     4
     5#from swb2_domain_ext import rotate
    46
    57class Transmissive_boundary(Boundary):
    … …  
    3234            q = self.domain.get_evolved_quantities(vol_id, edge = edge_id)
    3335        return q
     36
     37#####
     38
     39class  zero_mass_flux_transmissive_momentum_flux_boundary(Boundary):
     40    """ Boundary which operates directly on the fluxes
     41        CAN BE UNSTABLE -- CAREFUL
     42        Sets boundary values of h, uh, vh as in transmissive boundary, but
     43        also ensures:
     44        flux[0] = 0.0,
     45    """
     46    def __init__(self, domain = None):
     47        Boundary.__init__(self)
     48
     49        if domain is None:
     50            msg = 'Domain must be specified for zero_mass_flux_transmissive_momentum_flux boundary'
     51            raise Exception, msg
     52
     53        self.domain = domain
     54
     55    def __repr__(self):
     56        return 'zero_mass_flux_transmissive_momentum_flux_boundary(%s)' %self.domain
     57
     58    def evaluate(self, vol_id, edge_id):
     59        """Transmissive boundaries return the edge values
     60        of the volume they serve.
     61        """
     62
     63        if self.domain.get_centroid_transmissive_bc() :
     64            q = self.domain.get_evolved_quantities(vol_id)
     65        else:
     66            q = self.domain.get_evolved_quantities(vol_id, edge = edge_id)
     67        return q
     68
     69#####
     70
     71class  zero_mass_flux_zero_momentum_boundary(Boundary):
     72    """ Boundary which operates directly on the fluxes
     73        Sets boundary values of h=neighbour_h, uh=0, vh=0.
     74        and ensure that:
     75        flux[0] = 0.0
     76    """
     77    def __init__(self, domain = None):
     78        Boundary.__init__(self)
     79
     80        if domain is None:
     81            msg = 'Domain must be specified for zero_mass_flux_zero_momentum_boundary'
     82            raise Exception, msg
     83
     84        self.domain = domain
     85
     86    def __repr__(self):
     87        return 'zero_mass_flux_zero_momentum_boundary(%s)' %self.domain
     88
     89    def evaluate(self, vol_id, edge_id):
     90        """ Apply chosen values to q[1], q[2]
     91        """
     92
     93        if self.domain.get_centroid_transmissive_bc() :
     94            q = self.domain.get_evolved_quantities(vol_id)
     95        else:
     96            q = self.domain.get_evolved_quantities(vol_id, edge = edge_id)
     97
     98        normal = self.domain.normals[vol_id, 2*edge_id:2*edge_id+2]
     99
     100        r = rotate(q, normal, direction = 1)
     101        #r[1] = -r[1]
     102        tdamp=0.0
     103        ndamp=0.0
     104        r[1]=r[1]*ndamp
     105        r[2]=r[2]*tdamp
     106        q = rotate(r, normal, direction = -1)
     107
     108        # q[1]=0.
     109        # q[2]=0.
     110
     111        return q
     112
     113
     114class  zero_mass_flux_zero_n_transmissive_t_momentum_boundary(Boundary):
     115    """ Boundary which operates directly on the fluxes
     116        Sets boundary values of h=neighbour_h, uh=0, vh=vh
     117        and ensure that:
     118        flux[0] = 0.0
     119    """
     120    def __init__(self, domain = None):
     121        Boundary.__init__(self)
     122
     123        if domain is None:
     124            msg = 'Domain must be specified for zero_mass_flux_zero_n_transmissive_t_momentum_boundary'
     125            raise Exception, msg
     126
     127        self.domain = domain
     128
     129    def __repr__(self):
     130        return 'zero_mass_flux_zero_n_transmissive_t_momentum_boundary(%s)' %self.domain
     131
     132    def evaluate(self, vol_id, edge_id):
     133        """ Apply chosen values to q[1], q[2]
     134        """
     135
     136        if self.domain.get_centroid_transmissive_bc() :
     137            q = self.domain.get_evolved_quantities(vol_id)
     138        else:
     139            q = self.domain.get_evolved_quantities(vol_id, edge = edge_id)
     140
     141        normal = self.domain.normals[vol_id, 2*edge_id:2*edge_id+2]
     142
     143        r = rotate(q, normal, direction = 1)
     144        #r[1] = -r[1]
     145        tdamp=1.0
     146        ndamp=0.0
     147        r[1]=r[1]*ndamp
     148        r[2]=r[2]*tdamp
     149        q = rotate(r, normal, direction = -1)
     150
     151        # q[1]=0.
     152        # q[2]=0.
     153
     154        return q
     155###
     156def rotate(q,normal,direction=1):
     157    #Function for rotation -- better to mimic the shallow water one
     158    # FIXME: Use C version
     159    out=q*1.     
     160
     161    if(direction==1):
     162        out[1] = q[1]*normal[0] + q[2]*normal[1]   
     163        out[2] = -q[1]*normal[1] + q[2]*normal[0]
     164    elif(direction==-1):
     165        out[1] = q[1]*normal[0] - q[2]*normal[1]   
     166        out[2] = +q[1]*normal[1] + q[2]*normal[0]
     167    else:
     168        print 'ERROR IN rotate'
     169        assert 0==1
     170
     171    return out
     172   

  • trunk/anuga_work/development/gareth/experimental/balanced_dev/swb2_domain.py

    r8384 r8385  
    5858                            number_of_full_nodes = number_of_full_nodes,
    5959                            number_of_full_triangles = number_of_full_triangles)
    60         
     60       
    6161        #------------------------------------------------
    6262        # set some defaults
    … …  
    102102
    103103        self.maximum_allowed_speed=0.0
     104
    104105
    105106        print '##########################################################################'
    … …  
    242243                                           Xmom.boundary_values,
    243244                                           Ymom.boundary_values,
     245                                           domain.boundary_flux_type,
    244246                                           Stage.explicit_update,
    245247                                           Xmom.explicit_update,
    … …  
    519521                  int(self.optimise_dry_cells),
    520522                  int(self.extrapolate_velocity_second_order))
     523
     524    def set_boundary(self, boundary_map):
     525        ## Specialisation of set_boundary, which also updates the 'boundary_flux_type'
     526        Sww_domain.set_boundary(self,boundary_map)
     527       
     528        # Add a flag which can be used to distinguish flux boundaries within
     529        # compute_fluxes_central
     530        # Initialise to zero (which means 'not a flux_boundary')
     531        self.boundary_flux_type = self.boundary_edges*0
     532       
     533        # HACK to set the values of domain.boundary_flux
     534        for i in range(len(self.boundary_objects)):
     535            # Record the first 10 characters of the name of the boundary.
     536            # FIXME: There must be a better way than this!
     537            bndry_name = self.boundary_objects[i][1].__repr__()[0:10]
     538
     539            if(bndry_name=='zero_mass_'):
     540                # Create flag 'boundary_flux_type' that can be read in
     541                # compute_fluxes_central
     542                self.boundary_flux_type[i]=1
     543
     544

  • trunk/anuga_work/development/gareth/experimental/balanced_dev/swb2_domain_ext.c

    r8384 r8385  
    280280        double* xmom_boundary_values,
    281281        double* ymom_boundary_values,
     282        long*   boundary_flux_type,
    282283        double* stage_explicit_update,
    283284        double* xmom_explicit_update,
    … …  
    441442            edgeflux[1] *= length;
    442443            edgeflux[2] *= length;
     444
     445            if(n<0){
     446                if(boundary_flux_type[m]>0){
     447                    // IMPOSE BOUNDARY CONDITIONS DIRECTLY ON THE FLUX
     448                    if(boundary_flux_type[m]==1){
     449                        // Zero_mass_flux_transmissive_momentum_flux boundary, or
     450                        // zero_mass_flux_zero_momentum_boundary
     451                        // Just need to set the mass flux to zero, as the
     452                        // transmissive momentum is ensured by the values of
     453                        // qr[0], qr[1], qr[2]
     454                        edgeflux[0] = 0.0;
     455                    }
     456                }
     457            }
    443458
    444459            // Update triangle k with flux from edge i
    … …  
    14851500
    14861501
     1502//PyObject *rotate(PyObject *self, PyObject *args, PyObject *kwargs) {
     1503//  //
     1504//  // r = rotate(q, normal, direction=1)
     1505//  //
     1506//  // Where q is assumed to be a Float numeric array of length 3 and
     1507//  // normal a Float numeric array of length 2.
     1508//
     1509//  // FIXME(Ole): I don't think this is used anymore
     1510//
     1511//  PyObject *Q, *Normal;
     1512//  PyArrayObject *q, *r, *normal;
     1513//
     1514//  static char *argnames[] = {"q", "normal", "direction", NULL};
     1515//  int dimensions[1], i, direction=1;
     1516//  double n1, n2;
     1517//
     1518//  // Convert Python arguments to C
     1519//  if (!PyArg_ParseTupleAndKeywords(args, kwargs, "OO|i", argnames,
     1520//                   &Q, &Normal, &direction)) {
     1521//    report_python_error(AT, "could not parse input arguments");
     1522//    return NULL;
     1523//  } 
     1524//
     1525//  // Input checks (convert sequences into numeric arrays)
     1526//  q = (PyArrayObject *)
     1527//    PyArray_ContiguousFromObject(Q, PyArray_DOUBLE, 0, 0);
     1528//  normal = (PyArrayObject *)
     1529//    PyArray_ContiguousFromObject(Normal, PyArray_DOUBLE, 0, 0);
     1530//
     1531//
     1532//  if (normal -> dimensions[0] != 2) {
     1533//    report_python_error(AT, "Normal vector must have 2 components");
     1534//    return NULL;
     1535//  }
     1536//
     1537//  // Allocate space for return vector r (don't DECREF)
     1538//  dimensions[0] = 3;
     1539//  r = (PyArrayObject *) anuga_FromDims(1, dimensions, PyArray_DOUBLE);
     1540//
     1541//  // Copy
     1542//  for (i=0; i<3; i++) {
     1543//    ((double *) (r -> data))[i] = ((double *) (q -> data))[i];
     1544//  }
     1545//
     1546//  // Get normal and direction
     1547//  n1 = ((double *) normal -> data)[0];
     1548//  n2 = ((double *) normal -> data)[1];
     1549//  if (direction == -1) n2 = -n2;
     1550//
     1551//  // Rotate
     1552//  _rotate((double *) r -> data, n1, n2);
     1553//
     1554//  // Release numeric arrays
     1555//  Py_DECREF(q);
     1556//  Py_DECREF(normal);
     1557//
     1558//  // Return result using PyArray to avoid memory leak
     1559//  return PyArray_Return(r);
     1560//}
     1561
    14871562//========================================================================
    14881563// Compute fluxes
    … …  
    15521627    *xmom_boundary_values,
    15531628    *ymom_boundary_values,
     1629    *boundary_flux_type,
    15541630    *stage_explicit_update,
    15551631    *xmom_explicit_update,
    … …  
    15651641   
    15661642  // Convert Python arguments to C
    1567   if (!PyArg_ParseTuple(args, "ddddOOOOOOOOOOOOOOOOOOOiOOO",
     1643  if (!PyArg_ParseTuple(args, "ddddOOOOOOOOOOOOOOOOOOOOiOOO",
    15681644            &timestep,
    15691645            &epsilon,
    … …  
    15821658            &xmom_boundary_values,
    15831659            &ymom_boundary_values,
     1660            &boundary_flux_type,
    15841661            &stage_explicit_update,
    15851662            &xmom_explicit_update,
    … …  
    16101687  CHECK_C_CONTIG(xmom_boundary_values);
    16111688  CHECK_C_CONTIG(ymom_boundary_values);
     1689  CHECK_C_CONTIG(boundary_flux_type);
    16121690  CHECK_C_CONTIG(stage_explicit_update);
    16131691  CHECK_C_CONTIG(xmom_explicit_update);
    … …  
    16421720                     (double*) xmom_boundary_values -> data,
    16431721                     (double*) ymom_boundary_values -> data,
     1722                     (long*)   boundary_flux_type -> data,
    16441723                     (double*) stage_explicit_update -> data,
    16451724                     (double*) xmom_explicit_update -> data,
    … …  
    19862065   * three.
    19872066   */
     2067  //{"rotate", (PyCFunction)rotate, METH_VARARGS | METH_KEYWORDS, "Print out"},
    19882068  {"compute_fluxes_ext_central", compute_fluxes_ext_central, METH_VARARGS, "Print out"},
    19892069  {"gravity_c",        gravity,            METH_VARARGS, "Print out"},

  • trunk/anuga_work/development/gareth/tests/util.py

    r8375 r8385  
    1 """ Utilities for reading data / plotting of channel/floodplain case
     1""" Random utilities for reading sww file data and for plotting (in ipython, or
     2    in scripts)
     3
     4    Functionality of note:
     5
     6    util.get_outputs -- read the data from a single sww file into a single object
     7   
     8    util.combine_outputs -- read the data from a list of sww files into a single object
     9   
     10    util.near_transect -- for finding the indices of points 'near' to a given
     11                       line, and assigning these points a coordinate along that line.
     12                       This is useful for plotting outputs which are 'almost' along a
     13                       transect (e.g. a channel cross-section) -- see example below
     14
     15
     16
     17    Here is an example ipython session which uses some of these functions:
     18
     19    > import util
     20    > from matplotlib import pyplot as pyplot
     21    > p=util.get_output('myfile.sww',minimum_allowed_height=0.01)
     22    > p2=util.get_centroids(p,velocity_extrapolation=True)
     23    > xxx=util.near_transect(p,[95., 85.], [120.,68.],tol=2.) # Could equally well use p2
     24    > pyplot.ion() # Interactive plotting
     25    > pyplot.scatter(xxx[1],p.vel[140,xxx[0]],color='red') # Plot along the transect
     26
    227"""
     28from Scientific.IO.NetCDF import NetCDFFile
     29import numpy
     30
     31class combine_outputs:
     32    """
     33    Read in a list of filenames, and combine all their outputs into a single object.
     34    e.g.:
     35
     36    p = util.combine_outputs(['file1.sww', 'file1_time_10000.sww', 'file1_time_20000.sww'], 0.01)
     37   
     38    will make an object p which has components p.x,p.y,p.time,p.stage, .... etc,
     39    where the values of stage / momentum / velocity from the sww files are concatenated as appropriate.
     40
     41    This is nice for interactive interrogation of model outputs, or for sticking together outputs in scripts
     42   
     43    WARNING: It is easy to use lots of memory, if the sww files are large.
     44
     45    Note: If you want the centroid values, then you could subsequently use:
     46
     47    p2 = util.get_centroids(p,velocity_extrapolation=False)
     48
     49    which would make an object p2 that is like p, but holds information at centroids
     50    """
     51    def __init__(self, filename_list, minimum_allowed_height=1.0e-03):
     52        #
     53        # Go through the sww files in 'filename_list', and combine them into one object.
     54        #
     55
     56        for i, filename in enumerate(filename_list):
     57            print i, filename
     58            # Store output from filename
     59            p_tmp = get_output(filename, minimum_allowed_height)
     60            if(i==0):
     61                # Create self
     62                p1=p_tmp
     63            else:
     64                # Append extra data to self
     65                # Note that p1.x, p1.y, p1.vols, p1.elev should not change
     66                assert (p1.x == p_tmp.x).all()
     67                assert (p1.y == p_tmp.y).all()
     68                assert (p1.vols ==p_tmp.vols).all()
     69                p1.time = numpy.append(p1.time, p_tmp.time)
     70                p1.stage = numpy.append(p1.stage, p_tmp.stage, axis=0)
     71                p1.xmom = numpy.append(p1.xmom, p_tmp.xmom, axis=0)
     72                p1.ymom = numpy.append(p1.ymom, p_tmp.ymom, axis=0)
     73                p1.xvel = numpy.append(p1.xvel, p_tmp.xvel, axis=0)
     74                p1.yvel = numpy.append(p1.yvel, p_tmp.yvel, axis=0)
     75                p1.vel = numpy.append(p1.vel, p_tmp.vel, axis=0)
     76
     77        self.x, self.y, self.time, self.vols, self.elev, self.stage, self.xmom, \
     78                self.ymom, self.xvel, self.yvel, self.vel, self.minimum_allowed_height = \
     79                p1.x, p1.y, p1.time, p1.vols, p1.elev, p1.stage, p1.xmom, p1.ymom, p1.xvel,\
     80                p1.yvel, p1.vel, p1.minimum_allowed_height
     81
     82
    383
    484class get_output:
    … …  
    27107    # Import modules
    28108
    29     from Scientific.IO.NetCDF import NetCDFFile
    30109
    31110
    … …  
    92171    # Input: p is the result of e.g. p=util.get_output('mysww.sww'). See the get_output class defined above
    93172    # Output: Values of x, y, Stage, xmom, ymom, elev, xvel, yvel, vel at centroids
    94     import numpy
     173    #import numpy
    95174   
    96175    # Make 3 arrays, each containing one index of a vertex of every triangle.
    … …  
    207286    # e.g.
    208287    # import util
    209     # #import transect_interpolate
    210288    # from matplotlib import pyplot
    211289    # p=util.get_output('merewether_1m.sww',0.01)
Note: See TracChangeset for help on using the changeset viewer.