PoissonSolver

PinCFlow.PoissonSolver — Module

PoissonSolver

Module for solving the Poisson equation and performing a corrector step.

Provides the functions needed to solve the Poisson equation of the semi-implicit time scheme, using a preconditioned BiCGSTAB algorithm, and correcting the Exner-pressure, momentum and density fluctuations accordingly.

Horizontal

Singleton for dispatch to application of the horizontal part of the linear operator.

PinCFlow.PoissonSolver.Total — Type

Total

Singleton for dispatch to application of the full linear operator.

PinCFlow.PoissonSolver.apply_bicgstab! — Function

apply_bicgstab!(state::State, tolref::AbstractFloat)::Tuple{Bool, <:Integer}

Solve the Poisson equation using a preconditioned BiCGSTAB algorithm and return a tuple containing an error flag and the number of iterations.

Arguments

state: Model state.
tolref: Reference tolerance for convergence criterion.

See also

PinCFlow.PoissonSolver.apply_operator!
PinCFlow.PoissonSolver.apply_preconditioner!
PinCFlow.MPIOperations.compute_global_dot_product

PinCFlow.PoissonSolver.apply_corrector! — Function

apply_corrector!(
    state::State,
    dt::AbstractFloat,
    rayleigh_factor::AbstractFloat,
)::Tuple{Bool, <:Integer}

Perform the corrector step and return a tuple containing an error flag and the number of BiCGSTAB iterations.

The left-hand side and the linear operator of the discrete Poisson equation are calculated. The equation is then solved for Exner-pressure differences, using a preconditioned BiCGSTAB algorithm. Finally, the Exner-pressure, wind and density fluctuations are corrected accordingly.

Arguments

state: Model state.
dt: Time step.
rayleigh_factor: Factor by which the Rayleigh-damping coefficient is multiplied.

apply_operator!(
    sin::AbstractArray{<:AbstractFloat, 3},
    ls::AbstractArray{<:AbstractFloat, 3},
    hortot::Total,
    state::State,
)

Apply the total linear operator to the solution array sin.

Before the operator is applied, the boundary/halo values of sin are set, using set_zonal_boundaries_of_field!, set_meridional_boundaries_of_field! and set_vertical_halos_of_field!. Note that in the vertical, only halo values need to be set (if the domain is parallelized in that direction), due to the solid-wall boundaries.

apply_operator!(
    sin::AbstractArray{<:AbstractFloat, 3},
    ls::AbstractArray{<:AbstractFloat, 3},
    hortot::Horizontal,
    state::State,
)

Apply the "horizontal part" of the linear operator (excluding the lower, center and upper diagonal) to the solution array sin.

Before the operator is applied, the boundary/halo values of sin are set, in the same way as in the method applying the total operator.

Arguments

sin: Solution array.
ls: Result of applying the operator to the solution array.
hortot: Linear-operator mode.
state: Model state.

apply_preconditioner!(
    sin::AbstractArray{<:AbstractFloat, 3},
    sout::AbstractArray{<:AbstractFloat, 3},
    state::State,
)

Apply a preconditioner to the Poisson problem.

This preconditioner integrates the auxiliary equation

\[\frac{\mathrm{d} s}{\mathrm{d} \eta} = \mathcal{L}_\mathrm{h} \left(s\right) + \mathcal{L}_\mathrm{v} \left(s\right) - b,\]

where $s$ is the iterative solution, $\eta$ is a pseudo-time variable, $\mathcal{L}_\mathrm{v}$ contains the lower, center and upper diagonals of the linear operator, $\mathcal{L}_\mathrm{h}$ contains all remaining elements, and $b$ is the left-hand side. The integration is performed in a semi-implicit manner, following

\[\left(1 - \Delta \eta \mathcal{L}_\mathrm{v}\right) \left(s^{\left(m + 1\right)}\right) = \left(1 + \Delta \eta \mathcal{L}_\mathrm{h}\right) \left(s^{\left(m\right)}\right) - \Delta \eta b,\]

where $\Delta \eta = \Delta \tau / 2 \left[\left(\Delta \widehat{x}\right)^{- 2} + \left(\Delta \widehat{y}\right)^{- 2}\right]^{- 1}$, with $\Delta \tau$ being a namelist parameter (state.namelist.poisson.dtau). Therein, the implicit problem is solved with the Thomas algorithm for tridiagonal matrices. The number of iterations is given by state.namelist.poisson.preconditioner_iterations. Since the Thomas algorithm consists of an upward elimination sweep and a downward pass, this method performs sequential one-way MPI communication if the domain is parallelized in the vertical.

Arguments

sin: Residual array.
sout: Solution of the preconditioner.
state: Model state.

See also

PinCFlow.PoissonSolver.apply_operator!

PinCFlow.PoissonSolver.compute_lhs! — Function

compute_lhs!(state::State)::AbstractFloat

Compute the scaled left-hand side of the Poisson equation and return a reference tolerance for the convergence criterion by dispatching to a model-specific method.

compute_lhs!(
    state::State,
    model::Union{Boussinesq, PseudoIncompressible},
)::AbstractFloat

Compute the scaled left-hand side of the Poisson equation in pseudo-incompressible/Boussinesq mode and return a reference tolerance for the convergence criterion.

The scaled left-hand side is given by

\[\begin{align*} b & = \frac{\sqrt{\overline{\rho}}}{P} \frac{1}{J c_p} \left\{\frac{1}{\Delta \widehat{x}} \left[\left(J P\right)_{i + 1 / 2} u_{i + 1 / 2} - \left(J P\right)_{i - 1 / 2} u_{i - 1 / 2}\right]\right.\\ & \qquad \qquad \quad + \frac{1}{\Delta \widehat{y}} \left[\left(J P\right)_{j + 1 / 2} v_{j + 1 / 2} - \left(J P\right)_{j - 1 / 2} v_{j - 1 / 2}\right]\\ & \qquad \qquad \quad + \left.\frac{1}{\Delta \widehat{z}} \left[\left(J P\right)_{k + 1 / 2} \widehat{w}_{k + 1 / 2} - \left(J P\right)_{k - 1 / 2} \widehat{w}_{k - 1 / 2}\right]\right\} \end{align*}\]

and the reference tolerance is given by

\[\tau_\mathrm{ref} = \frac{\sum_{i, j, k} b_{i, j, k}^2}{\sum_{i, j, k} \left(b_{u, i, j, k}^2 + b_{v, i, j, k}^2 + b_{\widehat{w}, i, j, k}^2\right)},\]

where $b_u$, $b_v$ and $b_{\widehat{w}}$ are the zonal, meridional and vertical parts of $b$, respectively. Note that in Boussinesq mode, $P = P_0$ will cancel out, so that the appropriate divergence constraint remains.

compute_lhs!(state::State, model::Compressible)::AbstractFloat

Compute the scaled left-hand side of the Poisson equation in compressible mode and return a reference tolerance for the convergence criterion.

The scaled left-hand side is given by

\[\begin{align*} b & = \frac{\sqrt{\overline{\rho}}}{P} \frac{1}{J c_p} \left(\frac{U_{i + 1 / 2} - U_{i - 1 / 2}}{\Delta \widehat{x}} + \frac{V_{j + 1 / 2} - V_{j - 1 / 2}}{\Delta \widehat{y}} + \frac{\widehat{W}_{k + 1 / 2} - \widehat{W}_{k - 1 / 2}}{\Delta \widehat{z}}\right) - \frac{\sqrt{\overline{\rho}}}{P} F^P, \end{align*}\]

where $F^P$ is the diabatic heating, as computed by compute_volume_force, and the reference tolerance is computed in the same way as in the method for pseudo-incompressible/Boussinesq mode, with $b_{F, i, j, k}^2$ added to the sum in the denominator, representing the heating term.

Arguments

state: Model state.
model: Dynamic equations.

See also

PinCFlow.Update.compute_volume_force

PinCFlow.PoissonSolver.compute_operator! — Function

compute_operator!(
    state::State,
    dt::AbstractFloat,
    rayleigh_factor::AbstractFloat,
)

Compute the tensor elements of the linear operator on the right-hand side of the Poisson equation.

The operator is obtained by rewriting the scaled Poisson equation

\[\frac{\sqrt{\overline{\rho}}}{P} \mathrm{LHS} = \frac{\sqrt{\overline{\rho}}}{P} \mathrm{RHS} \left(\frac{\sqrt{\overline{\rho}}}{P} s\right)\]

\[\frac{\sqrt{\overline{\rho}}}{P} \mathrm{LHS} = \sum_{\lambda, \mu, \nu} A_{i + \lambda, j + \mu, k + \nu} s_{i + \lambda, j + \mu, k + \nu},\]

where the Exner-pressure differences are given by $\Delta \pi = \left(\sqrt{\overline{\rho}} / P\right) \left(s / \Delta t\right)$.

Arguments

state: Model state.
dt: Time step.
rayleigh_factor: Factor by which the Rayleigh-damping coefficient is multiplied.

See also

PinCFlow.Update.compute_buoyancy_factor

PinCFlow.PoissonSolver.correct! — Function

correct!(state::State, dt::AbstractFloat, rayleigh_factor::AbstractFloat)

Correct the Exner-pressure, wind and density fluctuations such that the divergence constraint is satisfied, using the Exner-pressure differences obtained from the solution to the Poisson problem.

This method calls specific methods for the correction of each individual variable.

correct!(
    state::State,
    dt::AbstractFloat,
    variable::U,
    rayleigh_factor::AbstractFloat,
)

Correct the zonal wind to account for the pressure differences obtained from the solution to the Poisson problem.

The correction is given by

\[u_{i + 1 / 2} \rightarrow u_{i + 1 / 2} - \mathcal{C}_{i + 1 / 2}^{\rho u}\]

in Boussinesq/pseudo-incompressible mode and

\[U_{i + 1 / 2} \rightarrow U_{i + 1 / 2} - \left(J P\right)_{i + 1 / 2} \mathcal{C}_{i + 1 / 2}^{\rho u}\]

in compressible mode, with

\[\mathcal{C}_{i + 1 / 2}^{\rho u} = \left(1 + \beta_{\mathrm{R}, i + 1 / 2} \Delta t\right)^{- 1} \Delta t c_p \frac{P_{i + 1 / 2}}{\rho_{i + 1 / 2}} \mathcal{D}_{i + 1 / 2}^{\rho u}.\]

Therein, $\Delta t$ is the fractional time step given as input to this method and $c_p \left(P_{i + 1 / 2} / \rho_{i + 1 / 2}\right) \mathcal{D}_{i + 1 / 2}^{\rho u}$ is computed with compute_pressure_gradient.

correct!(
    state::State,
    dt::AbstractFloat,
    variable::V,
    rayleigh_factor::AbstractFloat,
)

Correct the meridional wind to account for the pressure differences obtained from the solution to the Poisson problem.

The correction is given by

\[v_{j + 1 / 2} \rightarrow v_{j + 1 / 2} - \mathcal{C}_{j + 1 / 2}^{\rho v}\]

in Boussinesq/pseudo-incompressible mode and

\[V_{j + 1 / 2} \rightarrow V_{j + 1 / 2} - \left(J P\right)_{j + 1 / 2} \mathcal{C}_{j + 1 / 2}^{\rho v}\]

in compressible mode, with

\[\mathcal{C}_{j + 1 / 2}^{\rho v} = \left(1 + \beta_{\mathrm{R}, j + 1 / 2} \Delta t\right)^{- 1} \Delta t c_p \frac{P_{j + 1 / 2}}{\rho_{j + 1 / 2}} \mathcal{D}_{j + 1 / 2}^{\rho v},\]

where $c_p \left(P_{j + 1 / 2} / \rho_{j + 1 / 2}\right) \mathcal{D}_{j + 1 / 2}^{\rho v}$ is computed with compute_pressure_gradient.

correct!(
    state::State,
    dt::AbstractFloat,
    variable::W,
    rayleigh_factor::AbstractFloat,
)

Correct the transformed vertical wind to account for the pressure differences obtained from the solution to the Poisson problem.

The correction is given by

\[\begin{align*} \widehat{w}_{k + 1 / 2} & \rightarrow \widehat{w}_{k + 1 / 2} - \left[1 + \beta_{\mathrm{R}, k + 1 / 2} \Delta t + \frac{\overline{\rho}_{k + 1 / 2}}{\rho_{k + 1 / 2}} N_{k + 1 / 2}^2 \left(\Delta t\right)^2\right]^{- 1}\\ & \quad \times \left\{\Delta t c_p \frac{P_{k + 1 / 2}}{\rho_{k + 1 / 2}} \mathcal{D}_{k + 1 / 2}^{\rho \widehat{w}} + \frac{\overline{\rho}_{k + 1 / 2}}{\rho_{k + 1 / 2}} N_{k + 1 / 2}^2 \left(\Delta t\right)^2\right.\\ & \qquad \quad \times \left.\left[\left(G^{1 3} \mathcal{C}^{\rho u}\right)_{k + 1 / 2} + \left(G^{23} \mathcal{C}^{\rho v}\right)_{k + 1 / 2}\right]\vphantom{\frac{P_{k + 1 / 2}}{\rho_{k + 1 / 2}}}\right\}, \end{align*}\]

in Boussinesq/pseudo-incompressible mode and

\[\begin{align*} \widehat{W}_{k + 1 / 2} & \rightarrow \widehat{W}_{k + 1 / 2} - \left[1 + \beta_{\mathrm{R}, k + 1 / 2} \Delta t + \frac{\left(P / \overline{\theta}\right)_{k + 1 / 2}}{\rho_{k + 1 / 2}} N_{k + 1 / 2}^2 \left(\Delta t\right)^2\right]^{- 1}\\ & \quad \times \left\{\Delta t c_p \left(J P\right)_{k + 1 / 2} \frac{P_{k + 1 / 2}}{\rho_{k + 1 / 2}} \mathcal{D}_{k + 1 / 2}^{\rho \widehat{w}} \vphantom{\frac{\left(P / \overline{\theta}\right)_{k + 1 / 2}}{\rho_{k + 1 / 2}}} + \left(J P\right)_{k + 1 / 2} \frac{\left(P / \overline{\theta}\right)_{k + 1 / 2}}{\rho_{k + 1 / 2}} N_{k + 1 / 2}^2 \left(\Delta t\right)^2\right.\\ & \qquad \quad \times \left.\left[\left(G^{1 3} \mathcal{C}^{\rho u}\right)_{k + 1 / 2} + \left(G^{23} \mathcal{C}^{\rho v}\right)_{k + 1 / 2}\right]\vphantom{\frac{\left(P / \overline{\theta}\right)_{k + 1 / 2}}{\rho_{k + 1 / 2}}}\right\}, \end{align*}\]

in compressible mode, where $c_p \left(P_{k + 1 / 2} / \rho_{k + 1 / 2}\right) \mathcal{D}_{k + 1 / 2}^{\rho \widehat{w}}$ is computed with compute_pressure_gradient.

correct!(
    state::State,
    dt::AbstractFloat,
    variable::RhoP,
    rayleigh_factor::AbstractFloat,
)

Correct the density fluctuations to account for the pressure differences obtained from the solution to the Poisson problem.

The correction is given by

\[\begin{align*} \rho' & \rightarrow \rho' + \frac{\rho}{g} \left[1 + \beta_\mathrm{R} \Delta t + \frac{\overline{\rho}}{\rho} \left(N \Delta t\right)^2\right]^{- 1}\\ & \quad \times \left[- \frac{\overline{\rho}}{\rho} \left(N \Delta t\right)^2 J \left(c_p \frac{P_{k + 1 / 2}}{\rho_{k + 1 / 2}} \mathcal{D}_{k + 1 / 2}^{\rho \widehat{w}}\right)\right.\\ & \qquad \quad + \left.\frac{\overline{\rho}}{\rho} N^2 \Delta t J \left(1 + \beta_\mathrm{R} \Delta t\right) \left(G^{1 3} \mathcal{C}^{\rho u} + G^{2 3} \mathcal{C}^{\rho v}\right)\right], \end{align*}\]

in Boussinesq/pseudo-incompressible mode and

\[\begin{align*} \rho' & \rightarrow \rho' + \frac{\rho}{g} \left[1 + \beta_\mathrm{R} \Delta t + \frac{P / \overline{\theta}}{\rho} \left(N \Delta t\right)^2\right]^{- 1}\\ & \quad \times \left[- \frac{P / \overline{\theta}}{\rho} \left(N \Delta t\right)^2 J \left(c_p \frac{P_{k + 1 / 2}}{\rho_{k + 1 / 2}} \mathcal{D}_{k + 1 / 2}^{\rho \widehat{w}}\right)\right.\\ & \qquad \quad + \left.\frac{P / \overline{\theta}}{\rho} N^2 \Delta t J \left(1 + \beta_\mathrm{R} \Delta t\right) \left(G^{1 3} \mathcal{C}^{\rho u} + G^{2 3} \mathcal{C}^{\rho v}\right)\right], \end{align*}\]

in compressible mode, where $c_p \left(P_{k + 1 / 2} / \rho_{k + 1 / 2}\right) \mathcal{D}_{k + 1 / 2}^{\rho \widehat{w}}$ and $c_p \left(P_{k - 1 / 2} / \rho_{k - 1 / 2}\right) \mathcal{D}_{k - 1 / 2}^{\rho \widehat{w}}$ are computed with compute_pressure_gradient, and used to interpolate to $\left(i, j, k\right)$.

correct!(state::State, variable::PiP)

Update the Exner-pressure fluctuations with the differences obtained from the solution to the Poisson problem.

Arguments

state: Model state.
dt: Time step.
variable: Variable to correct.
rayleigh_factor: Factor by which the Rayleigh-damping coefficient is multiplied.

solve_poisson!(
    state::State,
    dt::AbstractFloat,
    rayleigh_factor::AbstractFloat,
    tolref::AbstractFloat,
)::Tuple{Bool, <:Integer}

Solve the Poisson equation and return a tuple containing an error flag and the number of iterations.

Given a left-hand side and reference tolerance, this method computes the elements of the linear operator and solves the Poisson equation, using a preconditioned BiCGSTAB algorithm. Both the Exner-pressure differences and the entire equation are scaled with $\sqrt{\overline{\rho}} / P$ in advance (the left-hand side has already been scaled at this point), so that the equation

\[\frac{\sqrt{\overline{\rho}}}{P} \mathrm{LHS} = \frac{\sqrt{\overline{\rho}}}{P} \mathrm{RHS} \left(\frac{\sqrt{\overline{\rho}}}{P} s\right)\]

is solved for $s$. The Exner-pressure differences are then given by $\Delta \pi = \left(\sqrt{\overline{\rho}} / P\right) \left(s / \Delta t\right)$.

Arguments

state: Model state.
dt: Time step.
rayleigh_factor: Factor by which the Rayleigh-damping coefficient is multiplied.
tolref: Reference tolerance for convergence criterion.

See also

PinCFlow.PoissonSolver.compute_operator!
PinCFlow.PoissonSolver.apply_bicgstab!