Changeset 4712


Ignore:
Timestamp:
Sep 7, 2007, 10:43:46 PM (15 years ago)
Author:
steve
Message:

Working on 2nd order time stepping

Location:
anuga_core/source/anuga
Files:
10 edited

Legend:

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

    r4711 r4712  
    3131
    3232import types
     33
     34
    3335
    3436class Domain(Mesh):
     
    176178        from anuga.config import max_smallsteps, beta_w, beta_h, epsilon
    177179        from anuga.config import CFL
     180        from anuga.config import timestepping_method
    178181        from anuga.config import protect_against_isolated_degenerate_timesteps
    179182        self.beta_w = beta_w
     
    182185        self.protect_against_isolated_degenerate_timesteps = protect_against_isolated_degenerate_timesteps
    183186       
     187       
    184188
    185189        # FIXME: Maybe have separate orders for h-limiter and w-limiter?
     
    194198        self.number_of_first_order_steps = 0
    195199        self.CFL = CFL
    196 
     200        self.set_timestepping_method(timestepping_method)
     201       
    197202        self.boundary_map = None  # Will be populated by set_boundary       
    198203       
     
    954959        return msg
    955960
     961    def get_timestepping_method(self):
     962        return self.timestepping_method
     963
     964    def set_timestepping_method(self,timestepping_method):
     965       
     966        if timestepping_method in ['euler', 'rk2']:
     967            self.timestepping_method = timestepping_method
     968            return
     969
     970        msg = '%s is an incorrect timestepping type'% timestepping_method
     971        raise Exception, msg
    956972
    957973    def get_name(self):
     
    973989    def set_datadir(self, name):
    974990        self.datadir = name
     991
     992    def get_starttime(self):
     993        return self.starttime
     994
     995    def set_starttime(self, time):
     996        self.starttime = float(time)       
    975997
    976998
     
    10851107
    10861108        while True:
    1087             # Compute fluxes across each element edge
    1088             self.compute_fluxes()
    1089 
    1090             # Update timestep to fit yieldstep and finaltime
    1091             self.update_timestep(yieldstep, finaltime)
    1092 
    1093             # Update conserved quantities
    1094             self.update_conserved_quantities()
    1095 
    1096             # Update ghosts
    1097             self.update_ghosts()
    1098 
    1099             # Update vertex and edge values
    1100             self.distribute_to_vertices_and_edges()
     1109
     1110            # Evolve One Step, using appropriate timestepping method
     1111            if self.get_timestepping_method() == 'euler':
     1112                self.evolve_one_euler_step(yieldstep,finaltime)
     1113               
     1114            elif self.get_timestepping_method() == 'rk2':
     1115                self.evolve_one_rk2_step(yieldstep,finaltime)               
    11011116
    11021117            # Update extrema if necessary (for reporting)
    11031118            self.update_extrema()           
    11041119
    1105             # Update boundary values
    1106             self.update_boundary()
    1107 
    1108             # Update time
    1109             self.time += self.timestep
     1120
    11101121            self.yieldtime += self.timestep
    11111122            self.number_of_steps += 1
     
    11451156
    11461157
     1158    def evolve_one_euler_step(self, yieldstep, finaltime):
     1159        """One Euler Time Step"""
     1160
     1161        #Compute fluxes across each element edge
     1162        self.compute_fluxes()
     1163
     1164        #Update timestep to fit yieldstep and finaltime
     1165        self.update_timestep(yieldstep, finaltime)
     1166
     1167        #Update conserved quantities
     1168        self.update_conserved_quantities()
     1169
     1170        #update ghosts
     1171        self.update_ghosts()
     1172
     1173        #Update vertex and edge values
     1174        self.distribute_to_vertices_and_edges()
     1175
     1176        #Update boundary values
     1177        self.update_boundary()
     1178
     1179        #Update time
     1180        self.time += self.timestep
     1181
     1182       
     1183
     1184
     1185    def evolve_one_rk2_step(self,yieldstep, finaltime):
     1186        """One 2nd order RK timestep"""
     1187
     1188        #Save initial initial conserved quantities values
     1189        self.backup_conserved_quantities()           
     1190
     1191        #--------------------------------------
     1192        #First euler step
     1193        #--------------------------------------
     1194
     1195        #Compute fluxes across each element edge
     1196        self.compute_fluxes()
     1197
     1198        #Update timestep to fit yieldstep and finaltime
     1199        self.update_timestep(yieldstep, finaltime)
     1200
     1201        #Update conserved quantities
     1202        self.update_conserved_quantities()
     1203
     1204        #update ghosts
     1205        self.update_ghosts()
     1206
     1207        #Update vertex and edge values
     1208        self.distribute_to_vertices_and_edges()
     1209
     1210        #Update boundary values
     1211        self.update_boundary()
     1212
     1213        #Update time
     1214        self.time += self.timestep
     1215
     1216        #------------------------------------
     1217        #Second Euler step
     1218        #------------------------------------
     1219           
     1220        #Compute fluxes across each element edge
     1221        self.compute_fluxes()
     1222
     1223        #Update conserved quantities
     1224        self.update_conserved_quantities()
     1225
     1226        #------------------------------------
     1227        #Combine initial and final values
     1228        #of conserved quantities
     1229        #------------------------------------
     1230
     1231        self.saxpy_conserved_quantities(0.5, 0.5)
     1232
     1233        #------------------------------------
     1234        #Clean up rk step
     1235        #------------------------------------
     1236           
     1237        #update ghosts
     1238        self.update_ghosts()
     1239
     1240        #Update vertex and edge values
     1241        self.distribute_to_vertices_and_edges()
     1242
     1243        #Update boundary values
     1244        self.update_boundary()
     1245
     1246
    11471247    def evolve_to_end(self, finaltime = 1.0):
    11481248        """Iterate evolve all the way to the end
     
    11531253
    11541254
     1255    def backup_conserved_quantities(self):
     1256        N = len(self) #number_of_triangles
     1257
     1258        #backup conserved_quantities centroid values
     1259        for name in self.conserved_quantities:
     1260            Q = self.quantities[name]
     1261            Q.backup_centroid_values()       
     1262
     1263    def saxpy_conserved_quantities(self,a,b):
     1264        N = len(self) #number_of_triangles
     1265
     1266        #backup conserved_quantities centroid values
     1267        for name in self.conserved_quantities:
     1268            Q = self.quantities[name]
     1269            Q.saxpy_centroid_values(a,b)       
     1270   
    11551271
    11561272    def update_boundary(self):
  • anuga_core/source/anuga/abstract_2d_finite_volumes/quantity.py

    r4704 r4712  
    13051305        self.explicit_update = zeros(N, Float )
    13061306        self.semi_implicit_update = zeros(N, Float )
     1307        self.centroid_backup_values = zeros(N, Float)       
     1308
     1309       
    13071310
    13081311
     
    13181321        return compute_gradients(self)
    13191322
    1320 
    13211323    def limit(self):
    13221324        #Call correct module function
     
    13241326        limit(self)
    13251327
    1326 
    13271328    def extrapolate_second_order(self):
    13281329        #Call correct module function
     
    13301331        extrapolate_second_order(self)
    13311332
     1333    def backup_centroid_values(self):
     1334        #Call correct module function
     1335        #(either from this module or C-extension)
     1336        backup_centroid_values(self)
     1337
     1338    def saxpy_centroid_values(self,a,b):
     1339        #Call correct module function
     1340        #(either from this module or C-extension)
     1341        saxpy_centroid_values(self,a,b)
     1342   
    13321343
    13331344def update(quantity, timestep):
     
    14171428
    14181429
    1419 def extrapolate_second_order(quantity):
    1420     """Extrapolate conserved quantities from centroid to
    1421     vertices for each volume using
    1422     second order scheme.
    1423     """
    1424 
    1425     a, b = quantity.compute_gradients()
    1426 
    1427     X = quantity.domain.get_vertex_coordinates()
     1430def backup_centroid_values(quantity):
     1431    """Copy centroid values to backup array"""
     1432
    14281433    qc = quantity.centroid_values
    1429     qv = quantity.vertex_values
     1434    qb = quantity.centroid_backup_values
    14301435
    14311436    #Check each triangle
    14321437    for k in range(len(quantity.domain)):
    1433         #Centroid coordinates
    1434         x, y = quantity.domain.centroid_coordinates[k]
    1435 
    1436         #vertex coordinates
    1437         x0, y0, x1, y1, x2, y2 = X[k,:]
    1438 
    1439         #Extrapolate
    1440         qv[k,0] = qc[k] + a[k]*(x0-x) + b[k]*(y0-y)
    1441         qv[k,1] = qc[k] + a[k]*(x1-x) + b[k]*(y1-y)
    1442         qv[k,2] = qc[k] + a[k]*(x2-x) + b[k]*(y2-y)
     1438        qb[k] = qc[k]
     1439
     1440
     1441def saxpy_centroid_values(quantity,a,b):
     1442    """saxpy operation between centroid value and backup"""
     1443
     1444    qc = quantity.centroid_values
     1445    qb = quantity.centroid_backup_values
     1446
     1447    #Check each triangle
     1448    for k in range(len(quantity.domain)):
     1449        qc[k] = a*qc[k]+b*qb[k]       
    14431450
    14441451
     
    15821589    #Replace python version with c implementations
    15831590
    1584     from quantity_ext import average_vertex_values
     1591    from quantity_ext import average_vertex_values, backup_centroid_values, \
     1592         saxpy_centroid_values
    15851593
    15861594    from quantity_ext import compute_gradients, limit,\
  • anuga_core/source/anuga/abstract_2d_finite_volumes/quantity_ext.c

    r4536 r4712  
    100100
    101101
     102
    102103int _extrapolate(int N,
    103104                 double* centroids,
     
    162163        return 0;
    163164}
     165
     166int _backup_centroid_values(int N,
     167                            double* centroid_values,
     168                            double* centroid_backup_values) {
     169    //Backup centroid values
     170
     171
     172    int k;
     173
     174    for (k=0; k<N; k++) {
     175        centroid_backup_values[k] = centroid_values[k];
     176    }
     177
     178
     179    return 0;
     180}
     181
     182
     183int _saxpy_centroid_values(int N,
     184                           double a,
     185                           double b,
     186                           double* centroid_values,
     187                           double* centroid_backup_values) {
     188    //saxby centroid values
     189
     190
     191    int k;
     192
     193
     194    for (k=0; k<N; k++) {
     195        centroid_values[k] = a*centroid_values[k] + b*centroid_backup_values[k];
     196    }
     197
     198
     199    return 0;
     200}
     201
    164202
    165203int _update(int N,
     
    305343
    306344
     345PyObject *backup_centroid_values(PyObject *self, PyObject *args) {
     346
     347        PyObject *quantity;
     348        PyArrayObject *centroid_values, *centroid_backup_values;
     349
     350        int N, err;
     351
     352
     353        // Convert Python arguments to C
     354        if (!PyArg_ParseTuple(args, "O", &quantity)) {
     355          PyErr_SetString(PyExc_RuntimeError,
     356                          "quantity_ext.c: backup_centroid_values could not parse input");
     357          return NULL;
     358        }
     359
     360        centroid_values        = get_consecutive_array(quantity, "centroid_values");
     361        centroid_backup_values = get_consecutive_array(quantity, "centroid_backup_values");
     362
     363        N = centroid_values -> dimensions[0];
     364
     365        err = _backup_centroid_values(N,
     366                      (double*) centroid_values -> data,
     367                      (double*) centroid_backup_values -> data);
     368
     369
     370        //Release and return
     371        Py_DECREF(centroid_values);
     372        Py_DECREF(centroid_backup_values);
     373
     374        return Py_BuildValue("");
     375}
     376
     377PyObject *saxpy_centroid_values(PyObject *self, PyObject *args) {
     378
     379        PyObject *quantity;
     380        PyArrayObject *centroid_values, *centroid_backup_values;
     381
     382        double a,b;
     383        int N, err;
     384
     385
     386        // Convert Python arguments to C
     387        if (!PyArg_ParseTuple(args, "Odd", &quantity, &a, &b)) {
     388          PyErr_SetString(PyExc_RuntimeError,
     389                          "quantity_ext.c: saxpy_centroid_values could not parse input");
     390          return NULL;
     391        }
     392
     393        centroid_values        = get_consecutive_array(quantity, "centroid_values");
     394        centroid_backup_values = get_consecutive_array(quantity, "centroid_backup_values");
     395
     396        N = centroid_values -> dimensions[0];
     397
     398        err = _saxpy_centroid_values(N,a,b,
     399                      (double*) centroid_values -> data,
     400                      (double*) centroid_backup_values -> data);
     401
     402
     403        //Release and return
     404        Py_DECREF(centroid_values);
     405        Py_DECREF(centroid_backup_values);
     406
     407        return Py_BuildValue("");
     408}
     409
     410
    307411PyObject *interpolate_from_vertices_to_edges(PyObject *self, PyObject *args) {
    308412
     
    645749        {"limit", limit, METH_VARARGS, "Print out"},
    646750        {"update", update, METH_VARARGS, "Print out"},
     751        {"backup_centroid_values", backup_centroid_values, METH_VARARGS, "Print out"},
     752        {"saxpy_centroid_values", saxpy_centroid_values, METH_VARARGS, "Print out"},
    647753        {"compute_gradients", compute_gradients, METH_VARARGS, "Print out"},
    648754        {"extrapolate_second_order", extrapolate_second_order,
  • anuga_core/source/anuga/abstract_2d_finite_volumes/test_domain.py

    r4704 r4712  
    319319
    320320
    321 
    322 
    323 
     321    def test_setting_timestepping_method(self):
     322        """test_set_quanitities_to_be_monitored
     323        """
     324
     325        a = [0.0, 0.0]
     326        b = [0.0, 2.0]
     327        c = [2.0,0.0]
     328        d = [0.0, 4.0]
     329        e = [2.0, 2.0]
     330        f = [4.0,0.0]
     331
     332        points = [a, b, c, d, e, f]
     333        #bac, bce, ecf, dbe, daf, dae
     334        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
     335
     336
     337        domain = Domain(points, vertices, boundary=None,
     338                        conserved_quantities =\
     339                        ['stage', 'xmomentum', 'ymomentum'],
     340                        other_quantities = ['elevation', 'friction', 'depth'])
     341
     342
     343        domain.timestepping_method = None
     344
     345
     346        # Check that invalid requests are dealt with
     347        try:
     348            domain.set_timestepping_method('eee')       
     349        except:
     350            pass
     351        else:
     352            msg = 'Should have caught illegal method'
     353            raise Exception, msg
     354
     355
     356        #Should have no trouble with euler or rk2
     357        domain.set_timestepping_method('euler')
     358        domain.set_timestepping_method('rk2')
     359
     360        #test get timestepping method
     361        assert domain.get_timestepping_method() == 'rk2'
    324362
    325363
  • anuga_core/source/anuga/abstract_2d_finite_volumes/test_quantity.py

    r4704 r4712  
    13681368
    13691369
     1370    def test_backup_saxpy_centroid_values(self):
     1371        quantity = Conserved_quantity(self.mesh4)
     1372
     1373        #Set up for a gradient of (3,1), f(x) = 3x+y
     1374        c_values = array([2.0+2.0/3, 4.0+4.0/3, 8.0+2.0/3, 2.0+8.0/3])
     1375        d_values = array([1.0, 2.0, 3.0, 4.0])
     1376        quantity.set_values(c_values, location = 'centroids')
     1377
     1378        #Backup
     1379        quantity.backup_centroid_values()
     1380
     1381        #print quantity.vertex_values
     1382        assert allclose(quantity.centroid_values, quantity.centroid_backup_values)
     1383
     1384
     1385        quantity.set_values(d_values, location = 'centroids')
     1386
     1387        quantity.saxpy_centroid_values(2.0, 3.0)
     1388
     1389        assert(quantity.centroid_values, 2.0*d_values + 3.0*c_values)
     1390
     1391
     1392
    13701393    def test_first_order_extrapolator(self):
    13711394        quantity = Conserved_quantity(self.mesh4)
  • anuga_core/source/anuga/abstract_2d_finite_volumes/util.py

    r4704 r4712  
    8585            quantities = domain.conserved_quantities
    8686           
    87         domain_starttime = domain.starttime
     87        domain_starttime = domain.get_starttime()
    8888    else:
    8989        domain_starttime = None
     
    140140           
    141141            if verbose: print msg
    142             domain.starttime = starttime #Modifying model time
     142            domain.set_starttime(starttime) #Modifying model time
    143143            if verbose: print 'Domain starttime is now set to %f'\
    144144               %domain.starttime
  • anuga_core/source/anuga/config.py

    r4699 r4712  
    104104           #make changes
    105105
     106# Choose type of timestepping,
     107timestepping_method = 'euler' # 1st order euler
     108#timestepping_method = 'rk2'   # 2nd Order TVD scheme
     109
    106110# Option to search for signatures where isolated triangles are
    107111# responsible for a small global timestep.
  • anuga_core/source/anuga/shallow_water/shallow_water_domain.py

    r4710 r4712  
    446446
    447447
    448         #Initial update of vertex and edge values before any storage
     448        #Initial update of vertex and edge values before any STORAGE
    449449        #and or visualisation
     450        #This is done again in the initialisation of the Generic_Domain
     451        #evolve loop but we do it here to ensure the values are ok for storage
    450452        self.distribute_to_vertices_and_edges()
    451453
  • anuga_core/source/anuga/shallow_water/test_system.py

    r4561 r4712  
    4747        domain.set_name(boundary_name)                 
    4848        domain.set_datadir(dir)         
    49         domain.starttime = boundary_starttime
     49        domain.set_starttime(boundary_starttime)
    5050       
    5151        # Setup initial conditions
     
    7171       
    7272        """
     73     
    7374        boundary_starttime = 500
    7475        boundary_filename = self.create_sww_boundary(boundary_starttime)
     
    9596        # Setup boundary conditions
    9697        domain.set_boundary({'exterior': Bf})
    97         for t in domain.evolve(yieldstep = 5, finaltime = 10.0):
     98
     99       
     100        for t in domain.evolve(yieldstep = 5.0, finaltime = 10.0):
    98101            pass
    99             #domain.write_time()
     102            #print domain.write_time()
    100103            #print "domain.time", domain.time
    101104
     
    117120        Test that starttime can be set in the middle of a boundary condition
    118121        """
     122       
    119123        boundary_starttime = 500
    120124        boundary_filename = self.create_sww_boundary(boundary_starttime)
     
    134138        domain.set_datadir(dir)
    135139        new_starttime = 510.
    136         domain.starttime = new_starttime
     140        domain.set_starttime(new_starttime)
    137141
    138142        # Setup initial conditions
     
    146150        for t in domain.evolve(yieldstep = 5, finaltime = 9.0):
    147151            pass
    148             #domain.write_time()
     152            #print domain.write_time()
    149153            #print "domain.time", domain.time
    150154
  • anuga_core/source/anuga/utilities/compile.py

    r4493 r4712  
    201201    # This will make compile work locally
    202202    utilities_include_dir = '.'
     203    utilities_include_dir = '../utilities'
    203204
    204205
Note: See TracChangeset for help on using the changeset viewer.