Peter - Interesting. I'd be happy to work with you to figure out this discretization in FiPy, either as a mesh or as a new term, if you find you'd like the benefits of FiPy's implicit solvers, parallelization, visualization, etc.
Thanks very much for your bug report on `UniformGrid3D._cellToCellIDs`. What we have is definitely wrong. Your fixes look right to me. It would be great if you submitted a pull request so that you get credit. We have tests of `_cellToCellIDs`, `_faceToCellDistanceRatio`, and `_cellToCellDistances` in that module, so it would be much appreciated if you could provide replacement tests that fail with the current code. - Jon On Oct 22, 2015, at 9:48 PM, Peter Williams <[email protected]> wrote: > 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 ] _______________________________________________ fipy mailing list [email protected] http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
