On Wed, Jun 13, 2012 at 11:36 PM, Allawala, Altan
<[email protected]> wrote:
> Hi Daniel,
>
> Thanks for the tip. The coupling functionality is working very nicely.
>
> I'm trying to solve the following differential equation under the condition
> that Psi=1 at x=0:
> 0 = \frac{\partial\phi}{\partial x} - \frac{\partial^{3}\phi}{\partial
> x^{3}} + x^{2}\phi,
>
> The results are qualitatively matching separate output from Mathematica.
> However, the graph's vertical scale is all off. In fact, it scales depending
> on the value of Lx.
That's weird. You seem to have set datamin and datamax, which fixes
the scaling. Is the scaling changing because the value of the solution
is changing. It's designed to rescale based on the solution if the
datamin and datamax arguments are not set.
> In addition, it isn't obeying the boundary on the left
> side of the mesh.
Am I right in saying that this equation can only needs three boundary
conditions to be well posed? That seems to make sense. The diffusion
equation needs a left and right boundary condition, but the convection
part only needs a left or right depending on the direction of the
velocity (-1 in this case) blowing from right to left. Do you have an
analytical result. That should only have three boundary conditions.
Due to backwards compatibility in FiPy, you need to set a constraint
on the outflow condition in order for the velocity not to be zero on
the outgoing face. The value of this boundary condition doesn't matter
since the flow is outgoing.
You have one of the phi boundary conditions commented out below. That
means you have no flux for the diffuision term
> I can't see why this is happening. Any idea what's going on? Here is the
> code: (I used yours as a template)
>
> from fipy import *
>
> xmax = 10.
> nx = 1000.
> mesh = Grid1D(nx=nx, Lx=xmax)
>
> x = mesh.getCellCenters()[0]
>
> psi = CellVariable(name=r'$\psi$', hasOld=True, mesh=mesh)#, value=0.)
> phi = CellVariable(name=r'$\phi$', hasOld=True, mesh=mesh)
>
> psi.constrain(1, mesh.facesLeft)
> psi.constrain(0, mesh.facesRight)
> #phi.constrain(0, mesh.facesLeft)
> phi.constrain(0, mesh.facesRight)
Keep all four boundary conditions for the reasons above.
> eqn1 = 0 == UpwindConvectionTerm(coeff = [[1]], var=psi) -
> DiffusionTerm(coeff = 1, var=phi) + ImplicitSourceTerm(coeff = x, var=psi)
> eqn2 = ImplicitSourceTerm(coeff = 1, var=phi) == UpwindConvectionTerm(coeff
> = [[1]], var = psi)
Have you swapped the meaning of phi and psi? What's written above
doesn't match the equation as far as I can tell. See the script below.
This should work better. Not sure if it is right, but if seems to obey
the boundary conditions.
from fipy import *
xmax = 10.
nx = 1000.
mesh = Grid1D(nx=nx, Lx=xmax)
x = mesh.getCellCenters()[0]
psi = CellVariable(name=r'$\psi$', hasOld=False, mesh=mesh)#, value=0.)
phi = CellVariable(name=r'$\phi$', hasOld=False, mesh=mesh)
psi.constrain(0, mesh.facesLeft)
psi.constrain(0, mesh.facesRight)
phi.constrain(1, mesh.facesLeft)
phi.constrain(0, mesh.facesRight)
eqn1 = (0 == UpwindConvectionTerm(coeff = [[1]], var=phi) -
UpwindConvectionTerm(coeff = [[1]], var=psi) +
ImplicitSourceTerm(x**2, var=phi))
eqn2 = ImplicitSourceTerm(coeff=1, var=psi) ==
DiffusionTerm(coeff=[[1]], var=phi)
eqn = eqn1 & eqn2
vi = Viewer(vars=(psi,phi))#, datamin=-1.5, datamax=1.5)
for i in range(5):
res = eqn.sweep()
vi.plot()
print res
raw_input('stopped')
>
> On a separate note, one of the things I'm going to do next is define the
> mesh from -xmax to +xmax. I know how to do this using the following line:
> mesh = Grid1D(nx=nx, Lx=2*xmax) + [[-xmax]]
>
> However, I would still need to set the value of $\psi$ at x=0. But since x=0
> is now no longer at the boundary of the mesh, it appears that one would
> somehow need to set a value of Psi right in the middle of the mesh - a
> "non-boundary condition". Can FiPy do this?
You can do this with a source term. You need to map the boundary
conditions to a volume integral using a Dirac-delta function and
smooth the Dirac-delta function. It can just be be smoothed over 1
cell so it is effectively an internal boundary condition.
--
Daniel Wheeler
_______________________________________________
fipy mailing list
[email protected]
http://www.ctcms.nist.gov/fipy
[ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
