SUPG weak form for a system of reaction-advection-diffusion PDEs

1

Hi,

I am trying to solve a system of 12 reaction-advection-diffusion PDEs given by:

\frac{D X_i}{D t}+ X_i \nabla \cdot{\bf v} = D_i \Delta X_i + f_i

where \frac{D }{D t} = \frac{\partial }{\partial t}+{\bf v} \cdot \nabla is the material derivative and the system is formulated in Eulerian. We know \nabla \cdot {\bf v} = \Psi (f) and {\bf v} is the solution of a simultaneous mechanical problem decoupled from the given PDE. Also, the diffusion rates are very small O(10^{-8}). Now, of course with such small diffusion rates, the direct Galerkin method doesn’t necessarily satisfy the nonnegativity criteria and it does show in my simulations. I used the following SUPG stabilization weak form:

\int \phi_i \frac{D X_i}{D t} dX + \int \phi_i X_i \nabla \cdot{\bf v} dX + \int D_i \nabla\phi_i \cdot \nabla X_i dX \\+ \int (\frac{D X_i}{D t} + X_i \nabla \cdot{\bf v} - D_i \Delta X_i - f_i ) (\tau \phi_i \nabla \cdot \mathbf{v}) dX = 0
where
\tau = \left( \left(2 \frac{||\mathbf{v}||}{h} \right)^2+ \left(4 \frac{{D_i}}{h^2} \right)^2 \right)^{-\frac{1}{2}}
and
h = 2*Circumradius(mesh) .

But I don’t get any improvements and the problem convergence rate decreases significantly. I am sure there is something wrong with my weak form and that I am not applying the method right. I’d really appreciate it if you could let me know how to do it properly or if there is an easier way to go about solving the nonnegativity issue.

1 Like
2

If the field \mathbf{v} is a solution to Stokes flow, then \nabla\cdot\mathbf{v} should be zero, so the term X_i\nabla\cdot\mathbf{v} can be neglected (unless it is meant to be read as (X_i\nabla)\cdot\mathbf{v}, where the \nabla acts to the left, but this term would be part of the material time derivative). This is not necessarily exactly true when \mathbf{v} is a discrete approximate solution to Stokes flow, but it seems to be common practice to remove the term when formulating the continuous problem, even if it will be solved with some small nonzero velocity divergence in the discrete setting.

Then you are left with the standard advection–diffusion problem. Since there appears to be no coupling between components of \mathbf{X}, I will simply drop the subscript i in the rest of my answer, with the understanding that it applies to each scalar component of \mathbf{X}. The stabilization term should be

\int_\Omega\left(\frac{\partial X}{\partial t} + \mathbf{v}\cdot\nabla X - D\Delta X - f\right)\left(\tau\mathbf{v}\cdot\nabla\phi\right)\,d\Omega\text{ ,}

which you can see will still be nonzero even when \nabla\cdot\mathbf{v} = 0. Note that I expanded out the definition of DX/Dt, to more clearly show how you get a term of the form

+\int_\Omega\left((\tau\mathbf{v}\otimes\mathbf{v})\nabla X\right)\cdot\nabla\phi\,d\Omega

which provides the desired “streamline diffusion”.

Lastly, if there is coupling between components, optimal stabilization becomes significantly more complicated, as discussed in this paper, although simplified versions may still work in practice.

3

Thanks for the reply. As for \nabla \cdot {\bf v} it is not zero. I don’t know why I forgot to mention that. It is a rather peculiar formulation. In fact, it is related to the reaction term in the mentioned PDEs. Basically, it says the movement of the cells is dictated by their reactions which causes some kind of pressure change. So \nabla \cdot {\bf v} = \Psi(f). In this case, do I simply add the term X_i \nabla \cdot {\bf v} in the integral you provided? Or is it more complicated?

4

If there is a significant reaction term, then that can also be a source of instabilities (even in the absence of advection). Take a look at Section 4 of this paper. The method spelled out in (34)–(40) of the linked paper should be applicable to your problem, with some changes in notation:

u\to X~,~v\to\phi~,~\sigma\to\nabla\cdot\mathbf{v}~,~\mathbf{a}\to\mathbf{v}~,~\kappa\to D\text{ .}
1 Like
5

Wonderful paper. I read all of it and was able to successfully implement it on FEniCS for the example given in the paper. However, my problem is time-dependent so I wanted to check something with you. Based on this paper the equation

X_i \nabla \cdot{\bf v}- D_i \Delta X_i = f_i

is stabilized by:

B(X_i,\phi_i) =(X_i \nabla \cdot{\bf v}, \phi_i)+ (D_i \nabla X_i,\nabla \phi_i) \\- \displaystyle \sum_{Triangles} (X_i \nabla \cdot{\bf v}- D_i \Delta X_i, \tau_k (\phi_i \nabla \cdot{\bf v}- D_i \Delta \phi_i))

F(X_i) = (f_i,\phi_i) - \displaystyle \sum_{Triangles} (f, \tau_k (\phi_i \nabla \cdot{\bf v}- D_i \Delta \phi_i)

Now I did some dig in the literature (see this for example) and I found out that in case of a time-dependent problem like mine:

\frac{D X_i}{D t}+ X_i \nabla \cdot{\bf v} = D_i \Delta X_i + f_i

You simply change the weak forms by:

B(X_i,\phi_i) =(\frac{D X_i}{D t},\phi_i)+(X_i \nabla \cdot{\bf v}, \phi_i)+ (D_i \nabla X_i,\nabla \phi_i) \\- \displaystyle \sum_{Triangles} (\frac{D X_i}{D t}+X_i \nabla \cdot{\bf v}- D_i \Delta X_i, \tau_k (\frac{D \phi_i}{D t}+\phi_i \nabla \cdot{\bf v}- D_i \Delta \phi_i))

F(X_i) = (f_i,\phi_i) - \displaystyle \sum_{Triangles} (f, \tau_k (\frac{D \phi_i}{D t}+\phi_i \nabla \cdot{\bf v}- D_i \Delta \phi_i)

Now this might be a stupid question but how do we implement the time derivative of a test function in Fenics?

6

That would only make sense in a space–time finite element method, which is somewhat involved to implement. In most work I’ve seen with SUPG and related VMS methods, people just leave the time derivative out from the operator acting on the test function. You can also view the midpoint finite difference scheme in time as a space–time finite element method where the test function is constant w.r.t. time in each time step, so the time derivative term would simply be zero. (See Remarks 4 and 5 here.)

I also recall seeing a paper once (I think by Franca and others, although I can’t find it now) suggesting that it makes more sense to discretize first in time, so that the time derivative’s finite difference contributes an additional reaction term and a source term. Then methods for stabilizing steady reaction–advection–diffusion can be applied to the spatial discretization of each time step.

7

I see. Thank you so much for your help. Even though I still have stabilization problems, I know much more than I knew two days ago. I will take a look at space-time methods to see if I can incorporate them into my problem.

8

You may also consider adding a shock(discontinuity)-capturing operator into your stabilized (SUPG) formulation if the instability problems you mentioned that you had still had are local instabilities, i.e., they occur near sharp gradients. Besides, you may also use the time-dependent version of the stabilization parameter \tau, by adding the term \left(\frac{2}{\varDelta t} \right)^2 into the definition of \tau. You can refer to this paper, Eq. (3.9).

9

This PDE model is considered in [1], see the scheme in Eq. 12.3.1, which shows the correct way to add stabilization. There if you take \rho=0, you get SUPG scheme.

Your PDE is

X_t + L(X) = f, \qquad L(X) = \nabla\cdot(Xv) - D \Delta X

SUPG stabilization is of the form

\int_\Omega (X_t + L(X) - f) \tau L_{SS}(\phi) d\Omega

where L_{SS} is the skew-symmetric part of L

L_{SS}(X) = \frac{1}{2}\nabla\cdot(v X) + \frac{1}{2}v \cdot \nabla X

If \nabla\cdot v=0 then you get the simpler form

L_{SS}(X) = v \cdot \nabla X

[1] Quarteroni and Valli, Numerical Approximation of PDEs