All Classes Files Functions Variables Typedefs Pages
Solver Howto
For proper usage of the Doxygen reference of the solver library, please make sure that EXTRACT_PRIVATE is set to YES in the Doxygen configuration. Otherwise, the detailed interface description is not generated!


The solver modules consists of several modules:

Solve function for preconditioned CG and Chebyshev

The easiest way to use the solver module is to setup all configuration and use the generic solve function to solve a stencil-based linear equation system A=b with the help of the CG method or Chebyshev iteration and a preconditioner.

In MPIOM, the solver configuration in generated in mo_solver_hlp.f90. There, the read_solver_config function reads the configuration from file SLVCTL and fills the type solver_config with it. SLVCTL is generated by the preparte_job_mpiom_hamocc script. It comes with detailed comments on its usage and the meaning of the various parameters. Additionally, this configuration is checked (check_solver_config), broadcast (broadcast_solver_config) and logged (logwrite_solver_config). With the help of the functions get_interior and set_extents the extents of the solution vector z1o is set.

During the startup of MPIOM in mpiom.f90, a solver_config is generated and passed to init_solver to initialize the solver library with the user-set configuration. The actual utilization of the solver takes place in mo_tro.f90 in the function solve_barotro_sys. There, first the calc_barotro_rhs function is called to calculate the right-hand-side of the barotropic system then stencil (i.e. the linear system) of the barotropic system is set via set_stencil. Next, the checkprep_params from mo_solver_hlp function checkes and prepares all parameters by looking at the parameters set by the solver initialization and complementing them. For instance the eigenvalues lambda_min and max can be set in SLVCTL or are retrieved automatically during this preparation step which is executed only once. After this a call to the solve function solves the system. As parameters this function needs the barotropic stencil operation as given by apply_stencil_dp , the right-hand-side b1o, the starting value z1o with its extents ext_z1o as well as a function for the boundary exchange (bounds_exchange_dp) which is defined in mo_solver_hlp in MPIOM. Since almost all functions are defined for single (suffix _sp) and double (suffix _dp) precision we pass the double precision variants of all funtions to solve. The corresponding code is:

! Read config and initialize
slv_config = read_solver_config()
CALL init_solver(slv_config)
! Set the matrix stencil
CALL set_stencil(uf, vf, ff, ext_z1o)
! Check and prepare all parameters
CALL checkprep_params(z1o, ext_z1o)
! Solve the barotropic subsystem
kiter = solve(apply_stencil_dp,b1o,z1o,ext_z1o,bounds_exchange_dp)

The variable z1o now holds the result and kiter the needed iterations. As stated before, the solve function is just a convenient interface to the precond_cg_method and the precond_chebyshev_method as well as all preconditioners in preconditioners module. Direct usage of these two solvers is pretty self-explanatory and can be looked up easily in this reference.

Additive Schwarz Method

The additive Schwarz method schwarz_method resides in the solvers module. This method uses a local method to solve the global problem first on each partition seperately with fixed Dirichlet boundaries. Then, the extended overlapping boundary halos are exchanged and another local problem is solved. This step is repeated until the solution of the global problem has converged.

To use this method, we first need to set up arrays with extended boundary halos. This is done for instance via the setup_halos function (which then uses sethaloN) in mo_tro.f90 of MPIOM. This function enlarges the defining arrays of the barotropic system, i.e., uf, vf, ff, as well as z1o and b1o. The next step is to use the enlarged fields zuf, zvf, zff to initialize the matrix stencil with set_stencil. Then all parameters are checked again via checkprep_params in MPIOM. After this we can directly call schwarz_method with the function reference apply_stencil. The following code snipped illustrates this.

! enlargen fields uf, vf, ff, z1o, b1o to zuf, zvf, zff, zz1o, zb1o
new_comm = setup_halos(zz1o, zuf, zvf, zff, zb1o, SLV_SCHWARZ_HALOS)
! Set up matrix stencil with extra halos
CALL set_stencil(zuf, zvf, zff, ext_zz1o)
! Check and prepare all parameters
CALL checkprep_params(zz1o, ext_zz1o)
! Solve the barotropic subsystem
kiter = schwarz_method(apply_stencil_dp, zb1o, zz1o, ext_zz1o, solve_dp, bounds_exchange_dp)
! Get result from field with extra halos
DO i=1,ie
DO j=1,je

Multi-Precision Iterative Refinement

It is also possible to use the CG method in single precision while obtaining double precision accuracy in the end. The hope is that solving in single precision (sp) is faster than in double precision (dp) because twice as many sp numbers can be transferred in the same time as db numbers. Also, if the VMC (on PowerPC) or SSE units (on Intel) is used, vector operations in sp can be twice as fast as in dp. On the POWER6 architecture the VMX units as well as the compiler is quite limited according to IBM but the future POWER7+ with its VSX units will drastically increase the sp performace compared to dp.

To make future use of this approach, nearly all functions and routines in the solver library have a single and double precision variant. This was accomplished by splitting each module into a declaration file, e.g. solvers.f90 and an implementation file, e.g. solvers_multi.f90. Preprocessor macros take care that the implementation file is included into the declaration file with the corresponding data kinds set accordingly. Also, each name of a multi-precision function gets the suffix _sp resp. _dp for single precision resp. double precision.

The iter_refinement implements the multi-precision iterative refinement. It takes as parameters the stencil operations in single and double precision. The code below illustrates the usage taken from solve_barotro_sys.

ALLOCATE(uf_sp(ie,je), vf_sp(ie,je), ff_sp(ie,je))
uf_sp = REAL(uf, sp)
vf_sp = REAL(vf, sp)
ff_sp = REAL(ff, sp)
! Setup single precision matrix stencil
CALL set_stencil(uf_sp, vf_sp, ff_sp, ext_z1o)
! Setup double precision matrix stencil
CALL set_stencil(uf, vf, ff, ext_z1o)
! Check and prepare all parameters
CALL checkprep_params(z1o, ext_z1o)
kiter = iter_refinement(apply_stencil_dp, apply_stencil_sp, b1o, z1o, ext_z1o, bounds_exchange_sp)
DEALLOCATE(uf_sp, vf_sp, ff_sp)

Das diesem Bericht zugrundeliegende Vorhaben wurde mit Mitteln des Bundesministeriums für Bildung, und Forschung unter dem Förderkennzeichen 01IH08004E gefördert. Die Verantwortung für den Inhalt dieser Veröffentlichung liegt beim Autor.