Hi Daniel,

Wonderful ! Indeed that works like you propose.
Many thanks.
Best,

Hervé

--- 
Hervé Turlier, PhD
European Molecular Biology Laboratory, 
Meyerhofstrasse 1 
69117 Heidelberg, Germany
[email protected]
[email protected]
Tel: +49 (0)6221 387-8182

> On Mar 10, 2016, at 16:24, 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 ]

Reply via email to