On Thu, Jun 14, 2012 at 10:55 AM, Daniel Wheeler <[email protected]>wrote:
> 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. > By scaling, I meant that the solutions themselves get scaled - not just the axis range. That is, increasing or decreasing Lx has the effect of a dilation in the y-direction. I checked with the code that you provided and Psi seems to be experiencing the same thing (ie: changes depending on the length of the mesh). > > > 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. > Spot on. There are exactly three conditions: that Phi goes to zero at -/+ infinity and that Phi=1 at x=0. I've taken on-board your explanation about why we still need to specify four. And this is where it gets interesting: In the code I sent you, I had broken up the triple derivative by placing a single inside a double (ie: convection term inside a diffusion term). In the code that you sent, it was the other way around (ie: a diffusion term inside a convection term). I would have thought that they would be identical. But strangely enough, they give different results - neither of which match the correct output from Mathematica. I think it might boil down to the boundary conditions again. In my code, it isn't obeying the condition at x=0 (ie that it must take the value of 1). In fact, changing this particular constraint doesn't change the output at all, suggesting that it's being treated as an outflow condition that you mentioned above. As for how to fix it and solve the PDE, I'm out of ideas. > > 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. > Yes I had swapped around Psi and Phi in the equation that I'd written. My bad. In addition, the coefficient in front of the third (last) term should have been only x, not x^2. I made the amendment to the code that you provided and it made no qualitative impact. > > > 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. > Makes sense - thanks. > > -- > 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 ]
