Merge pull request #185 from Da-Be-Ru/test/cachesolve

Add cache-based solve test for multiple right-hand sides

Author: pjaap

CHANGELOG.md

@@ -6,6 +6,10 @@
   - Try to remove type piracies
   - Remove `params` from edge, node structs (apparently never used)
 
+## v2.9 April 3, 2025
+  - Fixed bug in `evaluate_residual_and_jacobian`
+  - Added `init_dirichlet` flag to `evaluate_residual_and_jacobian!` to suppress initialization of Dirichlet values in argument vector `u`
+
 ## v2.8 February 23, 2025
   - Add `integrate` methods for transient solutions returning instances of `DiffEqArray`
 

Project.toml

@@ -1,7 +1,7 @@
 name = "VoronoiFVM"
 uuid = "82b139dc-5afc-11e9-35da-9b9bdfd336f3"
 authors = ["Jürgen Fuhrmann <juergen.fuhrmann@wias-berlin.de>", "Patrick Jaap", "Daniel Runge", "Dilara Abdel", "Jan Weidner", "Alexander Seiler", "Patricio Farrell", "Matthias Liero"]
-version = "2.8.0"
+version = "2.9.0"
 
 [deps]
 BandedMatrices = "aae01518-5342-5314-be14-df237901396f"

src/VoronoiFVM.jl

@@ -77,7 +77,7 @@ solutionarray(a::AbstractSolutionArray) = a
 export AbstractSolutionArray
 
 include("vfvm_physics.jl")
-# see https://discourse.julialang.org/t/is-compat-jl-worth-it-for-the-public-keyword/119041/34m
+# see https://discourse.julialang.org/t/is-compat-jl-worth-it-for-the-public-keyword/119041/34
 VERSION >= v"1.11.0-DEV.469" && eval(Meta.parse("public Physics, AbstractPhysics, AbstractData"))
 
 include("vfvm_functions.jl")
@@ -91,6 +91,7 @@ export NewtonSolverHistory, TransientSolverHistory, details
 
 include("vfvm_densesolution.jl")
 include("vfvm_sparsesolution.jl")
+VERSION >= v"1.11.0-DEV.469" && eval(Meta.parse("public dofs"))
 export num_dof
 export dof
 export getdof

src/vfvm_solver.jl

@@ -192,8 +192,10 @@ function solve_step!(
     return solution
 end
 
-function evaluate_residual_and_jacobian!(state, u; params = Float64[], t = 0.0, tstep = Inf, embed = 0.0)
-    _initialize_dirichlet!(u, state.system, state.data)
+function evaluate_residual_and_jacobian!(state, u; params = Float64[], t = 0.0, tstep = Inf, embed = 0.0, init_dirichlet = false)
+    if init_dirichlet
+        _initialize_dirichlet!(u, state.system, state.data)
+    end
     eval_and_assemble(
         state.system,
         u,
@@ -213,16 +215,17 @@ end
 
 """
     evaluate_residual_and_jacobian(system,u;
-                                   t=0.0, tstep=Inf,embed=0.0)
+                                   t=0.0, tstep=Inf, embed=0.0, init_dirichlet=false)
 
 Evaluate residual and jacobian at solution value u.
 Returns a solution vector containing a copy of  residual, and an ExendableSparseMatrix
-containing a copy of the linearization at u.
+containing a copy of the linearization at u. The flag `init_dirichlet` controls whether
+u should be adjusted to satisfy the Dirichlet boundary conditions specified by the system.
 
 """
 function evaluate_residual_and_jacobian(sys, u; kwargs...)
     state = SystemState(sys)
-    eval_and_linearize!(state, u; kwargs...)
+    evaluate_residual_and_jacobian!(state, u; kwargs...)
     return copy(state.residual), copy(state.matrix)
 end
 

test/alltests.jl

@@ -7,11 +7,17 @@ ExampleJuggler.verbose!(true)
 # Include all Julia files in `testdir` whose name starts with `prefix`,
 # Each file `prefixModName.jl` must contain a module named
 # `prefixModName` which has a method test() returning true
-# or false depending on success.
+# or false depending on success. All files with filenames contained in
+# the `ignore` list will not be included.
 #
-function run_tests_from_directory(testdir, prefix)
+function run_tests_from_directory(testdir, prefix; ignore = [])
     @info "Directory $(testdir):"
-    examples = filter(ex -> length(ex) >= length(prefix) && ex[1:length(prefix)] == prefix, basename.(readdir(testdir)))
+    examples = filter(
+        filename -> length(filename) >= length(prefix) &&
+            filename[1:length(prefix)] == prefix &&
+            filename ∉ ignore,
+        basename.(readdir(testdir))
+    )
     @info examples
     @testmodules(testdir, examples)
     return nothing

test/test_cachesol.jl

@@ -0,0 +1,53 @@
+module test_cachesol
+using ExtendableGrids
+using VoronoiFVM
+using LinearSolve
+using ExtendableSparse
+using LinearAlgebra
+using SparseArrays
+using Test
+
+function flux!(f, u, edge, data)
+    f[1] = u[1, 1] - u[1, 2]
+    return nothing
+end
+
+function bcondition!(f, u, bnode, data)
+    boundary_dirichlet!(f, u, bnode; species = 1, region = 1, value = 0.0)
+    boundary_dirichlet!(f, u, bnode; species = 1, region = 2, value = 1.0)
+    return nothing
+end
+
+function setup_linear_solve(n)
+    # define system, initialize system matrix
+    X = linspace(0.0, 1.0, n)
+    system = VoronoiFVM.System(X; flux = flux!, bcondition = bcondition!, species = [1])
+    F = unknowns(system; inival = 0.0)
+    # assemble system matrix A and the residual F of A*U=RHS in U=0, i.e.
+    # F = A*0 - RHS = -RHS
+    F, A = VoronoiFVM.evaluate_residual_and_jacobian(system, F)
+
+    p = LinearProblem(A, -VoronoiFVM.dofs(F))
+    linear_cache = init(p, LinearSolve.UMFPACKFactorization())
+
+    return X, linear_cache
+end
+
+function solve_problem(cache)
+    sol = LinearSolve.solve!(cache)
+    return sol
+end
+
+function main(; n = 10)
+    X, cache = setup_linear_solve(n)
+    approx_sol = solve_problem(cache)
+    true_sol = map((x) -> (x), X)
+    @test norm(approx_sol .- true_sol, Inf) ≤ 3.0e-16
+    return nothing
+end
+
+function runtests()
+    main()
+    return nothing
+end
+end