CODEWITHSUNDEEP
pythonpdefipy
kcast49
kcast49

Reputation: 1

FiPy Solving Couple PDEs Results in IndexError

I am new to FiPy and am trying to solve a pair of coupled PDEs describing adsorption of a material to an adsorbent packed in a cylinder. The material is flowing through the cylinder and its concentration is in either of 2 phases: a mobile phase or a stationary phase. I took the mesh for the mobile phase to be a 1D grid along the long axis pf the cylinder. I took the mesh for stationary phase to be a cylindrical 1D grid for the cross-sectional area of the cylinder. The 2 phases are coupled and I'm trying to solve for the concentration of the material in the mobile phase over time.

import numpy as np
import fipy as fp

# meshes and variables
nx = 50.
dx = 1.

# mobile phase
mesh_m = fp.Grid1D(nx=nx, dx=dx)
cm = fp.CellVariable(name="Concentration mobile", mesh=mesh_m, value=0.)

# stationary phase
mesh_s = fp.CylindricalGrid1D(nr=50, dr=dx)
cs = fp.CellVariable(name="Concentration stationary", mesh=mesh_s, value=0.)

# phase equations 
# mobile
eqCm = fp.TransientTerm(coeff = 0.5, var = cm) \
    == fp.ExplicitDiffusionTerm(coeff = 0.03, var = cm) \
    - fp.PowerLawConvectionTerm((1,), var = cm) \
    - fp.ImplicitSourceTerm(coeff= 3 * (cm - cs), var = cm)  

# stationary
eqCs = fp.TransientTerm(coeff = 0.4, var = cs) \
    == fp.ExplicitDiffusionTerm(coeff = 0.2, var = cs) 
   
# boundary conditions 
# mobile
cm_Left = 1.
cm_Right = 0.
cm_gradValueLeft = 30 * (cm_Left - 1)
cm_gradValueRight = 0.

cm.constrain(cm_Left, mesh_m.facesLeft)

cm.faceGrad.constrain(cm_gradValueLeft, mesh_m.facesLeft)
cm.faceGrad.constrain(cm_gradValueRight, mesh_m.facesRight)

# stationary
cs_Left = 0.
cs_Right = 1.
cs_gradValueLeft = 0.
cs_gradValueRight = 0.01 * (cm - cs)

cs.constrain(cs_Right, mesh_s.facesRight)

cs.faceGrad.constrain(cs_gradValueLeft, mesh_s.facesLeft)
cs.faceGrad.constrain(cs_gradValueRight, mesh_s.facesRight)

# **************** setup system and solve *********************************
eqn = eqCm and eqCs

T = 100
timeStep = 0.5
steps = int(T/timeStep)
viewer = fp.Viewer(vars=(cm, cs), datamin=0., datamax=4)

for step in range(steps):
    eqn.solve(var=cm, dt=timeStep)
    print(cm)
    viewer.plot()

I keep running into the following error during the eqn.solve(var=cm, dt=timeStep) command: IndexError: boolean index did not match indexed array along dimension 0; dimension is 50 but corresponding boolean dimension is 51. Is there something wrong with the way the cell variables or the equations are set up? Any help on this is appreciated.

Upvotes: 0

Views: 77

Answers (1)

jeguyer
jeguyer

Reputation: 2484

The error is very likely because cm and cs have 50 cells, leading to cs_gradValueRight having 50 cells. mesh_s.facesRight accesses 51 faces, though, so cs.faceGrad.constrain(cs_gradValueRight, mesh_s.facesRight) doesn't know what to do. To get rid of this error, you could write:

cs.faceGrad.constrain(cs_gradValueRight.faceValue, where=mesh_s.facesRight)

You're still going to have problems, though:

  • FiPy cannot couple equations between different meshes. In general, meshes can be concatenated, but I seriously doubt that a Grid1D can be concatenated to a CylindricalGrid1D. I'm pretty confident that it will lose the cylindrical geometry.
  • FiPy couples equations with "&", not "and".

If I understand what you're trying to do, I think you will have much better luck solving a single concentration field on a CylindricalGrid2D and changing the properties of the medium between the mobile region and the stationary region.

Not a problem, but not not ideal: Use DiffusionTerm rather than ExplicitDiffusionTerm. It will solve much more efficiently.

Upvotes: 0

Related Questions