> > I am trying to solve a pure convection problem with a dirac-shaped
> > source term depending both in size and position on the solution
> > variable. The problem is that the source term becomes a
> > OperatorVariable, while a CellVariable is required by the equation. If
> > I try to do something like this
> >
> >>>> import fipy as fp
> >>>> mesh = fp.Grid1D(dx=L/nr, nx=nr)
> >>>> b = mesh.getFaceCenters()[0]
> >>>> dr = mesh.getCellVolumes()
> >>>> j = some_scalar_function(phi)
> >>>> rbirth = another_scalar_function(phi) S = j / dr * (b[:-1] <=
> >>>> rbirth) * (rbirth < b[1:])
>
> The issue is not that S is an OperatorVariable (we do that all the time), but
> that that OperatorVariable is a FaceVariable, through b.
>
> There isn't really a FiPyish way to tell what cell a particular position lies
> in,
> although I think this is close:
>
> >>> x = mesh.getCellCenters()[0]
> >>> S = j / dr * (abs(rbirth - x) < dr)
Thank you very much for the explanation and solution. I now ran into another
problem. I have the following lines of code
>>> rcrit = 0.0
>>> f = 4.0 * pi / 3.0 * ((phi * r**3 * dr) * (r > rcrit)).sum() # integrate
>>> from rcrit to infinity
>>> C = (C0 - f * Cp) / (1.0 - f)
>>> rcrit = alpha / log(C / Ceq)
where C0, Cp, Ceq and alpha are scalar constants and phi is the solution
variable. The problem with this is that rcrit on the second line will still
refer to zero and not the new value calculated on the last line. Is it possible
to handle this kind of circular reference in fipy?
Inspired by your answer, I have tried the following
class MyVariable(fp.Variable):
def __init__(self, value=0.0, unit=None, array=None, name='', cached=0,
hookVariable=None):
# turn off cached by default, otherwise _calcValue() will not be called.
fp.Variable.__init__(self, value, unit, array, name, cached)
self._hookVariable = hookVariable
def _calcValue(self):
#print '*** _calcValue() called'
if self._hookVariable is not None:
return self._hookVariable
else:
return fp.Variable._calcValue(self)
>>> rcrit = MyVariable()
>>> f = 4.0 * pi / 3.0 * ((phi * r**3 * dr) * (r > rcrit)).sum() # integrate
>>> from rcrit to infinity
>>> C = (C0 - f * Cp) / (1.0 - f)
>>> rcrit._hookVariable = alpha / log(C / Ceq)
but that results in a maximum recursion depth RuntimeError. Do you have a
suggestion on how to solve this?
Regards
/Jesper
_______________________________________________
fipy mailing list
[email protected]
http://www.ctcms.nist.gov/fipy
[ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
