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 ]
