tColumnSystem.cc - 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
---
tColumnSystem.cc (12093B)
---
1 // Copyright (C) 2004-2019 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 #include <cassert>
20 #include <fstream>
21 #include <iostream>
22
23 #include "pism/util/pism_utilities.hh"
24 #include "pism/util/iceModelVec.hh"
25 #include "ColumnSystem.hh"
26
27 #include "pism/util/error_handling.hh"
28 #include "pism/util/ColumnInterpolation.hh"
29
30 namespace pism {
31
32
33 //! Allocate a tridiagonal system of maximum size max_system_size.
34 /*!
35 Let N = `max_system_size`. Then allocated locations are like this:
36 \verbatim
37 D[0] U[0] 0 0 0 ...
38 L[1] D[1] U[1] 0 0 ...
39 0 L[2] D[2] U[2] 0 ...
40 0 0 L[3] D[3] U[3] ...
41 \endverbatim
42 with the last row
43 \verbatim
44 0 0 ... 0 L[N-1] D[N-1]
45 \endverbatim
46 Thus the index into the arrays L, D, U is always the row number.
47 */
48 TridiagonalSystem::TridiagonalSystem(unsigned int max_size,
49 const std::string &prefix)
50 : m_max_system_size(max_size), m_prefix(prefix) {
51 assert(m_max_system_size >= 1 && m_max_system_size < 1e6);
52
53 m_L.resize(m_max_system_size);
54 m_D.resize(m_max_system_size);
55 m_U.resize(m_max_system_size);
56 m_rhs.resize(m_max_system_size);
57 m_work.resize(m_max_system_size);
58 }
59
60 //! Zero all entries.
61 void TridiagonalSystem::reset() {
62 #if Pism_DEBUG==1
63 memset(&m_L[0], 0, (m_max_system_size)*sizeof(double));
64 memset(&m_U[0], 0, (m_max_system_size)*sizeof(double));
65 memset(&m_D[0], 0, (m_max_system_size)*sizeof(double));
66 memset(&m_rhs[0], 0, (m_max_system_size)*sizeof(double));
67 memset(&m_work[0], 0, (m_max_system_size)*sizeof(double));
68 #endif
69 }
70
71 //! Compute 1-norm, which is max sum of absolute values of columns.
72 double TridiagonalSystem::norm1(unsigned int system_size) const {
73 assert(system_size <= m_max_system_size);
74 if (system_size == 1) {
75 return fabs(m_D[0]); // only 1x1 case is special
76 }
77 double z = fabs(m_D[0]) + fabs(m_L[1]);
78 for (unsigned int k = 1; k < system_size; k++) { // k is column index (zero-based)
79 z = std::max(z, fabs(m_U[k-1])) + fabs(m_D[k]) + fabs(m_L[k+1]);
80 }
81 z = std::max(z, fabs(m_U[system_size-2]) + fabs(m_D[system_size-1]));
82 return z;
83 }
84
85
86 //! Compute diagonal-dominance ratio. If this is less than one then the matrix is strictly diagonally-dominant.
87 /*!
88 Let \f$A = (a_{ij})\f$ be the tridiagonal matrix
89 described by L, D, U for row indices 0 through `n`. The computed ratio is
90 \f[ \max_{j=1, \dots, n} \frac{|a_{j, j-1}|+|a_{j, j+1}|}{|a_{jj}|}, \f]
91 where \f$a_{1, 0}\f$ and \f$a_{n, n+1}\f$ are interpreted as zero.
92
93 If this is smaller than one then it is a theorem that the tridiagonal solve will
94 succeed.
95
96 We return -1.0 if the absolute value of any diagonal element is less than
97 1e-12 of the 1-norm of the matrix.
98 */
99 double TridiagonalSystem::ddratio(unsigned int system_size) const {
100 assert(system_size <= m_max_system_size);
101
102 const double scale = norm1(system_size);
103
104 if ((fabs(m_D[0]) / scale) < 1.0e-12) {
105 return -1.0;
106 }
107 double z = fabs(m_U[0]) / fabs(m_D[0]);
108
109 for (unsigned int k = 1; k < system_size - 1; k++) { // k is row index (zero-based)
110 if ((fabs(m_D[k]) / scale) < 1.0e-12) {
111 return -1.0;
112 }
113 const double s = fabs(m_L[k]) + fabs(m_U[k]);
114 z = std::max(z, s / fabs(m_D[k]));
115 }
116
117 if ((fabs(m_D[system_size - 1]) / scale) < 1.0e-12) {
118 return -1.0;
119 }
120 z = std::max(z, fabs(m_L[system_size - 1]) / fabs(m_D[system_size - 1]));
121
122 return z;
123 }
124
125 //! Utility for simple ascii view of a vector (one-dimensional column) quantity.
126 /*!
127 Give first argument NULL to get standard out. No binary viewer.
128
129 Give description string as `info` argument.
130
131 Result should be executable as a Python (SciPy) script.
132
133 Does not stop on non-fatal errors.
134 */
135 void TridiagonalSystem::save_vector(std::ostream &output,
136 const std::vector<double> &v,
137 unsigned int system_size,
138 const std::string &variable) const {
139 assert(system_size >= 1);
140
141 output << variable << " = numpy.array([";
142 for (unsigned int k = 0; k < system_size; k++) {
143 output << v[k] << ", ";
144 }
145 output << "])" << std::endl;
146 }
147
148
149 //! View the tridiagonal matrix.
150 void TridiagonalSystem::save_matrix(std::ostream &output,
151 unsigned int system_size,
152 const std::string &variable) const {
153
154 if (system_size < 2) {
155 std::cout << "\n\n<nmax >= 2 required to view tri-diagonal matrix " << variable
156 << " ... skipping view" << std::endl;
157 return;
158 }
159
160 save_vector(output, m_U, system_size, variable + "_U");
161 save_vector(output, m_D, system_size, variable + "_D");
162 save_vector(output, m_L, system_size, variable + "_L");
163
164 // prepare to convert to a sparse matrix
165 output << variable << "_diagonals = ["
166 << variable << "_L[1:], "
167 << variable << "_D, "
168 << variable << "_U[:-1]]" << std::endl;
169 output << "import scipy.sparse" << std::endl;
170 output << variable << " = scipy.sparse.diags(" << variable << "_diagonals, [-1, 0, 1])" << std::endl;
171 }
172
173
174 //! View the tridiagonal system A x = b to an output stream, both A as a full matrix and b as a vector.
175 void TridiagonalSystem::save_system(std::ostream &output,
176 unsigned int system_size) const {
177 save_matrix(output, system_size, m_prefix + "_A");
178 save_vector(output, m_rhs, system_size, m_prefix + "_rhs");
179 }
180
181 //! Write system matrix, right-hand-side, and (provided) solution into an already-named Python
182 //! script.
183 void TridiagonalSystem::save_system_with_solution(const std::string &filename,
184 unsigned int M,
185 const std::vector<double> &x) {
186 std::ofstream output(filename.c_str());
187 output << "import numpy" << std::endl;
188 output << "# system has 1-norm = " << norm1(M)
189 << " and diagonal-dominance ratio = " << ddratio(M) << std::endl;
190
191 save_system(output, M);
192 save_vector(output, x, M, m_prefix + "_x");
193 }
194
195
196 void TridiagonalSystem::solve(unsigned int system_size, std::vector<double> &result) {
197 result.resize(m_max_system_size);
198
199 solve(system_size, result.data());
200 }
201
202
203 //! The actual code for solving a tridiagonal system.
204 /*!
205 This is modified slightly from a Numerical Recipes version.
206
207 Input size n is size of instance. Requires n <= TridiagonalSystem::m_max_system_size.
208
209 Solution of system in x.
210 */
211 void TridiagonalSystem::solve(unsigned int system_size, double *result) {
212 assert(system_size >= 1);
213 assert(system_size <= m_max_system_size);
214
215 if (m_D[0] == 0.0) {
216 throw RuntimeError(PISM_ERROR_LOCATION, "zero pivot at row 1");
217 }
218
219 double b = m_D[0];
220
221 result[0] = m_rhs[0] / b;
222 for (unsigned int k = 1; k < system_size; ++k) {
223 m_work[k] = m_U[k - 1] / b;
224
225 b = m_D[k] - m_L[k] * m_work[k];
226
227 if (b == 0.0) {
228 throw RuntimeError::formatted(PISM_ERROR_LOCATION, "zero pivot at row %d", k + 1);
229 }
230
231 result[k] = (m_rhs[k] - m_L[k] * result[k-1]) / b;
232 }
233
234 for (int k = system_size - 2; k >= 0; --k) {
235 result[k] -= m_work[k + 1] * result[k + 1];
236 }
237 }
238
239 std::string TridiagonalSystem::prefix() const {
240 return m_prefix;
241 }
242
243 //! A column system is a kind of a tridiagonal system.
244 columnSystemCtx::columnSystemCtx(const std::vector<double>& storage_grid,
245 const std::string &prefix,
246 double dx, double dy, double dt,
247 const IceModelVec3 &u3,
248 const IceModelVec3 &v3,
249 const IceModelVec3 &w3)
250 : m_dx(dx), m_dy(dy), m_dt(dt), m_u3(u3), m_v3(v3), m_w3(w3) {
251 assert(dx > 0.0);
252 assert(dy > 0.0);
253 assert(dt > 0.0);
254
255 init_fine_grid(storage_grid);
256
257 m_solver = new TridiagonalSystem(m_z.size(), prefix);
258
259 m_interp = new ColumnInterpolation(storage_grid, m_z);
260
261 m_u.resize(m_z.size());
262 m_v.resize(m_z.size());
263 m_w.resize(m_z.size());
264 }
265
266 columnSystemCtx::~columnSystemCtx() {
267 delete m_solver;
268 delete m_interp;
269 }
270
271 unsigned int columnSystemCtx::ks() const {
272 return m_ks;
273 }
274
275 double columnSystemCtx::dz() const {
276 return m_dz;
277 }
278
279 const std::vector<double>& columnSystemCtx::z() const {
280 return m_z;
281 }
282
283 void columnSystemCtx::fine_to_coarse(const std::vector<double> &fine, int i, int j,
284 IceModelVec3& coarse) const {
285 double *array = coarse.get_column(i, j);
286 m_interp->fine_to_coarse(&fine[0], array);
287 }
288
289 void columnSystemCtx::coarse_to_fine(const IceModelVec3 &coarse, int i, int j,
290 double* fine) const {
291 const double *array = coarse.get_column(i, j);
292 m_interp->coarse_to_fine(array, m_ks, fine);
293 }
294
295 void columnSystemCtx::init_fine_grid(const std::vector<double>& storage_grid) {
296 // Compute m_dz as the minimum vertical spacing in the coarse
297 // grid:
298 unsigned int Mz = storage_grid.size();
299 double Lz = storage_grid.back();
300 m_dz = Lz;
301 for (unsigned int k = 1; k < Mz; ++k) {
302 m_dz = std::min(m_dz, storage_grid[k] - storage_grid[k - 1]);
303 }
304
305 size_t Mz_fine = static_cast<size_t>(ceil(Lz / m_dz) + 1);
306 m_dz = Lz / (Mz_fine - 1);
307
308 m_z.resize(Mz_fine);
309 // compute levels of the fine grid:
310 for (unsigned int k = 0; k < Mz_fine; ++k) {
311 m_z[k] = storage_grid[0] + k * m_dz;
312 }
313 // Note that it *is* allowed to go over Lz.
314 }
315
316 void columnSystemCtx::init_column(int i, int j,
317 double ice_thickness) {
318 m_i = i;
319 m_j = j;
320 m_ks = static_cast<unsigned int>(floor(ice_thickness / m_dz));
321
322 // Force m_ks to be in the allowed range.
323 if (m_ks >= m_z.size()) {
324 m_ks = m_z.size() - 1;
325 }
326
327 m_solver->reset();
328
329 // post-condition
330 #if Pism_DEBUG==1
331 // check if m_ks is valid
332 if (m_ks >= m_z.size()) {
333 throw RuntimeError::formatted(PISM_ERROR_LOCATION, "ks = %d computed at i = %d, j = %d is invalid,\n"
334 "possibly because of invalid ice thickness (%f meters) or dz (%f meters).",
335 m_ks, m_i, m_j, ice_thickness, m_dz);
336 }
337 #endif
338 }
339
340 //! Write system matrix and right-hand-side into an Python script. The file name contains ZERO_PIVOT_ERROR.
341 void columnSystemCtx::reportColumnZeroPivotErrorMFile(unsigned int M) {
342
343 auto filename = pism::printf("%s_i%d_j%d_ZERO_PIVOT_ERROR.py",
344 m_solver->prefix().c_str(), m_i, m_j);
345
346 std::ofstream output(filename);
347 output << "# system has 1-norm = " << m_solver->norm1(M)
348 << " and diagonal-dominance ratio = " << m_solver->ddratio(M) << std::endl;
349
350 m_solver->save_system(output, M);
351 }
352
353
354 //! @brief Write system matrix, right-hand-side, and (provided)
355 //! solution into Python script. Constructs file name from m_prefix.
356 void columnSystemCtx::save_to_file(const std::vector<double> &x) {
357
358 auto filename = pism::printf("%s_i%d_j%d.py", m_solver->prefix().c_str(), m_i, m_j);
359
360 std::cout << "saving "
361 << m_solver->prefix() << " column system at (i,j)"
362 << " = (" << m_i << "," << m_j << ") to " << filename << " ...\n" << std::endl;
363
364 this->save_to_file(filename, x);
365 }
366
367 void columnSystemCtx::save_to_file(const std::string &filename,
368 const std::vector<double> &x) {
369 m_solver->save_system_with_solution(filename, m_z.size(), x);
370 }
371
372 } // end of namespace pism