time_integration Module

SSPRK = Strong Stability Preserving Runge-Kutta (Shu & Osher, 1988).

Supported schemes: 'euler' — Explicit Euler (1st order) 'ssprk22' — SSPRK(2,2) (Shu & Osher, 1988) 'rk3' — TVD-RK3 (Shu & Osher, 1988) [default] 'rk4' — Classic RK4 'ssprk54' — SSPRK(5,4) (Spiteri & Ruuth, 2002) 'beuler' — Backward Euler with Newton-Raphson (implicit, 1st order) 'bdf2' — BDF2 with Newton-Raphson (implicit, 2nd order)

Implicit banded solver selection is controlled per-instance via state%cfg%lapack_solver = .true. (default) — LAPACK dgbsv (pivoted, faster) state%cfg%lapack_solver = .false. — built-in Gaussian elimination (no pivoting) Abstract interface satisfied by every single-step stepper.


Uses


Variables

Type Visibility Attributes Name Initial
procedure(stepper_iface), public, pointer :: step => null()

Procedure pointer to the active time-stepping scheme. Initialised to null; set once by init_time_scheme() before the time loop.
integer, private, parameter :: n_newton = 3

Maximum Newton-Raphson iterations per time step (beuler and bdf2). 3 iterations is sufficient for BDF2 at CFL ≤ 10; increase for stiff problems.
real(kind=wp), private, parameter :: tol_newton = 1.0e-10_wp

Newton convergence tolerance: max-norm of correction ΔQ must fall below this. 1e-10 gives full convergence well below solver tolerance for double precision.
real(kind=wp), private, parameter :: eps_jac = 1.0e-7_wp

Finite-difference step size for Jacobian column approximation. The classic optimum for a one-sided FD balances truncation error O(h) against floating-point cancellation O(ε/h) at h ≈ √ε ≈ √(2.2e-16) ≈ 1.5e-8. We use 1e-7 (~6.5× larger): somewhat above the theoretical √ε optimum, chosen for robustness against cancellation in the column FD of the residual where the per-component magnitudes vary widely.

Abstract Interfaces

abstract interface

  • public subroutine stepper_iface(state)

    Arguments

    Type IntentOptional Attributes Name
    type(solver_state_t), intent(inout), target :: state

Derived Types

type, private ::  local_layout_save_t

Snapshot of a solver_state_t's per-rank local layout, used by the implicit gather-solve-scatter wrappers when rank 0 temporarily swaps into a single-rank decomp.

Components

Type Visibility Attributes Name Initial
type(decomp_t), public :: decomp
real(kind=wp), public, allocatable :: ub_interior(:,:)
real(kind=wp), public, allocatable :: bdf2_ub_prev(:,:)
logical, public :: bdf2_was_initialized = .false.

Subroutines

private subroutine save_local_layout(state, snap)

Capture the per-rank local layout of state so it can be restored after the implicit gather-solve-scatter window. Saves the decomp, the ub-interior slice, and bdf2_ub_prev if present.

Arguments

Type IntentOptional Attributes Name
type(solver_state_t), intent(in) :: state
type(local_layout_save_t), intent(out) :: snap

private subroutine restore_local_layout(state, snap)

Restore the per-rank local layout captured by save_local_layout. Releases all work arrays, swaps state%decomp back to the saved local decomp, reallocates work arrays at local size, and copies the saved ub-interior (and bdf2_ub_prev if present) back in. The bdf2_initialized flag is also restored.

Arguments

Type IntentOptional Attributes Name
type(solver_state_t), intent(inout) :: state
type(local_layout_save_t), intent(in) :: snap

public subroutine resolve_time_scheme(stepper, scheme)

Resolve a scheme name to a specific stepper procedure pointer.

Read more…

Arguments

Type IntentOptional Attributes Name
procedure(stepper_iface), intent(out), pointer :: stepper
character(len=*), intent(in) :: scheme

public subroutine init_time_scheme(scheme)

Bind the procedure pointer @p step to the requested scheme.

Read more…

Arguments

Type IntentOptional Attributes Name
character(len=*), intent(in) :: scheme

private subroutine rhs_1d(ctx)

stepper_kernels rhs callback: recompute the 1D residual (halo exchange + BCs + interior resid happen inside compute_resid).

Arguments

Type IntentOptional Attributes Name
class(*), intent(inout), target :: ctx

private subroutine tvd_rk3_step(state)

Advance the solution by one time step using the TVD-RK3 scheme.

Read more…

Arguments

Type IntentOptional Attributes Name
type(solver_state_t), intent(inout), target :: state

private subroutine ssprk54_step(state)

Advance the solution by one time step using the SSPRK(5,4) scheme.

Read more…

Arguments

Type IntentOptional Attributes Name
type(solver_state_t), intent(inout), target :: state

private subroutine euler_step(state)

Advance the solution by one time step using the explicit Euler scheme.

Read more…

Arguments

Type IntentOptional Attributes Name
type(solver_state_t), intent(inout), target :: state

private subroutine ssprk22_step(state)

Advance the solution by one time step using the SSPRK(2,2) scheme.

Read more…

Arguments

Type IntentOptional Attributes Name
type(solver_state_t), intent(inout), target :: state

private subroutine rk4_step(state)

Advance the solution by one time step using the classic RK4 scheme.

Read more…

Arguments

Type IntentOptional Attributes Name
type(solver_state_t), intent(inout), target :: state

private subroutine setup_band_storage(state, n_dof, kl, ku, ldab, diag_row)

Compute banded-storage dimensions for the Newton-step Jacobian.

Read more…

Arguments

Type IntentOptional Attributes Name
type(solver_state_t), intent(in) :: state

Solver state (reads n_pt and cfg%lapack_solver).
integer, intent(out) :: n_dof

Total degrees of freedom (neq * n_pt)
integer, intent(out) :: kl

Total degrees of freedom (neq * n_pt) Lower bandwidth
integer, intent(out) :: ku

Total degrees of freedom (neq * n_pt) Lower bandwidth Upper bandwidth
integer, intent(out) :: ldab

Total degrees of freedom (neq * n_pt) Lower bandwidth Upper bandwidth Leading dimension of band matrix
integer, intent(out) :: diag_row

Total degrees of freedom (neq * n_pt) Lower bandwidth Upper bandwidth Leading dimension of band matrix Band-storage row that holds the diagonal

private subroutine band_lapack_solve(ab, kl, ku, n, b, info)

Banded linear solver using LAPACK dgbsv.

Read more…

Arguments

Type IntentOptional Attributes Name
real(kind=wp), intent(inout) :: ab(2*kl+ku+1,n)
integer, intent(in) :: kl

Lower bandwidth
integer, intent(in) :: ku

Lower bandwidth Upper bandwidth
integer, intent(in) :: n

Lower bandwidth Upper bandwidth Matrix order (number of unknowns)
real(kind=wp), intent(inout) :: b(n)
integer, intent(out) :: info

LAPACK return code (0 = success)

private subroutine run_newton(st, ub_n_loc, ub_nm1_loc, c_resid, hist_c1, hist_c2, scheme_label, n_dof_loc, kl_loc, ku_loc, ldab_loc, diag_row_loc, ab_loc, rhs_loc, dq_loc, resid_base_loc, resid_pert_loc)

Shared Newton + banded-solve inner loop for the implicit steppers.

Read more…

Arguments

Type IntentOptional Attributes Name
type(solver_state_t), intent(inout) :: st
real(kind=wp), intent(in) :: ub_n_loc(neq,st%n_pt)
real(kind=wp), intent(in) :: ub_nm1_loc(neq,st%n_pt)
real(kind=wp), intent(in) :: c_resid
real(kind=wp), intent(in) :: hist_c1
real(kind=wp), intent(in) :: hist_c2
character(len=*), intent(in) :: scheme_label
integer, intent(in) :: n_dof_loc
integer, intent(in) :: kl_loc
integer, intent(in) :: ku_loc
integer, intent(in) :: ldab_loc
integer, intent(in) :: diag_row_loc
real(kind=wp), intent(inout) :: ab_loc(ldab_loc,n_dof_loc)
real(kind=wp), intent(inout) :: rhs_loc(n_dof_loc)
real(kind=wp), intent(inout) :: dq_loc(n_dof_loc)
real(kind=wp), intent(inout) :: resid_base_loc(n_dof_loc)
real(kind=wp), intent(inout) :: resid_pert_loc(n_dof_loc)

private subroutine beuler_step_serial_body(state)

Serial (rank-aware-decomp-agnostic) Newton + banded-solve body for backward Euler. Operates entirely on state%ub / state%resid; the caller is responsible for sizing the state correctly. Does NOT call compute_resid_glob — the wrapper does, after restore_local_layout brings every rank back to its local view so the par_sum_real call involves all ranks.

Read more…

Arguments

Type IntentOptional Attributes Name
type(solver_state_t), intent(inout) :: state

public subroutine beuler_step(state)

Advance the solution by one time step using backward (implicit) Euler.

Read more…

Arguments

Type IntentOptional Attributes Name
type(solver_state_t), intent(inout), target :: state

private subroutine bdf2_step_serial_body(state)

Serial Newton + banded-solve body for BDF2 (including the backward-Euler bootstrap path). Operates entirely on state%ub / state%bdf2_ub_prev / state%resid at whatever size state%decomp%n_local says. Does NOT call compute_resid_glob — the wrapper does, from every rank, after restore_local_layout.

Arguments

Type IntentOptional Attributes Name
type(solver_state_t), intent(inout) :: state

public subroutine bdf2_step(state)

Advance the solution by one time step using the BDF2 (Gear) scheme.

Read more…

Arguments

Type IntentOptional Attributes Name
type(solver_state_t), intent(inout), target :: state

private pure subroutine pack_field(src, dst, nq, np)

Flatten a (neq x n_pt) field into a length-neq*n_pt 1-D array. Ordering: (eq=1,cell=1), (eq=2,cell=1), ..., (eq=neq,cell=n_pt).

Arguments

Type IntentOptional Attributes Name
real(kind=wp), intent(in) :: src(nq,np)
real(kind=wp), intent(out) :: dst(nq*np)
integer, intent(in) :: nq
integer, intent(in) :: np

private subroutine unpack_add(src, dst, nq, np)

Add a flattened 1-D correction to a (neq x n_pt) 2-D field in place.

Arguments

Type IntentOptional Attributes Name
real(kind=wp), intent(in) :: src(nq*np)
real(kind=wp), intent(inout) :: dst(nq,np)
integer, intent(in) :: nq
integer, intent(in) :: np

private pure subroutine jac_store_col(ab_loc, j_col, col, kl_loc, ku_loc, diag_row_loc, n)

Store one Jacobian column into band storage AB(diag_row_loc+i-j, j) = J(i,j). Works for both the custom (diag_row_loc = ku+1) and LAPACK (diag_row_loc = kl+ku+1) storage layouts.

Arguments

Type IntentOptional Attributes Name
real(kind=wp), intent(inout) :: ab_loc(:,:)
integer, intent(in) :: j_col
real(kind=wp), intent(in) :: col(n)
integer, intent(in) :: kl_loc
integer, intent(in) :: ku_loc
integer, intent(in) :: diag_row_loc
integer, intent(in) :: n

private subroutine band_lu_solve(ab_loc, n, kl_loc, ku_loc, b, x)

Banded Gaussian elimination without pivoting. Band storage: AB(ku+1+i-j, j) = A(i,j). On exit x = A^{-1} b. The matrix ab_loc is overwritten with LU factors.

Arguments

Type IntentOptional Attributes Name
real(kind=wp), intent(inout) :: ab_loc(kl_loc+ku_loc+1,n)
integer, intent(in) :: n
integer, intent(in) :: kl_loc
integer, intent(in) :: ku_loc
real(kind=wp), intent(in) :: b(n)
real(kind=wp), intent(out) :: x(n)

private subroutine compute_resid_glob(state)

Accumulate the global L2 norm of the residual into state%resid_glob.

Read more…

Arguments

Type IntentOptional Attributes Name
type(solver_state_t), intent(inout) :: state