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
