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