No problem Daniel, I appreciate that you at least looked over it. Thanks. On Thu, Sep 27, 2012 at 11:05 AM, Daniel Wheeler <[email protected]>wrote:
> Altan, I've played with this a bit and I don't have enough experience > of these third order terms to really understand what's causing the > instability. Probably worth finding the literature on these types of > equations and looking at the schemes used. It obviously requires the > combination of diffusion and convection in the correct way. Sorry I > can't help more. > > On Tue, Sep 25, 2012 at 7:11 PM, Allawala, Altan > <[email protected]> wrote: > > Hi Jonathan and Daniel, > > > > Thank you both for your replies. Yes, I got myself muddled up with the > > coefficients. Let me take a step back. Say I want to find the time > evolution > > of the following 1D PDE: > > > > \frac{\partial (\psi/x) }{\partial t} = - \frac{\partial \psi }{\partial > x} > > + \frac{\partial^{^{3}} \psi }{\partial x^{^{3}}} - x \psi > > > > I do this by first splitting up this third order PDE into a second order > and > > first order PDE as follows: > > > > \frac{\partial (\psi/x) }{\partial t} = - \frac{\partial \psi }{\partial > x} > > + \frac{\partial^{^{2}} \phi }{\partial x^{^{3}}} - x \psi > > > > and > > > > \phi = \frac{\partial \psi }{\partial x} > > > > I recognize that I need to sweep the solution since my PDEs are coupled. > > (Ref: /FAQ.html#iterations-timesteps-and-sweeps-oh-my) > > > > So I have placed sweeps inside each time step. However the steady-state > > solution isn't matching the correct one from Mathematica. In fact it > isn't > > even stable. I thought the problem was that I was using an inappropriate > > convection term but I've tried all of them by trial and error. Can you > see > > where I went wrong? > > > > Below is my code. I've commented out my BCs because the initial setValue > > functions which I assign to psi and phi satisfy the boundary condition > and > > so the values at the boundaries shouldn't change as time elapses. (My BCs > > are that psi(x=0)=1 and d(psi(x=0))/dx=0.) > > > > > > Many thanks, > > > > Altan > > > > > > import math > > from numpy import * > > from pylab import * > > > > import matplotlib > > from matplotlib import pyplot > > > > from fipy import * > > from scipy import * > > > > xmax = 5. > > dx = 0.01 > > > > nx = xmax/dx > > mesh = Grid1D(dx=dx, Lx=xmax) > > > > phi = CellVariable(name="phi", mesh=mesh, hasOld=True) > > psi = CellVariable(name="psi", mesh=mesh, hasOld=True) > > > > #mask = CellVariable(mesh=mesh, value=0.) > > #mask[0] = 1e+10 > > > > x = mesh.getCellCenters()[0] > > > > eqn1 = (TransientTerm(coeff=1/x, var=psi) == > > - CentralDifferenceConvectionTerm(coeff = [[1]], var=psi) > > + DiffusionTerm(coeff = 1, var=phi) > > - ImplicitSourceTerm(coeff = x, var=psi)) > > # - ImplicitSourceTerm(mask, var=psi) + mask) > > eqn2 = (ImplicitSourceTerm(coeff = 1, var=phi) == > > CentralDifferenceConvectionTerm(coeff = [[1]], var=psi)) > > # - ImplicitSourceTerm(mask, var=phi)) > > > > eqn = eqn1 & eqn2 > > > > phi.setValue(-x*exp(-x**2/2.)) #ansatz satisfies BCs > > psi.setValue(exp(-x**2/2.)) #ansatz satisfies BCs > > > > vi = Viewer(vars=(psi,phi),datamin=-1,datamax=2) > > > > for i in range(25): > > # move forward in time by one time step > > psi.updateOld() > > phi.updateOld() > > > > # sweep multiple times each time step > > for j in range(3): > > res = eqn.sweep(dt=0.001) > > print res > > vi.plot() > > > > raw_input() > > > > > > On Thu, Sep 20, 2012 at 10:59 AM, Daniel Wheeler < > [email protected]> > > wrote: > >> > >> On Wed, Sep 19, 2012 at 1:46 PM, Allawala, Altan > >> <[email protected]> wrote: > >> > Hello, > >> > > >> > I'm trying to find the steady state solution to the following linear > >> > PDE: > >> > > >> > \frac{\partial \varphi}{\partial t} = - x * \frac{\partial > >> > \varphi}{\partial > >> > x} - x^2 \varphi > >> > > >> > Below is my code to solve this. However, the numerical solution is not > >> > matching the analytical one. Any idea where the problem is? > >> > >> As Jon pointed out the code doesn't match the equation. Make the > >> following changes and you get the desired answer. > >> > >> $ diff old.py new.py > >> 30,32c30,32 > >> < eqn = (TransientTerm(coeff=1, var=numerical) == > >> < - PowerLawConvectionTerm(coeff = x*[[1]], var=numerical) > >> < - ImplicitSourceTerm(coeff = x*x, var=numerical)) > >> --- > >> > eqn = (TransientTerm(coeff=1 / x, var=numerical) == > >> > - PowerLawConvectionTerm(coeff = [[1]], var=numerical) > >> > - ImplicitSourceTerm(coeff = x, var=numerical)) > >> > >> Cheers > >> > >> -- > >> Daniel Wheeler > >> _______________________________________________ > >> fipy mailing list > >> [email protected] > >> http://www.ctcms.nist.gov/fipy > >> [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ] > > > > > > > > _______________________________________________ > > fipy mailing list > > [email protected] > > http://www.ctcms.nist.gov/fipy > > [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ] > > > > > > -- > Daniel Wheeler > _______________________________________________ > fipy mailing list > [email protected] > http://www.ctcms.nist.gov/fipy > [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ] >
_______________________________________________ fipy mailing list [email protected] http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
