Porous Flow Tutorial Page 02. Numerical issues
This tutorial page discusses a number of numerical issues involved in the very simple model presented in Page 01. So, you're really only interested in modelling, and not numerical mathematics rubbish? Unfortunately, to run PorousFlow effectively, the user must be keenly aware of all these numerical issues, so it is useful to discuss them in full in such a simple setting.
Newton nonlinear solves
Many users of MOOSE use the PJFNK method. This is not recommended in PorousFlow. Instead, a full Newton solve is recommended. This is implemented by the solve_type = Newton in the Executioner block:
[Executioner]
type = Transient
solve_type = Newton
end_time = 1E6
dt = 1E5
nl_abs_tol = 1E-13
[]
This means that MOOSE uses the Newton-Raphson method (see Wikipedia) to solve the system of equations generated by the DE. The Jacobian used in this method is the derivative of the residual (the DE) with respect to the Variables (just porepressure in this example). Much of the PorousFlow code is devoted to computing this Jacobian precisely, so users are strongly encouraged to use it. It typically improves convergence speed by a few orders of magnitude (although in the simple example of this tutorial page, which is a linear problem, it actually makes little difference).
When using solve_type = Newton, the overall solve strategy is:
Form the residual, by computing the DEs using the current variables. The nonlinear tolerances are checked to see if the solution has been found. If not the following steps are taken.
Form the Jacobian matrix, by computing the derivatives of the DEs with respect to the variables, and evaluating using the current variables
Invert the Jacobian matrix. This is the "linear solve". It may be an iterative process, depending on the flags provided to PETSc. It is deemed successful when certain linear tolerances are met
Using the inverted Jacobian, update the current variables, employing the Newton-Raphson method
This completes one "nonlinear iteration". Another nonlinear iteration is started from the first step, above.
The two difficulties most commonly encountered are:
PETSc finds it difficult to invert the Jacobian, ie, the linear solve fails or takes thousands of iterations to finally converge. This almost definitely means that a stronger preconditioner is needed.
Many nonlinear iterations are needed at each timestep. This is usually because your simulation is highly nonlinear. Sometimes this is simply "too bad", but often it is due to difficult numerical situations that can be fixed with some careful analysis.
Nonlinear solves and the Dictator's variables
In this simple example, the PorousFlowDictator isn't built explicitly, since it is built internally by the PorousFlowBasicTHM Action. However, in many input files you will see something like
[dictator]
type = PorousFlowDictator
porous_flow_vars = 'temp pp disp_x disp_y disp_z'
number_fluid_phases = 1
number_fluid_components = 1
[]
The PorousFlowDictator is given the name dictator and the number of fluid phases and fluid components are specified. The porous_flow_vars specifies the PorousFlow variables. All derivatives with respect to these variables enter the Jacobian, while derivatives with respect to any other Variables are not computed. Therefore it is strongly recommended that all Variables are entered into the porous_flow_vars, unless there is a good reason to not compute the derivatives with respect to a particular variable.
Preconditioning and the linear solve
In each nonlinear iteration (as MOOSE slowly converges to the solution for that given timestep) a linear solve must be performed to invert the Jacobian. Often in PorousFlow simulations, the Jacobian is quite ill conditioned, so a strong preconditioner is needed. Typically, one of two choices are made:
[Preconditioning]
active = basic
[basic]
type = SMP
full = true
[]
[preferred_but_might_not_be_installed]
type = SMP
full = true
petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
petsc_options_value = ' lu mumps'
[]
[]
More remarks can be found in Preconditioning and solvers. The PorousFlow tests and examples use some alternative preconditioners, and users struggling with convergence of the linear solve may like to explore those, but the two methods shown above work nicely in most cases.
In the simple example studied in this page, the Preconditioning makes very little difference, but in full-scale nonlinear PorousFlow simulations it can decrease or increase solve times by orders of magnitude.
Linear and nonlinear tolerances
A detailed discussion of setting the nonlinear tolerances may be found in convergence criteria. For single-phase simulation of this tutorial chapter, the residual is (The missing factor of between this formula and that in convergence criteria is because the DE solved in this page does not multiply by density, which is unusual for PorousFlow.) Here is the volume of the "interesting region" in the model, and is the tolerable error in . Supposing that the "interesting region" is the m surrounding the borehole, and a tolerable error is 1Pa.m, this yields . This provides a rough estimate for the nonlinear residual in the Executioner block:
[Executioner]
type = Transient
solve_type = Newton
end_time = 1E6
dt = 1E5
nl_abs_tol = 1E-13
[]
Numerical stabilization
A very careful reader will have noticed that during the simulation the porepressure falls below its initial value of 0: look again at the scale on the resulting animation shown in the previous page. This is completely unphysical.
This is purely a numerical issue that occurs in all simple-minded finite-element (and finite-difference) simulations, and is discussed in the numerical stabilization page. In later tutorial pages, numerical stabilization will be employed, since it is the default in most of PorousFlow, but let us ignore this issue for now.