The following script runs without exploding. The results are not
interesting as the field seems to equilibrate to .5. If the variance
is increased the solution diverges. Do we know what the stable
parameter choices are?
#!/usr/bin/env python
from fipy import *
nx = ny = 100
dx = dy = 0.5
mesh = Grid2D(nx=nx, ny=ny, dx=dx, dy=dy)
var = CellVariable(name="variable", mesh=mesh)
var.setValue(GaussianNoiseVariable(mesh=mesh,
mean=0.5,
variance=0.01))
if __name__ == "__main__":
viewer = Viewer(vars=(var,), datamin=0., datamax=1.)
b = 0.5
h = 5.0
steps = 1000
eq = TransientTerm() == -DiffusionTerm((var**3, 1)) +
DiffusionTerm(-exp(-var) * var**3 + 3 * b / var)
dexp = -1
for step in range(steps):
dt = min(100, exp(dexp))
dexp += 0.01
eq.solve(var, dt = dt, solver=LinearLUSolver())
print
print 'max(var)',max(var),'min(var)',min(var)
print
if __name__ == '__main__':
viewer.plot(filename="output_%d.png" %step)
print 'step',step,'dt',dt
On Sun, Feb 8, 2009 at 4:17 PM, Ionut Vancea <[email protected]> wrote:
>
> Hello,
>
> On Mon, Feb 2, 2009 at 5:48 PM, <[email protected]> wrote:
>
>> OK, this is the problem. fDerivative is a FaceVariable. To get the proper
>> order for the discretization, you want the second derivative of f(\phi) to
>> be a FaceVariable, so fDerivative needs to be a CellVariable. So, either
>>
>> fDerivative = exp(-var) - (b / var**3)
>> (h**3 * fDerivative.getFaceGrad()).getDivergence()
>>
>> or, better,
>>
>> fSecondDerivative = -exp(-faceVar) + (3 * b / faceVar**4)
>> DiffusionTerm(coeff=fSecondDerivative)
>>
>>
>
>
> Thank you for your help, but I still have a "small" problem with my
> script. Here is the complete script:
>
> ======
> #!/usr/bin/env python
>
> from fipy import *
>
> nx = ny = 300
> mesh = Grid2D(nx=nx, ny=ny, dx=0.25, dy=0.25)
>
> var = CellVariable(name="variable", mesh=mesh)
> var.setValue(GaussianNoiseVariable(mesh=mesh,
> mean=0.5,
> variance=0.01))
>
> #faceVar = var.getArithmeticFaceValue()
>
> if __name__ == "__main__":
> viewer = Viewer(vars=(var,), datamin=0., datamax=1.)
>
> b = 0.5
> h = 5.0
> steps = 1000
>
> fDerivative = exp(-var)-(b / var**3)
> #fSecondDerivative = -exp(-faceVar) + (3 * b / faceVar**4)
>
> diffTerm4 = - DiffusionTerm((h**3, 1))
> eq = TransientTerm() == -DiffusionTerm((h**3, 1)) +
> (fDerivative.getFaceGrad() * h**3).getDivergence()
> #eq = TransientTerm() == -DiffusionTerm((h**3, 1)) +
> DiffusionTerm(coeff=fSecondDerivative)
>
> dexp = -5
> for step in range(steps):
> dt = min(100, exp(dexp))
> dexp += 0.01
> eq.solve(var, dt = dt)
>
> if __name__ == '__main__':
> viewer.plot(filename="output_%d.png" %step)
> print 'step',step,'dt',dt
>
> ======
> but when I am trying to run it, I obtain the following: (after 2 steps
> I can see something useful but nothing after step 3, 4, ...)
> ..
> /fipy/solvers/pysparse/linearPCGSolver.py:70: DeprecationWarning:
> PyArray_As1D: use PyArray_AsCArray.
> info, iter, relres = itsolvers.pcg(A, b, x, self.tolerance,
> self.iterations, Assor)
> /fipy/solvers/pysparse/linearPCGSolver.py:70: DeprecationWarning:
> PyArray_FromDimsAndDataAndDescr: use PyArray_NewFromDescr.
> info, iter, relres = itsolvers.pcg(A, b, x, self.tolerance,
> self.iterations, Assor)
> step 0 dt 0.00673794699909
> step 1 dt 0.00680566449223
> thinfilm3.py:33: StagnatedSolverWarning: The solver stagnated.
> Iterations: 1. Relative error: nan
> eq.solve(var, dt = dt)
> step 2 dt 0.0068740625575
> step 3 dt 0.00694314803475
> step 4 dt 0.00701292783259
> .......
>
> I saw in mail archives that some people obtained the same error and I
> tried to decrese/increase the time step, or to consider
> bigger/smaller lattice, but the error still persist.
> Maybe I did something wrong, but I can no spot the mistake.
>
> Thank you again,
> Ionut
>
>
--
Daniel Wheeler
