1 | """Example of shallow water wave equation. |
---|

2 | |
---|

3 | Two levels - testing that momentum is not genearted in the upward direction |
---|

4 | |
---|

5 | """ |
---|

6 | |
---|

7 | ###################### |
---|

8 | # Module imports |
---|

9 | # |
---|

10 | from os import sep, path |
---|

11 | from mesh_factory import rectangular |
---|

12 | from shallow_water import Domain, Reflective_boundary, Dirichlet_boundary,\ |
---|

13 | Constant_height |
---|

14 | from Numeric import array |
---|

15 | from anuga.pyvolution.util import Polygon_function, read_polygon |
---|

16 | |
---|

17 | |
---|

18 | #Create basic mesh |
---|

19 | N = 2 |
---|

20 | points, vertices, boundary = rectangular(3, N, 20, 10) |
---|

21 | |
---|

22 | #Create shallow water domain |
---|

23 | domain = Domain(points, vertices, boundary) |
---|

24 | domain.store = True |
---|

25 | domain.set_name('two_levels') |
---|

26 | print "Output being written to " + domain.get_datadir() + sep + \ |
---|

27 | domain.filename + "_size%d." %len(domain) + domain.format |
---|

28 | |
---|

29 | |
---|

30 | domain.default_order=2 |
---|

31 | domain.smooth = False |
---|

32 | domain.visualise = True |
---|

33 | |
---|

34 | #PLAY WITH THIS [0;1]: |
---|

35 | # |
---|

36 | # beta_h == 0.0 reveals the problem |
---|

37 | # beta_h > 0.2 alleviates it |
---|

38 | domain.beta_h = 0.2 |
---|

39 | |
---|

40 | #IC |
---|

41 | domain.set_quantity('friction', 0.0) |
---|

42 | domain.set_quantity('stage', 3.0) |
---|

43 | |
---|

44 | def twostage(x, y): |
---|

45 | z = 0*x |
---|

46 | for i in range(len(x)): |
---|

47 | if x[i]<10: |
---|

48 | z[i] = 4. |
---|

49 | |
---|

50 | return z |
---|

51 | #Set elevation |
---|

52 | domain.set_quantity('elevation', twostage) |
---|

53 | |
---|

54 | |
---|

55 | ###################### |
---|

56 | # Boundary conditions |
---|

57 | Br = Reflective_boundary(domain) |
---|

58 | |
---|

59 | domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br}) |
---|

60 | domain.check_integrity() |
---|

61 | |
---|

62 | |
---|

63 | #print domain.quantities['elevation'].vertex_values |
---|

64 | ###################### |
---|

65 | #Evolution |
---|

66 | for t in domain.evolve(yieldstep = 0.01, finaltime = 0.05): |
---|

67 | domain.write_time() |
---|

68 | print domain.quantities['xmomentum'].edge_values |
---|

69 | #print domain.quantities['stage'].centroid_values |
---|

70 | |
---|

71 | #print domain.quantities['stage'].vertex_values |
---|

72 | #print |
---|

73 | #print domain.quantities['xmomentum'].vertex_values |
---|

74 | #print |
---|

75 | #print domain.quantities['ymomentum'].vertex_values |
---|

76 | |
---|

77 | print 'Done' |
---|

78 | |
---|