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

2 | |
---|

3 | Island surrounded by water |
---|

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 util import Polygon_function, read_polygon |
---|

16 | from math import exp |
---|

17 | |
---|

18 | |
---|

19 | #Create basic mesh |
---|

20 | N = 20 |
---|

21 | points, vertices, boundary = rectangular(N, N, 100, 100) |
---|

22 | |
---|

23 | #Create shallow water domain |
---|

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

25 | domain.store = True |
---|

26 | domain.smooth = False |
---|

27 | domain.visualise = False |
---|

28 | domain.set_name('island') |
---|

29 | print 'Output being written to ' + domain.get_datadir() + sep + \ |
---|

30 | domain.get_name() + '.' + domain.format |
---|

31 | |
---|

32 | |
---|

33 | domain.default_order=2 |
---|

34 | |
---|

35 | |
---|

36 | #I tried to introduce this parameter top control the h-limiter, |
---|

37 | #but it doesn't remove the 'lapping water' |
---|

38 | #NEW (ON): I think it has been fixed (or at least reduced significantly) now |
---|

39 | # |
---|

40 | # beta_h == 1.0 means that the largest gradients (on h) are allowed |
---|

41 | # beta_h == 0.0 means that constant (1st order) gradients are introduced |
---|

42 | # on h. This is equivalent to the constant depth used previously. |
---|

43 | domain.beta_h = 0.2 |
---|

44 | |
---|

45 | |
---|

46 | #Initial conditions |
---|

47 | |
---|

48 | def island(x, y): |
---|

49 | z = 0*x |
---|

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

51 | z[i] = 8*exp( -((x[i]-50)**2 + (y[i]-50)**2)/100 ) |
---|

52 | return z |
---|

53 | |
---|

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

55 | domain.set_quantity('elevation', island) |
---|

56 | domain.set_quantity('stage', 1) |
---|

57 | |
---|

58 | |
---|

59 | |
---|

60 | ###################### |
---|

61 | # Boundary conditions |
---|

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

63 | |
---|

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

65 | domain.check_integrity() |
---|

66 | |
---|

67 | |
---|

68 | ###################### |
---|

69 | #Evolution |
---|

70 | import time |
---|

71 | for t in domain.evolve(yieldstep = 1, finaltime = 20): |
---|

72 | domain.write_time() |
---|

73 | |
---|

74 | print 'Done' |
---|

75 | |
---|