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?


1 Answer 1

f=v[x,y] c[x,y] phi[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.

  • $\begingroup$ I see, so I should find a way to do it easily in Julia. So you confirm me that my idea should work? $\endgroup$
    – Iddingsite
    Mar 31, 2022 at 14:00
  • $\begingroup$ 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. $\endgroup$
    – user20857
    Mar 31, 2022 at 18:11
  • $\begingroup$ 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}$ ? $\endgroup$
    – Iddingsite
    Apr 2, 2022 at 21:26
  • $\begingroup$ 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). $\endgroup$
    – user20857
    Apr 2, 2022 at 21:55
  • 1
    $\begingroup$ @Iddingsite my mistake, you are right. I updated the answer via an edit. $\endgroup$
    – user20857
    Apr 3, 2022 at 13:37

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.