MOLE.jl Software Development Plan

# Goal

Produce a native, well-tested, high-performance Julia package `MOLE.jl` that matches the functionality and numerical correctness of the MATLAB/Octave reference in `src/matlab_octave`, offers idiomatic Julia APIs (multiple dispatch, sparse arrays, types), integrates with Julia ecosystem tooling (Documenter, BenchmarkTools, Test, etc), and provides parity tests and examples.

We will use the MATLAB/Octave code in `src/matlab_octave` as the canonical reference to port from

# Architecture & directory layout

```
MOLE.jl/
├── Project.toml
├── src/
│ ├── MOLE.jl # top-level module (re-exports)
│ ├── Grid/ # grid types & helpers
│ │ ├── types.jl
│ │ ├── uniform.jl
│ │ └── curvilinear.jl
│ ├── Operators/ # core operators; each source file will contain the 1D, 2D and 3D implementation
│ │ ├── interpolation.jl
│ │ ├── gradient.jl
│ │ ├── divergence.jl
│ │ ├── curl.jl
│ │ └── laplacian.jl
│ ├── BCs/ # boundary condition helpers
│ ├── Utils/ # numerical utilities (staggering, indexing)
│ └── IO/ # load/save, sample meshes
├── test/
│ ├── runtests.jl
│ ├── Grid/
│ | ├── grid_tests.jl
| ├── Operators/
│ | ├── operator_tests.jl
│ ├── BCs/
│ | ├── bc_tests.jl
│ ├── Utils/
│ | ├── util_tests.jl
│ └── IO/
│ | ├── io_tests.jl
├── docs/ # Documenter.jl docs source
└── examples/

```

- Precision: default `Float64`, accept `T<:AbstractFloat`

# Interoperability and compatibility
- Consider compatibility with DifferentialEquations.jl time integrators.

# Performance considerations
- Julia has an interactive REPL, is Just-In-Time (JIT) compiled (triggered by first evaluation of function). It allows polymorphism via multiple dispatch (at compile or run time). We definitely want to take advantage of this.
- Can support CPUs and GPUs programming using a common open-source (single-source) code base
- We can write generic code and the compiler will specialize on types of calling arguments, e.g., `f(x::AbstractArray)` where `AbstractArray` can be `Array` of `Float32`, `Float64` or a `CuArray`.
- We want to use Julia broadcasting:
* apply a vectorized function point-wise to an array. Scalar values are "broadcast" over arrays; Fusion of multiple operations.
* Support user-extensible API: can be specialized for custom functions or argument types (e.g., `CuArray` compiles and applies a custom CUDA kernel).
* The current MATLAB/Octave implementation relies on sparse storage, but
- Do we want MOLE Operators (`grad`, `div`, `interpol`) as "pseudo-functions"? In that case, they need to act like functions when broadcasted over an array of values, but can’t be called on a single value; can be composed and fused w/ function calls. (Very much like performant methods such as matrix-free methods, i.e., no assembly required; they specify the action of the operator.)
- Use `@inbounds` for bounds-checking and `@view`/`@views` for slices

# Milestones and plan
Here is a tentative inventory for the porting

[MOLE_Inventory__MATLAB-_Julia__Preview.csv](https://github.com/user-attachments/files/25350240/MOLE_Inventory__MATLAB-_Julia__Preview.csv)

## Subtasks

- [x] #280
- [x] #299