19 #ifndef PISM_COLUMNSYSTEM_HH
20 #define PISM_COLUMNSYSTEM_HH
86 double norm1(
unsigned int system_size)
const;
87 double ddratio(
unsigned int system_size)
const;
92 void solve(
unsigned int system_size, std::vector<double> &result);
93 void solve(
unsigned int system_size,
double *result);
96 unsigned int system_size,
97 const std::vector<double> &solution);
102 unsigned int system_size)
const;
105 unsigned int system_size,
106 const std::string &variable)
const;
109 const std::vector<double> &v,
110 unsigned int system_size,
111 const std::string &variable);
113 std::string
prefix()
const;
115 double&
L(
size_t i) {
118 double&
D(
size_t i) {
121 double&
U(
size_t i) {
143 columnSystemCtx(
const std::vector<double>& storage_grid,
const std::string &prefix,
144 double dx,
double dy,
double dt,
149 void save_to_file(
const std::string &filename,
const std::vector<double> &x);
151 unsigned int ks()
const;
153 const std::vector<double>&
z()
const;
154 void fine_to_coarse(
const std::vector<double> &fine,
int i,
int j,
180 void init_column(
int i,
int j,
double ice_thickness);
A virtual class collecting methods common to ice and bedrock 3D fields.
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.
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.
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.
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.
void fine_to_coarse(const std::vector< double > &fine, int i, int j, IceModelVec3 &coarse) const
const std::vector< double > & z() const
const IceModelVec3 & m_u3
pointers to 3D velocity components
const IceModelVec3 & m_w3
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)
const IceModelVec3 & m_v3
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.
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
Base class for tridiagonal systems in the ice.