Fix type stabilities for FVMSystems

.github/workflows/CI.yml

@@ -62,11 +62,19 @@ jobs:
           files: lcov.info
   docs:
     name: Documentation
-    runs-on: ubuntu-latest
+    runs-on: ubuntu-20.04
     permissions:
       contents: write
       statuses: write
     steps:
+      - name: Increase swapfile
+        run: |
+          sudo swapoff -a
+          sudo fallocate -l 8G /swapfile
+          sudo chmod 600 /swapfile
+          sudo mkswap /swapfile
+          sudo swapon /swapfile
+          sudo swapon --show      
       - uses: actions/checkout@v4
       - uses: julia-actions/setup-julia@v1
         with:

.github/workflows/DocCleanup.yml

@@ -0,0 +1,28 @@
+name: Doc Preview Cleanup
+
+on:
+  pull_request:
+    types: [closed]
+
+jobs:
+  doc-preview-cleanup:
+    runs-on: ubuntu-latest
+    permissions:
+      contents: write
+    steps:
+      - name: Checkout gh-pages branch
+        uses: actions/checkout@v3
+        with:
+          ref: gh-pages
+      - name: Delete preview and history + push changes
+        run: |
+            if [ -d "previews/PR$PRNUM" ]; then
+              git config user.name "Documenter.jl"
+              git config user.email "[email protected]"
+              git rm -rf "previews/PR$PRNUM"
+              git commit -m "delete preview"
+              git branch gh-pages-new $(echo "delete history" | git commit-tree HEAD^{tree})
+              git push --force origin gh-pages-new:gh-pages
+            fi
+        env:
+            PRNUM: ${{ github.event.number }}

.gitignore

@@ -22,6 +22,7 @@ docs/src/tutorials/reaction_diffusion_equation_with_a_time_dependent_dirichlet_b
 docs/src/tutorials/diffusion_equation_on_a_square_plate.md
 docs/src/tutorials/equilibrium_temperature_distribution_with_mixed_boundary_conditions_and_using_ensembleproblems.md
 docs/src/tutorials/diffusion_equation_in_a_wedge_with_mixed_boundary_conditions.md
+docs/src/tutorials/gray_scott_model_turing_patterns_from_a_coupled_reaction_diffusion_system.md 
 docs/src/tutorials/solving_mazes_with_laplaces_equation.md
 docs/src/wyos/laplaces_equation.md
 docs/src/wyos/mean_exit_time.md

README.md

@@ -15,7 +15,7 @@ in two dimensions using the finite volume method, with support also provided for
 
 1. `DiffusionEquation`s: $\partial_tu = \boldsymbol\nabla\boldsymbol\cdot[D(\boldsymbol x)\boldsymbol\nabla u]$.
 2. `MeanExitTimeProblem`s: $\boldsymbol\nabla\boldsymbol\cdot[D(\boldsymbol x)\boldsymbol\nabla T(\boldsymbol x)] = -1$.
-3. `LinearReactionDiffusionEquation`s: $\partial_tu= \boldsymbol\nabla\boldsymbol\cdot[D(\boldsymbol x)\boldsymbol\nabla u] +f(\boldsymbol x)u$.
+3. `LinearReactionDiffusionEquation`s: $\partial_tu = \boldsymbol\nabla\boldsymbol\cdot[D(\boldsymbol x)\boldsymbol\nabla u] + f(\boldsymbol x)u$.
 4. `PoissonsEquation`: $\boldsymbol\nabla\boldsymbol\cdot[D(\boldsymbol x)\boldsymbol\nabla u] = f(\boldsymbol x)$.
 5. `LaplacesEquation`: $\boldsymbol\nabla\boldsymbol\cdot[D(\boldsymbol x)\boldsymbol\nabla u] = 0$.
 

docs/src/index.md

@@ -22,8 +22,8 @@ The [tutorials](tutorials/overview.md) in the sidebar demonstrate the many possi
 
 1. `DiffusionEquation`s: $\partial_tu = \div[D(\vb x)\grad u]$.
 2. `MeanExitTimeProblem`s: $\div[D(\vb x)\grad T(\vb x)] = -1$.
-3. `LinearReactionDiffusionEquation`s: $\partial_tu + \div[D(\vb x)\grad u] + f(\vb x)u$.
-4. `PoissonsEquation`: $\div[D(\vb x)\grad u]  f(\vb x)$.
+3. `LinearReactionDiffusionEquation`s: $\partial_tu = \div[D(\vb x)\grad u] + f(\vb x)u$.
+4. `PoissonsEquation`: $\div[D(\vb x)\grad u]  = f(\vb x)$.
 5. `LaplacesEquation`: $\div[D(\vb x)\grad u] = 0$.
 
 See the [Solvers for Specific Problems, and Writing Your Own](wyos/overview.md) section for more information on these templates.

docs/src/literate_tutorials/gray_scott_model_turing_patterns_from_a_coupled_reaction_diffusion_system.jl

@@ -0,0 +1,99 @@
+# # Gray-Scott Model: Turing Patterns from a Coupled Reaction-Diffusion System
+#
+# In this tutorial, we explore some pattern formation from the
+# Gray-Scott model:
+# ```math
+# \begin{equation}
+# \begin{aligned}
+# \pdv{u}{t} &= \varepsilon_1\grad^2u+b(1-u)-uv^2, \\
+# \pdv{v}{t} &= \varepsilon_2\grad^2 v - dv+uv^2,
+# \end{aligned}
+# \end{equation}
+# ```
+# where $u$ and $v$ are the concentrations of two chemical species. The
+# initial conditions we use are:
+# ```math
+# \begin{align*}
+# u(x, y, 0) &= 1 -\exp\left[-80\left(x^2 + y^2\right)\right], \\
+# v(x, y, 0) &= \exp\left[-80\left(x^2+y^2\right)\right].
+# \end{align*}
+# ```
+# The domain we use is $[-1, 1]^2$, and we use
+# zero flux boundary conditions.
+using FiniteVolumeMethod, DelaunayTriangulation
+tri = triangulate_rectangle(-1, 1, -1, 1, 200, 200, single_boundary=true)
+mesh = FVMGeometry(tri)
+
+#-
+bc = (x, y, t, (u, v), p) -> zero(u) * zero(v)
+u_BCs = BoundaryConditions(mesh, bc, Neumann)
+v_BCs = BoundaryConditions(mesh, bc, Neumann)
+
+#-
+ε₁ = 0.00002
+ε₂ = 0.00001
+b = 0.04
+d = 0.1
+u_q = (x, y, t, α, β, γ, _ε₁) -> (-α[1] * _ε₁, -β[1] * _ε₁)
+v_q = (x, y, t, α, β, γ, _ε₂) -> (-α[2] * _ε₂, -β[2] * _ε₂)
+u_S = (x, y, t, (u, v), _b) -> _b * (1 - u) - u * v^2
+v_S = (x, y, t, (u, v), _d) -> -_d * v + u * v^2
+u_qp = ε₁
+v_qp = ε₂
+u_Sp = b
+v_Sp = d
+u_icf = (x, y) -> 1 - exp(-80 * (x^2 + y^2))
+v_icf = (x, y) -> exp(-80 * (x^ 2 + y^2))
+u_ic = [u_icf(x, y) for (x, y) in each_point(tri)]
+v_ic = [v_icf(x, y) for (x, y) in each_point(tri)]
+u_prob = FVMProblem(mesh, u_BCs;
+    flux_function=u_q, flux_parameters=u_qp,
+    source_function=u_S, source_parameters=u_Sp,
+    initial_condition=u_ic, final_time=6000.0)
+v_prob = FVMProblem(mesh, v_BCs;
+    flux_function=v_q, flux_parameters=v_qp,
+    source_function=v_S, source_parameters=v_Sp,
+    initial_condition=v_ic, final_time=6000.0)
+prob = FVMSystem(u_prob, v_prob)
+
+# Now that we have our system, we can solve.
+using OrdinaryDiffEq, Sundials
+sol = solve(prob, CVODE_BDF(linear_solver=:GMRES), saveat=10.0, parallel=Val(false)) 
+
+# Here is an animation of the solution, looking only at the $v$ variable.
+using CairoMakie
+fig = Figure(fontsize=33)
+ax = Axis(fig[1, 1], xlabel=L"x", ylabel=L"y")
+tightlimits!(ax)
+i = Observable(1)
+u = map(i -> reshape(sol.u[i][2, :], 200, 200), i)
+x = LinRange(-1, 1, 200)
+y = LinRange(-1, 1, 200)
+heatmap!(ax, x, y, u, colorrange=(0.0, 0.4))
+hidedecorations!(ax)
+record(fig, joinpath(@__DIR__, "../figures", "gray_scott_patterns.mp4"), eachindex(sol);
+    framerate=60) do _i
+    i[] = _i
+end
+
+# ![Animation of the Gray-Scott model](../figures/gray_scott_patterns.mp4)
+
+using ReferenceTests #src
+fig = Figure(fontsize=66) #src
+times = [0,1000,2000,3000,4000,5000,6000] #src
+t_idx = [findlast(≤(τ), sol.t) for τ in times] #src
+plotij = [(1,1),(1,2),(1,3),(2,1),(2,2),(2,3),(3,1)] #src
+for (i, j) in enumerate(t_idx) #src
+    ax = Axis(fig[plotij[i]...], width=600,height=600, xlabel=L"x", ylabel=L"y", title="t = $(sol.t[j])") #src
+    heatmap!(ax, x, y, reshape(sol.u[j][2,:],200,200), colorrange=(0.0, 0.4)) #src
+    tightlimits!(ax) #src
+end  #src
+plotij = [(4,1),(4,2),(4,3),(5,1),(5,2),(5,3),(6,1)] #src
+for (i, j) in enumerate(t_idx) #src
+    ax = Axis(fig[plotij[i]...], width=600,height=600, xlabel=L"x", ylabel=L"y", title="t = $(sol.t[j])") #src
+    heatmap!(ax, x, y, reshape(sol.u[j][1,:],200,200), colorrange=(0.0, 1)) #src
+    tightlimits!(ax) #src
+end  #src
+resize_to_layout!(fig) #src
+fig #src
+@test_reference joinpath(@__DIR__, "../figures", "gray_scott_patterns.png") fig #src

docs/src/literate_wyos/diffusion_equations.jl

@@ -280,7 +280,7 @@ fvm_prob = FVMProblem(mesh, BCs;
     initial_condition,
     final_time)
 
-using LinearSolve #src
+using LinearSolve
 
 #-
 # ```@example