[2225] | 1 | """ |
---|
| 2 | Main meribula script using new interface |
---|
| 3 | """ |
---|
| 4 | |
---|
| 5 | #------------------------------- |
---|
| 6 | # Module imports |
---|
| 7 | #------------------------------- |
---|
| 8 | import sys, os |
---|
| 9 | from pyvolution.shallow_water import Domain, Reflective_boundary,\ |
---|
| 10 | File_boundary, Transmissive_Momentum_Set_Stage_boundary |
---|
| 11 | from pyvolution.mesh_factory import rectangular_cross |
---|
| 12 | from pyvolution.pmesh2domain import pmesh_to_domain_instance |
---|
| 13 | from Numeric import array, zeros, Float, allclose |
---|
| 14 | import project |
---|
| 15 | from caching import cache |
---|
| 16 | |
---|
| 17 | |
---|
| 18 | |
---|
| 19 | |
---|
| 20 | #------------------------------- |
---|
| 21 | # Domain |
---|
| 22 | #------------------------------- |
---|
| 23 | print 'Creating domain from', project.mesh_filename |
---|
| 24 | |
---|
| 25 | domain = cache(pmesh_to_domain_instance, |
---|
| 26 | (project.mesh_filename, Domain), |
---|
| 27 | dependencies = [project.mesh_filename]) |
---|
| 28 | |
---|
| 29 | #domain = pmesh_to_domain_instance(project.mesh_filename, Domain) |
---|
| 30 | |
---|
| 31 | domain.check_integrity() |
---|
| 32 | print 'Number of triangles = ', len(domain) |
---|
| 33 | print 'The extent is ', domain.get_extent() |
---|
| 34 | |
---|
| 35 | |
---|
| 36 | |
---|
| 37 | #------------------------------- |
---|
| 38 | # Initial Conditions |
---|
| 39 | #------------------------------- |
---|
| 40 | print 'Initial values' |
---|
| 41 | |
---|
| 42 | domain.set_quantity('elevation', |
---|
| 43 | filename = project.bathymetry_filename[:-4] + '.xya', |
---|
| 44 | alpha = 10.0, |
---|
| 45 | verbose = True, |
---|
| 46 | use_cache = True) |
---|
| 47 | |
---|
| 48 | domain.set_quantity('friction', 0.03) |
---|
| 49 | domain.set_quantity('stage', 0.0) |
---|
| 50 | |
---|
| 51 | #------------------------------- |
---|
| 52 | # Boundary conditions |
---|
| 53 | #------------------------------- |
---|
| 54 | print 'Boundaries' |
---|
| 55 | |
---|
| 56 | # Tidal cycle recorded at Eden as open |
---|
| 57 | print 'Open sea boundary condition from ',project.boundary_filename |
---|
| 58 | from pyvolution.util import file_function |
---|
| 59 | tide_function = file_function(project.boundary_filename[:-4] + '.tms', domain, |
---|
| 60 | verbose = True) |
---|
| 61 | Bts = Transmissive_Momentum_Set_Stage_boundary(domain, tide_function) |
---|
| 62 | |
---|
| 63 | # All other boundaries are reflective |
---|
| 64 | Br = Reflective_boundary(domain) |
---|
| 65 | |
---|
| 66 | domain.set_boundary({'exterior': Br, 'open': Bts}) |
---|
| 67 | |
---|
| 68 | #------------------------------- |
---|
| 69 | # Setup domain runtime parameters |
---|
| 70 | #------------------------------- |
---|
| 71 | domain.visualise = True |
---|
| 72 | domain.visualise_color_stage = True |
---|
| 73 | |
---|
| 74 | base = os.path.basename(sys.argv[0]) |
---|
| 75 | domain.filename, _ = os.path.splitext(base) |
---|
| 76 | domain.default_order = 2 |
---|
| 77 | domain.store = True #Store for visualisation purposes |
---|
| 78 | |
---|
| 79 | |
---|
| 80 | |
---|
| 81 | #------------------------------- |
---|
| 82 | # Evolve |
---|
| 83 | #------------------------------- |
---|
| 84 | import time |
---|
| 85 | t0 = time.time() |
---|
| 86 | yieldstep = 60 |
---|
| 87 | finaltime = 3600*24 |
---|
| 88 | |
---|
| 89 | |
---|
| 90 | for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime): |
---|
| 91 | domain.write_time() |
---|
| 92 | |
---|
| 93 | print 'That took %.2f seconds' %(time.time()-t0) |
---|