Douwe -
Sorry for the delay in answering. I had to think through this a bit.
The most immediate cause of the error you received is that you are attempting
to solve for u declared as a FaceVariable. FiPy can only solve equations for
CellVariables. This is easily remedied by changing the declaration in line 80
to:
u = CellVariable(name ="velocity",mesh=mesh,rank=1)
and line 102 to
ExponentialConvectionTerm(coeff=epsilon*u.faceValue, var=y) + \
to make the error go away.
There are deeper problems than this, though.
\epsilon \vec u \frac{\partial y}{\partial z}
is not represented by an ExponentialConvectionTerm. A FiPy ConvectionTerm
represents
\frac{\partial}{\partial z}\left(\epsilon u_z y\right)
or, more generally,
\nabla\cdot\left(\epsilon \vec u y\right)
(I don't like working out problems in 1D with partial derivatives because it's
too easy to morph vectors into scalars and back. At a minimum, you need to
rigorously track your \hat{k} unit vectors).
Looking at your equations, the combination of (1) and (2) gives a FiPy
ConvectionTerm: [*]
1')
epsilon_t \frac{\partial y}{\partial t}+ \epsilon \frac{\partial u_z
y}{\partial z} +A \frac{\partial q_1}{\partial t} = 0
You can still use (2) to determine your velocity u, but \epsilon \frac{\partial
u}{\partial z} is not a ConvectionTerm for u. In 1D, you can integrate this
equation to obtain
u_z = \int -\frac{A}{\epsilon} \left(\frac{\partial q_1}{\partial
t}+\frac{\partial q_2}{\partial t} \right ) dz
and use uH to establish the constant of integration, rather than as a FiPy
constraint.
I've only skimmed the paper you cite in the code by Sun, le Queré & Levan, but
I note that it says in the "Numerical Method" section:
The variation of the velocity is caused by mass transfer and adsorption
inside the adsorbent particles and can be obtained by integration of eq.
(3).
In higher dimensions, (2) is not directly integrable and, as I understand it,
you need to use some form of the momentum equation to obtain solution for \vec
u.
Integration is not something FiPy does on its own, but you can use the
scipy.integrate.cumtrapz() function to achieve this. The best way to implement
this (IMO) is to create a new FiPy Variable class that does the integration.
I've shown my approach starting at line #93 of http://pastebin.com/QN1xWyFa.
I can't say that this implementation is *right*, but it at least solves.
- Jon
[*] In general, FiPy prefers that you solve Eqs. (2) & (3) of Sun, le Queré &
Levan, rather than (3) & (4).
On Aug 28, 2015, at 11:48 AM, Douwe Bruinsma <[email protected]> wrote:
> Hi all,
>
> I have a complete installation of fipy that is working correctly but I'm
> having trouble getting my model to work.
>
> I'm modeling langmuir adsorption of a binary gas flowing through a 1D space.
> There are four equations that need to be solved
>
> 1) gas concentration of component 1
> \epsilon_t \frac{\partial y}{\partial t}+\epsilon \vec u \frac{\partial
> y}{\partial z}+A \frac{\partial q_1}{\partial t}-A\left( \frac{\partial
> q_1}{\partial t}+\frac{\partial q_2}{\partial t}\right)y = 0
>
> 2) gas velocity and total mass balance:
> \epsilon \frac{\partial u}{\partial z} =-A\left(\frac{\partial q_1}{\partial
> t}+\frac{\partial q_2}{\partial t} \right )
>
> 3) langmuir adsorption for both species:
> q_1=\frac{D_1y}{1+B_1y+C_1\left(1-y \right )}
> q_2=\frac{D_2(1-y)}{1+B_2y+C_2\left(1-y \right )}
>
> The code that I'm trying to execute is available at the following link:
> http://pastebin.com/CYJkxUd1
>
> Just as a side note, the constants in the expressions above have different
> values in the code.
>
> When I try to run the code I get the following error trail:
>
> Traceback (most recent call last):
> File "/home/douwe/adsorption model/build-up.py", line 126, in <module>
> equ.solve(var=u,dt=timeStep)
> File "/usr/local/lib/python2.7/dist-packages/fipy/terms/term.py", line 211,
> in solve
> solver = self._prepareLinearSystem(var, solver, boundaryConditions, dt)
> File "/usr/local/lib/python2.7/dist-packages/fipy/terms/term.py", line 170,
> in _prepareLinearSystem
> buildExplicitIfOther=self._buildExplcitIfOther)
> File "/usr/local/lib/python2.7/dist-packages/fipy/terms/binaryTerm.py",
> line 68, in _buildAndAddMatrices
> buildExplicitIfOther=buildExplicitIfOther)
> File "/usr/local/lib/python2.7/dist-packages/fipy/terms/unaryTerm.py", line
> 99, in _buildAndAddMatrices
> diffusionGeomCoeff=diffusionGeomCoeff)
> File
> "/usr/local/lib/python2.7/dist-packages/fipy/terms/abstractConvectionTerm.py",
> line 199, in _buildMatrix
> constraintMask = var.faceGrad.constraintMask |
> var.arithmeticFaceValue.constraintMask
> AttributeError: 'FaceVariable' object has no attribute 'faceGrad'
>
> I think the problem lies with the expression for du/dz, but I'm not sure how
> to fix it after trying several approaches. The problem may also be the
> linking of the transient terms for q1 and q2 between the expressions. Should
> I use a transientTerm for those or should I use a discrete difference using
> q1 and q1.old?
>
> I appreciate any input you have.
>
> Kind regards,
>
> Douwe
>
> _______________________________________________
> fipy mailing list
> [email protected]
> http://www.ctcms.nist.gov/fipy
> [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
_______________________________________________
fipy mailing list
[email protected]
http://www.ctcms.nist.gov/fipy
[ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
