Hello all,
I am happily discovering the power of Fipy, I think it is a tool
adaptable to many problems.
I am currently trying to code a system of 3 equations which are the
following :
[1] dn/dt + nabla . (nv) = 0
[2] dv/dt + v . nabla v = mu * nabla c
[3] dc/dt = D nabla ^2 c + alpha n - c/tau
where the variables are n, v and c, and they depend on time t, x and y
(2D problem). v is a vector field whereas n and c are scalar.
mu, D and tau are numerical parameters. nabla is the gradient operator.
[1] is a mass conserveration equation, [2] a Burger's one, and [3] a
reaction diffusion.
Here is my code (the interesting parts) :
________________________________________________________________________________________________
# (...) imports
#$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
# Preparation
#$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
nx=100;
ny=100;
m = PeriodicGrid2D(nx=nx, ny=ny);
x,y=m.cellCenters;
#==== Variables
===============================================================================
n = CellVariable(mesh=m)
c = CellVariable(mesh=m)
v = CellVariable(mesh=m,rank=1) #v is a vector field
#==== Initial Conditions
======================================================================
# (...) I set initial conditions with n randomly dispersed on the grid
#==== Equations
================================================================================
#parameters
D=1
tau=64*60.;
mu=1;
alpha=1;
#equations
eq1 = TransientTerm(var=n) == PowerLawConvectionTerm(coeff=v,var=n);
eq2 = TransientTerm(var=v) == -v*v.getGrad()+mu*c.getGrad();
eq3 = TransientTerm(var=c) == DiffusionTerm(coeff=D, var=c) +alpha * n -
c/tau;
#equation coupling
eqn = eq1 & eq2 & eq3
#$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
# Computing
#$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
time_step=1e-5
tmax=2000;
for t in range(tmax):
eqn.solve(dt=time_step);
________________________________________________________________________________________________
I am somewhat uncertain about this code. I have trouble running it and
think that the problem comes from the way I set the problem and
equations.
For starters, I am unsure on how to write the equations. I *think* this
is the right way to do so, but am not that sure especially for equation
[2]. I could not manage to get that "v*v.getGrad()" in a ConvectionTerm
form. When the code works, it gives results which seem to slowly
approach what I want (I know what shape n is supposed to have at the
end).
But code as I gave it to you does not work. I get the following error :
File
"/home/douarre/Downloads/FiPy-3.1/myvirtualenv/lib/python2.7/site-packages/FiPy-3.1-py2.7.egg/fipy/terms/term.py",
line 151, in _prepareLinearSystem
self._checkVar(var)
File
"/home/douarre/Downloads/FiPy-3.1/myvirtualenv/lib/python2.7/site-packages/FiPy-3.1-py2.7.egg/fipy/terms/abstractBinaryTerm.py",
line 108, in _checkVar
self.term._checkVar(var)
File
"/home/douarre/Downloads/FiPy-3.1/myvirtualenv/lib/python2.7/site-packages/FiPy-3.1-py2.7.egg/fipy/terms/abstractBinaryTerm.py",
line 108, in _checkVar
self.term._checkVar(var)
File
"/home/douarre/Downloads/FiPy-3.1/myvirtualenv/lib/python2.7/site-packages/FiPy-3.1-py2.7.egg/fipy/terms/abstractBinaryTerm.py",
line 108, in _checkVar
self.term._checkVar(var)
File
"/home/douarre/Downloads/FiPy-3.1/myvirtualenv/lib/python2.7/site-packages/FiPy-3.1-py2.7.egg/fipy/terms/unaryTerm.py",
line 141, in _checkVar
and (numerix.sctype2char(var.getsctype()) not in
numerix.typecodes['Float'])):
File
"/home/douarre/Downloads/FiPy-3.1/myvirtualenv/lib/python2.7/site-packages/FiPy-3.1-py2.7.egg/fipy/variables/coupledCellVariable.py",
line 157, in getsctype
self.typecode = numerix.obj2sctype(rep=self.numericValue,
default=default)
File
"/home/douarre/Downloads/FiPy-3.1/myvirtualenv/lib/python2.7/site-packages/FiPy-3.1-py2.7.egg/fipy/variables/coupledCellVariable.py",
line 113, in numericValue
return numerix.concatenate([var.numericValue for var in self.vars])
ValueError: all the input arrays must have same number of dimensions
This is due I think to the coupling of [2] with [1] and [3].
If I don't couple the equations and replace
eqn.solve(dt=time_step);
by
eq1.solve(dt=dtt);
eq2.solve(dt=dtt);
eq3.solve(dt=dtt);
or
eqn.solve(dt=dtt); #where eqn = eq1 & eq3
eq2.solve(dt=dtt);
(and add has_Old(), update_Old() at the right places)
then the code works. But after some time computing, most of the time I
get the error :
/home/douarre/Downloads/FiPy-3.1/myvirtualenv/lib/python2.7/site-packages/FiPy-3.1-py2.7.egg/fipy/solvers/scipy/linearLUSolver.py:66:
RuntimeWarning: overflow encountered in square
error0 = numerix.sqrt(numerix.sum((L * x - b)**2))
/home/douarre/Downloads/FiPy-3.1/myvirtualenv/lib/python2.7/site-packages/FiPy-3.1-py2.7.egg/fipy/solvers/scipy/linearLUSolver.py:71:
RuntimeWarning: overflow encountered in square
if (numerix.sqrt(numerix.sum(errorVector**2)) / error0) <=
self.tolerance:
/home/douarre/Downloads/FiPy-3.1/myvirtualenv/lib/python2.7/site-packages/FiPy-3.1-py2.7.egg/fipy/solvers/scipy/linearLUSolver.py:71:
RuntimeWarning: invalid value encountered in double_scalars
if (numerix.sqrt(numerix.sum(errorVector**2)) / error0) <=
self.tolerance:
/home/douarre/Downloads/FiPy-3.1/myvirtualenv/lib/python2.7/site-packages/FiPy-3.1-py2.7.egg/fipy/variables/variable.py:1151:
RuntimeWarning: overflow encountered in multiply
return self._BinaryOperatorVariable(lambda a,b: a*b, other)
/home/douarre/Downloads/FiPy-3.1/myvirtualenv/lib/python2.7/site-packages/FiPy-3.1-py2.7.egg/fipy/variables/variable.py:1151:
RuntimeWarning: invalid value encountered in multiply
return self._BinaryOperatorVariable(lambda a,b: a*b, other)
Traceback (most recent call last):
File "test_fipy.py", line 200, in <module>
eq2.solve(dt=dtt);
File
"/home/douarre/Downloads/FiPy-3.1/myvirtualenv/lib/python2.7/site-packages/FiPy-3.1-py2.7.egg/fipy/terms/term.py",
line 213, in solve
solver._solve()
File
"/home/douarre/Downloads/FiPy-3.1/myvirtualenv/lib/python2.7/site-packages/FiPy-3.1-py2.7.egg/fipy/solvers/scipy/scipySolver.py",
line 61, in _solve
self.var[:] = numerix.reshape(self._solve_(self.matrix,
self.var.ravel(), numerix.array(self.RHSvector)), self.var.shape)
File
"/home/douarre/Downloads/FiPy-3.1/myvirtualenv/lib/python2.7/site-packages/FiPy-3.1-py2.7.egg/fipy/solvers/scipy/linearLUSolver.py",
line 64, in _solve_
permc_spec=3)
File
"/home/douarre/Downloads/FiPy-3.1/myvirtualenv/lib/python2.7/site-packages/scipy/sparse/linalg/dsolve/linsolve.py",
line 257, in splu
ilu=False, options=_options)
RuntimeError: Factor is exactly singular
which by my little bit of knowledge on numerical schemes mean that the
scheme diverged at some point, probably because of a too big time step
(error happens almost all the time if time_step=10e-3).
I have tried using different solvers (pysparse and scipy).
I have also tried different types of Convection Terms but am not sure
which one to use.
I am sorry if I made beginner mistakes while using Fipy, coding
vectors/variables, etc. I went through the examples and Theoretical and
Numerical Background but am still having trouble.
Do you think there is maybe a better way to code this so that it runs
more smoothly and for bigger time steps ?
Thanks in advance,
Clément
_______________________________________________
fipy mailing list
[email protected]
http://www.ctcms.nist.gov/fipy
[ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
