boutproject / zoidberg

A Flux-Coordinate Independent (FCI) Grid Generator for BOUT++
GNU Lesser General Public License v3.0
2 stars 2 forks source link

Dommasck potentials bugfixes and improvements #11

Closed nnikulshin closed 2 years ago

nnikulshin commented 2 years ago

I have found and fixed some bugs in the Dommaschk potential implementation, namely: 1) The summation in the delta coefficients (see functions CD and CN and eq 33 in Dommaschk's paper) should go up to n, not n+1. Note that the Python range builtin does not include the endpoint, so range(1,n+1) gives 1,2,...,n; however the original implementation, which used SymPy Piecewise functions instead of Python lambdas (see point 4) had a SymPy summation that included n+1. 2) In the V_hat function, I rewrote the expression to avoid having m in the denominator, which would be a problem if m is zero. 3) In the U_hat function, I check if i is zero before dividing by it. Even though the piecewise function is only supposed to return 1/i when i > 0, passing 1/i to Piecewise as a potential value will throw an exception if i is zero. 4) Not a bugfix, but I noticed that replacing the SymPy Piecewise functions by lambdas in the CD and CN functions noticeably speeds up computation.

ZedThree commented 2 years ago

LGTM, but I'll let @bshanahan check the maths and approve it