> 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 ]
