Hi Joakim -
On Mar 23, 2014, at 7:16 PM, Joakim Odqvist <odqv...@kth.se> wrote:
> All these samples produce different results (where code1.py is the correct).
> In order to get the coupling with the external thermodynamics package we need
> to use an approach that is shown in code3.py.
Are you sure code1.py is correct? The `[i]` in
DiffusionTerm(coeff=D*(13077.922*(1/(PHI[i]-PHI[i]**2))-96000.0))
mean that the diffusion coefficient depends only on the value of PHI at the
i'th position.
If I correct this to
DiffusionTerm(coeff=D*(13077.922*(1/(PHI-PHI**2))-96000.0))
then code1.py and code2.py give the same results.
code3.py gives a different answer because function() is only evaluated when eq
is first declared and the coefficient is unchanging after that.
> How can we define the problem in such a way so we can obtain the same results
> in code1.py and code3.py?
The attached script is how I'd do it by declaring a specialized FaceVariable
subclass, with the caveat that it's exceptionally slow because of the Python
loop, so that's *not* how I'd do it.
Ideally, you want to query TC/TQ (or whatever) for the full array of values at
one time and let the loop happen in C. If you can't do that, then you'll want
to look into declaring _calcValue using Cython to do the loop for you in C.
from fipy import *
L=10.0e-9
nx=400
dx=L/nx
mesh = Grid1D(nx=nx, dx=dx)
phi = CellVariable(name=r"$\phi$", mesh=mesh)
class CalculatedDiffusionCoefficient(FaceVariable):
def __init__(self, D, PHI):
FaceVariable.__init__(self, mesh=PHI.mesh)
self.PHI = self._requires(PHI)
self.D = self._requires(D)
def _calcValue(self):
arr = numerix.empty(self.PHI.shape)
for i in xrange(self.PHI.mesh.numberOfFaces):
arr[i]=self.D*(13077.922*(1/(self.PHI[i]-self.PHI[i]**2))-96000.0)
return arr
noice = GaussianNoiseVariable(mesh=mesh,mean=0.5,variance=0.001)
pi=3.14159265359
for i in range(nx):
z=i*dx/L*2*pi*30
y=0.5+0.1*numerix.sin(z)
noice[i] = y
phi.setValue(noice)
PHI = phi.getArithmeticFaceValue()
D = 1e-23
viewer = Viewer(vars=(phi,), datamin=0., datamax=1.)
dexp = -5
elapsed = 0.
dt=1.0e-8
x_new=[]
x_old=[]
eq = TransientTerm() == DiffusionTerm(coeff=CalculatedDiffusionCoefficient(D,PHI))- DiffusionTerm(coeff=(D, 7.8e-16))
while elapsed < 1.2:
dt = min(100,numerix.exp(dexp))
x_old=PHI.all()
elapsed += dt
dexp += 0.01
eq.solve(phi, dt=dt)
x_new = PHI.all()
if ((x_new-x_old)/x_old)>0.01:
dt=dt/2
viewer.plot()
print 'dt: ',dt,elapsed
raw_input('Done!')
_______________________________________________
fipy mailing list
fipy@nist.gov
http://www.ctcms.nist.gov/fipy
[ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
