Changeset 8944


Ignore:
Timestamp:
Jul 5, 2013, 10:04:02 PM (12 years ago)
Author:
steve
Message:

Mostly working on operators

Location:
trunk/anuga_core/source
Files:
11 edited
1 moved

Legend:

Unmodified
Added
Removed

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

    r8828 r8944  
    173173
    174174        for name in self.evolved_quantities:
    175             self.quantities[name] = Quantity(self)
     175            self.quantities[name] = Quantity(self, name=name)
    176176        for name in self.other_quantities:
    177             self.quantities[name] = Quantity(self)
     177            self.quantities[name] = Quantity(self, name=name)
    178178
    179179        # Create an empty list for forcing terms
    … …  
    317317        self.timestep = 0.0
    318318        self.flux_timestep = 0.0
     319        self.evolved_called = False
    319320
    320321        self.last_walltime = walltime()
    … …  
    13491350        return self.simulation_name
    13501351
    1351     def set_name(self, name):
     1352    def set_name(self, name=None, timestamp=False):
    13521353        """Assign a name to this simulation.
    13531354        This will be used to identify the output sww file.
    1354         """
    1355 
     1355        Without parameters the name will be derived from the script file,
     1356        ie run_simulation.py -> output_run_simulation.sww
     1357        """
     1358
     1359        import os
     1360        import inspect
     1361
     1362        if name is None:
     1363            frame = inspect.currentframe()
     1364            script_name = inspect.getouterframes(frame)[1][1]
     1365            name = 'output_'+os.path.splitext(script_name)[0]
     1366       
    13561367        # remove any '.sww' end
    13571368        if name.endswith('.sww'):
    13581369            name = name[:-4]
     1370
     1371
     1372        from time import localtime, strftime, gmtime
     1373        time = strftime('%Y%m%d_%H%M%S',localtime())
     1374
     1375        if timestamp:
     1376            name = name+'_'+time
    13591377
    13601378        self.simulation_name = name
    … …  
    13961414
    13971415
    1398         # Proc 0 gathers full and ghost nodes from self and other processors
     1416        # Gather full and ghost nodes
    13991417        fx = vertices[full_mask,0]
    14001418        fy = vertices[full_mask,1]
    … …  
    14661484                   % self.get_boundary_tags())
    14671485        assert hasattr(self, 'boundary_objects'), msg
    1468        
     1486
     1487
     1488        if self.evolved_called:
     1489            skip_initial_step = True
     1490
     1491        self.evolved_called = True
     1492
     1493        # We assume evolve has already been called so we should now
     1494        # set starttime to match actual time
     1495        if skip_initial_step:
     1496            self.set_starttime(self.get_time())
     1497
     1498
     1499        # This can happen on the first call to evolve
    14691500        if self.get_time() != self.get_starttime():
    14701501            self.set_time(self.get_starttime())
    … …  
    14771508        self._order_ = self.default_order
    14781509
    1479         assert finaltime >= self.get_starttime(), 'finaltime is less than starttime!'
     1510
    14801511       
    14811512        if finaltime is not None and duration is not None:
    … …  
    14871518            if duration is not None:
    14881519                self.finaltime = self.starttime + float(duration)
     1520
     1521        assert self.finaltime >= self.get_starttime(), 'finaltime is less than starttime!'
    14891522
    14901523        N = len(self)                             # Number of triangles

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

    r8930 r8944  
    3333class Quantity:
    3434
    35     def __init__(self, domain, vertex_values=None):
     35
     36    counter = 0
     37
     38    def __init__(self, domain, vertex_values=None, name=None):
    3639        from anuga.abstract_2d_finite_volumes.generic_domain \
    3740                            import Generic_Domain
    … …  
    8588
    8689        self.set_beta(1.0)
     90
     91        # Keep count of quantities
     92        Quantity.counter += 1
     93
     94        self.set_name(name)
     95
     96
     97
     98
     99
    87100
    88101    ############################################################################
    … …  
    248261        self.centroid_values[:] = num.maximum(self.centroid_values, Q.centroid_values)
    249262
     263        return self
     264
    250265       
    251266
    … …  
    273288        self.centroid_values[:] = num.minimum(self.centroid_values, Q.centroid_values)
    274289
    275 
     290        return self
    276291
    277292
    … …  
    279294    # Setters/Getters
    280295    ############################################################################
    281 
     296   
     297    def save_centroid_data_to_csv(self,filename=None):
     298
     299
     300        #FIXME SR: Should add code to deal with parallel
     301       
     302        c_v = self.centroid_values.reshape((-1,1))
     303        c_x = self.domain.centroid_coordinates[:,0].reshape((-1,1))
     304        c_y = self.domain.centroid_coordinates[:,1].reshape((-1,1))
     305
     306        import numpy
     307
     308        c_xyv = numpy.hstack((c_x, c_y, c_v))
     309
     310        if filename is None:
     311            filename = self.name+'_centroid_data.csv'
     312
     313        numpy.savetxt(filename, c_xyv, delimiter=',', fmt =  ['%.15e', '%.15e', '%.15e' ])
     314
     315
     316
     317    def plot_quantity(self, filename=None, show=True):
     318
     319        X, Y, A, V = self.get_vertex_values(smooth=True)
     320
     321        import matplotlib.pyplot as plt
     322        import numpy as np
     323
     324
     325        plt.tripcolor(X, Y, V, A)
     326        plt.colorbar()
     327
     328        if filename is None:
     329            filename = self.name+'.png'
     330
     331        plt.savefig(filename)
     332       
     333        if show:
     334            plt.show()
     335       
     336
     337
     338
     339    def set_name(self, name=None):
     340       
     341        if name is not None:
     342            self.name = name
     343        else:
     344            self.name = 'quantity_%g' % Quantity.counter
     345
     346
     347
     348    def get_name(self):
     349
     350        return self.name
     351
     352   
    282353    def set_beta(self, beta):
    283354        """Set default beta value for limiting """
    … …  
    13001371        self._set_vertex_values(vertex_list, A)
    13011372
     1373
     1374
    13021375    # Note Padarn 27/11/12:
    13031376    # This function has been changed and now uses an external c function

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

    r8939 r8944  
    3535        # Setup a quantity to store max_stage
    3636        #------------------------------------------
    37         self.max_stage = Quantity(domain)
     37        self.max_stage = Quantity(domain, name = 'max_stage')
    3838        self.max_stage.set_values(-1.0e+100)
    3939
    … …  
    5050        """
    5151
    52         self.max_stage.maximum(self.stage) #=self.max_stage.maximum(self.stage)
    53         #self.max_stage=self.max_stage.maximum(self.stage)
     52        self.max_stage.maximum(self.stage)
    5453
    5554
    … …  
    7372
    7473
     74    def save_centroid_data_to_csv(self, filename=None):
     75
     76        self.max_stage.save_centroid_data_to_csv(filename)
     77
     78
     79    def plot_quantity(self, filename=None, show=True):
     80
     81        self.max_stage.plot_quantity(filename=filename, show=show)
     82
     83
     84
     85

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

    r8574 r8944  
    1212import numpy as num
    1313from anuga.geometry.polygon import inside_polygon
     14from anuga import indent
     15
     16from anuga import Domain
     17from anuga import Quantity
    1418from anuga.operators.base_operator import Operator
    15 from anuga import indent
    16 
    17 
    18 class Erosion_operator(Operator):
     19from anuga.operators.region import Region
     20
     21
     22class Erosion_operator(Operator, Region):
    1923    """
    2024    Simple erosion operator in a region (careful to maintain continuitiy of elevation)
    … …  
    3236                 base=0.0,
    3337                 indices=None,
     38                 polygon=None,
     39                 center=None,
     40                 radius=None,
    3441                 description = None,
    3542                 label = None,
    … …  
    4047        Operator.__init__(self, domain, description, label, logging, verbose)
    4148
     49
     50
     51        Region.__init__(self, domain,
     52                        indices=indices,
     53                        polygon=polygon,
     54                        center=center,
     55                        radius=radius,
     56                        verbose=verbose)
     57
    4258        #------------------------------------------
    4359        # Local variables
    … …  
    4561        self.threshold = threshold
    4662        self.base = base
    47         self.indices = indices
     63
    4864
    4965        #------------------------------------------
    … …  
    409425                 verbose=False):
    410426
    411         assert center is not None
    412         assert radius is not None
    413 
    414 
    415         # Determine indices in update region
    416         N = domain.get_number_of_triangles()
    417         points = domain.get_centroid_coordinates(absolute=True)
    418 
    419 
    420         indices = []
    421 
    422         c = center
    423         r = radius
    424 
    425         self.center = center
    426         self.radius = radius
    427 
    428         intersect = False
    429         for k in range(N):
    430             x, y = points[k,:]    # Centroid
    431 
    432             if ((x-c[0])**2+(y-c[1])**2) < r**2:
    433                 intersect = True
    434                 indices.append(k)
    435 
    436 
    437         msg = 'No centroids intersect circle center'+str(center)+' radius '+str(radius)
    438         assert intersect, msg
    439 
    440 
    441 
    442 
    443         # It is possible that circle doesn't intersect with mesh (as can happen
    444         # for parallel runs
    445427
    446428
    … …  
    449431                                  threshold,
    450432                                  base,
    451                                   indices=indices,
     433                                  center=center,
     434                                  radius=radius,
    452435                                  verbose=verbose)
    453436
    … …  
    472455
    473456
    474         # Determine indices in update region
    475         N = domain.get_number_of_triangles()
    476         points = domain.get_centroid_coordinates(absolute=True)
    477 
    478 
    479         indices = inside_polygon(points, polygon)
    480         self.polygon = polygon
    481 
    482         # It is possible that circle doesn't intersect with mesh (as can happen
    483         # for parallel runs
    484 
    485 
    486457        Erosion_operator.__init__(self,
    487458                                  domain,
    488459                                  threshold=threshold,
    489460                                  base=base,
    490                                   indices=indices,
     461                                  polygon=polygon,
    491462                                  verbose=verbose)
     463
     464
     465
    492466
    493467

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

    r8879 r8944  
    1010
    1111
    12 from anuga import Domain
    13 from anuga import Quantity
     12
    1413from anuga import indent
    1514import numpy as num
    1615from warnings import warn
    1716import anuga.utilities.log as log
    18 
    1917from anuga.geometry.polygon import inside_polygon
    2018
    21 from anuga.operators.base_operator import Operator
     19
    2220from anuga.fit_interpolate.interpolate import Modeltime_too_early, \
    2321                                              Modeltime_too_late
    2422from anuga.utilities.function_utils import evaluate_temporal_function
    2523
    26 
    27 class Rate_operator(Operator):
     24from anuga import Domain
     25from anuga import Quantity
     26from anuga.operators.base_operator import Operator
     27from anuga.operators.region import Region
     28
     29class Rate_operator(Operator,Region):
    2830    """
    2931    Add water at certain rate (ms^{-1} = vol/Area/sec) over a
    … …  
    4345                 factor=1.0,
    4446                 indices=None,
     47                 polygon=None,
     48                 center=None,
     49                 radius=None,
    4550                 default_rate=None,
    4651                 description = None,
    … …  
    5257        Operator.__init__(self, domain, description, label, logging, verbose)
    5358
     59
     60        Region.__init__(self, domain,
     61                        indices=indices,
     62                        polygon=polygon,
     63                        center=center,
     64                        radius=radius,
     65                        verbose=verbose)
     66
     67
    5468        #------------------------------------------
    5569        # Local variables
    5670        #------------------------------------------
    57         self.indices = indices
    5871        self.factor = factor
    5972
    … …  
    131144
    132145        rate = evaluate_temporal_function(self.rate, t, default_right_value=self.default_rate)
    133 #        if  self.rate_callable:
    134 #            try:
    135 #                rate = self.rate(t)
    136 #            except Modeltime_too_early, e:
    137 #                raise Modeltime_too_early(e)
    138 #            except Modeltime_too_late, e:
    139 #                if self.default_rate is None:
    140 #                    msg = '%s: ANUGA is trying to run longer than specified data.\n' %str(e)
    141 #                    msg += 'You can specify keyword argument default_rate in the '
    142 #                    msg += 'rate operator to tell it what to do in the absence of time data.'
    143 #                    raise Modeltime_too_late(msg)
    144 #                else:
    145 #                    # Pass control to default rate function
    146 #                    rate = self.default_rate(t)
    147 #
    148 #                    if self.default_rate_invoked is False:
    149 #                        # Issue warning the first time
    150 #                        msg = ('\n%s'
    151 #                           'Instead I will use the default rate: %s\n'
    152 #                           'Note: Further warnings will be supressed'
    153 #                           % (str(e), self.default_rate(t)))
    154 #                        warn(msg)
    155 #
    156 #                        # FIXME (Ole): Replace this crude flag with
    157 #                        # Python's ability to print warnings only once.
    158 #                        # See http://docs.python.org/lib/warning-filter.html
    159 #                        self.default_rate_invoked = True
    160 #        else:
    161 #            rate = self.rate
    162 
    163146
    164147        if rate is None:
    … …  
    190173        assert x is not None
    191174        assert y is not None
    192         #print x.shape,y.shape
     175
    193176        assert isinstance(t, (int, float))
    194177        assert len(x) == len(y)
    … …  
    362345                 verbose=False):
    363346
    364         assert center is not None
    365         assert radius is not None
    366 
    367 
    368         # Determine indices in update region
    369         N = domain.get_number_of_triangles()
    370         points = domain.get_centroid_coordinates(absolute=True)
    371 
    372 
    373         indices = []
    374 
    375         c = center
    376         r = radius
    377 
    378         for k in range(N):
    379             x, y = points[k,:]    # Centroid
    380 
    381             if ((x-c[0])**2+(y-c[1])**2) < r**2:
    382                 indices.append(k)
    383 
    384 
    385         # It is possible that circle doesn't intersect with mesh (as can happen
    386         # for parallel runs
    387 
    388347
    389348        Rate_operator.__init__(self,
    … …  
    391350                               rate=rate,
    392351                               factor=factor,
    393                                indices=indices,
     352                               center=center,
     353                               radius=radius,
    394354                               default_rate=default_rate,
    395355                               verbose=verbose)
    396356
    397 
    398         self.center = center
    399         self.radius = radius
    400357
    401358
    … …  
    421378                 verbose=False):
    422379
    423         assert polygon is not None
    424 
    425 
    426         # Determine indices in update region
    427         N = domain.get_number_of_triangles()
    428         points = domain.get_centroid_coordinates(absolute=True)
    429 
    430 
    431         indices = inside_polygon(points, polygon)
    432         self.polygon = polygon
    433 
    434 
    435         # It is possible that circle doesn't intersect with mesh (as can happen
    436         # for parallel runs
    437 
    438380
    439381        Rate_operator.__init__(self,
    … …  
    441383                               rate=rate,
    442384                               factor=factor,
    443                                indices=indices,
     385                               polygon=polygon,
    444386                               default_rate=default_rate,
    445387                               verbose=verbose)
     388                               
    446389
    447390

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

    r8480 r8944  
    77# Import necessary modules
    88#------------------------------------------------------------------------------
     9import numpy
    910from anuga import rectangular_cross
    1011from anuga import Domain
    … …  
    2324                                               len1=length, len2=width)
    2425domain = Domain(points, vertices, boundary)
    25 domain.set_name('change_elevation') # Output name
     26domain.set_name()
    2627print domain.statistics()
    2728domain.set_quantities_to_be_stored({'elevation': 2,
    28                                     'stage': 2})
     29                                    'stage': 2,
     30                                    'xmomentum': 2,
     31                                    'ymomentum': 2})
    2932
    3033#------------------------------------------------------------------------------
    … …  
    3639    z = -x/100
    3740
    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
     41    # Step
     42    id = (2 < x) & (x < 4)
     43    z[id] += 0.4 - 0.05*y[id]
    4744
    4845
    49         # Dam
    50         if 12 < x[i] < 13:
    51             z[i] += 0.4
     46    # Permanent pole
     47    #id = (x - 8)**2 + (y - 2)**2 < 0.4**2
     48    #z[id] += 1
     49
     50
     51    # Dam
     52    id = (12 < x) & (x  < 13)
     53    z[id] += 0.4
     54
    5255
    5356           
    … …  
    5861
    5962    z = -x/100
    60 
    61     N = len(x)
    62     for i in range(N):
    63         # Step
    64         if 2 < x[i] < 4:
    65             z[i] += 0.4 - 0.05*y[i]
    66 
    67         # Permanent pole
    68         #if (x[i] - 8)**2 + (y[i] - 2)**2 < 0.4**2:
    69         #    z[i] += 1
     63   
     64    # Step
     65    id = (2 < x) & (x < 4)
     66    z[id] += 0.4 - 0.05*y[id]
    7067
    7168
    72         # Dam
    73         if 12 < x[i] < 13 and y[i] > 3.0:
    74             z[i] += 0.4
     69    # Permanent pole
     70    #id = (x - 8)**2 + (y - 2)**2 < 0.4**2
     71    #z[id] += 1
    7572
    76         if 12 < x[i] < 13 and y[i] < 2.0:
    77             z[i] += 0.4
     73    # Dam with hole
     74    id = (12 < x) & (x < 13) and (y > 3.0)
     75    z[id] += 0.4
     76
     77    id = (12 < x) & (x < 13) and (y < 2.0)
     78    z[id] += 0.4
     79
    7880
    7981    return z
    … …  
    9799#------------------------------------------------------------------------------
    98100
    99 #from anuga.operators.set_elevation_operators import Circular_set_elevation_operator
    100 
    101 #op1 = Circular_set_elevation_operator(domain, elevation=pole, radius=0.5, center = (12.0,3.0))
     101from anuga.operators.set_elevation import Set_elevation
     102op1 = Set_elevation(domain, elevation=lambda x,y : -x/100, radius=1.0, center = (12.5,3.0))
    102103
    103104
    104 dam_break = False
     105#dam_break = False
    105106
    106107for t in domain.evolve(yieldstep=0.1, finaltime=30.0):
    … …  
    108109    domain.print_operator_timestepping_statistics()
    109110
    110     if t >= 10 and not dam_break:
     111    if numpy.allclose(t, 10.0):
    111112        print 'changing elevation'
    112         domain.set_quantity('elevation', topography_dam_break)
    113         dam_break = True
     113        op1()
    114114
    115115
    … …  
    117117
    118118
     119

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

    r8936 r8944  
    1919
    2020
    21 
    22 #-------------------------------------------------------------------------------
    23 # Copy scripts to time stamped output directory and capture screen
    24 # output to file
    25 #-------------------------------------------------------------------------------
    26 #time = strftime('%Y%m%d_%H%M%S',localtime())
    27 
    28 #output_dir = 'dam_break_'+time
    29 output_file = 'dam_break'
    30 
    31 #anuga.copy_code_files(output_dir,__file__)
    32 #start_screen_catcher(output_dir+'_')
    33 
    34 
    3521#------------------------------------------------------------------------------
    3622# Setup domain
    … …  
    4531domain = anuga.rectangular_cross_domain(int(L/dx), int(W/dy), L, W, (0.0, -W/2))
    4632
    47 domain.set_name(output_file)               
     33domain.set_name() # based on script name
    4834
    4935#------------------------------------------------------------------------------
    5036# Setup Algorithm
    5137#------------------------------------------------------------------------------
    52 domain.set_timestepping_method('rk2')
    53 domain.set_default_order(2)
    54 domain.set_beta(1.7)
    55 
    56 #------------------------------------------------------------------------------
    57 # Setup Kinematic Viscosity Operator
    58 #------------------------------------------------------------------------------
    59 domain.set_use_kinematic_viscosity(True)
     38domain.set_flow_algorithm('2_0')
    6039
    6140#------------------------------------------------------------------------------
    … …  
    7655domain.set_quantity('friction', 0.0)
    7756domain.add_quantity('stage', height)
     57
    7858#-----------------------------------------------------------------------------
    7959# Setup boundary conditions
    … …  
    8969
    9070
    91 ##===============================================================================
    92 #from anuga.visualiser import RealtimeVisualiser
    93 #vis = RealtimeVisualiser(domain)
    94 #vis.render_quantity_height("stage", zScale = 50000/(h0-h1), dynamic=True)
    95 #vis.colour_height_quantity("stage",
    96 ##        (lambda q: numpy.sqrt( (q["xvelocity"]**2 ) + (q['yvelocity']**2 )), 0.0, 10.0) )
    97 #        (lambda q: q["yvelocity"], -10.0, 10.0) )
    98 #vis.start()
    99 ##===============================================================================
    10071
    101 
     72#-----------------------------------------------------------------------------
     73# Setup operators that will be applied each inner timestep
     74#------------------------------------------------------------------------------
    10275from anuga.operators.collect_max_stage_operator import Collect_max_stage_operator
    10376max_operator = Collect_max_stage_operator(domain)
    … …  
    10679# Evolve system through time
    10780#------------------------------------------------------------------------------
    108 
    10981for t in domain.evolve(yieldstep = 100.0, finaltime = 60*60.):
    11082    #print domain.timestepping_statistics(track_speeds=True)
    … …  
    11284    domain.print_operator_timestepping_statistics()
    11385
    114     #vis.update()
     86
     87
     88# save the max_stage centroid data to a text file
     89max_operator.save_centroid_data_to_csv()
     90
     91max_operator.plot_quantity()
    11592
    11693
    11794
    118 #print max_operator.max_stage.centroid_values
    11995
    120 c_v = max_operator.max_stage.centroid_values.flatten()
    121 c_x = domain.centroid_coordinates[:,0].flatten()
    122 c_y = domain.centroid_coordinates[:,1].flatten()
    123 
    124 import numpy
    125 
    126 numpy.savetxt('max_stage.txt', c_v)
    127 
    128 #test against know data
    129    
    130 #vis.evolveFinished()
    131 

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

    r8506 r8944  
    112112length = 24.
    113113width = 5.
    114 dx = dy = 0.1 #.1           # Resolution: Length of subdivisions on both axes
     114dx = dy = 0.2 #.1           # Resolution: Length of subdivisions on both axes
    115115
    116116points, vertices, boundary = rectangular_cross(int(length/dx), int(width/dy),
    117117                                               len1=length, len2=width)
    118118domain = Domain(points, vertices, boundary)
    119 domain.set_name('polygon_erosion') # Output name
     119domain.set_name() # Output name based on script
    120120print domain.statistics()
    121121domain.set_quantities_to_be_stored({'elevation': 2,

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

    r8477 r8944  
    1212from anuga import Dirichlet_boundary
    1313from anuga import Time_boundary
     14import os
    1415
    1516
    … …  
    2425                                               len1=length, len2=width)
    2526domain = Domain(points, vertices, boundary)
    26 domain.set_name('rate_polygon') # Output name
     27domain.set_name() # Output name based on script name. You can add timestamp=True
    2728print domain.statistics()
    2829
    … …  
    3637    z = -x/100
    3738
    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
     39    # Step
     40    id = (2 < x) & (x < 4)
     41    z[id] += 0.4 - 0.05*y[id]
    4742
    4843
    49 #        # Dam
    50 #        if 12 < x[i] < 13 and y[i] > 3.0:
    51 #            z[i] += 0.4
    52 #
    53 #        if 12 < x[i] < 13 and y[i] < 2.0:
    54 #            z[i] += 0.4
     44    # Permanent pole
     45    #id = (x - 8)**2 + (y - 2)**2 < 0.4**2
     46    #z[id] += 1
    5547
    5648
    57 #        # Dam
    58 #        if 12 < x[i] < 13:
    59 #            z[i] += 0.4
     49    # Dam
     50    #id = (12 < x) & (x  < 13)
     51    #z[id] += 0.4
     52
     53
    6054
    6155           
    … …  
    6862    """
    6963
    70 
    71     z = 0.0*x
    72    
     64    z = 0.0*x   
    7365
    7466    if t<10.0:
    … …  
    7668   
    7769
    78     N = len(x)
    79     for i in range(N):
    80         # Pole 1
    81         if (x[i] - 12)**2 + (y[i] - 3)**2 < 0.4**2:
    82             z[i] += 0.1
     70    # Pole 1
     71    id = (x - 12)**2 + (y - 3)**2 < 0.4**2
     72    z[id] += 0.1
    8373
    84     for i in range(N):
    85         # Pole 2
    86         if (x[i] - 14)**2 + (y[i] - 2)**2 < 0.4**2:
    87             z[i] += 0.05
     74
     75    # Pole 2
     76    id = (x - 14)**2 + (y - 2)**2 < 0.4**2
     77    z[id] += 0.05
     78
    8879
    8980    return z
    … …  
    119110polygon2 = [ [12.0, 2.0], [13.0, 2.0], [13.0, 3.0], [12.0, 3.0] ]
    120111
    121 from anuga.operators.rate_operators import Polygonal_rate_operator
    122 from anuga.operators.rate_operators import Circular_rate_operator
     112from anuga.operators.rate_operators import Rate_operator
    123113
    124 op1 = Polygonal_rate_operator(domain, rate=10.0, polygon=polygon2)
    125 op2 = Circular_rate_operator(domain, rate=10.0, radius=0.5, center=(10.0, 3.0))
     114op1 = Rate_operator(domain, rate=lambda t: 10.0 if (t>=0.0) else 0.0, polygon=polygon2)
     115op2 = Rate_operator(domain, rate=lambda t: 10.0 if (t>=0.0) else 0.0, radius=0.5, center=(10.0, 3.0))
    126116
    127117
    128 for t in domain.evolve(yieldstep=0.1, finaltime=40.0):
     118domain.set_starttime(-0.1)
     119for t in domain.evolve(yieldstep=0.01, finaltime=0.0):
     120    domain.print_timestepping_statistics()
     121    domain.print_operator_timestepping_statistics()
     122
     123    stage = domain.get_quantity('stage')
     124    elev  = domain.get_quantity('elevation')
     125    height = stage - elev
     126
     127    print 'integral = ', height.get_integral()
     128
     129
     130for t in domain.evolve(yieldstep=0.1, duration=5.0):
     131
    129132    domain.print_timestepping_statistics()
    130133    domain.print_operator_timestepping_statistics()

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

    r8871 r8944  
    2727                                               len1=length, len2=width)
    2828domain = Domain(points, vertices, boundary)
    29 domain.set_name('output_rate_spatial_operator') # Output name
     29domain.set_name() # Output name based on script
    3030print domain.statistics()
    3131
    … …  
    119119polygon2 = [ [12.0, 2.0], [13.0, 2.0], [13.0, 3.0], [12.0, 3.0] ]
    120120
    121 from anuga.operators.rate_operators import Polygonal_rate_operator
    122 from anuga.operators.rate_operators import Circular_rate_operator
     121
    123122from anuga.operators.rate_operators import Rate_operator
    124123
    125 op1 = Polygonal_rate_operator(domain, rate=10.0, polygon=polygon2)
     124op1 = Rate_operator(domain, rate=10.0, polygon=polygon2)
    126125
    127126area1 = numpy.sum(domain.areas[op1.indices])
    … …  
    129128print 'op1 Q ',Q1
    130129
    131 op2 = Circular_rate_operator(domain, rate=10.0, radius=0.5, center=(10.0, 3.0))
     130op2 = Rate_operator(domain, rate=10.0, radius=0.5, center=(10.0, 3.0))
    132131
    133132area2 = numpy.sum(domain.areas[op2.indices])
    … …  
    141140    abd t a scalar
    142141    """
    143     if t<=5.0:
     142    if t<=4.0:
    144143        return (x+y)
    145144    else:
    … …  
    160159accum = 0.0
    161160yieldstep = 0.1
    162 finaltime = 40.0
     161finaltime = 5.0
    163162for t in domain.evolve(yieldstep=yieldstep, finaltime=finaltime):
    164163    domain.print_timestepping_statistics()

  • trunk/anuga_core/source/anuga/utilities/log.py

    r8427 r8944  
    6464
    6565# The default name of the file to log to.
     66
     67
    6668log_filename = os.path.join('.', 'anuga.log')
    6769

  • trunk/anuga_core/source/anuga_validation_tests/case_studies/patong/run_model.py

    r8922 r8944  
    4343
    4444# Tell log module to store log file in output dir
    45 log.log_filename = os.path.join(project.output_run, 'anuga.log')
     45log.log_filename = os.path.join(project.output_run, 'anuga_P_%g.log'%myid)
    4646
    4747if(myid==0):
    … …  
    221221## conflated in the code
    222222##
    223 ## Skip over the first 6000 seconds
    224 #for t in domain.evolve(yieldstep=2000,
    225 #                       finaltime=6000):
    226 #    log.critical(domain.timestepping_statistics())
    227 #    log.critical(domain.boundary_statistics(tags='ocean'))
    228 #
    229 ## Start detailed model
    230 #for t in domain.evolve(yieldstep=project.yieldstep,
    231 #                       finaltime=project.finaltime,
    232 #                       skip_initial_step=True):
    233 #    log.critical(domain.timestepping_statistics())
    234 #    log.critical(domain.boundary_statistics(tags='ocean'))
     223## Skip over the first 60 seconds
     224for t in domain.evolve(yieldstep=2,
     225                       duration=60):
     226    log.critical(domain.timestepping_statistics())
     227    log.critical(domain.boundary_statistics(tags='ocean'))
    235228
    236229
    … …  
    239232import time
    240233time.sleep(myid)
    241 print 'Processor ', myid , ' : ', project.yieldstep, domain.flow_algorithm 
     234print 'Processor ', myid , ' : ', project.yieldstep, domain.flow_algorithm
    242235barrier()
    243236print 'project.yieldstep', project.yieldstep
     237# Start detailed model
    244238for t in domain.evolve(yieldstep=project.yieldstep,
    245239                       finaltime=project.finaltime):
    … …  
    250244    print 'Processor ', myid, ' project.yieldstep ', project.yieldstep
    251245
     246
    252247## Final print out
    253248barrier()
Note: See TracChangeset for help on using the changeset viewer.