Roadmap to v1.0¶

This chapter was written by Claude at the request of the author (see Claude’s contributions to Janus, task The road to v1.0). It analyses the evolutions that would turn Janus into a flexible, easily maintainable and extensible library, and proposes an ordered sequence of releases leading to version 1.0. It is a working document: none of the evolutions described below has been implemented yet.

Goals¶

Janus was initially designed for two purposes: rapid prototyping, and large (distributed) simulations. The second purpose is now obsolete, because new developments in the community will soon offer excellent alternatives. Janus should therefore focus on prototyping:

installation should be as easy as possible (ideally pip install, without compiler);
new discretizations and new physics (conductivity, Darcy flow, finite strain solid mechanics…) should be implementable in pure Python, without compilation, even if the resulting code is sub-optimal;
GPU execution and automatic differentiation are not required, but the architecture should not rule them out;
any dependency that is distributed as binary wheels is acceptable;
breaking changes are allowed.

The speed estimates given below are not backed by benchmarks (at the author’s request, the analysis relies on the code and on the known properties of the tools); benchmarks are part of the proposed milestones.

Where does the friction come from?¶

The following observations are based on a complete reading of the current code base (about 2300 lines of Cython, in janus/operators.pyx, janus/green.pyx, janus/material/elastic/linear/isotropic.pyx and janus/fft/).

Extensions must be written in Cython. All the methods that do the actual work (c_apply, c_set_frequency, c_apply_by_freq, c_to_memoryview) are cdef methods, which cannot be overridden in Python. The docstrings of init_sizes and init_shapes suggest that pure Python subclasses are supported, but this only holds for direct calls to apply. Inside the library, the Cython methods are called, and Python overrides are silently ignored. For example, a Python subclass of AbstractOperator which doubles its input returns the expected result when called directly, but when used as a local operator of a BlockDiagonalOperator2D, the output array is left untouched (the no-op AbstractOperator.c_apply is called instead); this was checked with the current build. Likewise, a continuum Green operator (for new physics) must be a Cython subclass of AbstractGreenOperator.
Operators are scalar and stateful. A continuum Green operator is evaluated one wave-vector at a time (set_frequency(k), then apply(tau, eta)), and discrete Green operators loop over all frequencies in Cython. This design is efficient in Cython, but it cannot be transposed to Python: a Python implementation would require one Python call per frequency (e.g. 8 million calls for a 256³ grid). Any Python-level extension mechanism requires operators that act on whole arrays of wave-vectors or cells at once.
Everything is duplicated for 2D and 3D. Almost every class exists in two versions (AbstractStructuredOperator2D/3D, BlockDiagonalOperator2D/3D, TruncatedGreenOperator2D/3D, FilteredGreenOperator2D/3D, FiniteDifferences2D/3D, _RealFFT2D/3D, _GreenOperatorForStrains2D/3D, FourthRankIsotropicTensor2D/3D…), and the 3D versions of some operators are missing (FourthRankCubicTensor3D).
Discretizations are tied to the physics. The filtered Green operators hard-code symmetric 3×3 (2D) and 6×6 (3D) matrices, i.e. strain-based elasticity; their code is generated by scripts/gencode.py. The filtered discretization cannot be reused for, say, conductivity (2×2 or 3×3 matrices).
The data model is restricted. Local data are vectors of fixed size, symmetric tensors are stored in Mandel–Voigt form, only float64 is supported, and complex Fourier coefficients are stored as interleaved real numbers, the Green operator being applied separately to the real and imaginary parts (which assumes a real Fourier symbol). Finite strain mechanics requires non-symmetric second-rank tensors (9 components in 3D), and some discretizations have complex symbols.
Distributed memory leaks into the serial code. Discrete Green operators and FFT objects distinguish local and global shapes (offset0, global_shape0, global_ishape…), although only the serial version is used in practice.
The build chain is fragile. Compilation requires a C compiler, a hand-written setup.cfg with machine-specific paths to FFTW, and a setup.py that relies on distutils (see Installation and the task Installation on a windows machine in Claude’s contributions to Janus).

Minor issues found along the way illustrate the maintenance burden of the current code: TruncatedGreenOperator3D.c_apply calls the Python-level set_frequency (with argument checks) in its innermost loop, whereas the 2D version calls c_set_frequency; some variables in FiniteDifferences3D.c_set_frequency are untyped (hence Python objects); the error messages of _RealFFT3D report wrong shapes (TODO in the code).

Analysis of the proposed evolutions¶

Each evolution is analysed in terms of gains, losses, consequences on the code and implementation problems, and ends with a recommendation. Evolutions E1 to E3 were proposed by the author; E4 to E12 are additional suggestions.

E1. Removal of the MPI dependency¶

Gains.

Installation: no mpi4py, no MPI-enabled FFTW, no parsing of mpicc -show in setup.py. This also removes a fragile step on Windows (setup.py fails if mpi4py is installed but mpicc is missing).
Code: the distinction between local and global shapes disappears from the FFT objects and the discrete Green operators (offset0, global_shape0, global_ishape, global_oshape, n0_loc). janus/fft/parallel/ (about 110 lines), tests/parallel/, sphinx/parallel_fft_tutorial.py and the [fftw_mpi] configuration can be deleted.
Documentation: the API of janus.fft.parallel is currently missing from the documentation built on Windows; this problem disappears.

Losses. Distributed-memory simulations are no longer possible, which is consistent with the new goals. Shared-memory parallelism remains available (multithreaded FFTs, Numba prange, JAX on GPU): 3D grids of 256³ to 512³ cells remain within reach of a workstation (a single 512³ field of symmetric tensors, 6 components in double precision, requires 6.4 GB).

Consequences on the code. Breaking change for users of janus.fft.parallel (presumably few). The constructors of the discrete Green operators are simplified.

Implementation problems. None: this is mostly code deletion. If distributed computing were needed again, it could be reintroduced at the level of the FFT backend (e.g. mpi4py-fft, or JAX sharding) without affecting the rest of the architecture, provided that the FFT is kept behind a small interface (see E2).

Recommendation. Do it first (milestone 0.2): it is cheap, and it reduces the amount of code to be ported by the other evolutions.

E2. Interfacing with FFTW¶

The hand-written wrapper (janus/fft/serial/, about 370 lines of Cython) could be replaced by pyFFTW. Another candidate should be considered: scipy.fft.

pyFFTW.

Gains: removes the Cython wrapper and, above all, the need to locate and link FFTW at build time (setup.cfg), which is the most fragile step of the current installation. pyFFTW is distributed as binary wheels (version 0.15.1 provides wheels for Python 3.14, including Windows), and exposes the planner flags, wisdom, multithreading and aligned arrays of FFTW.
Losses: an additional compiled dependency, with a slow release cadence (the latest release, 0.15.1, dates from October 2025). FFTW is licensed under the GPL, which constrains the distribution of software that links it; note that Janus (BSD-3-Clause) is already in this situation.

scipy.fft.

Gains: no additional dependency (SciPy is needed anyway for iterative solvers), BSD license, real-to-complex transforms along arbitrary axes of an array (rfftn(x, axes=...), hence directly applicable to fields of tensors), multithreading (workers argument), no planning step. Moreover, pyFFTW can be plugged in as a backend of scipy.fft (scipy.fft.set_backend(pyfftw.interfaces.scipy_fft)), without any change to the code that calls scipy.fft.
Losses: typically somewhat slower than FFTW with FFTW_MEASURE (to be measured); no wisdom.

Consequences on the code. The FFT objects (_RealFFT2D/3D, with their explicit copies to and from an internal buffer) and the planner flags (janus.fft.FFTW_*) disappear. The current wrapper copies the data to and from an FFTW buffer at each call anyway, so no zero-copy optimization is lost. Complex Fourier coefficients become genuine complex arrays, instead of interleaved real arrays.

Implementation problems. The layout of the output of real-to-complex transforms (last axis of length n//2 + 1) and the handling of the Nyquist frequency for even grid sizes must be reproduced exactly, so that the new implementation matches the reference data in tests/data.

Recommendation. Use scipy.fft by default (through the array namespace, see E7), and document pyFFTW as an optional accelerator, via the scipy.fft backend mechanism. Do not make pyFFTW a hard dependency.

E3. Replacing Cython with other (dynamic) tools¶

Two kinds of loops are compiled in the current code: the frequency-wise application of the Green operator (a small dense matrix per frequency), and the cell-wise application of local operators. Both are “embarrassingly vectorizable”. Three options are considered.

Pure NumPy (vectorized array expressions).

Gains: no compilation at all, runs everywhere, easiest to read, write and debug; new physics can be prototyped in a few lines (see E4).
Losses: array expressions create temporaries, and are limited by memory bandwidth. The order of magnitude of the cost should remain comparable to the current code, which itself is not optimal (memoryview slicing and virtual calls at each frequency), but this must be measured. Storing the Green operator for all frequencies is only affordable in 2D: in 3D, a 6×6 matrix per frequency of a 256³ grid requires 2.4 GB, so the operator must be applied on the fly (e.g. through a closed-form expression, or by blocks of frequencies).

Numba.

Gains: loops written in Python syntax and compiled at run time (@njit, parallel=True, cache=True), with a performance similar to Cython, and no compilation at installation (Numba 0.67 provides wheels for Python 3.14, including Windows).
Losses: only a subset of Python and NumPy is supported; polymorphism (user-defined physics passed to generic kernels) is awkward (jitted functions passed as arguments, experimental jitclass); compilation latency at the first call; harder debugging; no automatic differentiation; GPU support is limited to NVIDIA hardware, through a separate package.

JAX.

Gains: NumPy-like API (jax.numpy, including jax.numpy.fft); jit compiles and fuses array expressions (which removes the temporaries of the pure NumPy approach); vmap turns a function of one cell into a function of the whole grid; automatic differentiation (grad, jacfwd), which provides consistent tangent operators for nonlinear and finite strain constitutive laws; GPU execution; matrix-free solvers (jax.scipy.sparse.linalg.cg, gmres). jaxlib 0.11 provides wheels for Python 3.14, including Windows (CPU only; GPU support on Windows requires WSL2).
Losses: arrays are immutable, so that the current in-place API (apply(x, y), apply(x, x)) is not possible; single precision by default (double precision must be enabled explicitly); tracing constraints (static shapes, no Python control flow on array values inside jitted functions); compilation latency; a large dependency.

Consequences on the code. Whatever the tool, the per-frequency/per-cell design must be replaced by operations on whole arrays (E4). With Numba, the kernels would remain explicit loops, specific to each physics; with NumPy or JAX, the same array code serves both purposes.

Implementation problems. Porting about 2300 lines of Cython, most of which are duplicated or generated code, is less work than it seems: the dimension-generic, vectorized equivalent should be much shorter. The reference data in tests/data provide a regression oracle.

Recommendation. Replace Cython entirely, in two steps. First, write a pure NumPy/SciPy implementation, in functional style (no in-place updates), so that the same code can later run on JAX arrays (see E7). Second, add JAX as an optional backend, which brings speed (fusion), GPU and automatic differentiation. Numba is not recommended as the foundation of the architecture, because it does not solve the extensibility problem (user-defined physics would still have to be written as specialized kernels); it may be used locally, as an optional accelerator of the NumPy backend, if benchmarks show a bottleneck.

E4. Vectorized, stateless operators¶

Proposal. A continuum Green operator becomes a function of an array of wave-vectors: symbol(k) maps an array of shape (..., dim) to an array of shape (..., m, n) (or applies directly to an array of Fourier coefficients, without forming the matrices). Local operators act on whole fields. There is no set_frequency state.

Gains. Pure Python extensions become possible and efficient (one call per field, not per frequency); the code is simpler and thread-safe; the same function can be evaluated at modified wave-vectors (E6). For example, the Green operator of isotropic conductivity (conductivity k0) is written as follows, without compilation:

def green_symbol(k, k0):
    k2 = (k * k).sum(axis=-1)[..., None, None]
    kk = k[..., :, None] * k[..., None, :]
    return np.where(k2 > 0., kk / (k0 * np.where(k2 > 0., k2, 1.)), 0.)

Losses. Memory: evaluating a symbol on all frequencies requires a temporary array of matrices (see E3), which calls for evaluation by blocks or for closed-form application (e.g. for isotropic elasticity, the image of a polarization can be computed from the normalized wave-vector and a few scalars, without forming the 6×6 matrix).

Consequences. Complete rewrite of janus/green.pyx and janus/material/. The API changes from green.set_frequency(k); green.apply(tau) to green(k) @ tau (or green.apply(k, tau)).

Implementation problems. Treatment of the zero frequency (the symbol is not defined at k = 0) and of the Nyquist frequencies must be explicit and documented.

Recommendation. Core of the new architecture (milestone 0.3).

E5. Dimension-generic implementation¶

Proposal. A single implementation for 2D and 3D (and possibly 1D, useful for tests), based on the number of spatial axes of the arrays.

Gains. Roughly halves the code, and removes the risk of inconsistencies between the 2D and 3D versions (see the minor issues listed above); missing 3D variants (FourthRankCubicTensor3D) come for free when the underlying formula is dimension-independent.

Losses. Hand-unrolled 3×3/6×6 expressions are replaced by generic array operations, which are slightly slower in NumPy (not in JAX, where jit specializes the code for the actual shapes).

Consequences. Classes suffixed by 2D/3D disappear from the API.

Recommendation. Adopt together with E4.

E6. Separation of discretization and physics¶

Proposal. A discretization scheme is described independently of the physics, as a list of weighted modified wave-vectors associated with each discrete frequency; the discrete Green operator is then the weighted sum of the continuum symbol evaluated at these wave-vectors:

truncated scheme: one term (the wave-vector of the discrete frequency itself, with weight 1);
filtered scheme: 4 terms in 2D and 8 terms in 3D (neighbouring wave-vectors, with cos² weights), as currently implemented in FilteredGreenOperator2D/3D;
finite differences (willot2015): one term, with a modified wave-vector.

Gains. Any scheme works with any physics (e.g. the filtered scheme for conductivity), and new schemes can be added in a few lines of Python; scripts/gencode.py becomes useless.

Losses. The filtered scheme evaluates the continuum symbol 4 or 8 times per frequency, instead of combining precomputed matrices; this is also the case in the current code, so there is no actual loss.

Consequences. New public concept (Scheme), and new signature of the discrete Green operators (continuum operator + scheme + grid).

Implementation problems. Some schemes involve complex modified wave-vectors (the current implementation of willot2015 uses an equivalent real form); symbols should therefore accept complex wave-vectors and use conjugates where required (e.g. k ⊗ conj(k)). Schemes that are not of the “weighted sum” type (e.g. staggered grids) may require a more general interface; this should be kept in mind, but not designed for upfront.

Recommendation. Adopt in milestone 0.3, and validate it by reproducing the three existing schemes against the reference data.

E7. Array API and functional interface¶

Proposal. Write the library against the Python array API standard (via array-api-compat), in functional style: y = op(x) rather than op.apply(x, y). NumPy (version 2 and later) and JAX both implement this standard, including its FFT extension.

Gains. The same code runs on NumPy (CPU, default), JAX (jit, GPU, autodiff), and possibly CuPy or PyTorch, which keeps the “GPU and autodiff” option open at low cost.

Losses. In-place operations and output arguments (apply(x, y), apply(x, x)) disappear from the core API, which slightly increases memory allocations with NumPy (an optional out argument may be kept for the NumPy backend only). Only the functions of the standard can be used in the core (np.einsum is not part of it, for instance, although most libraries provide it).

Consequences. Breaking change for all operators (the documented in-place semantics, see Performing in-place operations, is dropped).

Implementation problems. Testing on several backends (at least NumPy and JAX) is required to guarantee compliance; double precision must be enabled explicitly with JAX.

Recommendation. Decide the functional interface as early as milestone 0.3, even if the JAX backend comes later (milestone 0.6): retrofitting it afterwards would require a second rewrite.

E8. Local operators as fields¶

Proposal. Replace block-diagonal operators made of arrays of Python objects (BlockDiagonalOperator2D of FourthRankIsotropicTensor, as in the tutorial) with array-based descriptions of the microstructure: a phase map (array of integers) and the stiffness of each phase, or more generally a constitutive function of one cell, applied to the whole field (vectorized in NumPy, vmap in JAX).

Gains. Pure Python, vectorized, and extensible to nonlinear laws; with JAX, the tangent operator of a nonlinear law is obtained by automatic differentiation, which is exactly what Newton–Krylov solvers for finite strain require.

Losses. Arbitrary heterogeneous Python objects per cell are no longer supported (they were not usable from Python anyway, see the first friction point).

Consequences. Rewrite of janus/operators.pyx; the fourth-rank tensor classes become small functions returning arrays (Mandel–Voigt matrices).

Recommendation. Milestone 0.5 (phase maps and linear laws may be needed as early as 0.3, to port the tutorial).

E9. Generic tensor representations¶

Proposal. Describe the local data of a field by a small object (vector of size dim, symmetric second-rank tensor in Mandel form, full second-rank tensor…), instead of assuming a vector of Mandel–Voigt components. The layout (*grid_shape, *local_shape) (grid axes first, local components last) should be kept: it is the current layout, and it is natural for NumPy broadcasting (matmul acts on the last two axes).

Gains. Conductivity and Darcy flow (vector fields), and finite strain (non-symmetric deformation gradients, 9 components in 3D) become possible. The Mandel form (janus/mandelvoigt.py) remains available for symmetric tensors; being orthonormal, it keeps transposes and scalar products simple.

Losses. Some added complexity in the API.

Recommendation. Milestone 0.5 (vectors), and 0.6 (full tensors, for the finite strain prototype).

E10. Solvers and interoperability¶

Proposal. Janus does not provide solvers, on purpose. For prototyping, however, a thin layer would reduce friction: adapters to scipy.sparse.linalg.LinearOperator (flattened vectors) and to the JAX solvers, and reference implementations of the basic scheme of Moulinec & Suquet and of a Newton–Krylov loop for nonlinear problems.

Gains. The tutorial becomes shorter; users compare schemes on a common basis.

Losses. Scope creep: this layer must remain minimal.

Recommendation. Adapters in milestone 0.4; reference schemes in 0.5.

E11. Packaging, continuous integration and distribution¶

Proposal. Once Cython is gone, Janus becomes a pure Python package: configuration in pyproject.toml only (no setup.py, no setup.cfg), a universal wheel, publication on PyPI (and conda-forge), and continuous integration on GitHub Actions (tests on Linux, Windows and macOS, build of the documentation, possibly its deployment to GitHub Pages).

Gains. pip install in seconds, without compiler; all the fragile steps identified in Claude’s contributions to Janus disappear.

Implementation problems. The name janus is already taken on PyPI (by an unrelated package), so that a distribution name must be chosen (the import name can remain janus, although a distinct import name avoids clashes in environments that contain both packages). Minimum supported versions of Python and NumPy should follow a policy such as SPEC 0.

Recommendation. Continuous integration as early as milestone 0.2 (with the current Cython code); pure Python packaging and publication in milestone 0.4.

E12. Testing strategy for the rewrite¶

Proposal. Keep the reference data (tests/data/*.npz) as a regression oracle; keep the current Cython implementation available (in a tagged release) to generate additional references if needed; add tests of mathematical properties, independent of the implementation: the discrete Green operator is a projector (Γ₀ C₀ Γ₀ = Γ₀), it is symmetric, it annihilates uniform fields, its Fourier symbol is homogeneous of degree zero, and the results do not depend on the backend.

Gains. Confidence in the rewrite, and tests that remain valid for new physics.

Recommendation. From milestone 0.3 on.

Milestones¶

The evolutions above are ordered as follows. Each milestone leaves the library in a usable, tested state.

0.2 — Clean-up. Remove MPI (E1). Set up continuous integration with the current Cython code (E11). Tag the last Cython-based release, so that it remains available as a reference.

0.3 — New core. New pure NumPy/SciPy implementation, alongside the Cython code (e.g. in a new subpackage): FFT through scipy.fft (E2), vectorized and stateless operators (E4), dimension-generic code (E5), separation of discretization and physics (E6), functional interface (E7). Isotropic linear elasticity with the three existing schemes, validated against the reference data (E12). Benchmark against the Cython implementation, on 2D and 3D grids: the acceptable slowdown is to be decided by the author.

0.4 — Switch. The new core becomes the default; Cython, FFTW, setup.py and setup.cfg are deleted. The tutorials and the documentation are ported; adapters to SciPy solvers (E10). Pure Python packaging and first publication (E11).

0.5 — Extensibility. A second physics (conductivity or Darcy flow, i.e. vector fields, E9) implemented in pure Python, as a validation of the extension mechanism, with a user guide on adding physics and discretizations. Local operators as fields, including nonlinear laws (E8); reference iterative schemes (E10).

0.6 — Backends. Array API compliance tested on NumPy and JAX (E7); JAX backend (jit, optional GPU, double precision); automatic differentiation of local laws; finite strain prototype with full second-rank tensors (E9).

1.0 — Stabilization. API freeze, complete documentation, removal of deprecated code, release on PyPI and conda-forge.

The order matters in two places. MPI is removed first, because it reduces the amount of code to port. The functional interface is adopted in 0.3, although JAX comes in 0.6, because changing the interface afterwards would require a second rewrite.

Decisions required from the author¶

Distribution name on PyPI (janus is taken), and import name.
Data layout: keep grid axes first and local components last (recommended, see E9)?
Default backend: NumPy with optional JAX (recommended), or JAX only?
Precision: float64 only, or float32 as well (relevant on GPU)?
Acceptable slowdown of the pure Python implementation with respect to the current Cython code (acceptance criterion of milestone 0.3).
Minimum supported versions of Python and NumPy.