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 ]
