Solvers

The algorithm choices for the simulation environment are set using the chidg%set method with a selection parameter and a string option:

call chidg%set('Time Integrator', String)
call chidg%set('Nonlinear Solver',String)
call chidg%set('Linear Solver',   String)
call chidg%set('Preconditioner',  String)

Time Integrators

Algorithm	String Selector
Steady	Steady
Backward Euler	Backward Euler
Diagonally-Implicit Runge-Kutta	DIRK
Harmonic Balance	Harmonic Balance
Forward Euler	Forward Euler
Explicit Runge Kutta	Runge-Kutta Method

ChiDG Representation

Integration of the temporal derivative term

\[\int \psi \frac{\partial Q}{\partial t} d\Omega\]

in

\[\int_{\Omega_e} \boldsymbol{\psi} \frac{\partial Q}{\partial t} d\Omega + \int_{\Gamma_e} \boldsymbol{\psi} \vec{F} \cdot \vec{n} d\Gamma - \int_{\Omega_e} \nabla \boldsymbol{\psi} \cdot \vec{F} d\Omega + \int_{\Omega_e} \boldsymbol{\psi} S d\Omega = 0\]

is handled by the time_integrator_t class.

type, abstract, public :: time_integrator_t

contains

    procedure(data_interface), deferred :: step     ! Must define this procedure in the extended type

end type time_integrator_t

Note that in the terminology used in this section, \(M\) is the Mass Matrix

\[\boldsymbol{M} = \int_{\Omega} \boldsymbol{\psi} \boldsymbol{\psi} d\Omega\]

and \(R(Q)\) is the right-hand-side vector of spatial operators

\[R(\hat{Q}) = \int_{\Gamma_e} \boldsymbol{\psi} \vec{F} \cdot \vec{n} d\Gamma - \int_{\Omega_e} \nabla \boldsymbol{\psi} \cdot \vec{F} d\Omega + \int_{\Omega_e} \boldsymbol{\psi} S d\Omega\]

Implicit integrators

Steady

\[\int \psi \frac{\partial Q}{\partial t} d\Omega = 0\]

Nonlinear problem:

\[R(Q) = 0\]

Newton Linearization:

\[\frac{\partial R(\hat{Q}^{n})}{\partial Q} \Delta Q = -R(\hat{Q}^{n})\]

Backward Euler

Solution advancement via first-order backward difference of the time derivative:

\[\boldsymbol{M} \frac{\hat{Q}^{n+1} - \hat{Q}^{n}}{\Delta t} = R(\hat{Q}^{n+1})\]

Nonlinear problem:

\[\frac{\Delta \hat{Q}}{\Delta t}\boldsymbol{M} + R(\hat{Q}^{n+1}) = 0\]

Newton Linearization:

\[\begin{split}\bigg(\frac{\boldsymbol{M}}{\Delta t} + \frac{\partial R(\hat{Q}^{m})}{\partial Q}\bigg) \delta Q^{m} & = -\frac{\boldsymbol{M}}{\Delta t}\Delta Q^{m} -R(\hat{Q}^{m})\\ \hat{Q}^{m} & = \hat{Q}^{n} + \Delta Q^{m}\end{split}\]

where, \(\delta Q^{m} = \Delta Q^{m + 1} -\Delta Q^{m}\) for the \(m^{th}\) Newton iteration.

Note

• The ChiDG nonlinear solver requires the assembly of the linearized system outlined above
• The nonlinear solver computes \(\hat{Q}^{m}\) instead of \(\Delta \hat{Q}^{m}\)
• For the assembly of the rhs of the linearized system, \(\Delta \hat{Q}^{m} = \hat{Q}^{m} - \hat{Q}^{n}\)

Diagonally-Implicit Runge-Kutta

\[\boldsymbol{M}\frac{\partial Q}{\partial t} + R(Q) = 0\]

With the coefficient arrays associated with the diagonally implicit Runge-Kutta method:

\[\begin{split}\boldsymbol{A} = \left[\begin{array}{ccc} \gamma & 0 & 0 \\ \tau_{2} - \gamma & \gamma & 0 \\ b_{1} & b_{2} & \gamma \end{array} \right]\end{split}\]

\[\boldsymbol{b} = \left[\begin{array}{ccc} b_{1} & b_{2} & \gamma \end{array} \right]\]

where \(\gamma\) is the root of \(x^{3} - 3x^{2} + \frac{3}{2}x - \frac{1}{6} = 0 \in \left(\frac{1}{6},\frac{1}{2}\right)\) and

\[\begin{split}\tau_{2} & = (1 + \gamma)/2\\ b_{1} & = -(6\gamma^{2} - 16\gamma + 1)/4\\ b_{2} & = (6\gamma^{2} - 20\gamma + 5)/4\end{split}\]

The solution is advanced in time as:

\[\hat{Q}^{n + 1} = \hat{Q}^{n} + b_{1}\Delta \hat{Q}_{1} + b_{2}\Delta \hat{Q}_{2} + b_{3} \Delta \hat{Q}_{3}\]

Implicit system:

\[\frac{\Delta \hat{Q}_{i}}{\Delta t}\boldsymbol{M} = -R\left(\hat{Q}^{n} + \sum_{j = 1}^{i}A_{ij}\Delta \hat{Q}_{i}\right)\;\;\;\text{for}\;i = 1,3\]

Newton linearization:

\[\left(\boldsymbol{M} + \gamma \Delta t \frac{\partial R(\hat{Q}^{m}_{i})}{\partial Q}\right)\delta \hat{Q}^{m}_{i} = -\boldsymbol{M}\Delta \hat{Q}^{m}_{i} - \Delta t R\left(\hat{Q}^{m}_{i}\right)\;\;\;\text{for}\;i = 1,3\]

with

\[\begin{split}\hat{Q}^{m}_{i} & = \hat{Q}^{n} + \sum_{j = 1}^{i - 1}A_{ij}\Delta \hat{Q}_{i} + \gamma \Delta \hat{Q}^{m}_{i}\\ \delta \hat{Q}^{m}_{i} & = \Delta \hat{Q}^{m + 1}_{i} - \Delta \hat{Q}^{m}_{i}\end{split}\]

Note

• The ChiDG nonlinear solver requires the assembly of the stagewise linearized systems
• The nonlinear solver computes \(\hat{Q}^{m}_{i}\) instead of \(\Delta \hat{Q}^{m}_{i}\)
• For the assembly of the rhs of the linearized system, \(\Delta \hat{Q}^{m}_{i} = (\hat{Q}^{m}_{i} - \hat{Q}^{n} - \sum_{j = 1}^{i - 1}A_{ij}\Delta \hat{Q}_{i})/\gamma\)

Harmonic Balance

Consider a set of \(N\) independent equations:

\[\frac{\partial \hat{\boldsymbol{Q}}^{*}}{\partial t} + \nabla \cdot \boldsymbol{F}^{*} + \boldsymbol{S}^{*} = 0\]

where

\[\begin{split}\hat{\boldsymbol{Q}^{*}} & = \left[\hat{Q}_{1}, \hat{Q}_{2}, \cdots, \hat{Q}^{N}\right]^{T}\\ \boldsymbol{F}^{*} & = \left[\vec{F}(\hat{Q}_{1}, \nabla \hat{Q}_{1}), \vec{F}(\hat{Q}_{2}, \nabla \hat{Q}_{2}), \cdots, \vec{F}(\hat{Q}_{N}, \nabla \hat{Q}_{N})\right]^{T}\\ \boldsymbol{S}^{*} & = \left[S(\hat{Q}_{1}, \nabla \hat{Q}_{1}), S(\hat{Q}_{2}, \nabla \hat{Q}_{2}), \cdots, S(\hat{Q}_{N}, \nabla \hat{Q}_{N})\right]^{T}\end{split}\]

In the harmonic balance method, a conservative solution vector at any instant of time is represented as a Fourier series in time as:

\[\hat{Q}_{n} = A_{0} + \sum_{k = 1}^{K}\left[A_{k}\text{sin}(\omega_{k}t_{n}) + B_{k}\text{cos}(\omega_{k}t_{n})\right]\]

with \(K\) frequencies, \(\boldsymbol{\omega} = [\omega_{1}, \omega_{2}, \cdots, \omega_{K}]\) and the instant of time \(t_{n}\) belongs to the set of time levels, \(\boldsymbol{t} = [t_{1}, t_{2}, \cdots, t_{N}]\) with \(N = 2K + 1\). Thus, the series of conservative solution vectors can be related to the Fourier coefficients vectors, \(\hat{\boldsymbol{Q}}_{F}\) as:

\[\hat{\boldsymbol{Q}}^{*} = E^{-1}\hat{\boldsymbol{Q}}_{F}\]

Defining the pseudo spectral operator as,

\[D = \frac{\partial E^{-1}}{\partial t}E\]

which couples \(\hat{\boldsymbol{Q}}^{*}\) such that the conservative solutions satisfy time-varying sinusoidal functions according to their Fourier representation, the governing equation can be rewritten as the Harmonic Balance equation:

\[D\hat{\boldsymbol{Q}}^{*} + \nabla \boldsymbol{F}^{*} + S^{*} = 0\]

Multiplying with a column of test functions, \(\psi\) and applying Gauss’ divergence theorem provides the working form of the Harmonic Balance equation:

\[\int_{\Omega_{e}}\psi D\hat{\boldsymbol{Q}}^{*}d\Omega + \int_{\Gamma_{e}}\boldsymbol{F}^{*} \cdot \vec{n}d\Gamma - \int_{\Omega_{e}}\nabla \psi \cdot \boldsymbol{F}^{*}d\Omega + \int_{\Omega_{e}}\psi \boldsymbol{S}^{*}d\Omega = 0\]

Newton Linearization:

Consider,

\[\begin{split}\mathscr{R}^{*} & = \int_{\Gamma_{e}}\boldsymbol{F}^{*} \cdot \vec{n}d\Gamma - \int_{\Omega_{e}}\nabla \psi \cdot \boldsymbol{F}^{*}d\Omega + \int_{\Omega_{e}}\psi \boldsymbol{S}^{*}d\Omega\\ \mathscr{D}^{*} & = \int_{\Omega_{e}}\psi D\hat{\boldsymbol{Q}}^{*}d\Omega\end{split}\]

Then, Newton linearization of the Harmonic Balance system of equations is:

\[\left(\frac{\partial \mathscr{D}^{*}}{\partial \hat{\boldsymbol{Q}}^{*}} + \frac{\partial \mathscr{R}^{*}}{\partial \hat{\boldsymbol{Q}}^{*}}\right)\Delta \hat{\boldsymbol{Q}^{*}} = -(\mathscr{D}^{*} + \mathscr{R}^{*})\]

Explicit integrators

Forward Euler

Solution advancement via a first-order forward difference of the time derivative:

\[\boldsymbol{M} \frac{\hat{Q}^{n+1} - \hat{Q}^{n}}{\Delta t} = R(\hat{Q}^{n})\]

Algebraic problem:

\[\frac{\Delta \hat{Q}}{\Delta t}\boldsymbol{M} + R(\hat{Q}^{n}) = 0\]

Solution iterated in time as:

\[\hat{Q}^{n+1} = \hat{Q}^n - {\Delta t} \boldsymbol{M}^{-1}R(\hat{Q}^{n})\]

Explicit Runge Kutta

For a general explicit runge Kutta method with \(s\) stages:

\[\boldsymbol{M}\frac{\hat{Q}^{n + 1} - \hat{Q}^{n}}{\Delta t} = \sum_{i = 1}^{s}b_{i}\Delta \hat{Q}_{i}\]

where

\[\Delta \hat{Q}_{i} = -R\left(\hat{Q}^{n} + \sum_{j = 1}^{i - 1}a_{ij}\Delta \hat{Q}_{j}\right)\]

Algebraic problem:

\[\frac{\Delta \hat{Q}}{\Delta t}\boldsymbol{M} - \sum_{i = 1}^{s}b_{i}\Delta \hat{Q}_{i} = 0\]

Solution iterated in time as:

\[\hat{Q}^{n + 1} = \hat{Q}^{n} + \Delta t \boldsymbol{M}^{-1}\sum_{i = 1}^{s}b_{i}\Delta \hat{Q}_{i}\]

Nonlinear solvers

Algorithm	String Selector
Newton	Newton
Quasi-Newton	Quasi-Newton

ChiDG Representation

ChiDG includes nonlinear solvers for solving the nonlinear sets of partial differential equations. In general, the implicit problem statement here is:

• Find \(Q\), such that \(\mathcal{R}(Q) = 0\)

type, abstract, public :: nonlinear_solver_t

contains

    procedure(data_interface), deferred :: solve   ! Must define this procedure in the extended type

end type nonlinear_solver_t

Newton

The Full-Newton solver solves the equation

\[\mathcal{R}(Q) = 0\]

where \(\mathcal{R}(Q)\) is some potentially nonlinear function of the solution. This depends on the discretization, the equation set, the solution order, and the time-integration scheme. The Newton solver linearizes the problem and computes an update of the solution by solving

\[\frac{\partial \mathcal{R}}{\partial Q} \Delta Q = -\mathcal{R}\]

So, at each Newton step, a linear system of equations is being solved for \(\Delta Q\). Once the update is solved for, the solution vector is updated as

\[Q^{n+1} = Q^{n} + \Delta Q\]

Considerations:

One item to consider when using the Full-Newton solver is that the Newton linearization(direction) is dependent on the current solution. Without a reasonable initial guess, Newton’s method can diverge by sending the solution too far in the wrong direction.

Quasi-Newton

The Quasi-Newton solver solves a modified set of equations

\[\int_{\Omega_e} \psi \frac{\partial Q}{\partial \tau} d\Omega + \mathcal{R}(Q) = 0\]

Note the addition of a pseudo-time term to the nonlinear system of equations. This is an effort increase robustness of the nonlinear solver by limiting the size of the solution update in a single Newton step. This is accomplished by adding the time-scaling to the diagonal of the Jacobian matrix, increasing the diagonal dominance of the matrix, and limiting the size of the soltion update. As the solution progresses, the pseudo-timestep is increased the pseudo-time derivative goes to zero and the original system of equations is recovered.

The Quasi-Newton solver linearizes the problem including the pseudo-time scaling of the solution update to the system of equations as

\[\int_{\Omega_e} \psi \frac{\Delta Q}{\Delta \tau} d\Omega + \frac{\partial \mathcal{R}}{\partial Q} \Delta Q = -\mathcal{R}\]

So, at each Quasi-Newton step, a linear system of equations is being solved for \(\Delta Q\). Once the update is solved for, the solution vector is updated as

\[Q^{n+1} = Q^{n} + \Delta Q\]

At each Quasi-Newton step, the pseudo-time step is updated as

\[d\tau = \frac{CFL^n h_e}{\bar{\lambda_e}}\]

where \(h_e = \sqrt[3]{\Omega_e}\) and \(\bar{\lambda_e} = |\bar{V_e}| + \bar{c}\) is a mean characteristic speed. The CFL term is computed from the ratio of the initial and current residual norms as

\[CFL^n = CFL^0 \frac{||\mathcal{R}(Q^0)||_2}{||\mathcal{R}(Q^n)||_2}\]

Options:

• CFL0: The initial CFL factor

Linear Solvers

Algorithm	String Selector
Flexible Generalized Minimum Residual	FGMRES

ChiDG Representation

Flexible Generalized Minimum Residual

A flexible version of the Generalized Minimum Residual(FGMRES) algorithm, which is an iterative method for solving linear systems of equations. The FGMRES algorithm allows the GMRES algorithm to be preconditioned in a flexible way such that the solution can be easily reconstructed.

Options:

• m: Number of iterations before the algorithm is restarted

Preconditioners

Algorithm	String Selector
block-Jacobi	Jacobi
block-ILU0	ILU0
Restricted Additive Schwarz + block-ILU0	RAS+ILU0

ChiDG Representation

block-Jacobi

block-ILU0

Restricted Additive Schwarz + block-ILU0