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 ]