If I replace the var**3 with h**3 the system seems stable, but uninteresting.
##eq = TransientTerm() == -DiffusionTerm((var**3, 1)) +
DiffusionTerm(-exp(-var) * var**3 + 3 * b / var)
eq = TransientTerm() == -DiffusionTerm((h**3, 1)) +
DiffusionTerm(-exp(-var) * h**3 + 3 * b * h**3 / var**4)
In the email you wrote above you have:
\partial_t \phi=-\nabla\{\phi^3\nabla[\triangle \phi
-\partial_\phi f(\phi)]\}
where f(\phi) = \frac{b}{2\phi^2}-e^{-\phi}
Why did you replace the phi^3 with a constant?
On Mon, Feb 9, 2009 at 1:00 PM, Daniel Wheeler
<[email protected]> wrote:
>
> 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
>
>
--
Daniel Wheeler
