 |
TRACE User Guide
TRACE Version 9.6.1
|
|
TRACE User Guide
TRACE Version 9.6.1
|
This page summarizes the implicit time-integration algorithms available in TRACE. The algorithms are derived for the following generic scalar equation
\begin{equation}\label{Advec1D} \frac{d u}{d t} = R(u) \end{equation}
where \(R(u)\) may be a linear or non-linear function of the scalar quantity \(u\).
A stationary solution to Eqn. Advec1D corresponds to \(u(x)\) such that \(R(u) = 0\). A simple iterative method for finding such roots is Newton's method. Newton's method can be derived by expanding the residual \(R\) in a Taylor series and neglecting terms of order \(\Delta u^2\) and higher
\begin{equation}\label{NewtonMethod01} R(u^n + \Delta u) \approx R^n + \Delta u \left.\frac{dR}{du}\right|_{n} + O(\Delta u^2) \end{equation}
Setting \(R(u^n + \Delta u) = 0\), Eqn. NewtonMethod01 may be rearranged for \(\Delta u\) to obtain
\begin{equation}\label{NewtonMethod021} \Delta u = -\frac{R^n}{\left.\frac{dR}{du}\right|_{n}} \end{equation}
where \(\Delta u = u^{n+1} - u^n\). If \(u^0\) corresponds to an initial guess, Eqn. NewtonMethod021 can be used to obtain an improved solution \(u^{n+1}\) that lies more closely to the true solution \(u\) corresponding to R(u) = 0. For additional robustness a relaxation parameter \(\alpha \geq 1\) may be included in Eqn. NewtonMethod021 as follows
\begin{equation}\label{NewtonMethod03} \Delta u = -\frac{1}{\alpha}\frac{R^n}{\left.\frac{dR}{du}\right|_{n}} \end{equation}
Higher-order methods can be derived by retaining more terms in the Taylor series expansion of the residual \(R\). The method corresponding to the Predictor-Corrector (PC) scheme in TRACE can be derived by retaining the second-order term in the Taylor series expansion and neglecting only those terms of order \(\Delta u^3\) and higher
\begin{eqnarray}\label{ModifiedNewtonMethod01} R(u^n + \Delta u) &=& R^n + \Delta u \left.\frac{dR}{du}\right|_{n} + \frac{\Delta u^2}{2!} \left.\frac{d^2R}{du^2}\right|_{n}+ O(\Delta u^3) \nonumber \\ &=& R^n + \Delta u \left.\frac{dR}{du}\right|_{n} + \frac{\Delta u^2}{2!\Delta u} \left(\left.\frac{dR}{du}\right|_{n+1} - \left.\frac{dR}{du}\right|_{n}\right)+ O(\Delta u^3) \nonumber \\ &=& R^n + \Delta u \left.\frac{dR}{du}\right|_{n} + \frac{\Delta u^2}{2!\Delta u} \left(\frac{R^{n+1} - R^n}{\Delta u} - \left.\frac{dR}{du}\right|_{n}\right)+ O(\Delta u^2) \nonumber \\ &=& \frac{1}{2}\left(R^{n+1} + R^n\right) + \frac{\Delta u}{2} \left.\frac{dR}{du}\right|_{n} + O(\Delta u^2) \end{eqnarray}
Again, setting \(R(u^n + \Delta u) = 0\), Eqn. ModifiedNewtonMethod01 may be rearranged for \(\Delta u\) to obtain
\begin{equation}\label{ModifiedNewtonMethod02} \Delta u = -\frac{1}{\alpha}\frac{\frac{1}{2}\left(R^{n+1} + R^n\right)}{\frac{1}{2}\left.\frac{dR}{du}\right|_{n}} \end{equation}
where \(\alpha\) is a relaxation parameter. Of course, \(R^{n+1}\) is unknown and must be approximated. This can be achieved using a predictor step based on the standard Newton's method. That is, using the standard Newton method we approximate \(u^{n+1}\) as \(\hat{u}^{n+1} = u^n + \Delta \hat{u}\) where \(\Delta \hat{u}\) is computed by solving
\begin{equation}\label{ModifiedNewtonMethod03a} \Delta \hat{u} = -\frac{R^n}{\left.\frac{dR}{du}\right|_{n}} \end{equation}
The value of \(u^{n+1}\) is then computed by solving
\begin{equation}\label{ModifiedNewtonMethod02a} \Delta u = -\frac{1}{\alpha}\frac{\frac{1}{2}(\hat{R}^{n+1} + R^n)}{\frac{1}{2}\left.\frac{dR}{du}\right|_{n}} \end{equation}
where \(\hat{R}^{n+1} = R(\hat{u}^{n+1})\).
An alternative approach to Newton's method is the pseudo-time or pseudo-transient method. In this approach the actual steady problem is treated as an unsteady problem and integrated in time until a steady state solution is obtained. As we are only interested in the steady solution time accuracy is unimportant and a first-order accurate time-integration scheme is sufficient. The method used in TRACE is based on the Euler-backwards method:
\begin{equation}\label{EulerBackwards} u^{n+1} = u^n + \Delta t R^{n+1} \end{equation}
As \(R^{n+1}\) is unknown it must be approximated. This is accomplished using a Taylor series expansion, for which we obtain
\begin{equation}\label{TaylorExp} R^{n+1} \approx R^n + \frac{\partial R}{\partial t}\Delta t + \frac{\partial^2 R}{\partial t^2}\frac{\Delta t^2}{2!} + \frac{\partial^3 R}{\partial t^3}\frac{\Delta t^3}{3!} + O(\Delta t^4) \end{equation}
Conventionally, those terms of order \(\Delta t^2\) and higher are dropped in the expansion and the first-order derivative is re-written
\begin{eqnarray} \frac{\partial R}{\partial t}\Delta t &=& \frac{\partial R}{\partial u}\frac{\partial u}{\partial t}\Delta t \\ &=& \frac{\partial R}{\partial u}\frac{\Delta u}{\Delta t}\Delta t \\ &=& \frac{\partial R}{\partial u}\Delta u \end{eqnarray}
where \(\Delta u = u^{n+1} - u^{n}\). The expression for the residual at time-level \(n+1\) therefore becomes
\begin{equation}\label{TaylorExpApprox} R^{n+1} \approx R^n + \frac{\partial R}{\partial u}\Delta u + O(\Delta t^2) \end{equation}
Substituting Eqn. TaylorExpApprox in Eqn. EulerBackwards the following expression is obtained for \(\Delta u\)
\begin{equation}\label{PseudoTime01} \left(\frac{1}{\Delta t} - \left.\frac{\partial R}{\partial u}\right|^n\right)\Delta u = R^{n} \end{equation}
Clearly as \(\Delta t \rightarrow \infty\), Eqn. PseudoTime01 tends to the standard Newton method. \(\Delta t\) acts as an under-relaxation parameter and greatly expands the radius of convergence compared to the standard Newton's method, providing greater flexibility in selecting the starting solution.
The predictor-corrector method can be derived by retaining additional terms in the Taylor series expansion of \(R^{n+1}\). Specifically, an approximation for the second-derivative is retained. The Taylor series expansion is
\begin{eqnarray}\label{TaylorExpApproxPC} R^{n+1} &\approx & R^n + \left.\frac{\partial R}{\partial t}\right|^n\Delta t + \left.\frac{\partial^2 R}{\partial t^2}\right|^n\frac{\Delta t^2}{2!} + O(\Delta t^3) \nonumber \\ R^{n+1} &\approx & R^n + \left.\frac{\partial R}{\partial t}\right|^n\Delta t + \left(\left.\frac{\partial R}{\partial t}\right|^{n+1} - \left.\frac{\partial R}{\partial \tau}\right|^n\right)\frac{\Delta t^2}{2!\Delta t} \nonumber \\ R^{n+1} &\approx & R^n + \left.\frac{\partial R}{\partial t}\right|^n\Delta t + \left(\frac{R^{n+1} - R^n}{\Delta t} - \left.\frac{\partial R}{\partial t}\right|^n\right)\frac{\Delta t^2}{2!\Delta t} \nonumber \\ R^{n+1} &\approx & \frac{1}{2}\left(R^{n+1} + R^n\right) + \left.\frac{\partial R}{\partial \tau}\right|^n\frac{\Delta t}{2} \nonumber \\ R^{n+1} &\approx & \frac{1}{2}\left(R^{n+1} + R^n\right) + \left.\frac{\partial R}{\partial u}\right|^n\frac{\Delta u}{2} \end{eqnarray}
Substituting Eqn. TaylorExpApproxPC in Eqn. EulerBackwards the following expression is obtained for \(\Delta u\)
\begin{equation}\label{PC01} \left(\frac{1}{\Delta t} - \frac{1}{2}\left.\frac{\partial R}{\partial u}\right|^n\right)\Delta u = \frac{1}{2}\left(R^{n+1} + R^n\right) \end{equation}
Unfortunately the value \(R^{n+1}\) is unknown and must be approximated. This is achieved using a predictor step identical to the standard method. That is \(R^{n+1}\) is approximated by solving
\begin{equation} \left(\frac{1}{\Delta t} - \left.\frac{\partial R}{\partial u}\right|^n\right)\Delta \hat{u} = R^{n} \end{equation}
where \(\Delta \hat{u} = \tilde{u}^{n+1} - u^n\) and \(\tilde{u}^{n+1}\) is a first approximation for \(u^{n+1}\). Using \(\tilde{u}^{n+1}\) an estimate of \(R^{n+1}\), e.g \(\tilde{R}^{n+1}\), is obtained. Eqn. PC01 can then be re-written as
\begin{equation}\label{PC02} \left(\frac{1}{\Delta t} - \frac{\alpha}{2}\left.\frac{\partial R}{\partial u}\right|^n\right)\Delta u = \frac{1}{2}\left(\tilde{R}^{n+1} + R^n\right) \end{equation}
where, as in TRACE, a relaxation parameter \(\alpha\) has been incorporated. Setting \(\alpha = 2\) we find the left-hand-sides of the predictor and corrector stages are identical. This is advantageous as the costly inversion of the left-hand-side matrix must be performed only once. The effect of including the relaxation parameter is twofold. To see this we rewrite Eqn. PC02 as follows
\begin{equation}\label{PC03} \left(\frac{1}{\alpha\Delta t} - \frac{1}{2}\left.\frac{\partial R}{\partial u}\right|^n\right)\Delta u = \frac{1}{2\alpha}\left(\tilde{R}^{n+1} + R^n\right) \end{equation}
We see that for \(\alpha > 1\) the right-hand-side is reduced as in the case of the standard or modified Newton methods. Therefore, in this sense the parameter acts to under-relax the system. However, we also see that the pseudo-time term \(\Delta t\) is increased by the same factor. Conversely, in this sense, the parameter acts to over-relax the system. Whether we obtain a net under- or over-relaxation is most probably dependent on the CFL number. For large \(\Delta t\) the PC scheme reduces to the modified Newton's method and \(\alpha\) acts as a relaxation parameter. For small \(\Delta t\), \(\alpha\) effectively increases the pseudo-time step-size thereby reducing the under-relaxation, or over-relaxing the system.
To solve Eqn. Advec1D time-accurately using Newton's method we define the modified residual \(R^*\) as follows
\begin{eqnarray}\label{Advec1DUnsteadyNewton} R^*(u) &=& \frac{\partial u}{\partial t} - \frac{\partial u}{\partial x} \nonumber \\ &=& \frac{\partial u}{\partial t} - R(u) \end{eqnarray}
To obtain a second-order accurate solution to Eqn. Advec1DUnsteadyNewton we can discretize the temporal operator using the Crank-Nicolson method. Eqn. Advec1DUnsteadyNewton then becomes:
\begin{eqnarray}\label{Advec1DUnsteadyNewtonCN} R^*(u) &=& \frac{\partial u}{\partial t} - \frac{\partial u}{\partial x} \nonumber \\ &=& \frac{u^{m+1} - u^n}{\Delta t} - \frac{1}{2}\left[R^{m+1} + R^n\right] \nonumber \\ &=& \frac{u^{m+1} - u^m + u^m - u^n}{\Delta t} - \frac{1}{2}\left[R^{m+1} + R^n\right] \nonumber \\ &=& \frac{\Delta u + u^m - u^n}{\Delta t} - \frac{1}{2}\left[R^{m+1} + R^n\right] \end{eqnarray}
A second order accurate solution to Eqn. Advec1D corresponds to \(u(x)\) such that \(R^*(u) = 0\). As before Newton's method is derived by expanding the residual \(R^*\) in a Taylor series and neglecting terms of order \(\Delta u^2\) and higher
\begin{equation}\label{NewtonMethodTimeAccurate} R(u^m + \Delta u) = R^m + \Delta u \left.\frac{dR}{du}\right|_{m} + O(\Delta u^2) \end{equation}
Setting \(R(u^m + \Delta u) = 0\), Eqn. NewtonMethodTimeAccurate may be rearranged for \(\Delta u\) to obtain
\begin{equation}\label{NewtonMethod02} \Delta u = -\frac{R^n}{\left.\frac{dR}{du}\right|_{n}} \end{equation}
where \(\Delta u = u^{n+1} - u^n\). If \(u^0\) corresponds to an initial guess, Eqn. NewtonMethod02 can be used to obtain an improved solution \(u^{n+1}\) that lies more closely to the true solution \(u\) corresponding to R(u) = 0. For additional robustness a relaxation parameter \(\alpha \geq 1\) may be included in Eqn. NewtonMethod02 as follows
\begin{equation}\label{NewtonMethod04} \Delta u = -\frac{1}{\alpha}\frac{R^n}{\left.\frac{dR}{du}\right|_{n}} \end{equation}
To treat Eqn. Advec1D time-accurately we introduce a pseudo-time operator and define the modified residual \(R^*\) as follows
\begin{eqnarray}\label{Advec1DPseudoTime} \frac{\partial u}{\partial \tau} &=& \frac{\partial u}{\partial t} - \frac{\partial u}{\partial x} \nonumber \\ &=& \frac{\partial u}{\partial t} - R(u) \nonumber \\ &=& R^*(u) \end{eqnarray}
The time-accurate problem can now be treated as a stationary problem in pseudo-time. In pseudo-time accuracy is unimportant and a first-order accurate time-integration scheme is sufficient.The method used in TRACE is based on the Euler-backwards method:
\begin{equation}\label{EulerBackwardsPseudoTime} u^{m+1} = u^m + \Delta\tau {R^*}^{m+1} \end{equation}
The precise form of \(R^*\) is dependent on the time-integration scheme used to discretize the physical-time operator. The standard approach employed in TRACE is to the use the second-order accurate Crank-Nicolson method. Using this method, the modified residual may be written as follows
\begin{eqnarray} {R^*}^{m+1} &=& \frac{u^{m+1} - u^{n}}{\Delta t} - \frac{1}{2}\left(R^{m+1} + R^n\right) \nonumber \\ &=& \frac{u^{m+1} - u^m + u^m - u^{n}}{\Delta t} - \frac{1}{2}\left(R^{m+1} + R^n\right) \nonumber \\ &=& \frac{\Delta u}{\Delta t} + \frac{u^m - u^{n}}{\Delta t} - \frac{1}{2}\left(R^{m+1} + R^n\right) \end{eqnarray}
where \(\Delta u = u^{m+1} - u^{m}\). As \({R}^{m+1}\) is unknown it must be approximated. This, as before, is accomplished using a Taylor series expansion.
\begin{equation}\label{TaylorExpPseudoTime} R^{m+1} \approx R^m + \frac{\partial R}{\partial \tau}\Delta \tau + \frac{\partial^2 R}{\partial \tau^2}\frac{\Delta \tau^2}{2!} + \frac{\partial^3 R}{\partial \tau^3}\frac{\Delta \tau^3}{3!} + O(\Delta\tau^4) \end{equation}
Dropping those terms of order \(\Delta\tau^2\) and rewriting the first-order derivative as follows
\begin{eqnarray} \frac{\partial R}{\partial \tau}\Delta \tau &=& \frac{\partial R}{\partial u}\frac{\partial u}{\partial \tau}\Delta \tau \\ &=& \frac{\partial R}{\partial u}\frac{\Delta u}{\Delta \tau}\Delta \tau \\ &=& \frac{\partial R}{\partial u}\Delta u \end{eqnarray}
The expression for the residual at time-level \(m+1\) becomes
\begin{equation}\label{TaylorExpApproxPseudoTime} R^{m+1} \approx R^m + \frac{\partial R}{\partial u}\Delta u + O(\Delta\tau^2) \end{equation}
Substituting Eqn. TaylorExpApprox in Eqn. EulerBackwards the following expression is obtained for \(\Delta u\)
\begin{equation}\label{PseudoTime01PseudoTime} \left[\left(\frac{1}{\Delta \tau} - \frac{1}{\Delta t}\right) + \frac{1}{2}\left.\frac{\partial R}{\partial u}\right|^m\right]\Delta u = \frac{u^m - u^{n}}{\Delta t} - \frac{1}{2}\left(R^{m} + R^n\right) \end{equation}
This is the basic form of the dual-time algorithm in TRACE, as solved by the IMPLICIT-SGS method.
The predictor-corrector method can be derived by retaining additional terms in the Taylor series expansion of \(R^{m+1}\). Specifically, an approximation for the second-derivative is retained. The Taylor series expansion is
\begin{eqnarray}\label{TaylorExpApproxPCPseudoTime} R^{m+1} &\approx & R^m + \left.\frac{\partial R}{\partial \tau}\right|^m\Delta \tau + \left.\frac{\partial^2 R}{\partial \tau^2}\right|^m\frac{\Delta \tau^2}{2!} + O(\Delta\tau^3) \nonumber \\ R^{m+1} &\approx & R^m + \left.\frac{\partial R}{\partial \tau}\right|^m\Delta \tau + \left(\left.\frac{\partial R}{\partial \tau}\right|^{m+1} - \left.\frac{\partial R}{\partial \tau}\right|^m\right)\frac{\Delta \tau^2}{2!\Delta \tau} \nonumber \\ R^{m+1} &\approx & R^m + \left.\frac{\partial R}{\partial \tau}\right|^m\Delta \tau + \left(\frac{R^{m+1} - R^m}{\Delta \tau} - \left.\frac{\partial R}{\partial \tau}\right|^m\right)\frac{\Delta \tau^2}{2!\Delta \tau} \nonumber \\ R^{m+1} &\approx & \frac{1}{2}\left(R^{m+1} + R^m\right) + \left.\frac{\partial R}{\partial \tau}\right|^m\frac{\Delta \tau}{2} \nonumber \\ R^{m+1} &\approx & \frac{1}{2}\left(R^{m+1} + R^m\right) + \left.\frac{\partial R}{\partial u}\right|^m\frac{\Delta u}{2} \end{eqnarray}
Substituting Eqn. TaylorExpApproxPCPseudoTime in Eqn. EulerBackwardsPseudoTime the following expression is obtained for \(\Delta u\)
\begin{equation}\label{PC01PseudoTime} \left[\left(\frac{1}{\Delta \tau} - \frac{1}{\Delta t}\right) + \frac{\alpha}{4}\left.\frac{\partial R}{\partial u}\right|^m\right]\Delta u = \frac{u^m - u^{n}}{\Delta t} - \frac{1}{2}\left[\frac{1}{2}\left(R^{m+1} + R^m\right) + R^n\right] \end{equation}
Where \(\alpha\) is a relaxation parameter. Unfortunately the value \(R^{m+1}\) is unknown and must be approximated. This is achieved using a predictor step identical to the standard method. That is, \(R^{m+1}\) is approximated by solving
\begin{equation} \left[\left(\frac{1}{\Delta \tau} - \frac{1}{\Delta t}\right) + \frac{1}{2}\left.\frac{\partial R}{\partial u}\right|^m\right]\Delta \hat{u} = \frac{u^m - u^{n}}{\Delta t} - \frac{1}{2}\left(R^{m} + R^n\right) \end{equation}
where \(\Delta \hat{u} = \tilde{u}^{m+1} - u^m\) and \(\tilde{u}^{m+1}\) is a first approximation for \(u^{m+1}\). Using \(\tilde{u}^{m+1}\) an estimate of \(R^{m+1}\), e.g \(\tilde{R}^{m+1}\), is obtained. Eqn. PC01PseudoTime can then be re-written as
\begin{equation}\label{PC02PseudoTime} \left[\left(\frac{1}{\Delta \tau} - \frac{1}{\Delta t}\right) + \frac{\alpha}{4}\left.\frac{\partial R}{\partial u}\right|^m\right]\Delta u = \frac{u^m - u^{n}}{\Delta t} - \frac{1}{2}\left[\frac{1}{2}(\hat{R}^{m+1} + R^m) + R^n\right] \end{equation}
Setting the relaxation parameter \(\alpha = 2\) it may be seen that the left-hand-sides of the predictor and corrector stages are identical. This is the basic form of the dual-time algorithm in TRACE, as solved by the IMPLICIT-PC method.
1.9.1