# Computing Jacobian in WENO scheme for advection in a porous media

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?

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.

• I see, so I should find a way to do it easily in Julia. So you confirm me that my idea should work? Mar 31, 2022 at 14:00
• The Jacobian is trivial in this case (correct me if I'm wrong). It doesn't matter what language you use to compute it, or just do it by hand. Formulate once and reuse.
– user20857
Mar 31, 2022 at 18:11
• I have to admit that my knowledge of Jacobians is quite shallow and there is something that is confusing for me. I understand how you came up with the 3 sums for each dimension, and yes, the Jacobian seems quite trivial. But I don't really understand how to do it by hand with finite differences. If I take the term with the derivative of the velocity for example, in the x direction: $\frac{\partial v}{\partial \tilde{C}} C \phi$. Should I just implement $\frac{v_{i+\frac{1}{2},j} - v_{i-\frac{1}{2},j}}{\tilde{C}_{i+\frac{1}{2},j} - \tilde{C}_{i-\frac{1}{2},j}} C_{i,j} \phi_{i,j}$ ? Apr 2, 2022 at 21:26
• One of the velocity derivative terms looks something like $C(x) \phi(x) \partial_x v(x)$ in 1D. You just take the finite difference approximation of the $\partial_x v(x)$ derivative w.r.t. the $x$-direction in the usual way (there is no $C$ or $\tilde{C}$ involved).
– user20857
Apr 2, 2022 at 21:55
• @Iddingsite my mistake, you are right. I updated the answer via an edit.
– user20857
Apr 3, 2022 at 13:37