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
<altan_allaw...@brown.edu> 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 <daniel.wheel...@gmail.com>
> wrote:
>>
>> On Wed, Sep 19, 2012 at 1:46 PM, Allawala, Altan
>> <altan_allaw...@brown.edu> 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
>> fipy@nist.gov
>> http://www.ctcms.nist.gov/fipy
>> [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
>
>
>
> _______________________________________________
> fipy mailing list
> fipy@nist.gov
> http://www.ctcms.nist.gov/fipy
> [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
>
--
Daniel Wheeler
_______________________________________________
fipy mailing list
fipy@nist.gov
http://www.ctcms.nist.gov/fipy
[ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
