Hi Daniel, Thank you very much for your valuable advices. I think that procedure which you suggested in your first mail is a viable approach. I already tried multiplying complete equation with gradient \phi but this didn't help a lot (this because when grad \phi is zero, \phi_x \phi_t is also zero and equation can't be satisfied for \beta not equal 0, which is the case when driving force is not zero).
from fipy import *
nx = 400
Lx = 25 * nx / 100.
dx = Lx / nx
mesh = Grid1D(dx=dx, nx=nx)
timeStepDuration = 1e-6
phaseTransientCoeff = 1
thetaSmallValue = 1e-6
desiredResidual=1e-4
beta = 0.001
alpha = 1
mu = 1e3
# we take a equal 1 so it doesn't influence solution
a=1
gamma = 1000
x0=40
#phaseOld = CellVariable(name='phase field', mesh=mesh, hasOld=1)
phase = CellVariable(name='phase field',mesh=mesh,value=1.,hasOld=1)
phaseOld = CellVariable(name='phase field old',mesh=mesh, value=1, hasOld=1)
phaseEx = CellVariable(name='exact solution', mesh=mesh, value=1, hasOld=1)
x = mesh.getCellCenters()[0]
#phase.setValue(-tanh(x),where=mesh.getCellCenters()[0] > 0 )
bcs=(FixedValue(mesh.getFacesLeft(),1),
FixedValue(mesh.getFacesRight(),-1))
#phase.setValue(-1,where=mesh.getCellCenters()[0] > Lx /2.)
phase.setValue(tanh(-(x-x0)*alpha/4),where=mesh.getCellCenters()[0] > 0 )
# set exact solution
phaseEx.setValue(tanh(-(x-x0)*alpha/2),where=mesh.getCellCenters()[0] > 0)
#phase.setValue(sin(-(x-40)/4),where=mesh.getCellCenters()[0] > 0 )
print phase
convectionCoeff=beta-alpha/2*phase.getOld().getFaceGrad()
phaseEq = TransientTerm(phase.getOld().getGrad()[0]) ==
ImplicitDiffusionTerm()+VanLeerConvectionTerm(convectionCoeff) \
+ImplicitSourceTerm((phaseOld.getGrad()[0]-phase.getOld().getGrad()[0])/timeStepDuration)
#print diffusionCoeff
if __name__ == '__main__':
phaseViewer = Viewer(vars=(phase,phaseEx), datamin=-2., datamax=2)
phaseViewer.plot()
steps = 100
for i in range(steps):
phaseOld = phase.copy()
#print "Old old time step"
#print phaseOld.getOld()
phase.updateOld()
#print "After update old"
#print phase.getOld()
while residual > desiredResidual:
residual=phaseEq.sweep(phase,dt=timeStepDuration,solver=LinearLUSolver()
,boundaryConditions=bcs)
if __name__ =='__main__':
phaseViewer.plot()
I don't know if this implementation is correct, I have tried to approximate
\phi_{xt}^* \phi term
as ImplicitSourceTerm with coefficient \phi_{xt}^* =
phaseOld.getGrad()[0]-phase.getOld().getGrad()[0]/ timeStepDuration
which represents time derivative of the \phi_x. The solution is still unstable
as you can see from picture (alpha=1 beta=0.1).
About the question which I raised in my previous email, about handling 1/grad
Phi as a very high diffusion coefficient:
In their paper Kobayashi and Giga (Journal of Statistical Physics, Vol. 95,
Nos. 5/6, 1999) said quote:
"If u_x never vanishes or vanishes only at the exceptional points,
it might be no problem to consider these equations. However, as it will be
shown in the following sections, the situation is quite the opposite. In fact,
the region where u_x vanishes will increase and finally cover the whole
region in the homogeneous case, while in the inhomogeneous case, we are
particularly interested in the solutions whose spatial derivative vanishes
almost everywhere."
This is exactly the case for the angular equation in the phase field approach,
as
function is always more or less step like or "staircase" like. But in the
case of the previous equation there a regions were spatial derivative vanishes
(solid and liquid bulk regions) and
also the regions in which is smoothly varies like tanh solution (solid-liquid
interface).
I am laymen to understand completely the paper by Kobayashi and Giga, but in my
opinion
there is no proof that changing "infinite" diffusive constant with tanh \gamma
p /p will work for "every" singular diffusion
equation (SDE). It would be helpful to have certain mathematical criteria which
would give a estimation for how big must be maximum diffusion constant \gamma,
which gives the correct solution for the SDE.
There are some papers (PHYSICAL REVIEW E 72, 046136 2005, Phys. Rev.
Lett. 64, 1361 1990, Phys. Rev. Lett. 73, 1311 1994; Phys. Rev. E54, 376
1996.) based
on the Exact Renormalization Group (ERG) which handle nonlinear diffusion
equation (Barenblatt equation).
Giving asymptotic behavior of the Barenblatt equation for a very large times.
In my opinion the SDE equations can
also profit from such approach. This is probably a difficult mathematical
problem:)
Thank you, sorry for long mail
Toni
-----Original Message-----
From: [email protected] [mailto:[email protected]] On Behalf Of Daniel Wheeler
Sent: Freitag, 19. März 2010 5:02
To: Multiple recipients of list
Subject: Re: How to solve infinite diffusion PDE with Fipy
If alpha and beta are zero, which they seem to be in the equation that
you wrote down. Then you just have the three terms, which is easy to
set up. The source term will just have phi.getGrad()[0] -
phi.getOld().getGrad()[0] as the ImplictSourceTerm coefficient. It
doesn't have to be implicit either, maybe just start with an explicit
source at first, you'll have to sweep to convergence at each time step
and think carefully about the boundary conditions, which I haven't
discussed.
On Fri, Mar 19, 2010 at 11:48 AM, Daniel Wheeler
<[email protected]> wrote:
> Hi Toni,
>
> I'm not sure about this problem. One thought is that you could try solving
>
> \left( \phi_x^* \phi \right)_t = \phi_{xx} + \left( \left[ \beta -
> \frac{\alpha}{2} \phi^*\right] \phi \right)_x + \phi_{xt}^* \phi
>
> where the * values are old sweeps or old time steps. The non-starred
> values are implicit. Basically multiply through by \phi_x. Get an
> awkward \phi_{xt} term, but maybe that will be okay. The value
> \phi_x^* will be part of the transient coefficient. I assume you
> understand that it is best to get coefficients inside the derivatives
> for FV schemes in general.
>
> Cheers
>
> On Thu, Mar 18, 2010 at 12:58 PM, Ivas Toni <[email protected]> wrote:
>> Hi,
>>
>>
>>
>> I would like to solve "infinite" diffusion PDE (3d):
>>
>> 1/(a\nabla\phi) a \otimes a \frac{\partial^2}{\partial x_i \partial x_j}
>> \phi -\alpha \phi + \beta = \partial \phi /\partial t
>>
>>
>>
>> where "a" is a vector and a \otimes a is outer product of this vector
>> (symmetric matrix a_ij). This equation reduces to:
>>
>> \frac{\partial^2 \phi / \partial x^2 }{\partial \phi / \partial x} -\alpha
>> \phi + \beta = \partial \phi /\partial t
>>
>> for one dimensional case and a scalar.
>>
>>
>>
>> Equilibrium solution for 1d case is reducing to Burgers equation which can
>> be solved analytically(Hopf-Cole transformation). The solution is
>>
>> -tanh(2\alpha/a(x-\delta)) where \delta is position of the interface
>> (boundary condition are +1 for x=-\inf and -1 for x=+\inf) . This represents
>>
>> a traveling wave solution with certain velocity \delta=velocity * time.
>>
>>
>>
>> I tried to solve 1d equation with Fipy something like:
>>
>>
>>
>> from fipy import *
>>
>>
>>
>> nx = 40
>>
>> Lx = 2.5 * nx / 100.
>>
>> dx = Lx / nx
>>
>> mesh = Grid1D(dx=dx, nx=nx)
>>
>>
>>
>> timeStepDuration = 0.0002
>>
>> phaseTransientCoeff = 1
>>
>> thetaSmallValue = 1e-6
>>
>> beta = 1e5
>>
>> mu = 1e3
>>
>> gamma = 1000
>>
>> s = 1
>>
>>
>>
>> phase = CellVariable(name='phase field',mesh=mesh, value=1.)
>>
>> x = mesh.getCellCenters()[0]
>>
>> phase.setValue(x*(1-x),where=mesh.getCellCenters()[0] < Lx )
>>
>> #phase.setValue(1./3,where=mesh.getCellCenters()[0] > Lx /3.)
>>
>> #phase.setValue(1-x,where=mesh.getCellCenters()[0] > 2* Lx / 3.)
>>
>> print phase
>>
>> gradMag = phase.getFaceGrad().getMag()
>>
>> eps = 1./gamma/10.
>>
>> gradMag+=(gradMag < eps) * eps
>>
>> IGamma = (gradMag > 1./gamma) * (1/gradMag-gamma) + gamma
>>
>> #print IGamma
>>
>> diffusionCoeff = s*IGamma
>>
>>
>>
>> phaseEq = TransientTerm(phaseTransientCoeff) ==
>> ImplicitDiffusionTerm(diffusionCoeff)
>>
>> #print diffusionCoeff
>>
>> if __name__ == '__main__':
>>
>> phaseViewer = Viewer(vars=phase, datamin=0., datamax=2)
>>
>> phaseViewer.plot()
>>
>> steps = 100
>>
>> for i in range(steps):
>>
>> phase.updateOld();
>>
>> phaseEq.solve(phase,dt=timeStepDuration)
>>
>> if __name__ =='__main__':
>>
>> phaseViewer.plot()
>>
>>
>>
>> If I fix the boundary condition in fipy to 1 an -1 at the domain boundaries,
>> the numerical solution converges to step like function.
>>
>> It seems very difficult to control the terms involving the 1/grad(Phi). Is
>> it risky to handle this term as a very high diffusion coefficient ???
>>
>> The 1d solution of the previous equation has exact equilibrium solution as
>> tangent hyperbolic and not step function.
>>
>>
>>
>> I would be grateful if somebody can help me with this problem.
>>
>>
>>
>> Thank you,
>>
>> Toni Ivas
>>
>>
>>
>>
>
>
>
> --
> Daniel Wheeler
>
--
Daniel Wheeler
<<attachment: myeq1.png>>
