Dear FiPy users and developers,
first of all I would like to thank the developers for this great
software package.
I am faced with a convection-diffusion equation
$$
\frac{\partial \phi}{\partial t} = D \nabla^2 \phi + \nabla (\phi
\nabla V_{ext})
$$
with a strongly inhomogeneous concentration profile \phi in an external
potential V_ext. I would like to determine the flux. Mathematically, I
would have to calculate
$$
j = -\nabla \phi - \phi \nabla V_{ext}
$$
However, a naive calculation of the flux in FiPy ("-phi.grad - V.grad *
phi" or "-phi.faceGrad - V.faceGrad * phi.faceValue") is strongly
deviating from what is expected, especially where the gradient in phi is
steep. See below for a minimal example.
Considering the numerical finite volume scheme, I guess what I would
want to calculate is
$$
(\vec{n} \cdot \vec{u})_f (\alpha_f \phi_P - (1-\alpha_f) \phi_A) +
\Gamma_f \frac{\phi_A - \phi_P}{d_{AP}},
$$
where \vec{u} is the gradient of V_{ext}. Summation of this field gives
the change of \phi in the cells.
Indeed, if I calculate the flux this way I get what I would expect.
However, to get the values for \alpha_f I extracted them from the
stencil in the ConvectionTerm. This is fine in the case of the simple
example given below but it becomes error-prone in a more complicated
context.
So my questions are:
- Do these thoughts make sense? Is the physical flux well approximated
by the flux from the numerical scheme?
- Am I missing some API to get the flux directly? Is there a more
convenient way to get the flux as calculated from the numerical scheme?
Thanks in advance and best wishes,
Urs
##### Example ###################
# Equilibrium concentration profile
# in a harmonic potential
##############################
from fipy import *
nx = 100
dx = 1.0
mesh = Grid1D(nx=nx, dx=dx)
x, = mesh.cellCenters()
D = 1.0
L = nx*dx
phi = CellVariable(name="solution variable", mesh=mesh, value=1.0,
hasOld=True)
V_ext = CellVariable(name="external potential", mesh=mesh, value=(x -
L/2.0)**2)
eq = TransientTerm() == DiffusionTerm(coeff=D) +
ExponentialConvectionTerm(coeff=V_ext.faceGrad)
# Solve for large timesteps to get to equilibrium
for step in range(10):
res = 10
phi.updateOld()
while res > 1e-10:
res = eq.sweep(var=phi, dt=1000.0)
print(res)
print("next step")
# Naive flux calculation
simple_flux1 = -phi.grad[0] - phi*(V_ext.grad[0])
simple_flux2 = (-phi.faceGrad - phi.faceValue * V_ext.faceGrad)[0]
# Calculate according to numeric CC-FVM scheme
# get stencil from ConvectionTerm to extract \alpha_f weight factor
stencil = eq.other.other.stencil['implicit']
diffusive_flux = -phi.faceGrad[0]
conv_flux = V_ext.faceGrad[0][1:-1] * (stencil['cell 1
diag'][1:-1]*phi.value[:-1] + stencil['cell 1 offdiag'][1:-1]*phi.value[1:])
CCFVM_flux = diffusive_flux[1:-1] - conv_flux
# Since this is the equilibrium configuration flux should be 0
print("-------------------------------------------------")
print("Cell centered naive flux calculation:
max(|j|)={0}".format(abs(simple_flux1).max()))
print("Face centered naive flux calculation:
max(|j|)={0}".format(abs(simple_flux2).max()))
print("Face centered CC-FVM flux calculation:
max(|j|)={0}".format(abs(CCFVM_flux).max()))
_______________________________________________
fipy mailing list
[email protected]
http://www.ctcms.nist.gov/fipy
[ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
