Tighter coupling of source terms with Corner Transport Upwind Scheme
Consider wave equation with a 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}
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'