Tighter coupling of source terms with Corner Transport Upwind Scheme
1D wave equation with source term ==
\frac{\partial q}{\partial t} = -\frac{\partial f(q)}{\partial x} + s(q)
Exact integral finite volume solution to this problem is
Q_i^{n+1} = Q_i^n + \frac{\Delta t}{\Delta x} \left ( F_{i-1/2} - F_{i+1/2} \right ) + S_i \Delta t
Where
- Q_i^n is the volume average of q over cell i at time t^n
- F_{i-1/2} is the time and area averaged flux f(q) at the boundary between cell i and cell i-1
- S_i is the time and volume averaged source term s(q) over cell i
In Riemann methods, in 1D, we can write
F_{i+1/2} = R(Q_{L_{i+1/2}},Q_{R_{i+1/2}})
where Q_{L_{i+1/2}} is the time averaged area average of q on the left side of the interface and R is the Riemann solver for f(q)
Now the trick is how to calculate Q_{L_{i+1/2}} = \displaystyle \frac{\int_{t^n}^{t^{n+1}} q_L(t) dt}{\Delta t} = \displaystyle \frac{\int_{0}^{\Delta t} q_L(t) dt}{\Delta t}
Method of Characteristics
If we operator split the source term from the hyperbolic term, and if the hyperbolic term is trivial
f(q) = aq
where a is a constant and via the method of characteristics
q(x,t)=q(x(t),0) \rightarrow \frac{\partial q}{\partial t}=\frac{\partial q}{\partial x} \frac{\partial x}{\partial t} =-a \frac{\partial q}{\partial x} \rightarrow \frac{\partial x}{\partial t}=-a
which gives x(t)=x-at
and q(x,t) = q(x-at,0)
and the integral for Q_{L_{i+1/2}} becomes
Q_{L_{i+1/2}} = \displaystyle \frac{\int_{0}^{\Delta t} q_L(t) dt}{\Delta t} = \displaystyle \frac{\int_{0}^{\Delta t} q(x_L,t) dt}{\Delta t}= \displaystyle \frac{\int_{0}^{\Delta t} q(x_L-at,0) dt}{\Delta t}
Q_{L_{i+1/2}} = \displaystyle \frac{\int_{x_L-a\Delta t}^{x_L} q(x',0) dx'}{a \Delta t}
Source term contribution
For the source term contribution to the time averaged interface state, things are a little tricky. Imaging a constant source term s(q)=b
Then the contribution to q(x,t) = q(x,0)+bt and the time averaged value for q would be
\frac{\displaystyle \int_0^{\Delta t} q(x_L,t') dt'}{\Delta t} = \frac{\displaystyle \int_0^{\Delta t} q(x_L,0)+bt' dt'}{\Delta t} = q(x_L,0) + \frac{1}{2} b \Delta t
Now if we include the solution from characteristic tracing in our source term calculation, the first term just gives us the same spatial integration, however the second term is a double integral.
\frac{\displaystyle \int_0^{\Delta t} q(x_L,t') dt'}{\Delta t} = \frac{\displaystyle \int_0^{\Delta t} q(x_L-at',0) dt'+\int_0^{\Delta t} \int_0^{t'} s(q(x_L-at'',0)) dt'' dt'}{\Delta t}
There are lots of ways to reduce this double integral into a single (or finite number) of source evaluations using any combination of midpoint rules or trapezoid rules. Ideally we would like to use our already calculated spatial integral which we can do by using the trapezoid rule on the outer integral and exchanging the source evaluation with the integral average. This gives
\frac{\int_0^{\Delta t} \int_0^{t'} s(q(x_L-at'',0)) dt'' dt'}{\Delta t} \approx \frac{1}{2} \int_0^{\Delta t} s(q(x_L-at'',0)) dt'' \approx \frac{1}{2} s \left ( \int_0^{\Delta t} q(x_L-at'',0) dt'' \right ) = \frac{1}{2} s \left( Q_{L_{i+1/2}}^* \right ) \Delta t
where Q_{L_{i+1/2}}^* = \displaystyle \frac{\int_{x_L-a\Delta t}^{x_L} q(x',0) dx'}{a \Delta t}
Summary of 1D
- Spatial reconstruction and averaging of solution over domain \lambda \Delta t to get interface states
- Application of source term to interface states over \Delta t/2
- Calculate Fluxes from interface states
- Update Cells using those fluxes
- Calculate source term using time averaged cell center
- Update final states using source term for \Delta t
CTU
PPM details
In more general terms, we can write \frac{\partial f}{\partial x} = \frac{\partial f}{\partial q} \frac {\partial q}{\partial x}
and we can estimate \left . \frac{\partial f}{\partial q} \right | _0 = A
Now A = R \Lambda L and we can also add and subtract \lambda_M I so that \frac{\partial f}{\partial x} = \lambda_M \frac{\partial q}{\partial x} + R (\Lambda - \lambda_M I) L \frac{\partial q}{\partial x}
and try to use characteristic tracing with q(x,t)=q(x-\lambda_M t,0) + q'(x,t)
We can use the midpoint rule to estimate
Q_{L_{i+1/2}} = \displaystyle \frac{\int_{0}^{\Delta t} q_L(t) dt}{\Delta t} \approx q_L \left ( \frac{\Delta t}{2} \right ) + O(\Delta t^3)
Usually some form of spatial reconstruction is used to calculate q_L(0) and then the evolution equation for q is used to approximate time derivatives using spatial derivatives.
q_L(\Delta t/2) = q(x_L,0) + \displaystyle \int_0^{t/2} \frac{\partial q(x_L,t')}{\partial t}
which we then use the evolution equation to replace time derivatives with spatial ones
q_L(\Delta t/2) = q_L(x_L,0) + \displaystyle \int_0^{t/2} \left [ -\frac{\partial f(x_L,t'))}{\partial x} + s(q(x_L,t')) \right ]
and then we can linearize the flux function about q
\frac{\partial f(q(x,t))}{\partial x} = \frac{\partial f}{\partial q} \frac{\partial q(x,t)}{\partial x} = \lambda \frac{\partial q(x,t)}{\partial x}
q_L(\Delta t/2) = q_L(0) + \displaystyle \int_0^{t/2} \lambda \frac{\partial (q(x_L,t'))}{\partial x} dt'+ \int_0^{t/2} s(q(x_L,t')) dt'
Now if we ignore the source term, we can directly solve for
q(x,t)=q(x-\lambda t, 0)
and under the change of variables x'=x_L-\lambda t' we can rewrite the integrals as
q_L(\Delta t/2) = q_L(0) + \displaystyle \int_{x_L}^{x_L-\lambda t/2} \frac{\partial (q(x',0))}{\partial x} dx' + \int_{x_L}^{x_L-\lambda t/2} s(q(x',0)) dx'
which just gives us
q_L(\Delta t/2) = q_L(0) - (q(x_L,0) - q(x_L-\lambda \Delta t/2)) + s(\bar{q})\Delta t/2
where \bar{q} = \displaystyle \int_{x_L-\lambda \Delta t/2}^{x_L}q(x',0) dx'
\frac{\partial f(q(x,t))}{\partial x} = \frac{\partial f}{\partial q} \frac{\partial q(x,t)}{\partial x} = A \frac{\partial q(x,t)}{\partial x} = R \Lambda L \frac{\partial q(x,t)}{\partial x} = \displaystyle \sum_i R_i \lambda_i L_i \frac{\partial q(x,t)}{\partial x} \approx \displaystyle \sum_i R_i \lambda_i L_i \frac{\partial q(x-\lambda_i t,0)}{\partial x}
and use the characteristics (eigenvectors of the Jacobian of the flux) to linearize the flux
q_L(\Delta t/2) = q_L(0) + \displaystyle \int_0^{t/2} \left [ R \Lambda L \frac{\partial (q(x_L,t'))}{\partial x} + s(q(x_L,t')) \right ]
q_L(\Delta t/2) = q_L(0) + \displaystyle \int_0^{t/2} \left [ \sum R_i \lambda_i L_i \frac{\partial (q(x_L,t'))}{\partial x} + s(q(x_L,t')) \right ]
q_L(\Delta t/2) = q_L(0) + \displaystyle \int_0^{t/2} \left [ \sum R_i \lambda_i \delta l (L_i \frac{\partial (q(x_L,t'))}{\partial x} + s(q(x_L,t')) \right ]
q_L(\Delta t/2) = q_L(0) + \displaystyle \int_0^{t/2} \left [ \sum R_i \lambda_i \frac{\partial (q(x_L,t'))}{\partial x} + s(q(x_L,t')) \right ]
and then we can use the solution for the characteristics q(x_L,t')=q(x_L-\lambda t',0) to turn the time integral into a spatial integral
q_L(\Delta t/2) = q_L(0) + A \displaystyle \int_0^{t/2} \left [ \frac{\partial (q(x_L-\lambda t',0))}{\partial x} + s(q(x_L-\lambda t',0)) \right ]
which we can do a variable substitution x' = \lambda t'