Jon, Thank you so much for your thorough input on this model. I have implemented your recommendations and I'm able to get a solution. Now I need to just need to troubleshoot the model so that I can get the _correct_ solution.
Again, I really appreciate your thoughtfulness and recommendations! Kind regards, Douwe Bruinsma On Fri, Sep 4, 2015 at 10:45 AM, Guyer, Jonathan E. Dr. < [email protected]> wrote: > 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 ] >
_______________________________________________ fipy mailing list [email protected] http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
