storing value of ImplicitSourceTerm

Question

I am working with fipy and I wish to simulate a flow with a free flux BC on some selected faces. Following other examples, I tried 2 different technics:

from fipy import *
from fipy.meshes.nonUniformGrid1D import NonUniformGrid1D as Grid1D
from fipy.tools import numerix as np

dx = [1.,  1. , 0.78584047] 
coeff = 2.
dt = min(0.1, 0.9 * min(dx) ** 2 / (2 * coeff))
steps = 10

mesh = Grid1D(nx = len(dx), dx = dx)

#phi values 
phi = CellVariable(mesh = mesh, value=[10., 0., 0.], hasOld = True)
phi.constrain(0., mesh.facesRight) #dirichlet BC
phi2 = CellVariable(mesh = mesh, value=[10., 0., 0.], hasOld = True)
phi2.constrain(0., mesh.facesRight) #dirichlet BC

#outflow value
phi2Out = CellVariable(mesh = mesh, value=[0., 0., 0.], hasOld = True)

intCoef = FaceVariable(mesh, value=coeff)
extCoef = FaceVariable(mesh, value=coeff)
intCoef.constrain(0., where = mesh.facesRight)
extCoef.constrain(0., where = ~mesh.facesRight)

eq1= TransientTerm() == DiffusionTerm(coeff= intCoef * phi.faceValue) \
  + ImplicitSourceTerm(coeff= (extCoef * phi.faceGrad).divergence)
  
eq2 =TransientTerm() == DiffusionTerm(coeff= intCoef * phi2.faceValue) \
  + phi2 * (extCoef * phi2.faceGrad).divergence
  
eq3 =TransientTerm() == - phi2 * (extCoef * phi2.faceGrad).divergence

phi.updateOld()
phi2.updateOld()
phi2Out.updateOld()
    
for _ in range(steps):
    res =  1e+10
    loop = 0
    while res > 1e-4 and loop < 1000:
        res =  eq1.sweep(var= phi,dt= dt)
        res2 =  eq2.sweep(var=phi2 ,dt= dt)
        res3 =  eq3.sweep(var=phi2Out ,dt= dt)
        res = max(res, res2, res3)
        loop += 1
    if res > 1e-4:
        print('no convergence! res: ', res)
        break
        
    phi.updateOld()
    phi2.updateOld()
    phi2Out.updateOld()
    
    print('\nphi: ', phi * mesh.cellVolumes, sum(phi* mesh.cellVolumes))
    print('phi2: ', phi2* mesh.cellVolumes,sum(phi2* mesh.cellVolumes))
    print('phi2_outFlow: ', phi2Out* mesh.cellVolumes, sum(phi2Out* mesh.cellVolumes))
    print('phi2_total: ', sum(phi2* mesh.cellVolumes + phi2Out* mesh.cellVolumes))

Output after 10 steps:

phi: [2.21786728 1.82700906 0.63381385] 4.678690190704105

phi2: [2.21786872 1.82701185 0.6338187 ] 4.678699267591382

phi2_outFlow: [0. 0. 5.32132192] 5.321321918864122

phi2_total: 10.000021186455504

Ideally, I would like to extract the value of the outflow to check the mass balance and know the exact value of each term for each cell (other sink/source terms will be added later, as well as other free-flow boundaries).

My questions are:

1. why are the values of phi and phi2 slightly different?
2. how could I extract the outflow term for each cell (when a more complex grid will be used) while keeping 'ImplicitSourceTerm', which is more efficient?

Thank you for your help!

PS: I have been working with python3 so far, so, although I am interested in solutions involving other solvers, I am especially looking for something that can be implemented with scipy.

Answer 1

1. why are the values of phi and phi2 slightly different?

phi and phi2 are different because eq2 doesn't converge as rapidly as eq1. This is because eq1 is more implicit than eq2. If you change the tolerance for the residual, e.g., res > 1e-10, you'll see the two solutions are in much closer agreement.

2. how could I extract the outflow term for each cell (when a more complex grid will be used) while keeping 'ImplicitSourceTerm', which is more efficient?

You can still evaluate the flux phi2 * extCoef * phi2.faceGrad, even when you use the ImplicitSourceTerm.

In general, it's not easy to extract what each Term is doing physically (see issue #461). You can use the FIPY_DISPLAY_MATRIX environment variable to see how each Term contributes to the solution matrix, but this may or may not give you much physical intuition for what's going on.