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?