Virtual base class. Abstracts a tridiagonal system to solve in a column of ice and/or bedrock. More...

#include <columnSystem.hh>

Inheritance diagram for columnSystemCtx:

Public Member Functions
	columnSystemCtx (PetscInt my_nmax, string my_prefix)
	Allocate a tridiagonal system of maximum size my_nmax.
	~columnSystemCtx ()
PetscErrorCode	setIndicesAndClearThisColumn (PetscInt i, PetscInt j, PetscInt ks)
PetscScalar	norm1 (const PetscInt n) const
	Compute 1-norm, which is max sum of absolute values of columns.
PetscScalar	ddratio (const PetscInt n) const
	Compute diagonal-dominance ratio. If this is less than one then the matrix is strictly diagonally-dominant.
PetscErrorCode	viewVectorValues (PetscViewer viewer, PetscScalar v, PetscInt m, const char info) const
	Utility for simple ascii view of a vector (one-dimensional column) quantity.
PetscErrorCode	viewMatrix (PetscViewer viewer, const char *info) const
	View the tridiagonal matrix. Views as a full matrix if nmax <= 120, otherwise by listing diagonals.
PetscErrorCode	viewSystem (PetscViewer viewer) const
	View the tridiagonal system A x = b to a PETSc viewer, both A as a full matrix and b as a vector.
PetscErrorCode	reportColumnZeroPivotErrorMFile (const PetscErrorCode errindex)
	Write system matrix and right-hand-side into an m-file. The file name contains ZERO_PIVOT_ERROR.
PetscErrorCode	viewColumnInfoMFile (PetscScalar *x, PetscInt n)
	Write system matrix, right-hand-side, and (provided) solution into an m-file. Constructs file name from prefix.
PetscErrorCode	viewColumnInfoMFile (char filename, PetscScalar x, PetscInt n)
	Write system matrix, right-hand-side, and (provided) solution into an already-named m-file.
Protected Member Functions
PetscErrorCode	solveTridiagonalSystem (const PetscInt n, PetscScalar **x)
	The actual code for solving a tridiagonal system. Return code has diagnostic importance.
Protected Attributes
PetscInt	nmax
PetscScalar *	L
PetscScalar *	Lp
PetscScalar *	D
PetscScalar *	U
PetscScalar *	rhs
PetscScalar *	work
PetscInt	i
PetscInt	j
PetscInt	ks
Private Member Functions
PetscErrorCode	resetColumn ()
	Zero all entries.
Private Attributes
bool	indicesValid
string	prefix

Detailed Description

Virtual base class. Abstracts a tridiagonal system to solve in a column of ice and/or bedrock.

Because both the age evolution and conservation of energy equations require us to set up and solve a tridiagonal system of equations, this is structure is worth abstracting.

This base class just holds the tridiagonal system and the ability to solve it, but does not insert entries into the relevant matrix locations. Derived classes will actually set up instances of the system.

The sequence requires setting the column-independent (public) data members, calling the initAllColumns() routine, and then setting up and solving the system in each column.

Definition at line 39 of file columnSystem.hh.

Constructor & Destructor Documentation

columnSystemCtx	(	PetscInt	my_nmax,
		string	my_prefix
	)

Allocate a tridiagonal system of maximum size my_nmax.

Let N = nmax. Then allocated locations are like this:

D[0]   U[0]    0      0      0    ...
L[1]   D[1]   U[1]    0      0    ...
 0     L[2]   D[2]   U[2]    0    ...
 0      0     L[3]   D[3]   U[3]  ...

with the last row

0       0     ...     0  L[N-1]  D[N-1]

Thus the index into the arrays L, D, U is always the row number.

Note L[0] is not allocated and U[N-1] is not allocated.

Definition at line 42 of file columnSystem.cc.

References D, indicesValid, L, Lp, nmax, resetColumn(), rhs, U, and work.

~columnSystemCtx ( )

Definition at line 67 of file columnSystem.cc.

References D, Lp, rhs, U, and work.

Member Function Documentation

PetscScalar ddratio ( const PetscInt n ) const

Compute diagonal-dominance ratio. If this is less than one then the matrix is strictly diagonally-dominant.

Let $A = (a_{ij})$ be the tridiagonal matrix described by L,D,U for row indices 0 through n. The computed ratio is

$\max_{j=1,\dots,n} \frac{|a_{j,j-1}|+|a_{j,j+1}|}{|a_{jj}|},$

where $a_{1,0}$ and $a_{n,n+1}$ are interpreted as zero.

If this is smaller than one then it is a theorem that the tridiagonal solve will succeed.

We return -1.0 if the absolute value of any diagonal element is less than 1e-12 of the 1-norm of the matrix.

Definition at line 117 of file columnSystem.cc.

References D, e, L, nmax, norm1(), and U.

Referenced by viewColumnInfoMFile().

PetscScalar norm1 ( const PetscInt n ) const

Compute 1-norm, which is max sum of absolute values of columns.

Definition at line 89 of file columnSystem.cc.

References D, L, n, nmax, and U.

Referenced by ddratio(), and viewColumnInfoMFile().

PetscErrorCode reportColumnZeroPivotErrorMFile ( const PetscErrorCode errindex )

Write system matrix and right-hand-side into an m-file. The file name contains ZERO_PIVOT_ERROR.

Definition at line 332 of file columnSystem.cc.

References i, j, prefix, and viewColumnInfoMFile().

Referenced by IceModel::ageStep(), IceModel::enthalpyAndDrainageStep(), and IceModel::temperatureStep().

PetscErrorCode resetColumn ( ) [private]

Zero all entries.

Definition at line 77 of file columnSystem.cc.

References D, Lp, nmax, rhs, U, and work.

Referenced by columnSystemCtx(), and setIndicesAndClearThisColumn().

PetscErrorCode setIndicesAndClearThisColumn	(	PetscInt	i,
		PetscInt	j,
		PetscInt	ks
	)

Definition at line 140 of file columnSystem.cc.

References i, indicesValid, j, ks, and resetColumn().

Referenced by IceModel::ageStep(), IceModel::enthalpyAndDrainageStep(), and IceModel::temperatureStep().

PetscErrorCode solveTridiagonalSystem	(	const PetscInt	n,
		PetscScalar **	x
	)		`[protected]`

The actual code for solving a tridiagonal system. Return code has diagnostic importance.

This is modified slightly from a Numerical Recipes version.

Input size n is size of instance. Requires n <= columnSystemCtx::nmax.

Solution of system in x.

Success is return code zero. Positive return code gives location of zero pivot. Negative return code indicates a software problem.

Definition at line 302 of file columnSystem.cc.

References D, indicesValid, L, n, nmax, rhs, U, and work.

Referenced by ageSystemCtx::solveThisColumn(), tempSystemCtx::solveThisColumn(), and enthSystemCtx::solveThisColumn().

PetscErrorCode viewColumnInfoMFile	(	PetscScalar *	x,
		PetscInt	n
	)

Write system matrix, right-hand-side, and (provided) solution into an m-file. Constructs file name from prefix.

An example of the use of this procedure is from examples/searise-greenland/ running the enthalpy formulation. First run spinup.sh in that directory (FIXME: which was modified to have equal spacing in z, when I did this example) to generate g20km_steady.nc. Then:

  $ pismr -ocean_kill -e 3 -atmosphere searise_greenland -surface pdd -config_override  config_269.0_0.001_0.80_-0.500_9.7440.nc \
    -no_mass -y 1 -i g20km_steady.nc -view_sys -id 19 -jd 79

    ...

  $ octave -q
  >> enth_i19_j79
  >> whos

    ...

  >> A = enth_A; b = enth_rhs; x = enth_x;
  >> norm(x - (A\b))/norm(x)
  ans =  1.4823e-13
  >> cond(A)
  ans =  2.6190

Of course we can also do spy(A), eig(A), and look at individual entries, and row and column sums, and so on.

Definition at line 371 of file columnSystem.cc.

References i, j, and prefix.

Referenced by IceModel::ageStep(), IceModel::enthalpyAndDrainageStep(), reportColumnZeroPivotErrorMFile(), and IceModel::temperatureStep().