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 ]
