I am attempting to solve a diffusion equation for the concentration field u(r,t) in cylindrical (p=1) or spherical (p=2) coordinates,
After spatial discretization with the Galerkin finite element method, I end up with a system of ODEs of the form,
where C is also known as the capacity or mass matrix, A is the diffusivity/conductivity matrix,
\mathbf{b} is the concentration supply vector (it captures the effect of boundary conditions), \mathbf{u} is the unknown the concentration vector, and \dot{\mathbf{u}} is the time (t) derivative of \mathbf{u}. Due to the nonlinear property relations, such as the diffusivity D(u(x,t)), the mass matrix and supply vector change with time, but we can assume they have been linearized around the previous time step.
For the time discretization, introductory textbooks generally cover the classic methods such as,
- explicit Euler method
- implicit Euler method
- Crank-Nicolson method
These can be combined into a common form known as the generalized trapezoidal method or \theta-method. Since the mass matrix is not diagonal, a matrix problem must be solved even for explicit time discretization methods, like explicit Euler or the standard explicit Runge-Kutta methods.
Iāve also seen three-time-level methods in engineering textbooks such as,
- fully implicit three-time-level (TTL) scheme (other terms are evaluated at time t_{k + 1}; it may also be blended with the implicit Euler scheme), also known as BDF2:
- three-step ad hoc finite difference method using arithmetic mean:
- three-level ad hoc method with improved stability, combining levels t_{k-2}, t_k, and t_{k+1} (other terms evaluated at time t_k):
The three-level methods are not self starting, and will typically be started with a few steps of the implicit Euler method. In the form shown above, the methods assume a fixed time-step, meaning no possibility for temporal error control. They are at best second-order accurate, but supposed to have good stability.
The alternative path is to use the method of lines (MOL), and combine the semi-discretized system with a ready-to-use (Fortran) library routine. Since I know the problem becomes stiff for fine grid sizes, here are the integrators Iāve been considering:
| Integrator | Problem form | Jacobian | Download Link: |
|---|---|---|---|
| LSODA | \dot{y}= f(t,y) | Dense, banded | ODEPACK | Computing |
| LSODI | M(t,y) \dot{y}= f(t,y) | Dense, banded | ODEPACK | Computing |
| VODE | \dot{y} = f(t,y) | Dense, banded | ODEPACK | Computing |
| MEBDF | M \dot{y} = f(t,y) | Dense, banded | https://www.netlib.org/ode/mebdfdae.f |
| MEBDFI | g(t,y,\dot{y})=0 | Dense, banded | https://www.netlib.org/ode/mebdfi.f |
| DASSL | g(t,y,\dot{y})=0 | Dense, banded | https://www.netlib.org/ode/ddassl.f |
| M3RK | \dot{y} = f(t,y) | (spectral radius, ā ) | Algorithm 553: M3RK, An Explicit Time Integrator for Semidiscrete Parabolic Equations [D3] | ACM Transactions on Mathematical Software |
| RKC | \dot{y} = f(t,y) | (spectral radius, ā ) | https://www.netlib.org/ode/rkc.f |
| DDEABM | \dot{y} = f(t,y) | / | ddeabm.f |
| DDEBDF | \dot{y} = f(t,y) | Dense, banded | ddebdf.f |
| RADAU5 | M \dot{y} = f(t,y) | Dense, banded | Fortran Codes |
| RADAU | M \dot{y} = f(t,y) | Dense, banded | Fortran Codes |
| RODAS | M \dot{y} = f(t,y) | Dense, banded | Fortran Codes |
| SEULEX | M \dot{y} = f(t,y) | Dense, banded | Fortran Codes |
| CVODE | \dot{y} = f(t,y) | Generic | CVODE | Computing |
| ARKODE | M(t) \dot{y} = f(t,y) | Generic | ARKODE | Computing |
| IDA | g(t,y,\dot{y})=0 | Generic | IDA | Computing |
| SPRINT | M(t,y)\dot{y}=f(t,y) | Dense, banded, sparse | D02MN Sub-chapter Introduction : NAG Library, Mark 29 |
ā The M3RK routine of Verwer, and the RKC routine of Sommeijer and Verwer can use an estimate of the spectral radius (in case itās known), to determine the maximum stable step length.
Since the FEM discretization given above is linearly implicit, it can alway be refactored into an explicit form,
hence explicit solvers can be used in principle. In practice we can do this by factorizing the matrix C once at the beginning, and performing a back substitution during the right-hand side evaluation. This is inconvenient as the factorized matrix must be available within the ODE callback function. Iām also afraid that for solvers of the form \dot{y} = f(t,y), which internally use the Jacobian matrix \partial f / \partial y obtained via finite differencing, using a solver such as LSODA, VODE, CVODE, or DDEBDF, would be more expensive than really needed. Moreover, while the matrix C is generally sparse, the inverse will be dense, so the Jacobian will be dense too. This technique is therefore only applicable to strictly explicit solvers.
Further Fortran packages are available, such as FOODIE, rklib, FLINT, or FATODE, but these appear to be focused mainly on non-stiff or mildly stiff problems, explicit time-stepping, and high accuracy.
Which of the solvers above would you recommend?
I know this is not strictly a Fortran-related question, the choice of solver is mostly dictated by the properties of the system of ODEs. Nonetheless, I am hoping that some Discourse members may have previous experience with this type of problem.
Iād also appreciate any pointers to literature discussing the selection of time discretization methods. The most authoritative reference I am aware of is ā Verwer, J. G., & Hundsdorfer, W. H. (2003). Numerical solution of time-dependent advection-diffusion-reaction equations. Berlin: Springer. ā in which the authors state:
[ā¦] the conventional explicit Runge-Kutta methods are in general of little practical use for diffusion equations. Diffusion equations are often better solved with A-stable implicit methods which provide unconditional stability, or with special explicit methods which possess very long stability intervals along the negative real line.
Based on the recommendations in the Julia DifferentialEquations.jl package, I should pick a method like TRBDF2 (one-step ESDIRK method), Rodas4P (A Rosenbrock method designed for stiff parabolic PDEs), or radau (Implicit Runge Kutta). The MATLAB pdepe solver for parabolic PDEās in one spatial variable, uses the integrator ode15s, based on variable order Numerical Differentation Formulas (similar to the methods used in GEAR and LSODE). The NAG parabolic PDE solver - dim1_parab_fd - that uses the same semi-discretization algorithm as pdepe in MATLAB, uses backward differentiation formulas (BDF).