Staring at the code, it looks like you're defining initial values (C_a_0,
C_b_0, ..., V0) as CellVariables instead of scalars (or Variables), which is
probably not what you meant. Substituting the following seems to work better:
#%% - INITIATION VALUES
C_a_0 = 0
C_b_0 = 0
C_c_0 = 0
C_d_0 = 0
V0 = 409.5183981003883 #V_flow/A_crossSection #m/hr
#%% - VARIABLES
C_a = CellVariable(name="Concentration of A", mesh=mesh, hasOld=True)
C_b = CellVariable(name="Concentration of B", mesh=mesh, hasOld=True)
C_c = CellVariable(name="Concentration of C", mesh=mesh, hasOld=True)
C_d = CellVariable(name="Concentration of D", mesh=mesh, hasOld=True)
V = CellVariable(name="Velocity", mesh=mesh, hasOld=True)
C_a.setValue(C_a_0)
C_b.setValue(C_b_0)
C_c.setValue(C_c_0)
C_d.setValue(C_d_0)
V.setValue(V0)
But this still gives errors relating to division by zero and singular matrices.
I think this is because the last term in each of your equations (factor*r_a)
needs to be an appropriate SourceTerm, since
r_a = -k1 * C_a
and you're solving for C_a. I substituted
#%% - DEFINE EQUATIONS
EqC_a =
(TransientTerm(var=C_a)==-ExponentialConvectionTerm(coeff=V_vec,var=C_a) +
DiffusionTerm(var=C_b,coeff=Da) + ImplicitSourceTerm(var=C_a,coeff=1))
EqC_b =
(TransientTerm(var=C_b)==-ExponentialConvectionTerm(coeff=V_vec,var=C_b) +
DiffusionTerm(var=C_b,coeff=Da) + ImplicitSourceTerm(var=C_a,coeff=11))
EqC_c =
(TransientTerm(var=C_c)==-ExponentialConvectionTerm(coeff=V_vec,var=C_c) +
DiffusionTerm(var=C_c,coeff=Da) + ImplicitSourceTerm(var=C_a,coeff=-17))
EqC_d =
(TransientTerm(var=C_d)==-ExponentialConvectionTerm(coeff=V_vec,var=C_d) +
DiffusionTerm(var=C_d,coeff=Da) + ImplicitSourceTerm(var=C_a,coeff=-8))
Eq = EqC_a & EqC_b & EqC_c & EqC_d
This still returns an error,
/users/tnk10/Downloads/fipy/fipy/tools/numerix.py:966: RuntimeWarning: overflow
encountered in square
return sqrt(add.reduce(arr**2))
/users/tnk10/Downloads/fipy/fipy/solvers/pysparse/linearLUSolver.py:93:
RuntimeWarning: overflow encountered in square
error0 = numerix.sqrt(numerix.sum((L * x - b)**2))
/users/tnk10/Downloads/fipy/fipy/solvers/pysparse/linearLUSolver.py:98:
RuntimeWarning: overflow encountered in square
if (numerix.sqrt(numerix.sum(errorVector**2)) / error0) <= self.tolerance:
/users/tnk10/Downloads/fipy/fipy/solvers/pysparse/linearLUSolver.py:98:
RuntimeWarning: invalid value encountered in double_scalars
if (numerix.sqrt(numerix.sum(errorVector**2)) / error0) <= self.tolerance:
/users/tnk10/Downloads/fipy/fipy/variables/faceGradVariable.py:158:
RuntimeWarning: overflow encountered in divide
N = (N2 - numerix.take(self.var,id1, axis=-1)) / dAP
/users/tnk10/Downloads/fipy/fipy/variables/variable.py:1105: RuntimeWarning:
invalid value encountered in multiply
return self._BinaryOperatorVariable(lambda a,b: a*b, other)
Assuming (perhaps incorrectly) that component A should be diffusing, I changed
the DiffusionTerm
EqC_a =
(TransientTerm(var=C_a)==-ExponentialConvectionTerm(coeff=V_vec,var=C_a) +
DiffusionTerm(var=C_a,coeff=Da) + ImplicitSourceTerm(var=C_a,coeff=1))
That seems to work. The modified script is attached.
Trevor
________________________________
From: fipy-boun...@nist.gov <fipy-boun...@nist.gov> on behalf of Keller, Trevor
(Fed) <trevor.kel...@nist.gov>
Sent: Wednesday, July 13, 2016 4:26:38 PM
To: FIPY
Subject: Re: Diffusion-Convection-Reactive Chemisty
Is the definition of C_a_BC correct? For a 1D grid, is the behavior of
C_a.constrain(C_a_BC, where=mesh.facesRight)
with a CellVariable instead of a scalar value
C_a_BC= C_a_0*(1-X)
meaningful?
Trevor
________________________________
From: fipy-boun...@nist.gov <fipy-boun...@nist.gov> on behalf of Daniel
DeSantis <desan...@gmail.com>
Sent: Wednesday, July 13, 2016 4:08:50 PM
To: FIPY
Subject: Re: Diffusion-Convection-Reactive Chemisty
I'm sorry. I was trying to fix the problem, and forgot to put a line back in,
which was masked by me not clearing a variable value for V. My apologies.
Please try this code instead. It has the traceback I mentioned before.
Traceback (most recent call last):
File "<ipython-input-1-d209891ef431>", line 1, in <module>
runfile('C:/Users/ddesantis/Dropbox/PythonScripts/CFD
Models/ConversionModel.py', wdir='C:/Users/ddesantis/Dropbox/PythonScripts/CFD
Models')
File
"C:\Users\ddesantis\AppData\Local\Continuum\Anaconda2\lib\site-packages\spyderlib\widgets\externalshell\sitecustomize.py",
line 699, in runfile
execfile(filename, namespace)
File
"C:\Users\ddesantis\AppData\Local\Continuum\Anaconda2\lib\site-packages\spyderlib\widgets\externalshell\sitecustomize.py",
line 74, in execfile
exec(compile(scripttext, filename, 'exec'), glob, loc)
File "C:/Users/ddesantis/Dropbox/PythonScripts/CFD
Models/ConversionModel.py", line 110, in <module>
res = Eq.sweep(dt=dt, solver=solver)
File
"C:\Users\ddesantis\AppData\Local\Continuum\Anaconda2\lib\site-packages\fipy\terms\term.py",
line 236, in sweep
solver = self._prepareLinearSystem(var=var, solver=solver,
boundaryConditions=boundaryConditions, dt=dt)
File
"C:\Users\ddesantis\AppData\Local\Continuum\Anaconda2\lib\site-packages\fipy\terms\term.py",
line 170, in _prepareLinearSystem
buildExplicitIfOther=self._buildExplcitIfOther)
File
"C:\Users\ddesantis\AppData\Local\Continuum\Anaconda2\lib\site-packages\fipy\terms\coupledBinaryTerm.py",
line 122, in _buildAndAddMatrices
buildExplicitIfOther=buildExplicitIfOther)
File
"C:\Users\ddesantis\AppData\Local\Continuum\Anaconda2\lib\site-packages\fipy\terms\binaryTerm.py",
line 68, in _buildAndAddMatrices
buildExplicitIfOther=buildExplicitIfOther)
File
"C:\Users\ddesantis\AppData\Local\Continuum\Anaconda2\lib\site-packages\fipy\terms\binaryTerm.py",
line 68, in _buildAndAddMatrices
buildExplicitIfOther=buildExplicitIfOther)
File
"C:\Users\ddesantis\AppData\Local\Continuum\Anaconda2\lib\site-packages\fipy\terms\binaryTerm.py",
line 68, in _buildAndAddMatrices
buildExplicitIfOther=buildExplicitIfOther)
File
"C:\Users\ddesantis\AppData\Local\Continuum\Anaconda2\lib\site-packages\fipy\terms\unaryTerm.py",
line 99, in _buildAndAddMatrices
diffusionGeomCoeff=diffusionGeomCoeff)
File
"C:\Users\ddesantis\AppData\Local\Continuum\Anaconda2\lib\site-packages\fipy\terms\abstractConvectionTerm.py",
line 211, in _buildMatrix
self.constraintB = -((1 - alpha) * var.arithmeticFaceValue *
constraintMask * exteriorCoeff).divergence * mesh.cellVolumes
File
"C:\Users\ddesantis\AppData\Local\Continuum\Anaconda2\lib\site-packages\fipy\variables\variable.py",
line 1151, in __mul__
return self._BinaryOperatorVariable(lambda a,b: a*b, other)
File
"C:\Users\ddesantis\AppData\Local\Continuum\Anaconda2\lib\site-packages\fipy\variables\variable.py",
line 1116, in _BinaryOperatorVariable
if not v.unit.isDimensionless() or len(v.shape) > 3:
File
"C:\Users\ddesantis\AppData\Local\Continuum\Anaconda2\lib\site-packages\fipy\variables\variable.py",
line 255, in _getUnit
return self._extractUnit(self.value)
File
"C:\Users\ddesantis\AppData\Local\Continuum\Anaconda2\lib\site-packages\fipy\variables\variable.py",
line 561, in _getValue
value[..., mask] = numerix.array(constraint.value)[..., mask]
IndexError: index 100 is out of bounds for axis 0 with size 100
On Wed, Jul 13, 2016 at 3:59 PM, Daniel Wheeler
<daniel.wheel...@gmail.com<mailto:daniel.wheel...@gmail.com>> wrote:
Hi Daniel,
It is giving a different error for me
Traceback (most recent call last):
File "ConversionModel.py", line 78, in <module>
V.constrain(V0,mesh.facesLeft)
NameError: name 'V' is not defined
Is this the correct script?
On Wed, Jul 13, 2016 at 3:27 PM, Daniel DeSantis
<desan...@gmail.com<mailto:desan...@gmail.com>> wrote:
> Hello,
>
> I am having some trouble getting a workable convection coefficient on a
> reactive chemistry model with a vector velocity. I have reviewed several
> examples from the FiPy mailing list and tried several of the variations
> listed there. Originally, I was receiving the error that a coefficient must
> be a vector. Now I seem to be getting an error that says the index is out of
> bounds. Specific traceback is shown below, along with the code attachment.
>
> Could anyone suggest what I am doing wrong and how to fix it?
>
> I'm relatively new to Python and FiPy, specifically, so please bear with me
> as I am learning.
>
> Thank you all!
>
> --
> Daniel DeSantis
>
> Traceback (most recent call last):
>
> File "<ipython-input-22-d209891ef431>", line 1, in <module>
> runfile('C:/Users/ddesantis/Dropbox/PythonScripts/CFD
> Models/ConversionModel.py',
> wdir='C:/Users/ddesantis/Dropbox/PythonScripts/CFD Models')
>
> File
> "C:\Users\ddesantis\AppData\Local\Continuum\Anaconda2\lib\site-packages\spyderlib\widgets\externalshell\sitecustomize.py",
> line 699, in runfile
> execfile(filename, namespace)
>
> File
> "C:\Users\ddesantis\AppData\Local\Continuum\Anaconda2\lib\site-packages\spyderlib\widgets\externalshell\sitecustomize.py",
> line 74, in execfile
> exec(compile(scripttext, filename, 'exec'), glob, loc)
>
> File "C:/Users/ddesantis/Dropbox/PythonScripts/CFD
> Models/ConversionModel.py", line 107, in <module>
> res = Eq.sweep(dt=dt, solver=solver)
>
> File
> "C:\Users\ddesantis\AppData\Local\Continuum\Anaconda2\lib\site-packages\fipy\terms\term.py",
> line 236, in sweep
> solver = self._prepareLinearSystem(var=var, solver=solver,
> boundaryConditions=boundaryConditions, dt=dt)
>
> File
> "C:\Users\ddesantis\AppData\Local\Continuum\Anaconda2\lib\site-packages\fipy\terms\term.py",
> line 170, in _prepareLinearSystem
> buildExplicitIfOther=self._buildExplcitIfOther)
>
> File
> "C:\Users\ddesantis\AppData\Local\Continuum\Anaconda2\lib\site-packages\fipy\terms\coupledBinaryTerm.py",
> line 122, in _buildAndAddMatrices
> buildExplicitIfOther=buildExplicitIfOther)
>
> File
> "C:\Users\ddesantis\AppData\Local\Continuum\Anaconda2\lib\site-packages\fipy\terms\binaryTerm.py",
> line 68, in _buildAndAddMatrices
> buildExplicitIfOther=buildExplicitIfOther)
>
> File
> "C:\Users\ddesantis\AppData\Local\Continuum\Anaconda2\lib\site-packages\fipy\terms\binaryTerm.py",
> line 68, in _buildAndAddMatrices
> buildExplicitIfOther=buildExplicitIfOther)
>
> File
> "C:\Users\ddesantis\AppData\Local\Continuum\Anaconda2\lib\site-packages\fipy\terms\binaryTerm.py",
> line 68, in _buildAndAddMatrices
> buildExplicitIfOther=buildExplicitIfOther)
>
> File
> "C:\Users\ddesantis\AppData\Local\Continuum\Anaconda2\lib\site-packages\fipy\terms\unaryTerm.py",
> line 99, in _buildAndAddMatrices
> diffusionGeomCoeff=diffusionGeomCoeff)
>
> File
> "C:\Users\ddesantis\AppData\Local\Continuum\Anaconda2\lib\site-packages\fipy\terms\abstractConvectionTerm.py",
> line 211, in _buildMatrix
> self.constraintB = -((1 - alpha) * var.arithmeticFaceValue *
> constraintMask * exteriorCoeff).divergence * mesh.cellVolumes
>
> File
> "C:\Users\ddesantis\AppData\Local\Continuum\Anaconda2\lib\site-packages\fipy\variables\variable.py",
> line 1151, in __mul__
> return self._BinaryOperatorVariable(lambda a,b: a*b, other)
>
> File
> "C:\Users\ddesantis\AppData\Local\Continuum\Anaconda2\lib\site-packages\fipy\variables\variable.py",
> line 1116, in _BinaryOperatorVariable
> if not v.unit.isDimensionless() or len(v.shape) > 3:
>
> File
> "C:\Users\ddesantis\AppData\Local\Continuum\Anaconda2\lib\site-packages\fipy\variables\variable.py",
> line 255, in _getUnit
> return self._extractUnit(self.value)
>
> File
> "C:\Users\ddesantis\AppData\Local\Continuum\Anaconda2\lib\site-packages\fipy\variables\variable.py",
> line 561, in _getValue
> value[..., mask] = numerix.array(constraint.value)[..., mask]
>
> IndexError: index 100 is out of bounds for axis 0 with size 100
>
> _______________________________________________
> fipy mailing list
> fipy@nist.gov<mailto:fipy@nist.gov>
> http://www.ctcms.nist.gov/fipy
> [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
>
--
Daniel Wheeler
_______________________________________________
fipy mailing list
fipy@nist.gov<mailto:fipy@nist.gov>
http://www.ctcms.nist.gov/fipy
[ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
--
Daniel DeSantis
# -*- coding: utf-8 -*-
"""
Created on Wed Jul 13 14:05:25 2016
@author: ddesantis
"""
# -*- coding: utf-8 -*-
"""
Created on Tue Jul 12 13:32:09 2016
@author: ddesantis
"""
#%% - IMPORTS
import numpy as np
from fipy import *
#%% - REACTOR DIMENSIONS
d = 0.4188 #m, diameter
L = 1.2563 #m, length
r = d / 2 #m, radius
A_crossSection = np.pi * r**2 #m^2, cross section
#%% - MESH
nx = 100
dx = L / nx
mesh = Grid1D(nx=nx, dx=dx)
x = mesh.cellCenters
#%% - INITIATION VALUES
C_a_0 = 0 #0.5
C_b_0 = 0 #0.375
C_c_0 = 0 #0.12
C_d_0 = 0 #0.005
V0 = 409.5183981003883 #V_flow/A_crossSection #m/hr
#C_a_0.setValue(0.5)
#C_b_0.setValue(0)
#C_c_0.setValue(0)
#C_d_0.setValue(0)
#V0.setValue(409.5183981003883) #V_flow/A_crossSection #m/hr
#%% - VARIABLES
C_a = CellVariable(name="Concentration of A", mesh=mesh, hasOld=True)
C_b = CellVariable(name="Concentration of B", mesh=mesh, hasOld=True)
C_c = CellVariable(name="Concentration of C", mesh=mesh, hasOld=True)
C_d = CellVariable(name="Concentration of D", mesh=mesh, hasOld=True)
V = CellVariable(name="Velocity", mesh=mesh, hasOld=True)
C_a.setValue(C_a_0)
C_b.setValue(C_b_0)
C_c.setValue(C_c_0)
C_d.setValue(C_d_0)
V.setValue(V0)
#%% - REACTION PARAMETERS
Da = 1 #Assumed diffusion Coefficient for all chemicals
k1 = 2000 #WAG
X = 0.8 #Conversion
Y = 0.8 #Yield
r_a = -k1*C_a
#%% - SET UP THE VIEWER
#if __name__ == '__main__':
# viewer = Matplotlib1DViewer(vars=(C_a,C_b,C_c,C_d))
# viewer.plot()
#%% - BOUNDARY CONDITIONS
# Feeds (Left side BCs)
C_a_f = 3.19 #kmol/hr
C_b_f = (C_a_f*11)*1.15 #kmol/hr
C_c_f = 0 #kmol/hr
C_d_f = 0 #kmol/hr
C_a_BC= C_a_0*(1-X) #Assumed amount of A converted, reactor right BC
#FiPy Constraints
C_a.constrain(C_a_BC, where=mesh.facesRight)
C_a.constrain(C_a_f, where=mesh.facesLeft)
C_b.constrain(C_b_f, where=mesh.facesLeft)
C_c.constrain(C_c_f, where=mesh.facesLeft)
C_d.constrain(C_d_f, where=mesh.facesLeft)
V.constrain(V0, where=mesh.facesLeft)
#%% - DEFINE EQUATIONS
V_vec=V[None]
#EqC_a = (TransientTerm(var=C_a)==-ExponentialConvectionTerm(coeff=V_vec,var=C_a) + DiffusionTerm(var=C_b,coeff=Da) + r_a)
#EqC_b = (TransientTerm(var=C_b)==-ExponentialConvectionTerm(coeff=V_vec,var=C_b) + DiffusionTerm(var=C_b,coeff=Da) + 11*r_a)
#EqC_c = (TransientTerm(var=C_c)==-ExponentialConvectionTerm(coeff=V_vec,var=C_c) + DiffusionTerm(var=C_c,coeff=Da) + -17*r_a)
#EqC_d = (TransientTerm(var=C_d)==-ExponentialConvectionTerm(coeff=V_vec,var=C_d) + DiffusionTerm(var=C_d,coeff=Da) + -8*r_a)
EqC_a = (TransientTerm(var=C_a)==-ExponentialConvectionTerm(coeff=V_vec,var=C_a) + DiffusionTerm(var=C_a,coeff=Da) + ImplicitSourceTerm(var=C_a,coeff=1))
EqC_b = (TransientTerm(var=C_b)==-ExponentialConvectionTerm(coeff=V_vec,var=C_b) + DiffusionTerm(var=C_b,coeff=Da) + ImplicitSourceTerm(var=C_a,coeff=11))
EqC_c = (TransientTerm(var=C_c)==-ExponentialConvectionTerm(coeff=V_vec,var=C_c) + DiffusionTerm(var=C_c,coeff=Da) + ImplicitSourceTerm(var=C_a,coeff=-17))
EqC_d = (TransientTerm(var=C_d)==-ExponentialConvectionTerm(coeff=V_vec,var=C_d) + DiffusionTerm(var=C_d,coeff=Da) + ImplicitSourceTerm(var=C_a,coeff=-8))
Eq = EqC_a & EqC_b & EqC_c & EqC_d
# SOLVE EQUATIONS
elapsed = 0.
dt = 1.
solver = LinearLUSolver(tolerance=1e-10)
while elapsed < 10.:
elapsed += dt
C_a.updateOld()
C_b.updateOld()
C_c.updateOld()
C_d.updateOld()
res = 1.
while res > 1.e-3:
res = Eq.sweep(dt=dt, solver=solver)
if __name__ == '__main__':
#viewer.plot()
raw_input("Press <return> for next step...")
print C_a
print V
_______________________________________________
fipy mailing list
fipy@nist.gov
http://www.ctcms.nist.gov/fipy
[ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
