Computing Jacobian in WENO scheme for advection in a porous media

Question

I am trying to implement an advection equation for a coupled system of a two-phase flow in a porous media using a WENO scheme [1].

My equation is of the form:

\begin{align} \frac{\partial (\phi(x,t) C(x,t))}{\partial t} = - \nabla . (\vec{q(\phi, t, x)} C(x,t)) \end{align}

with $\phi$ beeing the porosity, which in this case is the proportion of fluid in the solide, C the concentration of an element, and $\vec{q}$ the fluid flux.

I am computing at each timestep $\phi$ and $\vec{q}$ in a separate solver.

As $\vec{q} = \vec{v} \times \phi$, I've decided to rewrite the advection equation in that form:

\begin{align} \frac{\partial (\tilde{C}(x,t))}{\partial t} = - \nabla . (\vec{v} \tilde{C}) = - \nabla . (f) \end{align}

with $\tilde{C} = \phi C$ (see Chemical advection of a fluid in a porous media) and $f$ the flux, which is a classical advection equation. I can then retrieve $C$ at the next timestep by dividing by the new porosity at the next timestep.

My problem right now is that I need the Jacobian of the flux to use a WENO scheme to compute the Lax-Friedrichs fluxes and I can't calculate it analytically because the link between $\vec{v}$ and $\tilde{C}$ is complexe as $\vec{v}$ depends on the proportion of fluid ($\phi$).

I guess my only solution would be to compute numerically the Jacobian. Is that right? My main problem is that I don't have a function that links $f$ and $\tilde{C}$, so it doesn't seem possible to use automatic differentiation. Should I use finite differences in this case? I am working in Julia, is there good tools for my case?

Answer 1

f=v[x,y] c[x,y] phi[x,y];
D[f,{{x,y}}]

Gives you the Jacobian (Mathematica). You just evaluate the variables at the relevant time step/instance.

Here she is:

$$\left\{c^{(1,0)}(x,y) v(x,y) \phi (x,y)+c(x,y) v^{(1,0)}(x,y) \phi (x,y)+c(x,y) v(x,y) \phi ^{(1,0)}(x,y), \\ c^{(0,1)}(x,y) v(x,y) \phi (x,y)+c(x,y) v^{(0,1)}(x,y) \phi (x,y)+c(x,y) v(x,y) \phi ^{(0,1)}(x,y)\right\}$$

Edit: For clarification of your question: derivatives wrt the state variable (in the case of a flux reconstruction method), you need $\frac{\partial f}{\partial \tilde{C}} = \frac{ \partial \vec{v} \tilde{C}}{\partial \tilde{C}} = \tilde{C} \frac{\partial \vec{v}}{ \partial \tilde{C}} + \vec{v}\frac{\partial{\tilde{C}}}{\partial \tilde{C}} = \tilde{C} \frac{\partial \vec{v}}{ \partial \tilde{C}} + \vec{v}$. Of course you are right, if you can't access the $\partial_{\tilde{C}} \vec{v}$ term then you need to approximate it. Finite differences should be fine.