Properly, I think it should be (A*Ca*zmu).grad.dot([[1]]), but in 1D it shouldn't be any different.
On Mar 10, 2016, at 3:46 PM, Hervé Turlier <[email protected]> wrote: > Dear Jonathan, > > Thanks for your inputs. > While I agree that the ConvectionTerm(var=A, coeff=Ca[None]*zmu) would not > work for higher dimensions, this should be okay for now in 1-dimension, right > ? > Also, replacing, in this term, A and Ca should lead to the same numerical > result ? > > For 2-dimensions, would you have another formulation to propose for this term > ? > > Another unrelated question: shall I preferentially make implicit my > “reaction” term (= source term with opposite sign in two PDEs) ? > > Many thanks again. > Best, > > Hervé > >> On Mar 10, 2016, at 17:49, Guyer, Jonathan E. Dr. (Fed) >> <[email protected]> wrote: >> >>> I'm not entirely sure why the tuple version breaks, but the convection >>> terms do store different variables. >> >> A tuple of Variable fields is never a higher rank field. It might work >> occasionally, but it really shouldn't be counted on. >> >> I have no idea why V[None] produces a vector from a scalar, but it does. >> >> I don't think ConvectionTerm(var=A, coeff=Ca[None] * emu) is conceptually >> right. Hervé's ODE can be expressed in arbitrary dimensions as: >> >> \gamma\vec{v} = 2\eta\nabla(a\nabla\cdot\vec{v}) + \nabla(\zeta a c_a) >> >> The middle term can also be written as 2\eta\nabla\cdot(a\nabla\vec{v}), but >> the end result seems to be the same. Either way, we have vector = vector + >> vector. Casting the concentration to a vector seems wrong to me. I got lost >> in a sea of dyadic notation trying to work it out, though. >> >> >> On Mar 10, 2016, at 10:24 AM, Daniel Wheeler <[email protected]> >> wrote: >> >>> Hi Hervé, >>> >>> As you correctly stated, the issue is the scalar/vector velocity and >>> the convection terms. Writing the equations using `V[None]` worked for >>> me rather than using a tuple, `(V,)`, to make the velocity a vector. I >>> don't know why exactly. Try this though, >>> >>> ~~~~ >>> V_vec = V[None] >>> >>> # DEFINE EQUATIONS >>> R = ka * Ci * A/A0 - ki * Ca # explicit reaction term >>> EqCa = ( TransientTerm(var=Ca) == - >>> ExponentialConvectionTerm(var=Ca, coeff=V_vec) + DiffusionTerm(var=Ca, >>> coeff=Da) + R ) >>> EqCi = ( TransientTerm(var=Ci) == - >>> ExponentialConvectionTerm(var=Ci, coeff=V_vec) + DiffusionTerm(var=Ci, >>> coeff=Di) - R ) >>> EqA = ( TransientTerm(var=A) == - >>> ConvectionTerm(var=A, coeff=V_vec) - kd*(A-A0) >>> ) >>> EqV = ( ImplicitSourceTerm(var=V, coeff=gam) == >>> DiffusionTerm(var=V, coeff=4*eta*A ) + ConvectionTerm(var=A, >>> coeff=Ca[None] * zmu) ) >>> Eq = EqCa & EqCi & EqA & EqV >>> ~~~~ >>> >>> The above worked for me and ran without any issues, I have no idea >>> whether the solution is correct or not. >>> >>> I'm not entirely sure why the tuple version breaks, but the convection >>> terms do store different variables. >>> >>> In [1]: import fipy as fp >>> >>> In [2]: m = fp.Grid1D(nx=4) >>> >>> In [3]: v = fp.CellVariable(mesh=m) >>> >>> In [4]: fp.ExponentialConvectionTerm(var=v, coeff=(v,)) >>> Out[4]: ExponentialConvectionTerm(coeff=((CellVariable(value=array([ >>> 0., 0., 0., 0.]), mesh=UniformGrid1D(dx=1.0, nx=4)) * [1]),), >>> var=CellVariable(value=array([ 0., 0., 0., 0.]), >>> mesh=UniformGrid1D(dx=1.0, nx=4))) >>> >>> In [12]: fp.ExponentialConvectionTerm(var=v, coeff=v[None]) >>> Out[12]: >>> ExponentialConvectionTerm(coeff=_ArithmeticCellToFaceVariable(value=array([[ >>> 0., 0., 0., 0., 0.]]), mesh=UniformGrid1D(dx=1.0, nx=4)), >>> var=CellVariable(value=array([ 0., 0., 0., 0.]), >>> mesh=UniformGrid1D(dx=1.0, nx=4))) >>> >>> The `v[None]` version has already been bad into a FaceVariable of rank >>> 1. Ideally, both of these versions would be the same. That's an issue >>> to resolve, but you should be able to make progress now. >>> >>> Cheers, >>> >>> Daniel >>> >>> On Thu, Mar 10, 2016 at 5:40 AM, Hervé Turlier <[email protected]> >>> wrote: >>>> Dear FiPy developers, >>>> >>>> Thanks a lot for sharing the very useful FiPy software and providing >>>> reactive support. >>>> >>>> I’ve got a set of equations for which FiPy should be the right tool to use >>>> I >>>> think. >>>> >>>> I have 3 simple convection-reaction-diffusion PDEs (convection-diffusion >>>> with reactive source term) in 1-dimenion, which are coupled together and >>>> with an additional ODE for the (convection) velocity (which comes from >>>> force >>>> balance for an active viscous layer). >>>> To make it more clear I attach a PDF with the governing equations. >>>> >>>> I manage to treat convection-reaction-diffusion correctly, but I block on >>>> the coupling to the ODE for the convective velocity with FiPy, which >>>> includes 1st-order spatial derivatives of the other coupled variables. I >>>> tried to treat these spatial derivatives as source terms, but .grad leads >>>> to >>>> rank-1 terms, while the ODE is written for the velocity defined as a rank-0 >>>> CellVariable. In the actual implementation, I try, on the contrary to treat >>>> these 1st-order spatial derivatives as a convection term, but I get the >>>> error below. >>>> I attach my current Python script to make it more clear. >>>> >>>> Traceback (most recent call last): >>>> File >>>> "/Users/Turlier/Dropbox/Research/Projects/AdvectionReactionDiffusion/FiPy/coupled-advection-reaction-diffusion-mechanics_1D.py", >>>> line 83, in <module> >>>> res = Eq.sweep(dt=dt, solver=solver) >>>> File >>>> "/usr/local/lib/python2.7/site-packages/FiPy-unknown-py2.7.egg/fipy/terms/term.py", >>>> line 223, in sweep >>>> solver = self._prepareLinearSystem(var=var, solver=solver, >>>> boundaryConditions=boundaryConditions, dt=dt) >>>> File >>>> "/usr/local/lib/python2.7/site-packages/FiPy-unknown-py2.7.egg/fipy/terms/term.py", >>>> line 157, in _prepareLinearSystem >>>> buildExplicitIfOther=self._buildExplcitIfOther) >>>> File >>>> "/usr/local/lib/python2.7/site-packages/FiPy-unknown-py2.7.egg/fipy/terms/coupledBinaryTerm.py", >>>> line 122, in _buildAndAddMatrices >>>> buildExplicitIfOther=buildExplicitIfOther) >>>> File >>>> "/usr/local/lib/python2.7/site-packages/FiPy-unknown-py2.7.egg/fipy/terms/binaryTerm.py", >>>> line 68, in _buildAndAddMatrices >>>> buildExplicitIfOther=buildExplicitIfOther) >>>> File >>>> "/usr/local/lib/python2.7/site-packages/FiPy-unknown-py2.7.egg/fipy/terms/binaryTerm.py", >>>> line 68, in _buildAndAddMatrices >>>> buildExplicitIfOther=buildExplicitIfOther) >>>> File >>>> "/usr/local/lib/python2.7/site-packages/FiPy-unknown-py2.7.egg/fipy/terms/binaryTerm.py", >>>> line 68, in _buildAndAddMatrices >>>> buildExplicitIfOther=buildExplicitIfOther) >>>> File >>>> "/usr/local/lib/python2.7/site-packages/FiPy-unknown-py2.7.egg/fipy/terms/unaryTerm.py", >>>> line 99, in _buildAndAddMatrices >>>> diffusionGeomCoeff=diffusionGeomCoeff) >>>> File >>>> "/usr/local/lib/python2.7/site-packages/FiPy-unknown-py2.7.egg/fipy/terms/abstractConvectionTerm.py", >>>> line 191, in _buildMatrix >>>> var, L, b = FaceTerm._buildMatrix(self, var, SparseMatrix, >>>> boundaryConditions=boundaryConditions, dt=dt, >>>> transientGeomCoeff=transientGeomCoeff, >>>> diffusionGeomCoeff=diffusionGeomCoeff) >>>> File >>>> "/usr/local/lib/python2.7/site-packages/FiPy-unknown-py2.7.egg/fipy/terms/faceTerm.py", >>>> line 178, in _buildMatrix >>>> self._implicitBuildMatrix_(SparseMatrix, L, id1, id2, b, >>>> weight['implicit'], var, boundaryConditions, interiorFaces, dt) >>>> File >>>> "/usr/local/lib/python2.7/site-packages/FiPy-unknown-py2.7.egg/fipy/terms/faceTerm.py", >>>> line 74, in _implicitBuildMatrix_ >>>> L.addAt(numerix.take(coeffMatrix['cell 1 diag'], interiorFaces, >>>> axis=-1).ravel(), id1.ravel(), id1.swapaxes(0,1).ravel()) >>>> File >>>> "/usr/local/lib/python2.7/site-packages/FiPy-unknown-py2.7.egg/fipy/matrices/offsetSparseMatrix.py", >>>> line 58, in addAt >>>> SparseMatrix.addAt(self, vector, id1 + self.mesh.numberOfCells * >>>> self.equationIndex, id2 + self.mesh.numberOfCells * self.varIndex) >>>> File >>>> "/usr/local/lib/python2.7/site-packages/FiPy-unknown-py2.7.egg/fipy/matrices/pysparseMatrix.py", >>>> line 246, in addAt >>>> self.matrix.update_add_at(vector, id1, id2) >>>> IndexError: id1 does not have the same size as b >>>> >>>> I went through the manual and spent several hours on past discussions on >>>> the >>>> mailing-list, but I could not find, yet a working solution. >>>> >>>> Your help would be greatly appreciated. >>>> Best, >>>> >>>> Hervé >>>> >>>> --- >>>> Hervé Turlier, PhD >>>> European Molecular Biology Laboratory, >>>> Meyerhofstrasse 1 >>>> 69117 Heidelberg, Germany >>>> [email protected] >>>> [email protected] >>>> Tel: +49 (0)6221 387-8182 >>>> >>>> >>>> >>>> >>>> _______________________________________________ >>>> fipy mailing list >>>> [email protected] >>>> http://www.ctcms.nist.gov/fipy >>>> [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ] >>>> >>> >>> >>> >>> -- >>> Daniel Wheeler >>> >>> _______________________________________________ >>> 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 ] _______________________________________________ fipy mailing list [email protected] http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
