[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 cairns_project |
---|
| 15 | from caching import cache |
---|
| 16 | |
---|
| 17 | |
---|
| 18 | |
---|
| 19 | |
---|
| 20 | #------------------------------- |
---|
| 21 | # Domain |
---|
| 22 | #------------------------------- |
---|
| 23 | #print 'Creating domain from', cairns_project.mesh_filename |
---|
| 24 | |
---|
| 25 | #domain = cache(pmesh_to_domain_instance, |
---|
| 26 | # (cairns_project.mesh_filename, Domain), |
---|
| 27 | # dependencies = [cairns_project.mesh_filename]) |
---|
| 28 | |
---|
| 29 | print 'Creating domain' |
---|
| 30 | N = 50 |
---|
| 31 | M = 60 |
---|
| 32 | #r_earth = 5e6 |
---|
| 33 | #from math import sin, pi |
---|
| 34 | #scale_y = r_earth*sin(74/180*pi) |
---|
| 35 | #print scale_y |
---|
| 36 | points, elements, boundary = rectangular_cross(N, M, len1=30.0, len2=16.0, |
---|
| 37 | origin = (145.0, -24.0)) |
---|
| 38 | domain = Domain(points, elements, boundary) |
---|
| 39 | |
---|
| 40 | domain.check_integrity() |
---|
| 41 | print 'Number of triangles = ', len(domain) |
---|
| 42 | print 'The extent is ', domain.get_extent() |
---|
| 43 | |
---|
| 44 | |
---|
| 45 | |
---|
| 46 | #------------------------------- |
---|
| 47 | # Initial Conditions |
---|
| 48 | #------------------------------- |
---|
| 49 | print 'Initial values' |
---|
| 50 | |
---|
| 51 | domain.set_quantity('elevation', |
---|
| 52 | filename = cairns_project.bathymetry_filename[:-4] + '.xya', |
---|
| 53 | alpha = 10.0, |
---|
| 54 | verbose = True, |
---|
| 55 | use_cache = True) |
---|
| 56 | |
---|
| 57 | domain.set_quantity('friction', 0.0) |
---|
| 58 | domain.set_quantity('stage', 0.0) |
---|
| 59 | |
---|
| 60 | #------------------------------- |
---|
| 61 | # Boundary conditions |
---|
| 62 | #------------------------------- |
---|
| 63 | print 'Boundaries' |
---|
| 64 | |
---|
| 65 | # Tidal cycle recorded at Eden as open |
---|
| 66 | #print 'Open sea boundary condition from ',cairns_project.boundary_filename |
---|
| 67 | #from pyvolution.util import file_function |
---|
| 68 | #tide_function = file_function(cairns_project.boundary_filename[:-4] + '.tms', domain, |
---|
| 69 | # verbose = True) |
---|
| 70 | #Bts = Transmissive_Momentum_Set_Stage_boundary(domain, tide_function) |
---|
| 71 | |
---|
| 72 | # All other boundaries are reflective |
---|
| 73 | Br = Reflective_boundary(domain) |
---|
| 74 | |
---|
| 75 | #domain.set_boundary({'exterior': Br, 'open': Br}) |
---|
| 76 | domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br}) |
---|
| 77 | |
---|
| 78 | #------------------------------- |
---|
| 79 | # Setup domain runtime parameters |
---|
| 80 | #------------------------------- |
---|
| 81 | domain.visualise = True |
---|
| 82 | domain.visualise_color_stage = True |
---|
| 83 | |
---|
| 84 | base = os.path.basename(sys.argv[0]) |
---|
| 85 | domain.filename, _ = os.path.splitext(base) |
---|
| 86 | domain.default_order = 2 |
---|
| 87 | domain.store = True #Store for visualisation purposes |
---|
| 88 | domain.smooth = False |
---|
| 89 | |
---|
| 90 | |
---|
| 91 | |
---|
| 92 | #------------------------------- |
---|
| 93 | # Evolve |
---|
| 94 | #------------------------------- |
---|
| 95 | import time |
---|
| 96 | t0 = time.time() |
---|
| 97 | yieldstep = 0.1 |
---|
| 98 | finaltime = 1 |
---|
| 99 | |
---|
| 100 | |
---|
| 101 | for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime): |
---|
| 102 | domain.write_time() |
---|
| 103 | |
---|
| 104 | print 'That took %.2f seconds' %(time.time()-t0) |
---|