FIPY Solving a PDE coupled with an ODE

Question

I have a coupled set of equations where the main PDE (function of time and position z) is given as:

The second equation is of the type:

where k_m = f(q) and q^* = f(c). As you can see the second equation is an ODE (no dependence of q on space directly). I am finding it hard to write the code to couple the two equations. As of now for a simplistic case where I neglect the second equation and take: q = A*c, where A is some constant, I was able to simplify and just solve for the following convection diffusion equation:

with the following code:

from fipy import Variable, FaceVariable, CellVariable, Grid1D, ExplicitDiffusionTerm, TransientTerm, DiffusionTerm, Viewer, AdvectionTerm, PowerLawConvectionTerm, VanLeerConvectionTerm
from fipy.tools import numerix

#define the grid
L = 3.
nx = L * 512
dx = L/nx
mesh = Grid1D(dx=dx, nx=nx)

# create the variable and initiate it's value on the mesh
conc = CellVariable(name="Conc", mesh=mesh, value=0.)

# physical parameters
Dapp = 1e-7 
u = 0.1
A = 0.85
e = 0.4
F = (1-e)/e

# provide the simplified coefficients
DiffCoeff = Dapp/(1+A*F)
ConvCoeff = ((u/(1+A*F)),)

#Boundary conditions
valueLeft = 1
valueRight = 0.
conc.constrain(valueLeft, mesh.facesLeft)
conc.faceGrad.constrain(valueRight, where=mesh.facesRight)

# define the equation
eqX = TransientTerm() == (DiffusionTerm(coeff=DiffCoeff) - VanLeerConvectionTerm(coeff=ConvCoeff)) 

# time stepping parameters
timeStepDuration = 0.001
steps = 50000

from tqdm import tqdm
for step in tqdm(range(steps), desc="Iterating..."):
    eqX.solve(var=conc,dt=timeStepDuration)
    # plot every 5000 iterations
    if step%5000 == 0: 
        viewer.plot()

Can someone help in coupling the convection diffusion equation with the ODE in the fipy framework. I am bit confused about how to take the right hand side which in the Finite Volume sense should be just a source term.

_{(https://www.codecogs.com/latex/eqneditor.php for generating the Latex equations)}

Answer 1

Here is your problem reworked with something that might give you a clue.

from fipy import Variable, FaceVariable, CellVariable, Grid1D, ExplicitDiffusionTerm, TransientTerm, DiffusionTerm, Viewer, AdvectionTerm, PowerLawConvectionTerm, VanLeerConvectionTerm
from fipy.tools import numerix

#define the grid
L = 3.
nx = L * 512
dx = L/nx
mesh = Grid1D(dx=dx, nx=nx)


# create the variable and initiate it's value on the mesh
conc = CellVariable(name="Conc", mesh=mesh, value=0.)

# physical parameters
Dapp = 1e-7
u = 0.1
A = 0.85
e = 0.4
F = (1-e)/e

# provide the simplified coefficients
DiffCoeff = Dapp
ConvCoeff = ((u),)

#Boundary conditions
valueLeft = 1
valueRight = 0.
conc.constrain(valueLeft, mesh.facesLeft)
conc.faceGrad.constrain(valueRight, where=mesh.facesRight)

def q_star_func(conc):
    return some_conc_expression

q_star = Variable(q_star_func(conc))
q_star_old = Variable(q_star_func(conc))

# define the equation
eqX = TransientTerm() + F * (q_star - q_star_old) / dt == (DiffusionTerm(coeff=DiffCoeff) - VanLeerConvectionTerm(coeff=ConvCoeff))

# time stepping parameters
timeStepDuration = 0.001
steps = 50000

q_value = 1.

def k_m(q):
    return some_q_expression

viewer = Viewer(conc)
from tqdm import tqdm
for step in tqdm(range(steps), desc="Iterating..."):
    q_star.set_value(q_star_func(conc))
    eqX.solve(var=conc,dt=timeStepDuration)

    q_star_old.set_value(q_star.value)
    q_value = (q_value + timeStepDuration * k_m(q) * q_star.value) / (1 + timeStepDuration * k_m(q))

    # plot every 5000 iterations
    if step%5000 == 0:
        viewer.plot()

You can add a source term for the transient term for q_star and add the ODE into the loop. q_star is required to be a variable as it needs to update inside the equation. The q_value doesn't seem to have any feedback into the main equation for conc, which is a little odd. The functions q_star_func and k_m need to be defined. You might also use Scipy's ODE solver if you want a more sophisticated method for solving for q. You'd need to collect the values of q_star through time from solving the main equation for conc. The equations don't need to necessarily be solved together as the main equation is independent of q.

Answer 2

It's not clear to me what the difference is between q and q* in your two equations.

You have two equations for unknowns q, q*, c, and u. I don't know what form F and D take.

It looks to me like you're one or two equations short.

Can you tell what the variables mean (e.g. u == velocity, c == concentration, q == heat flux)? Maybe an answer would be easier if you said what physics you were trying to solve.