Kernel formulation

.JuliaFormatter.toml

@@ -1,2 +1,5 @@
+whitespace_in_kwargs = true
 remove_extra_newlines = true
 trailing_comma = true
+separate_kwargs_with_semicolon = true
+long_to_short_function_def = true

.github/workflows/CI.yml

@@ -6,6 +6,7 @@ on:
       - main
     tags: ['*']
   pull_request:
+  workflow_dispatch:
 
 concurrency:
   # Skip intermediate builds: always.
@@ -51,10 +52,12 @@ jobs:
           version: '1'
       - uses: julia-actions/julia-buildpkg@v1
       - uses: julia-actions/julia-docdeploy@v1
+        with:
+          install-package: false
         env:
           GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }}
       - run: |
-          julia --project=docs -e '
+          julia --project=docs -t auto -e '
             using Documenter: DocMeta, doctest
             using IncompressibleNavierStokes
             DocMeta.setdocmeta!(

.github/workflows/FormatCheck.yml

@@ -5,6 +5,7 @@ on:
     branches:
       - "main"
   pull_request:
+  workflow_dispatch:
 
 jobs:
   build:

.github/workflows/SpellCheck.yml

@@ -0,0 +1,13 @@
+name: Spell Check
+
+on: [pull_request]
+
+jobs:
+  typos-check:
+    name: Spell Check with Typos
+    runs-on: ubuntu-latest
+    steps:
+      - name: Checkout Actions Repository
+        uses: actions/checkout@v4
+      - name: Check spelling
+        uses: crate-ci/typos@master

.github/workflows/TagBot.yml

@@ -5,7 +5,22 @@ on:
     types:
       - created
   workflow_dispatch:
-
+    inputs:
+      lookback:
+        default: 3
+permissions:
+  actions: read
+  checks: read
+  contents: write
+  deployments: read
+  issues: read
+  discussions: read
+  packages: read
+  pages: read
+  pull-requests: read
+  repository-projects: read
+  security-events: read
+  statuses: read
 jobs:
   TagBot:
     if: github.event_name == 'workflow_dispatch' || github.actor == 'JuliaTagBot'

.gitignore

@@ -1,12 +1,15 @@
+.vscode
 *.jl.*.cov
 *.jl.cov
 *.jl.mem
 *Manifest.toml
 /docs/build/
 /docs/src/generated/
 *output/
-*notebooks/
-examples/test.jl
+scratch/*
+test.jl
+toto.jl
+tata.jl
 
 # Vim swap
 *.swp

Project.toml

@@ -5,29 +5,50 @@ version = "0.4.1"
 
 [deps]
 Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
+ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
 FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
 IterativeSolvers = "42fd0dbc-a981-5370-80f2-aaf504508153"
+KernelAbstractions = "63c18a36-062a-441e-b654-da1e3ab1ce7c"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a"
+NNlib = "872c559c-99b0-510c-b3b7-b6c96a88d5cd"
+Observables = "510215fc-4207-5dde-b226-833fc4488ee2"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 WriteVTK = "64499a7a-5c06-52f2-abe2-ccb03c286192"
 
+[weakdeps]
+CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
+CUDSS = "45b445bb-4962-46a0-9369-b4df9d0f772e"
+
+[extensions]
+IncompressibleNavierStokesCUDAExt = ["CUDA", "CUDSS"]
+
 [compat]
 Adapt = "3, 4"
+CairoMakie = "0.11"
+ChainRulesCore = "1"
+CUDA = "5"
+CUDSS = "0.2"
 FFTW = "1"
+GLMakie = "0.9"
 IterativeSolvers = "0.9"
-Makie = "0.19, 0.20"
+KernelAbstractions = "0.9"
+LinearAlgebra = "1"
+Makie = "0.20"
+NNlib = "0.9"
+Observables = "0.5"
+Printf = "1"
+Random = "1"
+SparseArrays = "1"
+StaticArrays = "1"
 Statistics = "1"
 WriteVTK = "1"
-julia = "1.7"
+julia = "1.9"
 
 [extras]
-Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
 CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
 GLMakie = "e9467ef8-e4e7-5192-8a1a-b1aee30e663a"
-Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
-
-[targets]
-test = ["Aqua", "CairoMakie", "Test"]

README.md

@@ -9,26 +9,33 @@
 
 This package implements energy-conserving solvers for the incompressible Navier-Stokes
 equations on a staggered Cartesian grid. It is based on the Matlab package
-[INS2D](https://github.com/bsanderse/INS2D)/[INS3D](https://github.com/bsanderse/INS3D).
+[INS2D](https://github.com/bsanderse/INS2D)/[INS3D](https://github.com/bsanderse/INS3D). The simulations can be run on the single/multithreaded CPUs or Nvidia GPUs.
 
+This package also provides experimental support for neural closure models for
+large eddy simulation.
 
 ## Installation
 
 To install IncompressibleNavierStokes, open up a Julia-REPL, type `]` to get
 into Pkg-mode, and type:
 
-```
-add IncompressibleNavierStokes
+```julia-repl
+(v1.10) pkg> add IncompressibleNavierStokes
 ```
 
 which will install the package and all dependencies to your local environment.
-Note that IncompressibleNavierStokes requires Julia version `1.7` or above.
+Note that IncompressibleNavierStokes requires Julia version `1.9` or above.
 
 See the
 [Documentation](https://agdestein.github.io/IncompressibleNavierStokes.jl/dev/generated/LidDrivenCavity2D/)
 for examples of some typical workflows. More examples can be found in the
 [`examples`](examples) directory.
 
+## Source code for paper
+
+See [here](./lib/PaperDC) for the source code used in the paper
+[Discretize first, filter next: learning divergence-consistent closure models for large-eddy simulation](https://arxiv.org/abs/2403.18088).
+
 ## Gallery
 
 The velocity and pressure fields may be visualized in a live session using
@@ -37,97 +44,96 @@ The velocity and pressure fields may be visualized in a live session using
 snapshot files using the `save_vtk` function, or the full time series using the
 `VTKWriter` processor.
 
-| ![](assets/examples/Actuator2D.png)     | ![](assets/examples/BackwardFacingStep2D.png)                 | ![](assets/examples/DecayingTurbulence2D.png)                | ![](assets/examples/TaylorGreenVortex2D.png)                |
-|:---------------------------------------:|:-------------------------------------------------------------:|:------------------------------------------------------------:|:-----------------------------------------------------------:|
+| ![](assets/examples/Actuator2D.png) | ![](assets/examples/BackwardFacingStep2D.png) | ![](assets/examples/DecayingTurbulence2D.png) | ![](assets/examples/TaylorGreenVortex2D.png) |
+|:-:|:-:|:-:|:-:|
 | [Actuator (2D)](examples/Actuator2D.jl) | [Backward Facing Step (2D)](examples/BackwardFacingStep2D.jl) | [Decaying Turbulence (2D)](examples/DecayingTurbulence2D.jl) | [Taylor-Green Vortex (2D)](examples/TaylorGreenVortex2D.jl) |
-| ![](assets/examples/Actuator3D.png)     | ![](assets/examples/BackwardFacingStep3D.png)                 | ![](assets/examples/DecayingTurbulence3D.png)                | ![](assets/examples/TaylorGreenVortex3D.png)                |
+| ![](assets/examples/Actuator3D.png) | ![](assets/examples/BackwardFacingStep3D.png) | ![](assets/examples/DecayingTurbulence3D.png) | ![](assets/examples/TaylorGreenVortex3D.png) |
 | [Actuator (3D)](examples/Actuator3D.jl) | [Backward Facing Step (3D)](examples/BackwardFacingStep3D.jl) | [Decaying Turbulence (3D)](examples/DecayingTurbulence3D.jl) | [Taylor-Green Vortex (3D)](examples/TaylorGreenVortex3D.jl) |
 
+IncompressibleNavierStokes also supports adding a temperature equation.
+
+https://github.com/agdestein/IncompressibleNavierStokes.jl/assets/40632532/e4894d97-d890-4d4a-a43a-8eb2772d5130
 
 ## Demo
 
-The following example code  using a negative body force on a small rectangle
+The following example code using a negative body force on a small rectangle
 with an unsteady inflow. It simulates a wind turbine (actuator) under varying
 wind conditions.
 
+Make sure to have the `GLMakie` and `IncompressibleNavierStokes` installed:
+
+```julia
+using Pkg
+Pkg.add(["GLMakie", "IncompressibleNavierStokes"])
+```
+
+Then run run the following code to make a short animation:
+
 ```julia
 using GLMakie
 using IncompressibleNavierStokes
 
-# Viscosity model
-viscosity_model = LaminarModel(; Re = 100.0)
-
-# Boundary conditions: Unsteady BC requires time derivatives
-u_bc(x, y, t) = x ≈ 0.0 ? cos(π / 6 * sin(π / 6 * t)) : 0.0
-v_bc(x, y, t) = x ≈ 0.0 ? sin(π / 6 * sin(π / 6 * t)) : 0.0
-dudt_bc(x, y, t) = x ≈ 0.0 ? -(π / 6)^2 * cos(π / 6 * t) * sin(π / 6 * sin(π / 6 * t)) : 0.0
-dvdt_bc(x, y, t) = x ≈ 0.0 ? (π / 6)^2 * cos(π / 6 * t) * cos(π / 6 * sin(π / 6 * t)) : 0.0
-bc_type = (;
-    u = (; x = (:dirichlet, :pressure), y = (:symmetric, :symmetric)),
-    v = (; x = (:dirichlet, :symmetric), y = (:pressure, :pressure)),
+# A 2D grid is a Cartesian product of two vectors
+n = 40
+x = LinRange(0.0, 10.0, 5n + 1)
+y = LinRange(-2.0, 2.0, 2n + 1)
+
+# Boundary conditions
+boundary_conditions = (
+    # Inlet, outlet
+    (
+        # Unsteady BC requires time derivatives
+        DirichletBC(
+            (dim, x, y, t) -> sin(π / 6 * sin(π / 6 * t) + π / 2 * (dim() == 1)),
+            (dim, x, y, t) ->
+                (π / 6)^2 *
+                cos(π / 6 * t) *
+                cos(π / 6 * sin(π / 6 * t) + π / 2 * (dim() == 1)),
+        ),
+        PressureBC(),
+    ),
+
+    # Sides
+    (PressureBC(), PressureBC()),
 )
 
-# A 2D grid is a Cartesian product of two vectors
-x = LinRange(0.0, 10.0, 200)
-y = LinRange(-2.0, 2.0, 80)
-
-# Actuator body force: A thrust coefficient `Cₜ` distributed over a thin rectangle
-xc, yc = 2.0, 0.0 # Disk center
-D = 1.0           # Disk diameter
-δ = 0.11          # Disk thickness
-Cₜ = 5e-4         # Thrust coefficient
-cₜ = Cₜ / (D * δ)
-inside(x, y) = abs(x - xc) ≤ δ / 2 && abs(y - yc) ≤ D / 2
-bodyforce_u(x, y) = -cₜ * inside(x, y)
-bodyforce_v(x, y) = 0.0
+# Actuator body force: A thrust coefficient distributed over a thin rectangle
+inside(x, y) = abs(x - 2.0) ≤ 0.055 && abs(y) ≤ 0.5
+bodyforce(dim, x, y, t) = dim() == 1 && inside(x, y) ? -1.82 : 0.0
 
 # Build setup and assemble operators
-setup = Setup(
-    x,
-    y;
-    viscosity_model,
-    u_bc,
-    v_bc,
-    dudt_bc,
-    dvdt_bc,
-    bc_type,
-    bodyforce_u,
-    bodyforce_v,
-);
-
-# Time interval
-t_start, t_end = tlims = (0.0, 40.0)
+setup = Setup(x, y; Re = 100.0, boundary_conditions, bodyforce);
 
 # Initial conditions (extend inflow)
-initial_velocity_u(x, y) = 1.0
-initial_velocity_v(x, y) = 0.0
-initial_pressure(x, y) = 0.0
-V₀, p₀ = create_initial_conditions(
-    setup,
-    initial_velocity_u,
-    initial_velocity_v,
-    t_start;
-    initial_pressure,
-);
-
-# Time step processors
-processors = (
-    # Record solution every fourth time step
-    animator(setup, "vorticity.mp4"; nupdate = 4),
-
-    # Log time step information
-    step_logger(),
-)
+ustart = create_initial_conditions(setup, (dim, x, y) -> dim() == 1 ? 1.0 : 0.0);
 
 # Solve unsteady Navier-Stokes equations
-V, p, outputs = solve_unsteady(
-    setup, V₀, p₀, tlims;
-    method = RK44P2(),
-    Δt = 0.05,
-    processors,
+solve_unsteady(;
+    setup, ustart, tlims = (0.0, 48.0), Δt = 0.05,
+    processors = (
+        anim = animator(; setup, path = "vorticity.mp4", nupdate = 4),
+        log = timelogger(),
+    ),
 )
 ```
 
 The resulting animation is shown below.
 
 https://github.com/agdestein/IncompressibleNavierStokes.jl/assets/40632532/6ee09a03-1674-46e0-843c-000f0b9b9527
+
+## Similar projects
+
+- [WaterLily.jl](https://github.com/weymouth/WaterLily.jl/)
+  Incompressible solver with immersed boundaries
+- [Oceananigans.jl](https://github.com/CliMA/Oceananigans.jl):
+  Ocean simulations
+- [ClimaCore.jl](https://github.com/CliMA/ClimaCore.jl):
+  Atmospheric simulations
+- [Trixi.jl](https://github.com/trixi-framework/Trixi.jl):
+  High order solvers for various hyperbolic equations
+- [Ferrite.jl](https://github.com/Ferrite-FEM/Ferrite.jl):
+  Finite element discretizations
+- [Gridap.jl](https://github.com/gridap/Gridap.jl):
+  Finite element discretizations
+- [FourierFlows.jl](https://github.com/FourierFlows/FourierFlows.jl):
+  Pseudo-spectral discretizations

docs/Project.toml

@@ -1,9 +1,7 @@
 [deps]
-CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 DocumenterCitations = "daee34ce-89f3-4625-b898-19384cb65244"
-FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
-GLMakie = "e9467ef8-e4e7-5192-8a1a-b1aee30e663a"
+Examples = "318dbb63-4243-420f-99f2-d56058123f9d"
 IncompressibleNavierStokes = "5e318141-6589-402b-868d-77d7df8c442e"
-LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
 Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306"
+NeuralClosure = "099dac27-d7f2-4047-93d5-0baee36b9c25"

docs/liveserve.jl

@@ -0,0 +1,9 @@
+push!(LOAD_PATH, "@live-server")
+
+using LiveServer
+
+servedocs(;
+    foldername = @__DIR__,
+    literate = joinpath(@__DIR__, "..", "examples"),
+    skip_dir = joinpath(@__DIR__, "src", "generated"),
+)