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
