[1089] | 1 | """Simple example of shallow water wave equation using Pyvolution |
---|
| 2 | |
---|
| 3 | Water driven by linear slope and Dirichlet boundary |
---|
| 4 | |
---|
| 5 | """ |
---|
| 6 | |
---|
| 7 | ###################### |
---|
| 8 | # Module imports |
---|
| 9 | # |
---|
| 10 | from mesh_factory import rectangular |
---|
| 11 | from shallow_water import Domain, Reflective_boundary, Dirichlet_boundary,\ |
---|
| 12 | Constant_height, Time_boundary, Transmissive_boundary |
---|
| 13 | from Numeric import array |
---|
| 14 | |
---|
| 15 | #Create basic mesh |
---|
| 16 | points, vertices, boundary = rectangular(10, 10) |
---|
| 17 | |
---|
| 18 | #Create shallow water domain |
---|
| 19 | domain = Domain(points, vertices, boundary) |
---|
| 20 | domain.smooth = False |
---|
| 21 | domain.visualise = False |
---|
| 22 | domain.store = True |
---|
| 23 | domain.filename = 'bedslope' |
---|
| 24 | domain.default_order=2 |
---|
[1133] | 25 | domain.quantities_to_be_stored=['stage'] |
---|
[1089] | 26 | |
---|
| 27 | ####################### |
---|
| 28 | #Bed-slope and friction |
---|
| 29 | domain.set_quantity('elevation', lambda x,y: -x/3) |
---|
| 30 | domain.set_quantity('friction', 0.1) |
---|
| 31 | |
---|
| 32 | ###################### |
---|
| 33 | # Boundary conditions |
---|
| 34 | from math import sin, pi |
---|
| 35 | Br = Reflective_boundary(domain) |
---|
| 36 | Bt = Transmissive_boundary(domain) |
---|
| 37 | Bd = Dirichlet_boundary([0.2,0.,0.]) |
---|
| 38 | Bw = Time_boundary(domain=domain, |
---|
| 39 | f=lambda t: [(0.1*sin(t*2*pi)), 0.0, 0.0]) |
---|
| 40 | |
---|
| 41 | domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br}) |
---|
| 42 | |
---|
| 43 | ###################### |
---|
| 44 | #Initial condition |
---|
| 45 | h = 0.05 |
---|
| 46 | elevation = domain.quantities['elevation'].vertex_values |
---|
| 47 | domain.set_quantity('stage', elevation + h) |
---|
| 48 | #elevation = domain.get_quantity('elevation',location='unique vertices') |
---|
| 49 | #domain.set_quantity('stage', elevation + h,location='unique vertices') |
---|
| 50 | |
---|
| 51 | domain.check_integrity() |
---|
| 52 | ###################### |
---|
| 53 | #Evolution |
---|
| 54 | for t in domain.evolve(yieldstep = 1, finaltime = 2.0): |
---|
[1133] | 55 | pass |
---|
| 56 | #domain.write_time() |
---|
[1089] | 57 | |
---|
| 58 | |
---|
| 59 | ##NOW TEST IT!!! |
---|
| 60 | |
---|
| 61 | from data_manager import sww2domain |
---|
| 62 | from Numeric import allclose |
---|
| 63 | |
---|
[1133] | 64 | filename = domain.datadir+'\\'+domain.filename+'.sww' |
---|
[1089] | 65 | |
---|
| 66 | try: |
---|
[1133] | 67 | domain2 = sww2domain(filename,verbose=False) |
---|
[1089] | 68 | assert True == False |
---|
| 69 | except: |
---|
[1133] | 70 | filler = 0 |
---|
| 71 | domain2 = sww2domain(filename,fail_if_NaN=False,NaN_filler = filler,verbose=False) |
---|
[1089] | 72 | |
---|
| 73 | bits = ['xllcorner','yllcorner','vertex_coordinates','time','starttime'] |
---|
| 74 | |
---|
| 75 | for quantity in ['elevation']+domain.quantities_to_be_stored: |
---|
| 76 | bits.append('get_quantity("%s")'%quantity) |
---|
| 77 | |
---|
| 78 | for bit in bits: |
---|
[1133] | 79 | # print 'testing that domain.'+bit+' has been restored' |
---|
[1089] | 80 | assert allclose(eval('domain.'+bit),eval('domain2.'+bit)) |
---|
| 81 | |
---|
[1133] | 82 | #print max(max(domain2.get_quantity('xmomentum'))) |
---|
| 83 | #print min(min(domain2.get_quantity('xmomentum'))) |
---|
| 84 | #print max(max(domain2.get_quantity('ymomentum'))) |
---|
| 85 | #print min(min(domain2.get_quantity('ymomentum'))) |
---|
[1090] | 86 | |
---|
[1133] | 87 | assert max(max(domain2.get_quantity('xmomentum')))==filler |
---|
| 88 | assert min(min(domain2.get_quantity('xmomentum')))==filler |
---|
| 89 | assert max(max(domain2.get_quantity('ymomentum')))==filler |
---|
| 90 | assert min(min(domain2.get_quantity('ymomentum')))==filler |
---|
[1090] | 91 | |
---|
[1133] | 92 | #print 'passed' |
---|
| 93 | |
---|
| 94 | #cleanup |
---|
| 95 | #import os |
---|
| 96 | #os.remove(domain.datadir+'/'+domain.filename+'.sww') |
---|