Richards 1D integration

As shown in detail at the July summer school about landslides, one of the first approximations to the full Richards equation is the 1-D  Richards equation (please take care of the fact that vertical is not slope-normal). If we want to solve it, besides using the celebrated Hydrus-1D, we can use the novel algorithm developed by Vincenzo Casulli and Paola Zanolli (finite elements, the first, finite differences the second - the second Authors actually implemented a 3-D algorithm). The conceptual line followed by the Casulli and Zanolli is "simple" (but simplicity, as usual, is the product of long and deep thinking):

•  discretise the system to obtain a non-linear system
•  use a Newton method to solve the non linear system

In turn, the Newton's method^1 is implemented in a way that it iterates over a sequence of linear equation in which the relevant equation matrix is symmetric, and therefore the Conjugate Gradient^2 (CG) method can be used to solve it efficiently.
That's the theory. And, in practice ?

S

You need to choose a language for implementing the algorithm. Assume you use Java, because you know it (or you have learnt it, or you just want to learn it, and this can be a challenging but nice exercise).  Since you will be doing a lot of matrix algebra, you first have to choose a library for solving linear algebra (or to implement it by yourself), And if the library can, without you very much take care of it,  thread your algorithms, it is a bargain.  One of these libraries is, for instance, the Parallel Colt. Obviously you need to get used to the library (after having installed it).  some other library can be more efficient, but in Parallel Colt you can find the Conjugate Gradient method already implemented. Once you manage both (PC and CG), you are pretty close to the solution of the problem (easy to say ;-) ) but ...

_________________________________________________________________________________
^1 - A good concise reference is the Davis' book "Direct methods for sparse linear system", SIAM 2006, but the appendix of the paper on Boussinesq equation can clarify several things.
^2 - Give a look to  An Introduction to the Conjugate Gradient Method Without the Agonizing Pain by Jonathan Richard Shewchuk.