stuart-knock
commented
6 years ago Below is an out of band response to PR, regarding task and time guesstimates for the spherical coordinate case, posted here for tracking purposes.
Implementing Spherical Case in NFTsim
Note: this is not regarding the simpler case of spherical boundary conditions applied to a regular grid (ie, it is regarding point 2. in the initial post of #128).
Development Tasks
Note: A prerequisite for this work is that an explicit solution for the spherical grid case, equivalent to the regular grid case in Chris Rennie's thesis or the appendix to the NFTsim paper, has been derived.
Pre-Spherical Coordinates
Rough time estimates:
- Best case: 4 days.
- Realistic case: 8-12 days
High level tasks (Ordered roughly from "required" to "may need to be addressed if we run into issues"):
- Review code thoroughly to identify full details of what needs to be done [1-3 days] -- a day was required to collate the rough details that are presented here.
- Address known clean-up/bug issues affecting related code [1-2 days] -- avoids running into these issues when doing main work (eg. #153, #164).
- While not required to implement the spherical case, completing the following tasks will facilitate validation of the new implementation:
- Arbitrary initial conditions [1-2 days], (see, #126 and #127); and
- Restart a simulation [1-? days], (see #24) -- The ? in the time estimate here is because it depends on details of object initialisation that would need to be understood and considered before starting work.
- Issues associated with current parsing of .conf files (#160, see also #19 and #22):
- required if it causes a problem when adding necessary new configuration parameters;
- may also be required if it interferes with other tasks above.
- Reconsider where and how spatial domain information is stored (#151) -- partly falls under initial review task mentioned above.
Spherical Coordinates:
Rough time estimates:
- Best case: 6 days.
- Realistic case: 10-?? days
High level tasks:
- Stencil [2-3 days] -- currently based on square/regular grid with unique boundary connections defining topology.
- a completely new stencil class will most likely be needed, to handle the pole (single point) boundary of spherical coordinates.
- Propagator [1-2 days]
- a completely new propagator derived class will need to be written, using the new Stencil and with a step method implementing the new explicit solution.
- Solver [0.5-1 day]
- Updated to make use of the new Propagator
- White noise normalisation [0.5-1 day] -- currently assumes regular, evenly-spaced, grid points, will need to be checked.
- Population [0.5-1 day] -- review, currently sheet length information is loaded-into/stored-on population objects.
- New configuration parameters [0.5-??? days].
- update parsing of input file to handle two cases (potential interaction with known parsing issue mentioned in "Pre-Spherical").
- General deltax, regular grid, to spherical case [??? days]:
- there are a range of other places in the code where assumptions are made about the grid regularity, these would need to be identified and addressed.
- This partly falls under the 2-3 days of thorough review mentioned in the "Pre-Spherical" section.
- Validating resulting implementation against known analytic solutions [1-3 days] -- assumes solutions already exist.
Post-Spherical Coordinates
- An awareness of the new spatial domain will need to be built into the existing Matlab interface.
- Update documentation.
- Add example .conf files.
- Add numerical test .conf files.
PR request: implement spherical topology.
Alternatives discussed with @stuart-knock
1) Spherical boundary conditions (i.,e., a plane folded into an envelope) are tricky to implement and very likely to introduce boundary artifacts.
2) Using a spherical coordinates transformation, where the wave equation is solved on a grid representing phi, theta-domain instead of x, y. Singularities at the poles (near the origin and at the edges of the rectangular ) still need to be handled. Also, this approach results in a variable discretization in terms of delta x, delta y. A circle close to the equator will be undersampled (delta x is large), while a circle close to the poles will be oversampled.
Courant condition?