cuPSS: a package for the pseudo-spectral integration of stochastic PDEs

paper:
1. attempting to incorporate stochastic components.
2. good choice of PDE examples.

software:
1. ease of installation and use.
2. text-based equation specification.

paper:
1. lack of convergence and correctness tests.
2. limited overview of possible applications.

software:
1. no phenomena specific to SPDEs (rather than plain PDEs) were demonstrated.
2. only spatially homogeneous equations are demonstrated.

# Acceptance criteria:
## 1. Benchmarking tests must be provided.
The manuscript includes a section titled "Benchmarking," where an extensive comparison of the speed of operation is presented, along with an advantage over existing solvers. However, the information required to judge whether faster time-stepping actually constitutes a true speedup is incomplete. What should be compared is the time required to achieve a given level of solution precision, not just the number of steps taken.
The methodology used for solving the problem in this work is quite different from that of existing solvers, such as FEniCSx. FEniCSx employs the finite element method (FEM), while this work uses the cuPSS spectral method. It would be valuable to investigate whether the error scales in the same manner, or differently, with respect to the parameter $N$.

This leads to a more fundamental issue with the submitted software. The GitHub repository contains tests that I would categorize as 'unit tests,' which assess the implementation of small software components. While these are a good practice and demonstrate mature coding, what is missing are 'integration tests.' These tests would examine whether the various components of the software work together as intended, and whether the software as a whole performs the desired task with satisfactory quality.

For deterministic PDEs, Figure 2 shows progress in this direction, with a comparison to the analytical temperature distribution in steady state. This is a strong first step. However, since the software aims to resolve dynamics as well as statics, it would be useful to examine the time-dependent behavior of some solutions. For example, one could plot the gradient at the midpoint of the distribution as a function of time and compare this with the analytical result to test the time dependence.

Constructing such tests is a challenging but important task. Analytical solutions to nonlinear PDEs are rare, so a judicious selection of test cases is crucial. Well-chosen examples can serve a dual purpose: highlighting the strengths and weaknesses of the software and satisfying the testing criteria put forward by SciPost Physics Codebases. Furthermore good selection of test cases is of interest even without the software to guide development of similar packages.

Moreover, cuPSS is advertised as a solver for stochastic PDEs, not just deterministic ones. Testing software that emulates nondeterministic behavior presents its own set of challenges. However, some summary statistics can still be compared with known analytical results. First, families of nonlinear stochastic differential equations with exact solutions are well-known and could be used to test the solver's accuracy in various ways (e.g., in terms of expected values, moments, distributions, and invariants). Second, certain results for SPDEs, such as phase transitions in Landau-Ginsburg-type models, are also known and could be used for testing.

Another potential testing approach is to check for self-consistency. One of the most common tests of this type is a convergence test, which examines the solver's behavior as the spatial or temporal mesh is refined. Given the explicit comparison with FEniCSx, it would be informative to compare the solutions themselves, not just the speed of the solvers.

To summarize, excellent benchmarking should test the software's correctness in solving progressively more complex problems, ideally in cases where ground truth can be established independently of the software. In the case of approximate methods, it is essential to assess the rate of convergence to the true solution with respect to control parameters. Once this is established, the time complexity required to obtain a figure of merit (such as mean error) should be compared.

While the manuscript has partially met these criteria, I hope that the concerns outlined above can be addressed in a revision.

## 2. At least one example application must be presented in detail
The models presented strike a good balance between complexity and simplicity, effectively showcasing cuPSS's capability in solving time-dependent PDEs. However, in my opinion, there is a notable absence of an application-focused component. Merely producing visualizations that resemble other simulations of the Cahn-Hilliard equation does not fully demonstrate the potential of a tool like cuPSS. While the visualizations are aesthetically pleasing (and the code to generate them from simulation results should ideally be included in the GitHub repository), a more insightful form of post-processing would better highlight the scientific usefulness of the presented tool. Given that the author cites several publications utilizing the deterministic capabilities of cuPSS, I assume this type of application-focused work would be relatively straightforward to include.

Additionally, none of the examples engage with the 'stochastic' aspect of SPDEs, which I hope is an omission due to the inherent complexity of this subject matter. Nonetheless, this disconnect between the manuscript’s abstract and its content is somewhat jarring.

## 3. High-level programming standards must be followed throughout the source code
The code repository is well-organized, and the code itself is clean, with sensible naming conventions and adherence to best practices, such as unit testing. While it is difficult to assess the full coverage of unit tests within the limited review timespan, it is clear that good practices have been followed. However, there is still room for improvement, particularly in the organization of examples. I would recommend focusing on this aspect as a priority. For inspiration, consider looking at [scikit-fem](https://github.com/kinnala/scikit-fem). While it is a relatively immature package, I believe it serves as a fair comparison and shows how "examples" can fuel package use and development.

## 4. The userguide must properly contextualize the software, describe the logic of its workings and highlight its added value as compared to existing software
The added value of the tool is presented by the authors as ease of use, speed, and the incorporation of stochastic terms. The case for ease of use and speed is strong and convincing. However, the argument regarding the incorporation of stochastic terms is less compelling.

## 5. The software must address a demonstrable need for the scientific community
The need for a generic nonlinear SPDE solver capable of eventually becoming an industry standard is significant. However, experience from computational fluid mechanics suggests that a single solver capable of handling any highly nonlinear, very noisy equation is unlikely to be realized. On the other hand, the absence of a widely agreed-upon engine for fast prototyping hampers investigations into areas such as the stability of weakly nonlinear equations under stochastic perturbations, to name just one example.

The decision to specify equations using text is a welcome design choice, as it is likely to appeal to many physicists. This approach stands in contrast to the object/function-driven specification commonly found in many other integrators, offering a novel and user-friendly alternative.

## 6. The documentation must be complete, including detailed instructions on downloading, installing and running the software
The installation instructions are straightforward and easy to follow. I was able to reproduce the Cahn-Hilliard example with minimal difficulty. For a CUDA-based application, cuPSS ranks highly in terms of ease of startup.

## Summary
Acceptance criteria 3, 5 and 6 are met to a satisfactory degree (although there is still room for improvement). Acceptance criteria 1, 2 and 4 are satisfied only partially, manuscript looks promising and I am hopeful these issues can be resolved with a revision.

-------

# Manuscript feedback
## Abstract
I am unaware of any SPDE solver in the most recent version of Mathematica. Listing MATLAB, Mathematica, and OpenFOAM in the abstract without further comparison feels unnecessary, especially since the solver is compared to Dedalus and FEniCSx.

## Introduction
Equation 1: There appears to be an issue with this equation. Assuming that $W$ represents the standard Wiener process, the path should be space-independent. This assumption conflicts with Equation 3. Simply stating that "I assume that SPDEs are written in the Itô formulation" does not resolve this issue. Stochastic differential equation (SDE) calculus fundamentally differs from ordinary differential equation (ODE) calculus. The Itô formulation specifically addresses the "pathological" behavior of the Wiener process with respect to time. It would be imprudent to apply the same methods used in ODEs to eliminate problematic terms in the spatial domain. This is because the choice of the Itô formulation (as opposed to, for example, the Stratonovich formulation) assigns a specific direction along the time axis, representing the direction of non-anticipation. Removing spatial parity symmetry in this context would be highly unconventional.

## 1.1 pseudo-spectral integration
Equation 2: This choice of notation is really unfortunate, especially in the case of nonlinear equations. Notation used in equation 1. is used to emphasize that stochastic calculus (rather than the ordinary calculus) is used. Changing notation to that of Equation 2. is a significant step backward in terms of notation development and is detrimental to science communication. Consult, for example, section "3.1.2 Criticism of the Naive Approach" of Öttinger's "Stochastic processes in polymeric fluids: tools and examples for developing simulation algorithms" (2012) for a great overview of problems that arise from pretending that SDEs are just ODEs with unusual forcing.

Line 49: Whether or not the noise term can be written as $D_0 q^2$ depends on the problem at hand and is not just an algebraic observation. It is a statement about regularity of the solutions and should be discussed further.

# 2. Models that cuPSS can solve
Equation 10: At face value this choice looks rather constraining, the equations appear to be homogeneous with no explicit dependence on location allowed - this might be just an illusion which can be fixed by a clever use of the "callback functions" capability but I think it should be explicitly discussed. This is a serious feature request (rather than an exotic nitpick) -- being able to work in, for example, cylindrical coordinates or spherical coordinates is an essential capability for a physics oriented package, this can be accommodated by allowing position dependent operators.

Line 144: Milstein scheme is stochastic equivalent of midpoint scheme, it improves strong convergence but not weak convergence - compare naming scheme with Kloeden Platen "Numerical Solution of Stochastic Differential Equations" (reference [20]). For a physicist strong convergence order is rarely important, choice of Milstein (or implicit Milstein) should be discussed further. Some comparison (in terms of convergence or stability) for a real SPDE would make arguments for Milstein scheme much more convincing.

# 3. Writing Models
Line 223: Picking a higher power would be more informative. does q^4 correspond to $\partial_{xxxx} + \partial_{yyyy}$ or $\partial_{xxxx} + 2 \partial_{xxyy} + \partial_{yyyy}$?

# 4.1 Model B (Cahn-Hillard)
Figure 1: what are the spatial dimensions of the picture? is it really 256x256 or is it 256*dx by 256*dy? Either way $\Delta x$ should be specified. I would prefer graphs plotted in real axes not as a function of discretization index.

# 5. Callback functions
Figure 2: Please use the same vertical scale for both panels of figure 2 for better clarity

# Software
Dependencies: it looks like the software works just fine with cuda toolkit v10 - this is important because its the current default on older ubuntu distros. NVIDIA CUDA development toolkit:
1. focal (20.04LTS): 10.1.243-3
2. jammy (22.04LTS): 11.5.1-1ubuntu1
3. noble (24.04LTS): 12.0.140~12.0.1-4build4

# Nitpicks
- Line 124: "mathcal" is a typo
- Equation 16: brackets a typo
- Figure 1: there is a typo "t = ???"
- Consider changing GitHub banner to animated Cahn-Hillard, the picture of a noisy Gaussian doesn't highlight the capabilities of the package very well.

First, the manuscript should engage more deeply with the challenges inherent in stochastic calculus. It should clearly articulate the specific difficulties associated with the integration of stochastic differential equations, with appropriate emphasis on problems arising in the context of stochastic partial differential equations. This deeper engagement should be evident in the following sections:
# 1. The introductory section, where the mathematical, physical, and historical context of SPDEs should be established.
# 2. The examples section, which should highlight behaviours and phenomena unique to SPDEs, rather than those already well-characterised in deterministic PDEs.
# 3. The benchmarking section, which should include tests that directly address the non-deterministic nature of SPDE solutions.

Second, the testing methodology should prioritise correctness before evaluating computational efficiency. Specifically:
# 1. Correctness and convergence tests should be included that demonstrate appropriate behaviour in both deterministic and stochastic settings.
# 2. Speed comparisons should be conducted with greater methodological care. It must be recognised, for instance, that the parameter $N$ in a spectral method is not directly comparable to $N$ in a finite element or finite difference method. Similarly, while Euler time-stepping may be faster than Runge-Kutta, this does not automatically imply superior performance in terms of accuracy or stability.

Ask for major revision