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 ]
