tSSAFD_Regional.cc - pism - [fork] customized build of PISM, the parallel ice sheet model (tillflux branch)
HTML git clone git://src.adamsgaard.dk/pism
DIR Log
DIR Files
DIR Refs
DIR LICENSE
---
tSSAFD_Regional.cc (4278B)
---
1 /* Copyright (C) 2015, 2016, 2017 PISM Authors
2 *
3 * This file is part of PISM.
4 *
5 * PISM is free software; you can redistribute it and/or modify it under the
6 * terms of the GNU General Public License as published by the Free Software
7 * Foundation; either version 3 of the License, or (at your option) any later
8 * version.
9 *
10 * PISM is distributed in the hope that it will be useful, but WITHOUT ANY
11 * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
12 * FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
13 * details.
14 *
15 * You should have received a copy of the GNU General Public License
16 * along with PISM; if not, write to the Free Software
17 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
18 */
19
20 #include "SSAFD_Regional.hh"
21 #include "pism/util/Vars.hh"
22 #include "pism/stressbalance/StressBalance.hh"
23
24 namespace pism {
25
26 namespace stressbalance {
27
28 SSAFD_Regional::SSAFD_Regional(IceGrid::ConstPtr g)
29 : SSAFD(g) {
30
31 m_h_stored = NULL;
32 m_H_stored = NULL;
33 m_no_model_mask = NULL;
34 }
35
36 SSAFD_Regional::~SSAFD_Regional() {
37 // empty
38 }
39
40 void SSAFD_Regional::init() {
41
42 SSAFD::init();
43
44 m_log->message(2, " using the regional version of the SSA solver...\n");
45
46 if (m_config->get_flag("stress_balance.ssa.dirichlet_bc")) {
47 m_log->message(2, " using stored SSA velocities as Dirichlet B.C. in the no_model_strip...\n");
48 }
49 }
50
51 void SSAFD_Regional::update(const Inputs &inputs, bool full_update) {
52 m_h_stored = inputs.no_model_surface_elevation;
53 m_H_stored = inputs.no_model_ice_thickness;
54 m_no_model_mask = inputs.no_model_mask;
55
56 SSA::update(inputs, full_update);
57
58 m_h_stored = NULL;
59 m_H_stored = NULL;
60 m_no_model_mask = NULL;
61 }
62
63 static int weight(int M_ij, int M_n, double h_ij, double h_n) {
64 // fjord walls, nunataks, headwalls
65 if ((mask::icy(M_ij) and mask::ice_free(M_n) and h_n > h_ij) or
66 (mask::ice_free(M_ij) and mask::icy(M_n) and h_ij > h_n)) {
67 return 0;
68 }
69
70 return 1;
71 }
72
73 void SSAFD_Regional::compute_driving_stress(const IceModelVec2S &ice_thickness,
74 const IceModelVec2S &surface_elevation,
75 const IceModelVec2CellType &cell_type,
76 const IceModelVec2Int *no_model_mask,
77 IceModelVec2V &result) const {
78
79 SSAFD::compute_driving_stress(ice_thickness, surface_elevation, cell_type, no_model_mask, result);
80
81 double
82 dx = m_grid->dx(),
83 dy = m_grid->dy();
84
85 int
86 Mx = m_grid->Mx(),
87 My = m_grid->My();
88
89 IceModelVec::AccessList list{&result, &cell_type, no_model_mask, m_h_stored, m_H_stored};
90
91 for (Points p(*m_grid); p; p.next()) {
92 const int i = p.i(), j = p.j();
93
94 auto M = no_model_mask->int_star(i, j);
95
96 if (M.ij == 0) {
97 // this grid point is in the modeled area so we don't need to modify the driving
98 // stress
99 continue;
100 }
101
102 double pressure = m_EC->pressure((*m_H_stored)(i, j));
103 if (pressure <= 0) {
104 result(i, j) = 0.0;
105 continue;
106 }
107
108 auto h = m_h_stored->star(i, j);
109 auto CT = cell_type.int_star(i, j);
110
111 // x-derivative
112 double h_x = 0.0;
113 {
114 double
115 west = M.w == 1 and i > 0,
116 east = M.e == 1 and i < Mx - 1;
117
118 // don't use differences spanning "cliffs"
119 west *= weight(CT.ij, CT.w, h.ij, h.w);
120 east *= weight(CT.ij, CT.e, h.ij, h.e);
121
122 if (east + west > 0) {
123 h_x = 1.0 / ((west + east) * dx) * (west * (h.ij - h.w) + east * (h.e - h.ij));
124 } else {
125 h_x = 0.0;
126 }
127 }
128
129 // y-derivative
130 double h_y = 0.0;
131 {
132 double
133 south = M.s == 1 and j > 0,
134 north = M.n == 1 and j < My - 1;
135
136 // don't use differences spanning "cliffs"
137 south *= weight(CT.ij, CT.s, h.ij, h.s);
138 north *= weight(CT.ij, CT.n, h.ij, h.n);
139
140 if (north + south > 0) {
141 h_y = 1.0 / ((south + north) * dy) * (south * (h.ij - h.s) + north * (h.n - h.ij));
142 } else {
143 h_y = 0.0;
144 }
145 }
146
147 result(i, j) = - pressure * Vector2(h_x, h_y);
148 } // end of the loop over grid points
149 }
150
151 } // end of namespace stressbalance
152
153 } // end of namespace pism