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?
It seems that neither of the left hand side boundary conditions are being
obeyed. In addition, I've observed that adding a minus sign in front of
every term on both sides of the equality entirely changes the solution,
even though mathematically it shouldn't. Another strange observation is
that dividing the right hand side of the equation by x and then solving for
Phi gives the correct answer.
Any insight would be much appreciated.
Regards,
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.05
nx = xmax/dx
mesh = Grid1D(dx=dx, Lx=xmax)
x = mesh.getCellCenters()[0]
numerical = CellVariable(name="numerical", mesh=mesh)
numerical.constrain(1, mesh.facesLeft)
numerical.constrain(0, mesh.facesRight)
numerical.faceGrad.constrain([0], mesh.facesLeft)
analytical = CellVariable(name="analytical value", mesh=mesh)
analytical.setValue(exp(-x**2/2.))
mask = CellVariable(mesh=mesh, value=0.)
mask[0] = 1e+10
eqn = (TransientTerm(coeff=1, var=numerical) ==
- PowerLawConvectionTerm(coeff = x*[[1]], var=numerical)
- ImplicitSourceTerm(coeff = x*x, var=numerical))
# - ImplicitSourceTerm(mask, var=numerical) + mask)
numerical.setValue(exp(-x**2/2.))
vi = Viewer(vars=(numerical,analytical),datamin=0,datamax=2)
for i in range(10):
eqn.solve(var=numerical,dt=1.)
vi.plot()
_______________________________________________
fipy mailing list
[email protected]
http://www.ctcms.nist.gov/fipy
[ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
