Modules Description

The driver routines have to assure that the modules below are fed with all necessary arrays and variables set:

1. The problem dependent sizes have to be defined in the module problem_size, but can be variables that are readin anywhere. The following variables rely on these definititions.

2. An array y(nb_rows,nequ), which holds all of the dependent variables for each equation node. This has to be filled up with an estimate for the solution for each call of the Newton-Raphson solver. The array y maybe placed into a module storage together with an allocation routine (see also next).

3. An array u(ngrid) for the independent variable. This has to be placed into a module storage which is usd by the solver package and maybe in user-written routines.

4. The module problem, which contains a routine defining the equations and a routine perturbing the dependent variables.

5. Optionally, a module physics, which holds all physics routines and is used by the user-supplied equation-defining routine (see problem).

6. Several option switches in the module switches to control debugging, timing etc. may be redefined from their default values.

Modules can hold both data and routines. The module data is accessible to all module routines and maybe to the outside (depending on declaration as public or private). Any public variable or module routine may be accessed via a statement like:

USE module, only : variable, routine

The Newton-Raphson solver module consists of the routines solution, get_corrections, apply_corrections and calls the solver package. The first two are mostly driver routines. Only solution is directly accessible from the module. The routine get_corrections is the driver for the solver package. Apply_corrections implements the correction strategy which can affect the convergence of the iteration severely, therefore, it is a user-modifiable routine. The whole dataflow between the core modules is shown in the diagram to the right.

After setting up problem_size (and, maybe, switches), solution can be called with the following variables:

call solution (model_no, y, restart)

The solver package module consists of the routines init_solver, set_jacobi, preconditioner, [analyser], lin_solver, mul_m_v and relies on proper setup of the problem_size module. The init_solver routine must be called first and for every solution if changes of one of the sizes in problem_size are expected during runtime. Then automatic reallocation is initiated. Set_jacobi automatically sets the Jacobian matrix for the Newton_Raphson solver. It calls the user-supplied routines equations, perturbation from the module problem. Preconditioner balances the matrix for stable inversion. [Analyser analyses the matrix for possible stability problems and tries to hint at possible causes. It uses lin_solver and mul_m_v.] The lin_solver routine is a special block band diagonal solver for arbitrary bandwidth and a variety of boundary forms. It uses a variant of Gauss elimination with partial or full pivoting on successive frames. It is also available as a standalone. Mul_m_v is a utility routine to multiply the special matrix form with a vector and is used by analyser.

The solver package routines are normally called by get_corrections from the Newton-Raphson module. For reference, below are the calling parameters. For details on the arrays, consult the Lin_solver documentation.

call init_solver ([reallocate]) :

call set_jacobi (S, y, update_jacobi) :

call preconditioner (S) :

[call analyser (S, failed, full_piv)] (not fully operational yet) :

call lin_solver (S, sol, failed, debug, full_piv) :

call mul_m_v (S, v) :

This module holds all problem dependent sizes and may be extended by any related variables the user likes to place there. These can be variables or constants, depending on whether they are to be changed during runtime or not. All modules should adapt automatically when one of these sizes change. As most of these variables affect the solver package, the Lin_solver documentation should be consulted for further details. See below for correct choices for selected discretization schemes.

This module holds all problem specific routines and is supplied by the user. It must contain at least two routines: equations and perturbation. These are called from the solver package routine set_jacobi and are imported therein by use of the module problem. Below you find some (skeleton) sample code.

subroutine equations (u, var, compute_structure, compute_equations, En, unperturbed)

USE precision   , only : dp, zero, one
USE problem_size, only : nequ 
USE physics     , only : structure 

real (KIND=dp), intent (   IN) : u(:)
real (KIND=dp), intent (INOUT) : var(:,:)
real (KIND=dp), intent (  OUT) : En(0:,:)
logical       , intent (   IN) : compute_structure, compute_equations, unperturbed

real (KIND=dp) :: rhs(SIZE(En,DIM=1)), dydu(SIZE(En,DIM=1))

  if (compute_structure) then
     call structure (var)
  endif

  if (compute_equations) then
     rhs  = .....
     dydu = .....
     En(..,..) = rhs-dydu
  endif

end subroutine equations

module physics
CONTAINS
  subroutine structure (var)
  USE precision   , only : dp
  USE problem_size, only : nequ 
  implicit NONE
  real (KIND=dp), intent(INOUT) :: var(:,:)
  real (KIND=dp), pointer       :: rho(SIZE(var,DIM=1))
 
  rho => var(:,nequ+1)
  rho = ...

  end subroutine structure
end module physics

In equations, the user has to fill up an array En(0:nb_rows,1:nequ), which is supplied from the calling routine, with the nequ equations En = rhs - dy/du at each equation node. u(nb_rows) is the independent variable supplied from above and computed by the user somewhere else before calling the Newton-Raphson module. rhs stands for the right hand side of the differential equations system to be computed by the user. The nequ dependent variables y at the nb_rows nodes in their current state are supplied in the array var(1:nb_rows,1:nequ) (note that possible eigenvalues are supplied in the array locations var(nb_rows-npar+1:nb_rows,:)). The zero index in En has to be filled if the first boundary equations do not fill up a whole block jb1<nequ,jb1>0). In the other cases it's the users choice if to start in En(0,*)(jb1=nequ) or En(1,*)(jb1=0).

Additionally, there are two logical switches, compute_structure and compute_equations, which request action from the user. These serve for the special case when there are more than the original nequ dependent variables y entering the equations. These may be secondary variables like, e.g., the density, which depend on the basic dependent ones. In this case, the user must compute these secondary variables first (for compute_structure = .true.) for the given basic dependent ones and has to put them into the additional space supplied for that purpose in var(*,nequ+1:nvars) (nvars has to be specified in problem_size). Thereafter, they can be used to compute the equations, which has to be done only on request of compute_equations = .true.. Note that the two branches should not be separated by an 'else if' construct, as they may be called concurrently. Due to the automatic Jacobian derivative computation, equations is called nequ+1+npar times.

The third logical switch, unperturbed, indicates that an unperturbed set of variables is provided in var. This switch may be used to store the computed secondary variables into another array which is exported via a module to other routines, e.g., for output to a file.

Computation of, e.g., microphysics, as in a stellar evolution code, may be placed inside an additional module named,say, physics, which is then used inside equations and can simply be called by call structure(var).

subroutine perturbation (y, delta_y)

USE precision, only : dp, zero
  implicit NONE
  real (KIND=dp), intent (IN) :: y (:)
  real (KIND=dp), intent(OUT) :: delta_y (:)

  real (KIND=dp), parameter :: factor = 0.001_dp
  real (KIND=dp) :: ymin

  ymin    = MINVAL( ABS(y), MASK=(y /= zero) )
  delta_y = factor*(y + SIGN(ymin,y))

end subroutine PERTURBATION

This subroutine gets the y(nb_rows,nequ) array of solution estimates as input and should return an array delta_y(nb_rows,nequ) with the perturbations for y. How to perturb reasonably depends largely on the problem and can even vary during the course of the iteration. As a first try, a relative perturbation of 10^-3 to 10^-6 may suffices. At zeros, the perturbation should be based on the adjacent nonzero data points.

This module holds some option switches for informational output and have default settings:

This module contains constants for KIND - selection of variables and other useful ones. The KIND parameters should be used instead of the obsolete double precision or even the nonstandard real*8 declarations. This also applies to constants, e.g., 6.67E-8_dp instead of 6.67D-8. This ensures portability and easy maintenance.

Home