Fix docs

docs/make.jl

@@ -92,7 +92,7 @@ makedocs(;
     # draft = true,
     # clean = !localdev,
     modules = [IncompressibleNavierStokes, NeuralClosure],
-    warnonly = true, # Reexporting KernelAbstractions.CPU fails otherwise
+    # warnonly = true, # Reexporting KernelAbstractions.CPU fails otherwise
     plugins = [bib],
     authors = "Syver Døving Agdestein, Benjamin Sanderse, and contributors",
     repo = Remotes.GitHub("agdestein", "IncompressibleNavierStokes.jl"),

docs/references.bib

@@ -22,14 +22,19 @@ @misc{Sanderse2023
     Month         = {Jul},
     Url           = {http://arxiv.org/abs/2307.10874v1},
 }
-@article{Agdestein2024,
-  archiveprefix = {arXiv},
-  author        = {Syver Døving Agdestein and Benjamin Sanderse},
-  eprint        = {2403.18088},
-  primaryclass  = {math.NA},
-  title         = {Discretize first, filter next: learning divergence-consistent closure models for large-eddy simulation},
-  year          = {2024}
+@article{Agdestein2025,
+    title = {Discretize first, filter next: Learning divergence-consistent closure models for large-eddy simulation},
+    journal = {Journal of Computational Physics},
+    volume = {522},
+    pages = {113577},
+    year = {2025},
+    issn = {0021-9991},
+    doi = {https://doi.org/10.1016/j.jcp.2024.113577},
+    url = {https://www.sciencedirect.com/science/article/pii/S0021999124008258},
+    author = {Syver Døving Agdestein and Benjamin Sanderse},
+    keywords = {Discrete filtering, Closure modeling, Divergence-consistency, Large-eddy simulation, Neural ODE, A-posteriori training}
 }
+
 @article{Beck2021,
   author   = {Beck, Andrea and Kurz, Marius},
   doi      = {10.1002/gamm.202100002},

docs/src/.vitepress/config.mts

@@ -94,6 +94,7 @@ export default defineConfig({
               { text: 'Differentiating code', link: '/manual/differentiability' },
               { text: 'Temperature equation', link: '/manual/temperature' },
               { text: 'Large eddy simulation', link: '/manual/les' },
+              { text: 'SciML', link: '/manual/sciml' },
             ],
           },
         ],
@@ -189,6 +190,7 @@ export default defineConfig({
               { text: 'Differentiating code', link: '/manual/differentiability' },
               { text: 'Temperature equation', link: '/manual/temperature' },
               { text: 'Large eddy simulation', link: '/manual/les' },
+              { text: 'SciML', link: '/manual/sciml' },
             ],
           },
         ],

docs/src/manual/sciml.md

@@ -17,20 +17,20 @@ Using IncompressibleNavierStokes it is possible to write the momentum equations
 \end{align*}
 ```
 
-The derivation and the drawbacks of this approach are discussed in the [documentation](/docs/src/manual/spatial.md).
+The derivation and the drawbacks of this approach are discussed in the [documentation](spatial.md).
 
 This projected right-hand side can be used in the SciML solvers to solve the Navier-Stokes equations. The following example shows how to use the SciML solvers to solve the ODEs obtained from the Navier-Stokes equations.
 
 ```julia
-using DifferentialEquations 
-f(u, p, t) = create_right_hand_side(setup, psolver) 
+using DifferentialEquations
+f(u, p, t) = create_right_hand_side(setup, psolver)
 u0 = INITIAL_CONDITION
 tspan = (0.0, 1.0)     # time span where to solve.
 problem = ODEProblem(f, u0, tspan) #SciMLBase.ODEProblem
 sol = solve(problem, Tsit5(), reltol = 1e-8, abstol = 1e-8) # sol: SciMLBase.ODESolution
 ```
 
-Alternatively, it is also possible to use an [in-place formulation](https://docs.sciml.ai/DiffEqDocs/stable/basics/problem/#In-place-vs-Out-of-Place-Function-Definition-Forms) 
+Alternatively, it is also possible to use an [in-place formulation](https://docs.sciml.ai/DiffEqDocs/stable/basics/problem/#In-place-vs-Out-of-Place-Function-Definition-Forms)
 
 ```julia
 f(du,u,p,t) = right_hand_side!(du, u, Ref([setup, psolver]), t)
@@ -39,8 +39,8 @@ that is usually faster than the out-of-place formulation.
 
 You can look [here](https://docs.sciml.ai/DiffEqDocs/stable/basics/overview/) for more information on how to use the SciML solvers and all the options available.
 
-## API 
+## API
 ```@autodocs
 Modules = [IncompressibleNavierStokes]
 Pages = ["sciml.jl"]
-```
+```

examples/TaylorGreenVortex2D.jl

@@ -49,11 +49,8 @@ function compute_convergence(; D, nlist, lims, Re, tlims, Δt, uref, backend = C
         )
         (; u, t), outputs = solve_unsteady(; setup, ustart, tlims, Δt, psolver)
         (; Ip) = setup.grid
-        a, b = T(0), T(0)
-        for α = 1:D
-            a += sum(abs2, u[α][Ip] - ut[α][Ip])
-            b += sum(abs2, ut[α][Ip])
-        end
+        a = sum(abs2, u[Ip, :] - ut[Ip, :])
+        b = sum(abs2, ut[Ip, :])
         e[i] = sqrt(a) / sqrt(b)
     end
     e

ext/IncompressibleNavierStokesMakieExt.jl

@@ -70,19 +70,6 @@ function plotgrid(x, y, z)
     fig
 end
 
-"""
-Animate a plot of the solution every `update` iteration.
-The animation is saved to `path`, which should have one
-of the following extensions:
-
-- ".mkv"
-- ".mp4"
-- ".webm"
-- ".gif"
-
-The plot is determined by a `plotter` processor.
-Additional `kwargs` are passed to `plot`.
-"""
 animator(;
     setup,
     path,
@@ -108,23 +95,6 @@ animator(;
         stream
     end
 
-"""
-Processor for plotting the solution in real time.
-
-Keyword arguments:
-
-- `plot`: Plot function.
-- `nupdate`: Show solution every `nupdate` time step.
-- `displayfig`: Display the figure at the start.
-- `screen`: If `nothing`, use default display.
-    If `GLMakie.screen()` multiple plots can be displayed in separate
-    windows like in MATLAB (see also `GLMakie.closeall()`).
-- `displayupdates`: Display the figure at every update (if using CairoMakie).
-- `sleeptime`: The `sleeptime` is slept at every update, to give Makie
-    time to update the plot. Set this to `nothing` to skip sleeping.
-
-Additional `kwargs` are passed to the `plot` function.
-"""
 realtimeplotter(;
     setup,
     plot = fieldplot,
@@ -149,30 +119,6 @@ realtimeplotter(;
         fig
     end
 
-"""
-Plot `state` field in pressure points.
-If `state` is `Observable`, then the plot is interactive.
-
-Available fieldnames are:
-
-- `:velocity`,
-- `:vorticity`,
-- `:streamfunction`,
-- `:pressure`.
-
-Available plot `type`s for 2D are:
-
-- `heatmap` (default),
-- `image`,
-- `contour`,
-- `contourf`.
-
-Available plot `type`s for 3D are:
-
-- `contour` (default).
-
-The `alpha` value gets passed to `contour` in 3D.
-"""
 fieldplot(state; setup, kwargs...) = fieldplot(
     setup.grid.dimension,
     state isa Observable ? state : Observable(state);
@@ -332,9 +278,6 @@ function fieldplot(
     fig
 end
 
-"""
-Create energy history plot.
-"""
 function energy_history_plot(state; setup)
     @assert state isa Observable "Energy history requires observable state."
     (; Ip) = setup.grid
@@ -351,22 +294,6 @@ function energy_history_plot(state; setup)
     fig
 end
 
-"""
-Create energy spectrum plot.
-The energy at a scalar wavenumber level ``\\kappa \\in \\mathbb{N}`` is defined by
-
-```math
-\\hat{e}(\\kappa) = \\int_{\\kappa \\leq \\| k \\|_2 < \\kappa + 1} | \\hat{e}(k) | \\mathrm{d} k,
-```
-
-as in San and Staples [San2012](@cite).
-
-Keyword arguments:
-
-- `sloperange = [0.6, 0.9]`: Percentage (between 0 and 1) of x-axis where the slope is plotted.
-- `slopeoffset = 1.3`: How far above the energy spectrum the inertial slope is plotted.
-- `kwargs...`: They are passed to [`observespectrum`](@ref).
-"""
 function energy_spectrum_plot(
     state;
     setup,

lib/PaperDC/README.md

@@ -1,7 +1,7 @@
 # PaperDC
 
 Scripts for generating results of the paper
-[Discretize first, filter next: learning divergence-consistent closure models for large-eddy simulation](https://arxiv.org/abs/2403.18088).
+[Discretize first, filter next: Learning divergence-consistent closure models for large-eddy simulation](https://www.sciencedirect.com/science/article/pii/S0021999124008258).
 
 ## Set up environment
 
@@ -11,7 +11,9 @@ Run:
 julia --project=lib/PaperDC -e 'using Pkg; Pkg.instantiate()'
 ```
 
-Now you can run the scripts in this directory:
+Now you can run the scripts in this directory. They generate results for
 
-- `prioranalysis.jl`: Generate results for section 5.1 "Filtered DNS (2D and 3D)"
-- `postanalysis.jl`: Generate results for section 5.2 "LES (2D)"
+- `prioranalysis.jl`: Section 5.1 "Filtered DNS (2D and 3D)"
+- `postanalysis.jl`: Section 5.2 "LES (2D)"
+- `postanalysis3D.jl`: Appendix G.2. "LES of forced turbulence (3D)"
+- `transferfunctions.jl`: Appendix E. "Continuous filters and transfer functions"

lib/PaperDC/postanalysis.jl

@@ -1,6 +1,6 @@
 # # A-posteriori analysis: Large Eddy Simulation (2D)
 #
-# This script is used to generate results for the the paper [Agdestein2024](@citet).
+# This script is used to generate results for the the paper [Agdestein2025](@citet).
 #
 # - Generate filtered DNS data
 # - Train closure models

lib/PaperDC/postanalysis3D.jl

@@ -1,6 +1,6 @@
 # # A-posteriori analysis: Large Eddy Simulation (2D)
 #
-# This script is used to generate results for the the paper [Agdestein2024](@citet).
+# This script is used to generate results for the the paper [Agdestein2025](@citet).
 #
 # - Generate filtered DNS data
 # - Train closure models

lib/PaperDC/prioranalysis.jl

@@ -1,6 +1,6 @@
 # # A-priori analysis: Filtered DNS (2D or 3D)
 #
-# This script is used to generate results for the the paper [Agdestein2024](@citet).
+# This script is used to generate results for the the paper [Agdestein2025](@citet).
 #
 # - Generate filtered DNS data
 # - Compute quantities for different filters

src/IncompressibleNavierStokes.jl

@@ -55,6 +55,12 @@ using WriteVTK: CollectionFile, paraview_collection, vtk_grid, vtk_save
 "$LICENSE"
 license = "MIT"
 
+# We are reexporting KernelAbstractions.CPU, but
+# Documenter cannot find the docstring and complains.
+# Reapply the docstring here to keep Documenter happy.
+s = @doc(KernelAbstractions.CPU)
+@doc s.text[1] KernelAbstractions.CPU
+
 # # Easily retrieve value from Val
 # (::Val{x})() where {x} = x
 

src/operators.jl

@@ -197,7 +197,7 @@ end
     p[I] += adjoint
 end
 
-# "Subtract pressure gradient (differentiable version)."
+"Subtract pressure gradient (differentiable version)."
 applypressure(u, p, setup) = applypressure!(copy.(u), p, setup)
 
 ChainRulesCore.rrule(::typeof(applypressure), u, p, setup) = (

src/processors.jl

@@ -326,3 +326,86 @@ function realtimeplotter end
 function fieldplot end
 function energy_history_plot end
 function energy_spectrum_plot end
+
+# Add docstrings here, otherwise Documenter can't find them
+
+"""
+Animate a plot of the solution every `update` iteration.
+The animation is saved to `path`, which should have one
+of the following extensions:
+
+- ".mkv"
+- ".mp4"
+- ".webm"
+- ".gif"
+
+The plot is determined by a `plotter` processor.
+Additional `kwargs` are passed to `plot`.
+"""
+animator
+
+"""
+Processor for plotting the solution in real time.
+
+Keyword arguments:
+
+- `plot`: Plot function.
+- `nupdate`: Show solution every `nupdate` time step.
+- `displayfig`: Display the figure at the start.
+- `screen`: If `nothing`, use default display.
+    If `GLMakie.screen()` multiple plots can be displayed in separate
+    windows like in MATLAB (see also `GLMakie.closeall()`).
+- `displayupdates`: Display the figure at every update (if using CairoMakie).
+- `sleeptime`: The `sleeptime` is slept at every update, to give Makie
+    time to update the plot. Set this to `nothing` to skip sleeping.
+
+Additional `kwargs` are passed to the `plot` function.
+"""
+realtimeplotter
+
+"""
+Plot `state` field in pressure points.
+If `state` is `Observable`, then the plot is interactive.
+
+Available fieldnames are:
+
+- `:velocity`,
+- `:vorticity`,
+- `:streamfunction`,
+- `:pressure`.
+
+Available plot `type`s for 2D are:
+
+- `heatmap` (default),
+- `image`,
+- `contour`,
+- `contourf`.
+
+Available plot `type`s for 3D are:
+
+- `contour` (default).
+
+The `alpha` value gets passed to `contour` in 3D.
+"""
+fieldplot
+
+"Create energy history plot."
+energy_history_plot
+
+"""
+Create energy spectrum plot.
+The energy at a scalar wavenumber level ``\\kappa \\in \\mathbb{N}`` is defined by
+
+```math
+\\hat{e}(\\kappa) = \\int_{\\kappa \\leq \\| k \\|_2 < \\kappa + 1} | \\hat{e}(k) | \\mathrm{d} k,
+```
+
+as in San and Staples [San2012](@cite).
+
+Keyword arguments:
+
+- `sloperange = [0.6, 0.9]`: Percentage (between 0 and 1) of x-axis where the slope is plotted.
+- `slopeoffset = 1.3`: How far above the energy spectrum the inertial slope is plotted.
+- `kwargs...`: They are passed to [`observespectrum`](@ref).
+"""
+energy_spectrum_plot