Hi Daniel, Eureka. That certainly fixed it. In fact, I applied your changes to the other approach too (ie: my original code where I expressed the triple derivative by placing a convection term inside a diffusion term) and it's now giving the same correct answer as well. It is interesting that changing the Upwind to the CentralDifference convection term didn't make a difference in this case (probably because the second equation in this setup had only a convection term, no diffusion) but the actual change occurred when the left constraint was abolished in favor of a dirac-delta function approach. The FiPy works in mysterious ways!
Thank you enormously for looking over this, Daniel. One of the great things that I've noticed about using FiPy is the plentiful examples, documentation and support available from people such as yourself. Cheers, Altan On Fri, Jun 22, 2012 at 11:26 AM, Daniel Wheeler <[email protected]>wrote: > Altan, The future is here. Try this script, it seems quite close to the > mathematica solution. Not sure why. I removed the left hand side boundary > conditions and shifted the mesh alignment and used a source term to fix phi > in the left. I also used central difference rather than upwind. Upwind is > quite strange for coupled equations. upwind is okay for the first > convection term, but power law is required for the second convection term. > > Cheers > > > from fipy import * > > xmax = 10. > nx = 1000. > dx = xmax / nx > mesh = Grid1D(nx=nx + 1, Lx=xmax) - [[dx / 2.]] > > > 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.facesRight) > phi.constrain(0, mesh.facesRight) > > mask = CellVariable(mesh=mesh, value=0.) > mask[0] = 1e+10 > > eqn1 = (0 == CentralDifferenceConvectionTerm(coeff = [[1]], var=phi) - > CentralDifferenceConvectionTerm(coeff = [[1]], var=psi) + > ImplicitSourceTerm(x, var=phi)) + ImplicitSourceTerm(mask, var=phi) - mask > 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 Mon, Jun 18, 2012 at 11:27 AM, Allawala, Altan < > [email protected]> wrote: > >> Thanks for your response, Daniel. >> >> On Mon, Jun 18, 2012 at 10:18 AM, Daniel Wheeler < >> [email protected]> wrote: >> >>> On Fri, Jun 15, 2012 at 12:26 AM, Allawala, Altan < >>> [email protected]> wrote: >>> >>>> >>>> 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 >>>> >>> >>> Not really following you here. >>> >> >> In my code, the DEs are: >> >> 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)) >> >> >> In your code, the DEs are: >> >> >> eqn1 = (0 == UpwindConvectionTerm(coeff = [[1]], var=phi) >> - UpwindConvectionTerm(coeff = [[1]], var=psi) >> + ImplicitSourceTerm(coeff = x, var=phi)) >> >> eqn2 = ImplicitSourceTerm(coeff=1, var=psi) == DiffusionTerm(coeff=[[1]], >> var=phi) >> >> (Our phi's and psi's are switched of course.) >> >> Both your approach and mine should be completely equivalent ways of >> splitting up the triple derivative into convection and derivative terms. >> But they're giving different answers. >> >> For example, one could express the differential equation d(phi)^3/dx^3 = >> 1 in FiPy either as: >> >> eqn1 = (DiffusionTerm(coeff = 1, var=psi) == 1.) >> eqn2 = (UpwindConvectionTerm(coeff = [[1]], var = phi) == >> ImplicitSourceTerm(coeff = 1, var=psi)) >> >> or else as: >> >> eqn1 = (UpwindConvectionTerm(coeff = [[1]], var = psi) == 1.) >> eqn2 = (DiffusionTerm(coeff = 1, var=phi) == ImplicitSourceTerm(coeff = >> 1, var=psi)) >> >> In the first case, it amounts to placing the convection term inside of >> the diffusion term whereas in the second case it is the diffusion term >> being placed inside the convection term. Both are identical once the >> different boundary conditions have been accounted for. The same principle >> applies for the difference between our codes. >> >> >>> >>>> would have thought that they would be identical. But strangely enough, >>>> they give different results - neither of which match the correct output >>>> from Mathematica. >>>> >>>> >>> In the code I sent you, the three boundary conditions were satisfied. Is >>> the solution still wrong even though the boundary conditions are satisfied? >>> >> >> Yes, the solution is still wrong. I've attached the correct solution. >> >>> >>> >>>> >>>>> >>>> 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. >>>> >>> >>> I see a slight difference. The curve for psi does not have multiple >>> changes in gradient. How does the fipy solution differ from the mathematica >>> solution? The residual is tiny and the boundary conditions seem to be >>> satisfied at least for the code that I am running. >>> >> >> The Mathematica output is attached. (FiPy's output seems to perhaps >> resemble a translated form of this.) >> >>> >> >> -- >>> 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 ] >> >> > > > -- > 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 ]
