# Handling time derivative as source term in SDIRK methods

I am currently facing some challenge implementing a less traditional PDE which take a form similar to the Navier-Stokes equation, except that the continuity equation is modified such that: $$\epsilon \nabla \vec{u} + \vec{u} \cdot \nabla \epsilon + \frac{d\epsilon}{dt} = m$$

where $$\epsilon(t,\vec{x})$$ is known at all time $$t$$ and $$m$$ is a mass source term which is used to build an MMS.

The associated momentum equation is : $$\epsilon \frac{d\vec{u}}{dt} + m \vec{v} + \epsilon \vec{u} \cdot \nabla \vec{u} = -\nabla p + \nu \nabla^2 \vec{u} + \vec{f}$$

where $$\nu$$ is known and $$\vec{f}$$ is a forcing term used to impose the MMS.

In steady-state, our code converge quite well to the appropriate order of convergence in the L2 norm. In transient, our code also converge quite well when using BDF1 and BDF2 time-integrators. Additionally, when using SDIRK and using a manufactured solution (MMS) where $$\epsilon = \epsilon(\vec{x})$$ we face no issue with convergence in time when using an SDIRK2 or SDIRK3 scheme. The issue arises when we consider cases with $$\epsilon = \epsilon(\vec{x},t)$$ and SDIRK2 or SDIRK3 time integrators.

How does one discretize the time derivative of the known variable $$\epsilon$$ using SDIRK2 or SDIRK3? This is a bit confusing to me, since $$\epsilon$$ is not an unknown. For the BDFs I use the appropriate backward derivative, but I am unsure how to handle this term well using SDIRK approaches...

For your information, this is all solved in the context of FEM, but I don't think this has an impact on the question itself.

How do you evaluate the term $$\epsilon$$ and $$\dfrac{d\epsilon}{dt}$$ in the unsteady case ? If I understand correctly, you use the Runge-Kutta scheme with a prescribed evolution for $$\epsilon$$ and solve each stage for the unknown $$\vec{u}$$, is that right ?

Then the $$i$$-th stage of a DIRK scheme should read: $$y_{n,i} = y_n + \Delta t \sum\limits_{j=1}^{i} a_{ij} f(t_{n,j}, y_{n,j})$$ with $$a_{ij}$$ the relevant element of the Butcher tableau of the method, $$y=(\epsilon, u_x, u_y, u_z)$$ in 3D, and $$f(t,y) = \left( m-\epsilon \nabla \vec{u} - \vec{u} \cdot \nabla \epsilon,\quad (1/\epsilon)(-\nabla p + \nu \nabla^2 \vec{u} + \vec{f} - m\vec{v} - \epsilon \vec{u}\cdot \nabla \vec{u})^t \right)^t$$with all the spatial gradients evaluated at time $$t$$ with the scheme of your choice.

In principle, you just have to plug in the values of $$\epsilon$$ you want to enforce, and I guess you then use a Newton method on the vector of unknowns $$X=(m, \vec{u}^t)^t$$. Is that what you are doing ?

EDIT: If you are in the incompressible regime, the system for the complete flow is differential-algebraic (DAE) of index 2 when the space dimension is higher than 1, for which DIRK methods that are not stiffly accurate are not well suited. Petzold discusses some of these aspects in the book "Computer Methods for Ordinary Differential Equations and Differential-Algebraic Equations", page 238 for the incompressible Navier-Stokes equation.

EDIT2: In the Runge-Kutta formulation (first equation of this answer), the term $$\dfrac{d \epsilon}{dt}$$ does not appear. Pragmatically, it is "reconstructed" naturally by the method via the sum $$\sum a_{ij} f(t_j, y_{n,j})$$. I've dealt with a similar problem in one spatial dimension and this works great. The formulation I've written in my answer for $$\epsilon$$ can be seen as a "backward substitution" (not in the numerical sense) of the known evolution of $$\epsilon$$ into the RK equation that would otherwise define it (if we were not forcing the value of $$\epsilon$$). Hence the iterative resolution of this equation will adapt the other variables (here the pressure and the velocity field) so that all the equations of your PDE are satisfied. So in principle just plug the values of $$\epsilon(t_{n,j}, X)$$ in the RK formulation, with $$t_{n,j} = t_n + c_j \Delta t$$ the time of the $$j$$-th stage.

• The term $\epsilon$ is known analytically. So i just have a value for it. The term $\frac{d\epsilon}{dt}$ is what I am unsure of. For the BDF methods I just numerically differentiate it by applying a backward derivative (finite element approximation) using the previous value of $\epsilon$ that are stored. For SDIRK I am unsure of the right approach. In our case the SDIRK methods are able to recover their right convergence order (at least in test cases like a taylor-green vortex). I just don't know how to treat the $\frac{d\epsilon}{dt}$ term while remaining consistent.
– BlaB
Nov 30 '20 at 15:11
• Ok. Are you in the incompressible regime ? I've added some info in my answer + an edit at the end. Nov 30 '20 at 15:21
• Not in my opinon. Let's forget you want to enforce $\epsilon$. Then the ODE is $d_t\epsilon = f(\epsilon,\vec{u},m) = m - \epsilon \nabla \vec{u} - \vec{u} \cdot \nabla \epsilon$ (I drop the equation on $\vec{u}$ for simplicity). The $i$-th DIRK stage is obtained by solving $0=\epsilon_{n,i} - \epsilon_n - \Delta t\sum\limits_{j=1}^{j=i} a_{ij} f(t_{n,j}, y_{n,j})$. IN the normal case (no forcing of $\epsilon$), you iterate on $\epsilon_{n,j}$ to solve this equation. In your case, you already know what $\epsilon$ is, so you iterate on $m$ instead... (still dropping $\vec{u}$ for simplicity). Nov 30 '20 at 16:30
• So the term to add to my non-linear equation would be the analytically imposed value of $\epsilon$ at the intermediate steps weighted by their SDIRK coefficients? I have tried this, but when I do this I lose the underlying order of the scheme, which is not the case when I use the same approach for BDF2 for example...
– BlaB
Nov 30 '20 at 16:32
• I guess the error has to lie somewhere in my implementation. In all cases, thank you very much for the help.
– BlaB
Nov 30 '20 at 17:23