[4693] | 1 | """Simple water flow example using ANUGA |
---|
| 2 | |
---|
| 3 | Water driven up a linear slope and time varying boundary, |
---|
| 4 | similar to a beach environment |
---|
| 5 | """ |
---|
| 6 | |
---|
| 7 | #------------------------------------------------------------------------------ |
---|
| 8 | # Import necessary modules |
---|
| 9 | #------------------------------------------------------------------------------ |
---|
| 10 | |
---|
| 11 | import sys |
---|
| 12 | from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross |
---|
| 13 | from anuga.shallow_water import Domain |
---|
| 14 | from anuga.shallow_water import Reflective_boundary |
---|
| 15 | from anuga.shallow_water import Dirichlet_boundary |
---|
| 16 | from anuga.shallow_water import Time_boundary |
---|
| 17 | from anuga.shallow_water import Transmissive_boundary |
---|
| 18 | from anuga.shallow_water import Transmissive_Momentum_Set_Stage_boundary |
---|
| 19 | from anuga.geospatial_data.geospatial_data import * |
---|
| 20 | from math import cos |
---|
| 21 | |
---|
| 22 | #------------------------------------------------------------------------------ |
---|
| 23 | # Setup computational domain |
---|
| 24 | #------------------------------------------------------------------------------ |
---|
[4703] | 25 | dx = 10. |
---|
[4694] | 26 | dy = dx |
---|
| 27 | L = 100000. |
---|
| 28 | W = 3000. |
---|
| 29 | |
---|
[4703] | 30 | # structured mesh |
---|
| 31 | #points, vertices, boundary = rectangular_cross(int(L/dx), int(W/dy), |
---|
| 32 | # L, W, (0.0, -W/2)) # Basic mesh |
---|
[4694] | 33 | #points, vertices, boundary = rectangular_cross(666, 3, 100000, 3000, (0.0, -0.0)) # Basic mesh |
---|
[4693] | 34 | #points, vertices, boundary = rectangular_cross(530, 10, 5300, 100, (-5000.0, -50.0)) # Basic mesh |
---|
| 35 | #points, vertices, boundary = rectangular_cross(1000, 100, 20, 3) # Basic mesh |
---|
[4703] | 36 | #domain = Domain(points, vertices, boundary) |
---|
[4693] | 37 | |
---|
[4703] | 38 | # unstructured mesh |
---|
| 39 | poly_domain = [[0,-W],[0,W],[L,W],[L,-W]] |
---|
| 40 | meshname = 'test.msh' |
---|
| 41 | from anuga.pmesh.mesh_interface import create_mesh_from_regions |
---|
| 42 | # Create mesh |
---|
| 43 | create_mesh_from_regions(poly_domain, |
---|
| 44 | boundary_tags={'left': [0], 'top': [1], |
---|
| 45 | 'right': [2], 'bottom': [3]}, |
---|
[4715] | 46 | maximum_triangle_area = 100000, |
---|
[4703] | 47 | filename=meshname) |
---|
| 48 | |
---|
| 49 | # Create domain |
---|
| 50 | domain = Domain(meshname, use_cache=True, verbose = True) |
---|
| 51 | domain.set_name('myexample2') |
---|
[4715] | 52 | domain.set_default_order(2) # Use second order spatial scheme |
---|
[4693] | 53 | domain.set_datadir('.') # Use current directory for output |
---|
| 54 | |
---|
| 55 | #------------------------------------------------------------------------------ |
---|
| 56 | # Setup initial conditions |
---|
| 57 | #------------------------------------------------------------------------------ |
---|
| 58 | #domain.set_quantity('elevation', topography) # Use function for elevation |
---|
[4715] | 59 | domain.set_quantity('elevation', -100) |
---|
[4703] | 60 | domain.set_quantity('friction', 0.00) |
---|
| 61 | domain.set_quantity('stage', 0.0) |
---|
[4693] | 62 | |
---|
| 63 | #----------------------------------------------------------------------------- |
---|
| 64 | # Setup boundary conditions |
---|
| 65 | #------------------------------------------------------------------------------ |
---|
| 66 | from math import sin, pi, exp |
---|
| 67 | Br = Reflective_boundary(domain) # Solid reflective wall |
---|
| 68 | Bt = Transmissive_boundary(domain) # Continue all values on boundary |
---|
| 69 | Bd = Dirichlet_boundary([1,0.,0.]) # Constant boundary values |
---|
[4703] | 70 | amplitude = 1 |
---|
[4693] | 71 | Bw = Time_boundary(domain=domain, # Time dependent boundary |
---|
| 72 | ## Sine wave |
---|
| 73 | f=lambda t: [(-amplitude*sin((1./300.)*t*2*pi)), 0.0, 0.0]) |
---|
| 74 | ## Sawtooth? |
---|
| 75 | # f=lambda t: [(-8.0*(sin((1./180.)*t*2*pi))+(1./2.)*sin((2./180.)*t*2*pi)+(1./3.)*sin((3./180.)*t*2*pi)), 0.0, 0.0]) |
---|
| 76 | ## Sharp rise, linear fall |
---|
| 77 | # f=lambda t: [(5.0*(-((t-0.)/300.)*(t<300.)-cos((t-300.)*2.*pi*(1./240.))*(t>=300. and t<420.)+(1.-(t-420.)/300.)*(t>=420. and t <720.))), 0.0, 0.0]) |
---|
| 78 | # f=lambda t: [amplitude*(1.-2.*(pi*(1./720.)*(t-720.))**2)/exp((pi*(1./720.)*(t-720.))**2) , 0.0, 0.0]) |
---|
| 79 | # f=lambda t: [(-8.0*sin((1./720.)*t*2*pi))*((t<720.)-0.5*(t<360.)), 0.0, 0.0]) |
---|
| 80 | |
---|
| 81 | # Associate boundary tags with boundary objects |
---|
| 82 | domain.set_boundary({'left': Bw, 'right': Bt, 'top': Br, 'bottom': Br}) |
---|
| 83 | |
---|
| 84 | #------------------------------------------------------------------------------ |
---|
| 85 | # Evolve system through time |
---|
| 86 | #------------------------------------------------------------------------------ |
---|
| 87 | |
---|
[4709] | 88 | for t in domain.evolve(yieldstep = 20.0, finaltime = 40*60.): |
---|
[4693] | 89 | domain.write_time() |
---|
| 90 | |
---|