Diffusion Reaction pde with nonlinear source term / Potential Term

Question

my first attempt to write a couple diffusion reaction equations with FiPy.

My equations are for three different concentrations c_i, just with the diffusion part and source (LateX input):

\partial_t c_1 = \nabla \cdot (D_1 \nabla c_1) + A_1 c_{1} + B_1 c_2^4 + W_1 \\
\partial_t c_2 = \nabla \cdot (D_2 \nabla c_2) + c_{1}(A_2 + 2c_1) + c_2^4 + W_2\\
\partial_t c_3 = \nabla \cdot (D_3 \nabla c_3) + A_3 c_{1} - W_3 \\

The coefficients D_i, A_i, B_i, P_i and W_1 are constants. I have already wrote the code with just de Diffusion term and the nolinear source term. But for a Range of 100 i get some odd behavior for c2. Is maybe the nonlinear source term, that i wrote wrong? I use the ImplicitSourceTerm command and i thought that this would be linearized the term. Am i missing something? do i have to do the linearization by myslef? (like with Taylor?)

How to Couple Advection Diffusion Reaction PDEs with FiPy and https://www.ctcms.nist.gov/fipy/examples/diffusion/generated/examples.diffusion.coupled.html

from fipy import *
from fipy import CellVariable, Variable, Grid1D, TransientTerm, DiffusionTerm, ImplicitSourceTerm, LinearLUSolver, Viewer
import math
from fipy.tools import numerix

# Konstanten
F       = 96485.3399  #C /mol
R       = 8.314472    # J/(mol*K) 
T       = 273.18      # K
alpha   = 0.5
c_std   = 1.           #  mol/s
c1_sat  = 1.           #  mol/s
eta_Zn  = 1.           # V
kIc = 6.3      # mol/s Von Reaktion I an der Kathode
kIa = 5.       # mol/s Von Reaktion I an der Anode
kII = 0.25        # mol/s von Reaktion II
mu_1I   = 1.    # Von Zinkat in Reaktion I
mu_1II  = -2.   # Von Zinkat in Reaktion II
mu_2I   = -4.   # Von OH in Realtion I
mu_2II  = 2.    # Von OH in Reaktion II
mu_3II  = 1.   # Von H2O in Reaktion II
ep1 = 1.    # Von Zinkat
ep2 = 1.   # Von OH
ep3 = 1.   # Von H20
D1 = .5     # Von Zinkat
D2 = 1.     # Von OH
D3 = 0.1    # Von H20

## Mesh
L = 10.
nx = 1000
mesh = Grid1D(dx=L/1000, nx=nx)
x = mesh.cellCenters[0]

## Initial Conditions
c1 = CellVariable(name="c1", mesh=mesh, value=1., hasOld=True)
c2 = CellVariable(name="c2", mesh=mesh, value=1., hasOld=True)
c3 = CellVariable(name="c3", mesh=mesh, value=1., hasOld=True)

## Boundary Conditions
c1.constrain(2., mesh.facesLeft)
c1.constrain(0., mesh.facesRight)
c2.constrain(0., mesh.facesLeft)
c2.constrain(2., mesh.facesRight)
c3.constrain(0., mesh.facesLeft)
c3.constrain(2., mesh.facesRight)
c2.faceGrad.constrain(1., where=mesh.facesRight)

# Definition Konstanten fuer SourceTerm
Zn1 = (kIc/c_std)*math.exp(-(1-alpha)*(F/(R*T))*eta_Zn)+kII
Zn2 = (kIa/(c_std)**4)*math.exp(alpha*(F/(R*T))*eta_Zn)
Zn3 = kII*c1_sat
OH1 = 4*(kIc/c_std)*math.exp(-(1-alpha)*(F/(R*T))*eta_Zn)+2*kII
OH2 = 4*(kIa/(c_std)**4)*math.exp(alpha*(F/(R*T))*eta_Zn)
OH3 = 2*kII*c1_sat
H1 = kII
H2 = 0.
H3 = kII*c1_sat

sourceCoeff1 = (Zn1*c1) + (Zn1*c2**4) + Zn3
sourceCoeff2 = c1*(OH1 + 2*c1) + c2**4 + OH3
sourceCoeff3 = c1*H1 - H3

eq1 = (TransientTerm(var=c1) == DiffusionTerm(coeff=D1, var=c1) + ImplicitSourceTerm(coeff=sourceCoeff1, var=c1))
eq2 = (TransientTerm(var=c2) == DiffusionTerm(coeff=D2, var=c2) + ImplicitSourceTerm(coeff=sourceCoeff2, var=c2))
eq3 = (TransientTerm(var=c3) == DiffusionTerm(coeff=D3, var=c3) + ImplicitSourceTerm(coeff=sourceCoeff3, var=c3))

eqn = eq1 & eq2 & eq3


vi = Viewer((c1, c2, c3))

for t in range(50):
    c1.updateOld()
    c2.updateOld()
    c3.updateOld()
    eqn.solve(dt=1.e-3)
    vi.plot()

I am trying also to add the potential term but i do not understand if that is even possible. Something like that (same as above but with the potential part)

\partial_t c_1 = \nabla \cdot (D_1 \nabla c_1) + \nabla \cdot \biggl(P_1c_1\nabla\Phi\biggr) + A_1 c_{1} + B_1 c_2^4 + W_1 \\
\partial_t c_2 =  \nabla \cdot (D_2 \nabla c_2) + \nabla \cdot \biggl(P_2c_2\nabla\Phi\biggr) + c_{1}(A_2 + 2c_1) + c_2^4 + W_2\\
\partial_t c_3 = \nabla \cdot (D_3 \nabla c_3) + \nabla \cdot \biggl(P_3c_3\nabla\Phi\biggr) + A_3 c_{1} - W_3 \\

I analysed the example from https://www.ctcms.nist.gov/fipy/examples/phase/generated/examples.phase.binaryCoupled.html to try to add the potential term like a Diffusionterm but i always get the error:

SolutionVariableNumberError: Different number of solution variables and equations.

what i expected but maybe i am missing something.

Hier that part of the code:

eq1 = (TransientTerm(var=c1) == DiffusionTerm(coeff=D1, var=c1) + ImplicitSourceTerm(coeff=sourceCoeff1, var=c1)) + DiffusionTerm(coeff=D1phi, var=phi
eq2 = (TransientTerm(var=c2) == DiffusionTerm(coeff=D2, var=c2) + ImplicitSourceTerm(coeff=sourceCoeff2, var=c2)) + DiffusionTerm(coeff=D2phi, var=phi)
eq3 = (TransientTerm(var=c3) == DiffusionTerm(coeff=D3, var=c3) + DiffusionTerm(coeff=D3phi, var=phi) + ImplicitSourceTerm(coeff=sourceCoeff3, var=c3))

Also try to use the command sweep to do some numerical analysis like residual error and the norm. I saw a great explanation about the sweep command here: Solver tolerance and residual error when using sweep function in FiPy

That part of the code (just for diffusion term with source)

dt = 1.e5

solver = LinearLUSolver(tolerance=1e-10)

c1.updateOld()
c2.updateOld()
c3.updateOld()
res = 1.
initialRes = None

while res > 1e-4:
     res = eq.sweep(dt=dt, solver=solver)
     if initialRes is None:
         initialRes = res
     res = res / initialRes

Neverthless i do not get any result or error and have to stop the process manually.

In summary my questions are: it is possible to solve with FiPy couple pdes with diffusion term and a potential term? It is possible to implement the command sweep also for coupled pdes?

Or i am missing something?

I am very grateful for any help or advice. I hope, i wrote this clearly. Thank you a lot.

Answer 1

Here is a working version of your code. I'm not sure if it is correct, but it seems more sensible.

Use an explicit source

The major change made is to remove the implicit source terms and replace them with explicit sources. For example the first equation

eq1 = (TransientTerm(var=c1) == DiffusionTerm(coeff=D1, var=c1) + ImplicitSourceTerm(coeff=sourceCoeff1, var=c1))

is now

eq1 = (TransientTerm() == DiffusionTerm(coeff=D1) + sourceCoeff1)

The use of the ImplicitSourceTerm in all three equations in the OP's code is incorrect. To understand its use take a source of the form value * phi where phi is the variable being solved. In this case the correct source would be ImplicitSourceTerm(value). The multiplication by phi is implicitly assumed. Additionally, this only works correctly if value is negative [edit]. See this example for more explanation.

[edit]: actually, ImplicitSourceTerm handles both positive and negative values so that if the value is positive then the source is added explicitly.

Uncouple the equations

It's best to uncouple the equations at first to better understand what's going on. Coupling equation works best when there are multiple variables inside spatial derivatives that are coupled inside a single equation.

In this case, it might be useful to couple the equations due the the source terms, but only for the purposes of efficiency and only after the problem solution is validated. The new code now solves the equations in an uncoupled form.

for t in range(100):
    c1.updateOld()
    c2.updateOld()
    c3.updateOld()
    for sweep in range(5):
        res1 = eq1.sweep(c1, dt=dt, solver=solver)
        res2 = eq2.sweep(c2, dt=dt, solver=solver)
        res3 = eq3.sweep(c3, dt=dt, solver=solver)
        print('res1', res1)
        print('res2', res2)
        print('res3', res3)

    vi.plot()

Use a sweep loop

Since the equations are non-linear in the sources it is a good idea to use an inner loop and print the residuals to check convergence at each time step as in the code above.

Use a different solver

To make things work it seemed to require a different solver from the default chosen by FiPy. This can often depend on the solver suite employed. In the case, I was using Python 3 so I was probably using the Scipy solvers. To make it work I used the LinearLUSolver otherwise it ran incredibly slowly.

Future

I the long run it might be a good idea to linearize the sources and use ImplicitSourceTerms as well as coupling the equations for the purposes of efficiency, but only after the solution is fully understood and validated. See this example and this explanation for help.