Changeset 8477


Ignore:
Timestamp:
Jul 24, 2012, 4:02:41 PM (13 years ago)
Author:
steve
Message:

Adding in rate_operator

Location:
trunk/anuga_core/source/anuga/operators
Files:
1 added
5 edited

Legend:

Unmodified
Added
Removed

  • trunk/anuga_core/source/anuga/operators/base_operator.py

    r8454 r8477  
    6060        return self.domain.get_timestep()
    6161
     62
     63    def get_time(self):
     64
     65        return self.domain.get_time()
     66
    6267    def parallel_safe(self):
    6368        """By default an operator is not parallel safe

  • trunk/anuga_core/source/anuga/operators/kinematic_viscosity_operator_ext.c

    r8162 r8477  
    149149
    150150int update_elliptic_matrix(int n,
    151                           int tot_len,
    152                           long *geo_indices,
    153                           double *geo_values,
    154                           double *cell_data,
    155                           double *bdry_data,
    156                           double *data,
    157                           long *colind) {
    158         int i, k, edge, j[4], sorted_j[4], this_index;
    159         double h_j, v[3], v_i; //v[k] = value of the interaction of edge k in a given triangle, v_i = (i,i) entry
    160         for (i=0; i<n; i++) {
    161                 v_i = 0.0;
    162                 j[3] = i;
    163        
     151        int tot_len,
     152        long *geo_indices,
     153        double *geo_values,
     154        double *cell_data,
     155        double *bdry_data,
     156        double *data,
     157        long *colind) {
     158    int i, k, edge, j[4], sorted_j[4], this_index;
     159    double h_j, v[3], v_i; //v[k] = value of the interaction of edge k in a given triangle, v_i = (i,i) entry
     160    for (i = 0; i < n; i++) {
     161        v_i = 0.0;
     162        j[3] = i;
     163
    164164        //Get the values of each interaction, and the column index at which they occur
    165         for (edge=0; edge<3; edge++) {
    166             j[edge] = geo_indices[3*i+edge];
    167             if (j[edge]<n) { //interior
     165        for (edge = 0; edge < 3; edge++) {
     166            j[edge] = geo_indices[3 * i + edge];
     167            if (j[edge] < n) { //interior
    168168                h_j = cell_data[j[edge]];
    169169            } else { //boundary
    170                 h_j = bdry_data[j[edge]-n];
     170                h_j = bdry_data[j[edge] - n];
    171171            }
    172             v[edge] = -0.5*(cell_data[i] + h_j)*geo_values[3*i+edge]; //the negative of the individual interaction
    173             v_i += 0.5*(cell_data[i] + h_j)*geo_values[3*i+edge]; //sum the three interactions
    174         }
    175         if (cell_data[i]<=0.0) {
    176             v_i  = 0.0;
     172            v[edge] = -0.5 * (cell_data[i] + h_j) * geo_values[3 * i + edge]; //the negative of the individual interaction
     173            v_i += 0.5 * (cell_data[i] + h_j) * geo_values[3 * i + edge]; //sum the three interactions
     174        }
     175        if (cell_data[i] <= 0.0) {
     176            v_i = 0.0;
    177177            v[0] = 0.0;
    178178            v[1] = 0.0;
    179179            v[2] = 0.0;
    180180        }
    181                 //Organise the set of 4 values/indices into the data and colind arrays
    182         for (k=0; k<4; k++) sorted_j[k] = j[k];
    183                 quicksort(sorted_j,0,3);
    184                 for (k=0; k<4; k++) { //loop through the nonzero indices
    185                         this_index = sorted_j[k];
    186                         if (this_index == i) {
    187                                 data[4*i+k] = v_i;
    188                                 colind[4*i+k] = i;
    189                         } else if (this_index == j[0]) {
    190                                 data[4*i+k] = v[0];
    191                                 colind[4*i+k] = j[0];
    192                         } else if (this_index == j[1]) {
    193                                 data[4*i+k] = v[1];
    194                                 colind[4*i+k] = j[1];
    195                         } else { //this_index == j[2]
    196                                 data[4*i+k] = v[2];
    197                                 colind[4*i+k] = j[2];
    198                         }
    199                 }
    200         }
    201         return 0;
     181        //Organise the set of 4 values/indices into the data and colind arrays
     182        for (k = 0; k < 4; k++) sorted_j[k] = j[k];
     183        quicksort(sorted_j, 0, 3);
     184        for (k = 0; k < 4; k++) { //loop through the nonzero indices
     185            this_index = sorted_j[k];
     186            if (this_index == i) {
     187                data[4 * i + k] = v_i;
     188                colind[4 * i + k] = i;
     189            } else if (this_index == j[0]) {
     190                data[4 * i + k] = v[0];
     191                colind[4 * i + k] = j[0];
     192            } else if (this_index == j[1]) {
     193                data[4 * i + k] = v[1];
     194                colind[4 * i + k] = j[1];
     195            } else { //this_index == j[2]
     196                data[4 * i + k] = v[2];
     197                colind[4 * i + k] = j[2];
     198            }
     199        }
     200    }
     201    return 0;
    202202}
    203203

  • trunk/anuga_core/source/anuga/operators/rate_operators.py

    r8454 r8477  
    1212from anuga import Domain
    1313from anuga import Quantity
     14from anuga import indent
    1415import numpy as num
    1516import anuga.utilities.log as log
    … …  
    3637                 domain,
    3738                 rate=0.0,
     39                 factor=1.0,
    3840                 indices=None,
    3941                 default_rate=None,
    … …  
    5153        self.rate = rate
    5254        self.indices = indices
     55        self.factor = factor
    5356
    5457        #------------------------------------------
    … …  
    97100        t = self.domain.get_time()
    98101        timestep = self.domain.get_timestep()
     102        factor = self.factor
     103        indices = self.indices
    99104
    100105        rate = self.get_rate(t)
    … …  
    105110
    106111        if self.indices is None:
    107             self.stage_centroid_values[:] = self.stage_centroid_values[:]  \
    108                    + rate*timestep*self.areas[:]
     112            self.stage_c[:] = self.stage_c[:]  \
     113                   + factor*rate*timestep
    109114        else:
    110             self.stage_centroid_values[indices] = self.stage_centroid_values[indices] \
    111                    + rate*timestep*self.areas[indices]
    112 
    113 
    114     def get_rate(self, t):
     115            self.stage_c[indices] = self.stage_c[indices] \
     116                   + factor*rate*timestep
     117
     118
     119    def get_rate(self, t=None):
    115120        """Provide a rate to calculate added volume
    116121        """
    117122
     123        if t is None:
     124            t = self.get_time()
     125
     126           
    118127        if callable(self.rate):
    119128            try:
    … …  
    168177    def timestepping_statistics(self):
    169178
    170         message  = 'You need to implement timestepping statistics for your operator'
     179        message  = indent + self.label + ': Rate = ' + str(self.get_rate())
    171180        return message
     181
     182
    172183
    173184
    … …  
    188199    def __init__(self, domain,
    189200                 rate=0.0,
     201                 factor=1.0,
    190202                 center=None,
    191203                 radius=None,
    … …  
    221233                               domain,
    222234                               rate=rate,
     235                               factor=factor,
    223236                               indices=indices,
    224237                               default_rate=default_rate,
    … …  
    244257    def __init__(self, domain,
    245258                 rate=0.0,
     259                 factor=1.0,
    246260                 polygon=None,
    247261                 default_rate=None,
    248262                 verbose=False):
    249263
    250         assert center is not None
    251         assert radius is not None
     264        assert polygon is not None
    252265
    253266
    … …  
    258271
    259272        indices = inside_polygon(points, polygon)
     273        self.polygon = polygon
    260274
    261275
    … …  
    267281                               domain,
    268282                               rate=rate,
     283                               factor=factor,
    269284                               indices=indices,
    270285                               default_rate=default_rate,
    … …  
    272287
    273288
    274         self.polygon = polygon
    275        
    276 
    277 
    278 
    279 
     289               
     290
     291
     292
     293

  • trunk/anuga_core/source/anuga/operators/run_set_elevation.py

    r8475 r8477  
    1 """Simple test of setting stage in a circular region
     1"""Simple water flow example using ANUGA
     2
     3Water flowing down a channel with a topography that varies with time
    24"""
    35
    … …  
    57# Import necessary modules
    68#------------------------------------------------------------------------------
    7 
    8 import sys
    9 import anuga
    10 
    11 
    12 from math import cos
    13 from numpy import zeros, float, where
    14 import numpy
    15 from time import localtime, strftime, gmtime
    16 
    17 
     9from anuga import rectangular_cross
     10from anuga import Domain
     11from anuga import Reflective_boundary
     12from anuga import Dirichlet_boundary
     13from anuga import Time_boundary
    1814
    1915#------------------------------------------------------------------------------
    20 # Setup domain
     16# Setup computational domain
    2117#------------------------------------------------------------------------------
    22 dx = 100.
    23 dy = dx
    24 L = 10000.
    25 W = L
    26 #W = dx
     18length = 24.
     19width = 5.
     20dx = dy = 0.2 #.1           # Resolution: Length of subdivisions on both axes
    2721
    28 # structured mesh
    29 domain = anuga.rectangular_cross_domain(int(L/dx), int(W/dy), L, W, (-L/2.0, -W/2.0))
    30 
    31 
    32 output_file = 'set_bed'
    33 domain.set_name(output_file)               
    34 
    35 #------------------------------------------------------------------------------
    36 # Setup Algorithm
    37 #------------------------------------------------------------------------------
    38 domain.set_flow_algorithm(2.0)
     22points, vertices, boundary = rectangular_cross(int(length/dx), int(width/dy),
     23                                               len1=length, len2=width)
     24domain = Domain(points, vertices, boundary)
     25domain.set_name('set_elevation') # Output name
     26print domain.statistics()
     27domain.set_quantities_to_be_stored({'elevation': 2,
     28                                    'stage': 2})
    3929
    4030#------------------------------------------------------------------------------
    4131# Setup initial conditions
    4232#------------------------------------------------------------------------------
    43 h0 = 1000.0
     33def topography(x,y):
     34    """Complex topography defined by a function of vectors x and y."""
    4435
    45 domain.set_quantity('elevation',0.0)
    46 domain.set_quantity('friction', 0.0)
    47 domain.add_quantity('stage', h0)
     36    z = -x/100
    4837
    49 #-----------------------------------------------------------------------------
     38    N = len(x)
     39    for i in range(N):
     40        # Step
     41        if 2 < x[i] < 4:
     42            z[i] += 0.4 - 0.05*y[i]
     43
     44        # Permanent pole
     45        if (x[i] - 8)**2 + (y[i] - 2)**2 < 0.4**2:
     46            z[i] += 1
     47    return z
     48
     49
     50def pole_increment(x,y,t):
     51    """This provides a small increment to a pole located mid stream
     52    For use with variable elevation data
     53    """
     54
     55
     56    z = 0.0*x
     57   
     58
     59    if t<10.0:
     60        return z
     61   
     62
     63    N = len(x)
     64    for i in range(N):
     65        # Pole 1
     66        if (x[i] - 12)**2 + (y[i] - 3)**2 < 0.4**2:
     67            z[i] += 0.1
     68
     69    for i in range(N):
     70        # Pole 2
     71        if (x[i] - 14)**2 + (y[i] - 2)**2 < 0.4**2:
     72            z[i] += 0.05
     73
     74    return z
     75
     76
     77def pole(t):
     78
     79    if t<10:
     80        return 0.0
     81    elif t>15:
     82        return 0.0
     83    else:
     84        return 0.05
     85
     86
     87domain.set_quantity('elevation', topography)           # elevation is a function
     88domain.set_quantity('friction', 0.01)                  # Constant friction
     89domain.set_quantity('stage', expression='elevation')   # Dry initial condition
     90
     91#------------------------------------------------------------------------------
    5092# Setup boundary conditions
    5193#------------------------------------------------------------------------------
    52 Br = anuga.Reflective_boundary(domain)      # Solid reflective wall
     94Bi = Dirichlet_boundary([0.4, 0, 0])          # Inflow
     95Br = Reflective_boundary(domain)              # Solid reflective wall
     96Bo = Dirichlet_boundary([-5, 0, 0])           # Outflow
    5397
    54 # Associate boundary tags with boundary objects
    55 domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
    56 
    57 #------------------------------------------------------------------------------
    58 # Setup Operators
    59 #------------------------------------------------------------------------------
    60 from anuga.operators.set_bed_operators import Circular_set_bed_operator
    61 
    62 import math
    63 
    64 bed1 = lambda t: 20.0 * math.sin(t/3.0)
    65 cop1 = Circular_set_bed_operator(domain, bed=bed1, center=(0.0, 0.0), radius=100.0 )
    66 
    67 #stage2 = lambda t: h0 + 30.0 * math.sin(t/6.0)
    68 #cop2 = Circular_set_stage_operator(domain, stage=stage2, center=(2000.0, 1000.0), radius=100.0 )
    69 
    70 #print cop1.statistics()
    71 #print cop2.statistics()
     98domain.set_boundary({'left': Bi, 'right': Bo, 'top': Br, 'bottom': Br})
    7299
    73100#------------------------------------------------------------------------------
    … …  
    75102#------------------------------------------------------------------------------
    76103
    77 for t in domain.evolve(yieldstep = 1.0, finaltime = 20.0):
    78     #print domain.timestepping_statistics(track_speeds=True)
     104from anuga.operators.set_elevation_operators import Circular_set_elevation_operator
     105
     106op1 = Circular_set_elevation_operator(domain, elevation=pole, radius=0.5, center = (12.0,3.0))
     107
     108growing = False
     109shrinking = False
     110
     111done = False
     112for t in domain.evolve(yieldstep=0.1, finaltime=40.0):
    79113    domain.print_timestepping_statistics()
    80114    domain.print_operator_timestepping_statistics()
    81115
     116    #w = domain.get_quantity('stage').\
     117    #    get_values(interpolation_points=[[18, 2.5]])
     118    #print 'Level at gauge point = %.2fm' % w
    82119
     120    #z = domain.get_quantity('elevation').\
     121    #    get_values(interpolation_points=[[12, 3]])
     122    #print 'Elevation at pole location = %.2fm' % z
     123
     124#    # Start variable elevation after 10 seconds
     125#    if t > 10 and not (shrinking or growing or done):
     126#        growing = True
     127#
     128#    # Stop growing when pole has reached a certain height
     129#    if t > 16 and growing:
     130#        growing = False
     131#        shrinking = False
     132#
     133#    # Start shrinking
     134#    if t > 20:
     135#        shrinking = True
     136#        growing = False
     137#
     138#    # Stop changing when pole has shrunk to original level
     139#    if t > 25 and shrinking:
     140#        done = True
     141#        shrinking = growing = False
     142#        domain.set_quantity('elevation', topography)
     143#
     144#    # Grow or shrink
     145#    if growing:
     146#        domain.add_quantity('elevation', pole_increment)
     147#
     148#    if shrinking:
     149#        domain.add_quantity('elevation', lambda x,y: -2*pole_increment(x,y))
     150
     151
     152
     153
     154

  • trunk/anuga_core/source/anuga/operators/set_elevation_operators.py

    r8475 r8477  
    6161        """
    6262
    63         if self.indices is []:
    64             return
    65 
    66         elevation = self.get_elevation()
    67 
    68         if self.verbose is True:
    69             log.critical('Bed of %s at time = %.2f = %f'
    70                          % (self.quantity_name, domain.get_time(), elevation))
     63        #if self.indices is []:
     64        #    return
     65
     66        #elevation = self.get_elevation()
     67
     68#        if self.verbose is True:
     69#            log.critical('Bed of %s at time = %.2f = %f'
     70#                         % (self.quantity_name, domain.get_time(), elevation))
     71
     72        #if self.indices is None:
     73        #    self.elev_c[:] = elevation
     74        #else:
     75        #    self.elev_c[self.indices] = elevation
     76
     77        t = self.get_time()
     78        dt = self.get_timestep()
     79
     80        v_coors = self.domain.vertex_coordinates
     81        self.elev_v = self.domain.quantities['elevation'].vertex_values
     82
    7183
    7284        if self.indices is None:
    73             self.elevation_c[:] = elevation
     85            self.elev_v[:] = self.elev_v + 0.0
    7486        else:
    75             self.elevation_c[self.indices] = elevation
     87            self.elev_v[self.indices] += self.elevation(t)*dt
     88
     89        ### make sure centroid is correct as well
     90       
     91        #self.domain.add_quantity('elevation', lambda x,y: dt*self.elevation(x,y,t))
     92
     93
     94
     95        # clean up discontinuities for now
     96        #self.domain.quantities['elevation'].smooth_vertex_values()
     97
     98
    7699
    77100
    … …  
    122145    def timestepping_statistics(self):
    123146
    124         message  = indent + self.label + ': Set_elevation = ' + str(self.get_elevation())
    125         message  += ' at center '+str(self.center)
    126         return message
     147        #message  = indent + self.label + ': Set_elevation = ' + str('')
     148        #message  += ' at center '+str(self.center)
     149        return 'test'
    127150
    128151
Note: See TracChangeset for help on using the changeset viewer.