Merging 'master' into feature/extendablefembaseext

Project.toml

@@ -1,7 +1,7 @@
 name = "VoronoiFVM"
 uuid = "82b139dc-5afc-11e9-35da-9b9bdfd336f3"
 authors = ["Jürgen Fuhrmann <[email protected]>", "Patrick Jaap", "Dilara Abdel", "Jan Weidner", "Alexander Seiler", "Patricio Farrell", "Matthias Liero"]
-version = "2.0.0"
+version = "2.0.1"
 
 [deps]
 BandedMatrices = "aae01518-5342-5314-be14-df237901396f"

examples/Example440_ParallelState.jl

@@ -35,7 +35,7 @@ function main(;nref=5, Plotter=nothing)
     masses=similar(influxes)
 
     ## Split the index range in as many chunks as threads
-    Threads.@threads for indexes in chunks(influxes;n=Threads.nthreads())
+    Threads.@threads for indexes in chunks(1:length(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

src/vfvm_solver.jl

@@ -429,9 +429,15 @@ function solve_transient!(state,
                              Δλ * Δu_opt / (Δu + 1.0e-14),
                              λ_predict,
                              λend - λ)
+
+                    # adjust for roundoff error, so there will be no accidental
+                    # very tiny timestep
+                    if isapprox(λ+Δλ, λend; atol=1.0e-15, rtol=1.0e-15)
+                        Δλ=λend-λ
+                    end
                 end
             else
-                break # break out of inner loop overt timestep
+                break # break out of inner loop over timestep
             end # if solved            
         end # while λ<λ_end
 
@@ -480,9 +486,11 @@ function CommonSolve.solve!(state::VoronoiFVM.SystemState;
                             tstep = Inf,
                             kwargs...,)
     fix_deprecations!(control)
-
+    
     if isa(inival, Number) || isa(inival, Matrix) 
         inival = unknowns(state.system; inival = inival)
+    elseif !isdensesystem(state.system) && isa(inival, SparseMatrixCSC)
+        inival = SparseSolutionArray(inival)
     elseif !VoronoiFVM.isunknownsof(inival, state.system)
         @error "wrong type of inival: $(typeof(inival))"
     end

src/vfvm_sparsesolution.jl

@@ -54,6 +54,7 @@ $(SIGNATURES)
 Array of values in sparse solution array.
 """
 values(a::SparseSolutionArray) = a.u.nzval
+values(a::SparseMatrixCSC) = a.nzval
 
 Base.size(a::SparseSolutionArray)=size(a.u)
 ##################################################################

test/test_linsolve.jl

@@ -61,86 +61,59 @@ end
 
 checklum(n, T) = checklum0(Val{n}, T)
 
-@static if VERSION >= v"1.9-rc0"
-    function checklurm0(::Type{Val{N}}, ::Type{T}) where {N, T}
-        A = MMatrix{N, N, Float64}(undef)
-        x = MVector{N, Float64}(undef)
-        b = MVector{N, Float64}(undef)
-        ipiv = MVector{N, Int64}(undef)
-        for i = 1:N
-            for j = 1:N
-                A[i, j] = -rand()
-            end
-            A[i, i] += 100
-            x[i] = 1
-        end
-        @gc_preserve mul!(b, A, x)
-        @gc_preserve inplace_linsolve!(A, b, ipiv)
-
-        nm = 0
-        for i = 1:N
-            nm += (b[i] - x[i])^2
-        end
-
-        @assert sqrt(nm) / N < 100.0 * eps(Float64)
-        nothing
-    end
 
-    checklurm(n, T) = checklurm0(Val{n}, T)
-
-    function checklurx0(::Type{Val{N}}, ::Type{T}) where {N, T}
-        A = StrideArray{T}(undef, StaticInt(N), StaticInt(N))
-        x = StrideArray{T}(undef, StaticInt(N))
-        b = StrideArray{T}(undef, StaticInt(N))
-        ipiv = StrideArray{Int64}(undef, StaticInt(N))
-        for i = 1:N
-            for j = 1:N
-                A[i, j] = -rand()
+function checklurx0(::Type{Val{N}}, ::Type{T}) where {N, T}
+    A = StrideArray{T}(undef, StaticInt(N), StaticInt(N))
+    x = StrideArray{T}(undef, StaticInt(N))
+    b = StrideArray{T}(undef, StaticInt(N))
+    ipiv = StrideArray{Int64}(undef, StaticInt(N))
+    for i = 1:N
+        for j = 1:N
+            A[i, j] = -rand()
             end
-            A[i, i] += 100
-            x[i] = 1
-        end
-        @gc_preserve mul!(b, A, x)
-        @gc_preserve inplace_linsolve!(A, b, ipiv)
-
+        A[i, i] += 100
+        x[i] = 1
+    end
+    @gc_preserve mul!(b, A, x)
+    @gc_preserve inplace_linsolve!(A, b, ipiv)
+    
         nm = 0
-        for i = 1:N
-            nm += (b[i] - x[i])^2
-        end
-
-        @assert sqrt(nm) / N < 100.0 * eps(Float64)
-        nothing
+    for i = 1:N
+        nm += (b[i] - x[i])^2
     end
 
-    checklurx(n, T) = checklurx0(Val{n}, T)
+    @assert sqrt(nm) / N < 100.0 * eps(Float64)
+    nothing
+end
 
-    function checklurm0(::Type{Val{N}}, ::Type{T}) where {N, T}
-        A = MMatrix{N, N, Float64}(undef)
-        x = MVector{N, Float64}(undef)
-        b = MVector{N, Float64}(undef)
-        ipiv = MVector{N, Int64}(undef)
-        for i = 1:N
-            for j = 1:N
-                A[i, j] = -rand()
-            end
-            A[i, i] += 100
-            x[i] = 1
-        end
-        @gc_preserve mul!(b, A, x)
-        @gc_preserve inplace_linsolve!(A, b, ipiv)
+checklurx(n, T) = checklurx0(Val{n}, T)
 
-        nm = 0
-        for i = 1:N
-            nm += (b[i] - x[i])^2
+function checklurm0(::Type{Val{N}}, ::Type{T}) where {N, T}
+    A = MMatrix{N, N, Float64}(undef)
+    x = MVector{N, Float64}(undef)
+    b = MVector{N, Float64}(undef)
+    ipiv = MVector{N, Int64}(undef)
+    for i = 1:N
+        for j = 1:N
+            A[i, j] = -rand()
         end
-
-        @assert sqrt(nm) / N < 100.0 * eps(Float64)
-        nothing
+        A[i, i] += 100
+        x[i] = 1
     end
-
-    checklurm(n, T) = checklurm0(Val{n}, T)
+    @gc_preserve mul!(b, A, x)
+    @gc_preserve inplace_linsolve!(A, b, ipiv)
+
+    nm = 0
+    for i = 1:N
+        nm += (b[i] - x[i])^2
+    end
+
+    @assert sqrt(nm) / N < 100.0 * eps(Float64)
+    nothing
 end
 
+checklurm(n, T) = checklurm0(Val{n}, T)
+
 function runtests()
     checklum(10, Float64)
     checklum(10, Dual64)

test/test_state.jl

@@ -0,0 +1,52 @@
+module test_state
+using VoronoiFVM
+using VoronoiFVM: SystemState
+using ExtendableGrids
+using LinearAlgebra
+using Test
+
+flux(y,u,edge, data)= y[1]=u[1,1]-u[1,2]
+
+function bcondition(y,u,bnode,data)
+    boundary_robin!(y,u,bnode,region=1, value=0.0, factor=0.1)
+    boundary_robin!(y,u,bnode,region=2, value=1.0, factor=0.1)
+end
+
+function main(; unknown_storage=:dense)
+    g=simplexgrid(0:0.1:1)
+    sys=VoronoiFVM.System(g; flux, bcondition, species=[1], unknown_storage)
+    sol1=solve(sys)
+
+    # Solution and solution with state shall be the same
+    state=VoronoiFVM.SystemState(sys)
+    sol2=solve!(state)
+    @test sol1≈sol2
+
+    control=SolverControl()
+    fixed_timesteps!(control,0.025)
+
+
+    # Allow initial values as result of previous time evolution
+    tsol1=solve(sys; inival=0.0, times=(0,0.1), control)
+    tsol2=solve(sys; inival=tsol1.u[end], times=(0.1,0.2), control)
+    tsol3=solve(sys; inival=0.0, times=[0.0,0.1,0.2], control)
+    @test tsol3.u[end]≈tsol2.u[end]
+
+    # Solution and solution with state shall be the same
+    xsol1=solve!(state; inival=0.0, times=(0,0.1), control)
+    xsol2=solve!(state; inival=tsol1.u[end], times=(0.1,0.2), control)
+    xsol3=solve!(state; inival=0.0, times=[0.0,0.1,0.2], control)
+    @test xsol3.u[end]≈xsol2.u[end]
+    @test xsol3.u[end]≈tsol3.u[end]
+
+
+
+end
+
+function runtests()
+    main(;unknown_storage=:dense)
+    main(;unknown_storage=:sparse)
+
+end
+end
+