Solving coupled PDEs using FiPy in Python

[frandgr]

Hi, I'm trying to solve a system of coupled partial differential equations using FiPy, but I keep running into an ExplicitVariableError. The error message says "Terms with explicit Variables cannot mix with Terms with implicit Variables".

Here are the equations I'm trying to solve:

$$\frac{\partial (C_j)}{\partial t} = D_{az} \frac{\partial^2 C_j}{\partial z_2} - u_i \frac{\partial C_j}{\partial z} - \frac{(1-\varepsilon)}{\varepsilon} \rho_p f \frac{\partial q_j}{\partial t}$$

$$\frac{\partial q_j}{\partial t} = k_{p,j}* a_p*(q_j^* - q_j)$$

The variables are Cj and qj, the rest of terms are constants.

Here is the code I've written so far:

import numpy as np
from fipy import Grid1D, CellVariable, TransientTerm, DiffusionTerm, ConvectionTerm, Viewer

Parameters
D_az = 1e-9  # axial diffusion coefficient
u_i = [[1e-5]] # interstitial velocity
porosity = 0.4  # resin porosity
density = 1e3  # resin density
f = 0.5  # f coefficient
k_pj = 0.01  # exchange rate coefficient
a_p = 10  # a coefficient
q_j_star = 30  # maximum theoretical amount of ion that resin can adsorb under ideal conditions

Boundary and initial conditions
Cj_z0 = 10
Cj_zL = 3
qj_z0 = 1
qj_zL = 5
Cj_t0 = 10
qj_t0 = 1

grid = Grid1D(nx=50, dx=0.01)

Cj = CellVariable(name="Cj", hasOld=True, mesh=grid, value=Cj_t0)
qj = CellVariable(name="qj", hasOld=True, mesh=grid, value=qj_t0)

Define equations
eq1 = TransientTerm(var=qj) == k_pj * a_p * (q_j_star - qj)
eq2 = (TransientTerm() == DiffusionTerm(coeff=D_az) - ConvectionTerm(coeff=u_i)
       - ((1 - porosity) / porosity) * density * f * k_pj * a_p * (q_j_star - qj))
eqn = eq1 & eq2

Apply boundary conditions
Cj.constrain(Cj_z0, grid.facesLeft)
Cj.constrain(Cj_zL, grid.facesRight)
qj.constrain(qj_z0, grid.facesLeft)
qj.constrain(qj_zL, grid.facesRight)

vi = Viewer((Cj, qj))

from builtins import range
for t in range(1):
    qj.updateOld()
    Cj.updateOld()
    eqn.solve(dt=1.e-3)
    vi.plot()

I'm new at using fipy so there are likely several errors. I would greatly appreciate any advice or guidance on my issue. Additionally, I have attempted to resolve it in parallel using PETSc, but it's not working either.
Thank you very much in advance.

[guyer]

In

eq1 = TransientTerm(var=qj) == k_pj * a_p * (q_j_star - qj)
eq2 = (TransientTerm() == DiffusionTerm(coeff=D_az) - ConvectionTerm(coeff=u_i)
       - ((1 - porosity) / porosity) * density * f * k_pj * a_p * (q_j_star - qj))
eqn = eq1 & eq2

you specify which variable that the TransientTerm() in eq1 applies to, but you don't specify which variables any of the terms in eq2 apply to.

If you use var= for any Term, then you must use it for all of them. Further, var= is only optional for solving isolated equations; if you couple equations, you must specify which variable each Term applies to.

[frandgr]

Thank you very much for your response. I have changed that part of the script to this:

eq1 = TransientTerm(var=qj) == k_pj * a_p * (q_j_star - qj)
eq2 = (TransientTerm(var=Cj) == DiffusionTerm(coeff=D_az, var=Cj) - ConvectionTerm(coeff=u_i, var=Cj)
       - ((1 - porosity) / porosity) * density * f * k_pj * a_p * (q_j_star - qj))
eqn = eq1 & eq2 

But this is the plot that I obtain:

Of course, this is not the result I was expecting. There is no variation at all which is wrong. I'm still thinking that the problem is the definition of the equations. In eq2, does the script understands that qj is also varying with the time? Maybe I should define the equations in a different way.
Thank you again for your time.

[guyer]

There are a few issues:

• You should have received

  UserWarning: sweep() or solve() are likely to produce erroneous results when `var` does not contain floats.

  This is because you initialize Cj and qj with Cj_t0 and qj_t0, which are integers. Change these to

  diff --git a/coupled.py b/coupled.py
  index aac48c6..aa927a1 100644
  --- a/coupled.py
  +++ b/coupled.py
  @@ -16,8 +16,8 @@ Cj_z0 = 10
   Cj_zL = 3
   qj_z0 = 1
   qj_zL = 5
  -Cj_t0 = 10
  -qj_t0 = 1
  +Cj_t0 = 10.
  +qj_t0 = 1.
   
   grid = Grid1D(nx=50, dx=0.01) 

• No, in eq2, the equation does not understand that qj is varying in time. You have coupled the equations with &, but they don't actually share any variables, so they're not really coupled. This will help

  diff --git a/coupled.py b/coupled.py
  index aa927a1..ccaf35a 100644
  --- a/coupled.py
  +++ b/coupled.py
  @@ -27,7 +27,7 @@ qj = CellVariable(name="qj", hasOld=True, mesh=grid, value=qj_t0)
   # Define equations
   eq1 = TransientTerm(var=qj) == k_pj * a_p * (q_j_star - qj)
   eq2 = (TransientTerm(var=Cj) == DiffusionTerm(coeff=D_az, var=Cj) - ConvectionTerm(coeff=u_i, var=Cj)
  -       - ((1 - porosity) / porosity) * density * f * k_pj * a_p * (q_j_star - qj))
  +       - TransientTerm(coeff=((1 - porosity) / porosity) * density * f, var=qj))
   eqn = eq1 & eq2     

With those changes, qj evolves toward 30 and Cj evolves toward whatever it's doing. I suspect it's still not what you want, but it's what these equations with these coefficients want to do.