thydrology_steady_test.py - pism - [fork] customized build of PISM, the parallel ice sheet model (tillflux branch)
(HTM) git clone git://src.adamsgaard.dk/pism
(DIR) Log
(DIR) Files
(DIR) Refs
(DIR) LICENSE
---
thydrology_steady_test.py (5057B)
---
1 #!/usr/bin/env python3
2 from unittest import TestCase
3 import numpy as np
4 import os
5
6 import PISM
7 ctx = PISM.Context()
8 ctx.log.set_threshold(1)
9
10 class SteadyHydrology(TestCase):
11 def setUp(self):
12 # domain size
13 C = 1000.0
14 x0 = 0.0
15 y0 = 0.0
16 Lx = C*0.5
17 Ly = C*0.5
18 # grid size
19 Mx = 101
20 My = 101
21
22 grid = PISM.IceGrid.Shallow(ctx.ctx, Lx, Ly, x0, y0, Mx, My,
23 PISM.CELL_CENTER, PISM.NOT_PERIODIC)
24 self.grid = grid
25
26 geometry = PISM.Geometry(grid)
27 self.geometry = geometry
28
29 with PISM.vec.Access(nocomm=[geometry.bed_elevation, geometry.ice_thickness]):
30 for (i, j) in grid.points():
31 x = grid.x(i)
32 y = grid.y(j)
33 geometry.bed_elevation[i, j] = (x / C)**2 + y / C + 0.25
34
35 if geometry.bed_elevation[i, j] >= 0:
36 geometry.ice_thickness[i, j] = 1000.0
37 geometry.sea_level_elevation.set(0.0)
38 geometry.ensure_consistency(0.0)
39
40 surface_input_rate = PISM.IceModelVec2S(grid, "water_input_rate", PISM.WITHOUT_GHOSTS)
41 surface_input_rate.set_attrs("", "water input rate", "kg m-2 s-1", "kg m-2 s-1", "", 0)
42 self.surface_input_rate = surface_input_rate
43
44 self.surface_input_file = "hydrology_steady_surface_input.nc"
45 surface_input_rate.dump(self.surface_input_file)
46
47 # center of the patch of non-zero input
48 cx = 0.0
49 cy = 0.25 * C
50 # size of the patch
51 R = 0.125 * C
52 with PISM.vec.Access(nocomm=surface_input_rate):
53 for (i, j) in grid.points():
54 x = grid.x(i)
55 y = grid.y(j)
56 if abs(x - cx) < R and abs(y - cy) < R:
57 surface_input_rate[i, j] = 1000.0
58 else:
59 surface_input_rate[i, j] = 0
60
61 ctx.config.set_flag("hydrology.add_water_input_to_till_storage", False)
62 ctx.config.set_string("hydrology.surface_input.file", self.surface_input_file)
63
64 self.model = PISM.SteadyState(grid)
65
66 zero = PISM.IceModelVec2S(grid, "zero", PISM.WITHOUT_GHOSTS)
67 zero.set(0.0)
68 self.zero = zero
69
70 self.model.init(zero, zero, zero)
71
72 def tearDown(self):
73 os.remove(self.surface_input_file)
74
75 def divergence_theorem_test(self):
76 "Test that the total input equals the total flux through the boundary."
77 grid = self.grid
78
79 inputs = PISM.HydrologyInputs()
80
81 inputs.no_model_mask = None
82 inputs.geometry = self.geometry
83 inputs.basal_melt_rate = self.zero
84 inputs.ice_sliding_speed = self.zero
85 inputs.surface_input_rate = self.surface_input_rate
86
87 dt = self.model.max_timestep(0).value()
88 self.model.update(0, dt, inputs)
89
90 flux_magnitude = PISM.IceModelVec2S(grid, "flux_magnitude", PISM.WITHOUT_GHOSTS)
91
92 flux_magnitude.set_to_magnitude(self.model.flux())
93
94 # Compute the total input. This approximates a double integral, hence
95 # the "dx dy" factor.
96 total_input = self.model.surface_input_rate().sum() * (grid.dx() * grid.dy())
97
98 # Compute the total flux through the grounding line. This is not
99 # exactly what we want, but it's close. It would be better to use the
100 # flux on the staggered grid as computed internally, but the flux
101 # magnitude will do, especially in the dx=dy case. This approximates a
102 # line integral over the grounding line, hence the dx factor below.
103 total_flux = 0.0
104 cell_type = self.geometry.cell_type
105 with PISM.vec.Access(nocomm=[cell_type, flux_magnitude]):
106 for (i, j) in grid.points():
107 if cell_type.ice_free_ocean(i, j) and cell_type.next_to_grounded_ice(i, j):
108 total_flux += flux_magnitude[i, j] * grid.dx()
109
110 total_flux = PISM.GlobalSum(ctx.com, total_flux)
111
112 # This is the relative error. Note that it is not sensitive to the
113 # value of hydrology.steady.volume_ratio.
114 relative_error = np.fabs(total_input - total_flux) / total_input
115
116 assert relative_error < 1e-5
117 ctx.log.message(1, "relative error: {}\n".format(relative_error))
118
119 def write_results(self):
120 geometry = self.geometry
121 model = self.model
122
123 output_file = ctx.config.get_string("output.file_name")
124
125 f = PISM.util.prepare_output(output_file)
126
127 geometry.bed_elevation.write(f)
128 geometry.cell_type.write(f)
129 geometry.ice_thickness.write(f)
130
131 model.surface_input_rate().write(f)
132 model.flux().write(f)
133
134 Q = PISM.IceModelVec2S(self.grid, "water_flux_magnitude", PISM.WITHOUT_GHOSTS)
135 Q.set_attrs("", "magnitude of the water flux", "m2 s-1", "m2 s-1", "", 0)
136 Q.set_to_magnitude(model.flux())
137 Q.write(f)
138
139 f.close()
140
141 if __name__ == "__main__":
142
143 t = SteadyHydrology()
144 t.setUp()
145 t.divergence_theorem_test()
146 t.write_results()