diff --git a/Project.toml b/Project.toml index 078c16550..b819deed9 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "VoronoiFVM" uuid = "82b139dc-5afc-11e9-35da-9b9bdfd336f3" authors = ["Jürgen Fuhrmann ", "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" diff --git a/examples/Example440_ParallelState.jl b/examples/Example440_ParallelState.jl index c07877e39..274e45dc1 100644 --- a/examples/Example440_ParallelState.jl +++ b/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 diff --git a/src/vfvm_solver.jl b/src/vfvm_solver.jl index ff5cfc351..3d2b03bae 100644 --- a/src/vfvm_solver.jl +++ b/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 diff --git a/src/vfvm_sparsesolution.jl b/src/vfvm_sparsesolution.jl index e474dbc1e..cc51d8003 100644 --- a/src/vfvm_sparsesolution.jl +++ b/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) ################################################################## diff --git a/test/test_linsolve.jl b/test/test_linsolve.jl index fb57aedf3..16479e907 100644 --- a/test/test_linsolve.jl +++ b/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) diff --git a/test/test_state.jl b/test/test_state.jl new file mode 100644 index 000000000..06e17cbe6 --- /dev/null +++ b/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 +