Example for using multiple states of a system in parallel

.github/workflows/ci.yml

@@ -53,7 +53,9 @@ jobs:
             ${{ runner.os }}-
       - uses: julia-actions/julia-buildpkg@v1
       - uses: julia-actions/julia-runtest@v1
-      - uses: julia-actions/julia-processcoverage@v1
+         env:
+          JULIA_NUM_THREADS: 4
+     - uses: julia-actions/julia-processcoverage@v1
       - uses: codecov/codecov-action@v4
   docs:
     name: Documentation

docs/Project.toml

@@ -3,6 +3,7 @@ AMGCLWrap = "4f76b812-4ba5-496d-b042-d70715554288"
 AlgebraicMultigrid = "2169fc97-5a83-5252-b627-83903c6c433c"
 CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
 Catalyst = "479239e8-5488-4da2-87a7-35f2df7eef83"
+ChunkSplitters = "ae650224-84b6-46f8-82ea-d812ca08434e"
 DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
 DiffResults = "163ba53b-c6d8-5494-b064-1a9d43ac40c5"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"

docs/src/solver.md

@@ -38,7 +38,8 @@ fixed_timesteps!
 VoronoiFVM.SystemState
 VoronoiFVM.SystemState(::Type, system::VoronoiFVM.AbstractSystem; data)
 VoronoiFVM.SystemState(system::VoronoiFVM.AbstractSystem; data)
-VoronoiFVM.solve!(system::VoronoiFVM.SystemState; kwargs...)
+VoronoiFVM.solve!(state::VoronoiFVM.SystemState; kwargs...)
+Base.similar(state::VoronoiFVM.SystemState; kwargs...)
 ```
 
 ### Linear solver strategies

examples/Example440_ParallelState.jl

@@ -0,0 +1,57 @@
+#=
+
+# 440: Parallel solves
+ ([source code](@__SOURCE_URL__))
+
+Demonstrate how to solve one system with different data in parallel
+using the SystemState (new in v2.0).
+=#
+module Example440_ParallelState
+
+using VoronoiFVM, ExtendableGrids
+using GridVisualize
+using ChunkSplitters
+
+function flux(y,u,node,data)
+    y[1]=u[1,1]^2-u[1,2]^2
+end
+
+function bcondition(y,u,node,data)
+    boundary_dirichlet!(y,u,node, species=1, region=1, value=0.1)
+    boundary_neumann!(y,u,node,species=1,region=2,value=data.influx)
+end
+
+function main(;nref=5, Plotter=nothing)
+    grid=simplexgrid(range(0,1,length=10*2^nref+1))
+    sys=VoronoiFVM.System(grid;flux,bcondition, species=[1], data=(influx=0.0,))
+
+    ## Initial state. First solution creates the matrix
+    state0=VoronoiFVM.SystemState(sys)
+    sol=solve!(state0;inival=0.1)
+
+    ## Prepare parameter and result data
+    influxes=range(0.0,10.0, length=100)
+    masses=similar(influxes)
+
+    ## Split the index range in as many chunks as threads
+    Threads.@threads for indexes in chunks(influxes;n=Threads.nthreads())
+        ## Create a new state sharing the system - one for each chunk
+        state=similar(state0)
+        ## Solve for all data values in chunk
+        for iflux in indexes
+            data=(influx=influxes[iflux],)
+            sol=solve!(state;data, inival=0.1, verbose="")
+            masses[iflux]=integrate(sys,sol)[1,1]
+        end
+    end
+    scalarplot(influxes, masses;Plotter, xlabel="influx", ylabel="mass")
+    sum(masses)
+end
+
+using Test
+function runtests()
+    testval=140.79872772042577
+    @test main() ≈ testval
+end
+
+end

src/VoronoiFVM.jl

@@ -47,7 +47,7 @@ using RecursiveArrayTools: RecursiveArrayTools, AbstractDiffEqArray
 import RecursiveFactorization
 using SciMLBase: SciMLBase
 using SparseArrays: SparseArrays, SparseMatrixCSC, dropzeros!, nonzeros,
-                    nzrange, spzeros
+                    nzrange, spzeros, issparse
 using SparseDiffTools: SparseDiffTools, forwarddiff_color_jacobian!,
                        matrix_colors
 using StaticArrays: StaticArrays, @MVector, @SArray, @SMatrix

src/vfvm_state.jl

@@ -134,4 +134,27 @@ Shortcut for creating state with value type defined by `Tv` type parameter of sy
 """
 SystemState(system::AbstractSystem{Tv,Tc,Ti,Tm}; kwargs...) where {Tv,Tc,Ti,Tm} =SystemState(Tv, system; kwargs...) 
 
+"""
+    similar(state; data=state.data)
 
+Create a new state of with the same system, different work arrays, and possibly different data.
+The matrix of the new state initially shares the sparsity structure with `state`.
+"""
+function Base.similar(state::SystemState; data=state.data)
+    system=state.system
+    solution=similar(state.solution)
+    if issparse(state.matrix)
+        csc=SparseMatrixCSC(state.matrix)
+        cscnew=SparseMatrixCSC(csc.m, csc.n, csc.colptr, csc.rowval, similar(csc.nzval))
+        matrix=ExtendableSparseMatrix(cscnew)
+    else
+        matrix=similar(state.matrix)
+    end
+    dudp=similar(state.dudp)
+    residual=similar(state.residual)
+    update=similar(state.update)
+    linear_cache=nothing
+    uhash=zero(UInt64)
+    history=nothing
+    SystemState(system, data, solution, matrix, dudp, residual, update, linear_cache, uhash, history)
+end

test/Project.toml

@@ -4,6 +4,7 @@ AlgebraicMultigrid = "2169fc97-5a83-5252-b627-83903c6c433c"
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
 CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
 Catalyst = "479239e8-5488-4da2-87a7-35f2df7eef83"
+ChunkSplitters = "ae650224-84b6-46f8-82ea-d812ca08434e"
 DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
 DiffResults = "163ba53b-c6d8-5494-b064-1a9d43ac40c5"
 ExampleJuggler = "3bbe58f8-ed81-4c4e-a134-03e85fcf4a1a"

test/runtests.jl

@@ -89,9 +89,13 @@ function run_all_tests(; run_notebooks = false, notebooksonly = false)
 end
 
 # Don't run notebooks on 1.12: https://github.com/fonsp/Pluto.jl/issues/2939
-run_all_tests(; run_notebooks = VERSION < v"1.12.0-DEV.0" , notebooksonly = false)
-#run_all_tests(; run_notebooks=true, notebooksonly = true)
-
+if haskey(ENV,"EXAMPLES_ONLY")
+    run_all_tests(; run_notebooks = false, notebooksonly = false)
+elseif haskey(ENV,"NOTEBOOKS_ONLY")
+    run_all_tests(; run_notebooks=true, notebooksonly = true)
+else
+    run_all_tests(; run_notebooks = VERSION < v"1.12.0-DEV.0" , notebooksonly = false)
+end