Hi,
That particular L^2 Laplacian-like term is definitely correct. The L
variable is actually precisely a radius in many ways (it's the
"McIlwain L-shell number"), but it's also an inverse of a canonical
coordinate that the base equations are expressed in, and the operator
is essentially a diffusion term in that coordinate, called Phi. (Ie,
Phi ~ 1/L.) I'm still wrapping my head around why the references write
their equations the way they do; my best explanation right now is that
many of the other relevant operations are much more easily written in
terms of d/dL rather than d/dPhi.
Another weird attribute of the problem is that while the L dimension
has a clear spatial interpretation, the other 2 coordinates in the
problem definitely do not, so a spherical mesh would not be
appropriate.
Anyway, I have to say that for now I've gone and implemented my own
very simple-minded finite difference PDE solver for this problem, since
I had trouble with the 3D mesh (
https://github.com/usnistgov/fipy/issues/470) and all of my operators
are, in the end, pretty easy to implement. In the course of working on
this problem I've written a self-reference document that works out the
relevant equations pretty carefully, and I'm happy to share it if
anyone is curious about the aspects of this problem that may not make
it a great fit for fipy. (Though I'm definitely not familiar with
numerical PDE solvers in general so I'm sure it's mostly my lack of
facility with the tools.)
Cheers,
Peter
On Thu, 2015-10-22 at 15:39 +0000, Guyer, Jonathan E. Dr. wrote:
> Peter -
>
> I'm glad sweeping was able to get you a more reasonable solution
> (although I'm not seeing the non-linearity that would require it).
> `AdvectionTerm` only ever gets used in our level set examples, which
> I'm not intimately familiar with. I don't know of any fundamental
> reason it wouldn't work, but it's not something we've broadly tested.
>
> As far as how you've rearranged terms, I think you've taken a
> reasonable approach. You could alternatively do something like (2 *
> L**2 * f.grad.dot([[1]])) (not tested), but that's probably not any
> different and possibly less implicit than the AdvectionTerm.
>
> Generally, though, when terms arise that are not in one of FiPy's
> canonical forms, it's good to ask why. `L^2 d/dL (L df/dL)` and, even
> more so, `L^2 d/dL (D/L^2 df/dL)` in Eq (1) of the paper you cite in
> the code is strongly reminiscent of a laplacian operator in spherical
> coordinates: `(1/r^2) d/dr (r^2 df/dr)`, but the r^2 and (1/r^2) are
> swapped. Far be it from me to question an eminence like Van Allen,
> but are you certain that equation is right? Is L some sort of scaling
> that is inverse to radius by any chance?
>
> We provide 1D and 2D CylindricalGrid meshes and it would not be hard
> to write a SphericalGrid1D class, if that's called for. The reason to
> do that is that our CylindricalGrid1D class accounts for the
> curvature of the faces and the expanding volume of the cells with
> increasing radius. The laplacian operator automatically takes the
> form `(1/r) d/dr (r df/dr)` when you do this.
>
> Even if the operator in question is not in spherical coordinates, I
> suspect the correct answer is to build `L^2 d/dL (1/L^2 ...)` into
> the mesh somehow. Once we understand why it looks like that, it
> shouldn't be too hard to create the appropriate mesh class.
>
> - Jon
>
> P.S. Thanks very much for your pull request! I've merged it to the
> develop branch.
>
>
> On Oct 19, 2015, at 7:05 PM, Peter Williams <
> [email protected]> wrote:
>
> > Hi,
> >
> > I'm trying to use fipy to solve some differential equations, and
> > I'm not sure how to work with the particular equation I've got at
> > hand. I'm a total newbie at this, so I'm probably doing something
> > silly, but I can't figure out how to proceed so I thought I'd ask
> > for some help.
> >
> > The attached Python file demonstrates my problem and what I've done
> > so far. The key issue is that I have an equation with diffusion
> > terms that have spatially-dependent prefactors; in the simplest
> > case, something like
> >
> > > > L^2 d/dL (L df/dL)
> >
> > where L is a spatial variable and I'm trying to determine f(L). As
> > far as I can tell, this kind of term is not directly allowed by
> > fipy, and I get a TermMultiplyError if I try to construct it.
> >
> > In the simplest case, demonstrated in the attachment, I can divide
> > out the prefactor and fipy works great. However, in my full
> > problem, there will be multiple diffusion terms (in multiple
> > dimensions) with differing prefactors, so I won't be able to make
> > them all disappear.
> >
> > I'm pretty sure that I can generically re-express any prefactored
> > diffusion term as the sum of a prefactor-free diffusion term and an
> > advection term. However, I've tried to implement this with fipy,
> > and it doesn't work. The result from the solver is as if I didn't
> > include the advection term at all. Again this is demonstrated in
> > the attached file.
> >
> > So, I have two questions:
> >
> > (1) I'm probably doing something wrong with the advection term in
> > the attached sample. Can anyone spot the problem?
> >
> > (2) It'd be nice not to have to do the error-prone elimination of
> > my diffusion prefactors. Are there any ways to avoid this?
> >
> > Thanks for any tips,
> >
> > Peter
> > <testcase.py>_______________________________________________
> > 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
> ]
_______________________________________________
fipy mailing list
[email protected]
http://www.ctcms.nist.gov/fipy
[ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
