Merge pull request #138 from WIAS-PDELib/jf/fix-impedance-timeembed

Allow to use result of time embedding a steady state for impedance
WIAS-PDELib · Oct 30, 2024 · 68e11a5 · 68e11a5 · j-fu · Oct 30, 2024
2 parents 6dffb18 + 9bde30e
commit 68e11a5
Show file tree

Hide file tree

Showing 12 changed files with 92 additions and 26 deletions.
diff --git a/CHANGELOG.md b/CHANGELOG.md
@@ -1,12 +1,20 @@
 # Changes
-## 2024-10-28
-
-Moved repositiory from https://github.com/j-fu/VoronoiFVM.jl to https://github.com/WIAS-PDELib/VoronoiFVM.jl. [WIAS-PDELib](https://github.com/WIAS-PDELib/) is a github organization created to collectively manage the Julia packages developed under the lead of the [WIAS Numerical Mathematics and Scientific Computing](https://wias-berlin.de/research/rgs/fg3)  research group. According to the [github docs on repository transfer](https://docs.github.com/en/repositories/creating-and-managing-repositories/transferring-a-repository#whats-transferred-with-a-repository), all links to the previous repository location are automatically redirected to the new location, and all relationships with forks stay intact.
 
 ## v3.0 Planned  (pending some improvements in LinearSolve)
   - use `precs` based linear solver API, see https://github.com/SciML/LinearSolve.jl/pull/514
   - stop re-exporting ForwardDiff.value
   - try to remove type piracies
+  - remove `params` from edge, node structs (appearantly never used)
+
+## v2.2.0 October 30, 2024
+  - Add `params` to SystemState, allow to pass params to ODEProblem
+  - Fix use of end results of time evolution as steady state for impedance calculations
+  - new internal solutionarray method. 
+
+## October 28, 2024
+
+Moved repositiory from https://github.com/j-fu/VoronoiFVM.jl to https://github.com/WIAS-PDELib/VoronoiFVM.jl. [WIAS-PDELib](https://github.com/WIAS-PDELib/) is a github organization created to collectively manage the Julia packages developed under the lead of the [WIAS Numerical Mathematics and Scientific Computing](https://wias-berlin.de/research/rgs/fg3)  research group. According to the [github docs on repository transfer](https://docs.github.com/en/repositories/creating-and-managing-repositories/transferring-a-repository#whats-transferred-with-a-repository), all links to the previous repository location are automatically redirected to the new location, and all relationships with forks stay intact.
+
 ## v2.1.0 October 23, 2024
   - Features:
     - Add VoronoiFVMExtendableFEMBaseExt  extension allowing to use velocity fields caculated  with [ExtendableFEM](https://github.com/chmerdon/ExtendableFEM.jl) in convection-diffusion problems

diff --git a/Project.toml b/Project.toml
@@ -1,7 +1,7 @@
 name = "VoronoiFVM"
 uuid = "82b139dc-5afc-11e9-35da-9b9bdfd336f3"
 authors = ["Jürgen Fuhrmann <[email protected]>", "Patrick Jaap", "Daniel Runge", "Dilara Abdel", "Jan Weidner", "Alexander Seiler", "Patricio Farrell", "Matthias Liero"]
-version = "2.1"
+version = "2.2"
 
 [deps]
 BandedMatrices = "aae01518-5342-5314-be14-df237901396f"

diff --git a/docs/src/internal.md b/docs/src/internal.md
@@ -121,6 +121,7 @@ prepare_diffeq!
 
 ## Misc tools
 ```@docs
+VoronoiFVM.solutionarray
 VoronoiFVM.integrate(::Type{<:Cartesian2D}, coordl, coordr, hnormal, velofunc; kwargs...)
 VoronoiFVM.integrate(::Type{<:Cylindrical2D}, coordl, coordr, hnormal, velofunc; kwargs...)
 VoronoiFVM.doolittle_ludecomp!

diff --git a/examples/Example151_Impedance1D.jl b/examples/Example151_Impedance1D.jl
@@ -34,8 +34,11 @@ using Printf
 using VoronoiFVM
 using ExtendableGrids: geomspace, simplexgrid
 using GridVisualize
+using OrdinaryDiffEqSDIRK
 
-function main(; nref = 0, Plotter = nothing, verbose = false, unknown_storage = :sparse, assembly = :edgewise,
+function main(; nref = 0, Plotter = nothing, verbose = false,
+              unknown_storage = :sparse, assembly = :edgewise,
+              time_embedding = :none,
               L = 1.0, R = 1.0, D = 1.0, C = 1.0,
               ω0 = 1.0e-3, ω1 = 5.0e1)
 
@@ -64,7 +67,7 @@ function main(; nref = 0, Plotter = nothing, verbose = false, unknown_storage =
         f[1] = data.R * u[1]
     end
 
-    excited_bc = 1
+    excited_bc= 1
     excited_bcval = 1.0
     excited_spec = 1
     meas_bc = 2
@@ -90,7 +93,21 @@ function main(; nref = 0, Plotter = nothing, verbose = false, unknown_storage =
     factory = TestFunctionFactory(sys)
     measurement_testfunction = testfunction(factory, [excited_bc], [meas_bc])
 
-    steadystate = solve(sys; inival = 0.0, params = [1.0])
+    tend=1.0
+    if time_embedding == :builtin
+        tsol=solve(sys; inival = 0.0, params = [1.0], times=(0.0,tend),force_first_step=true)
+        steadystate=tsol.u[end]
+    elseif time_embedding == :ordinarydiffeq
+        inival=unknowns(sys,inival=0)
+        problem = ODEProblem(sys,inival,(0,tend); params = [1.0])
+        odesol = solve(problem,ImplicitEuler())
+        tsol=reshape(odesol,sys)
+        steadystate=tsol.u[end]
+    elseif time_embedding==:none
+        steadystate = solve(sys; inival = 0.0, params = [1.0])
+    else
+        error("time_embedding must be one of :builtin, :ordinarydiffeq, :none")
+    end
 
     function meas_stdy(meas, U)
         u = reshape(U, sys)
@@ -150,10 +167,10 @@ function main(; nref = 0, Plotter = nothing, verbose = false, unknown_storage =
         ω = ω * 1.1
     end
 
-    p = GridVisualizer(; Plotter = Plotter)
-    scalarplot!(p, real(allIxL), imag(allIxL); label = "exact", color = :red,
+    vis = GridVisualizer(; Plotter = Plotter)
+    scalarplot!(vis, real(allIxL), imag(allIxL); label = "exact", color = :red,
                 linestyle = :dot)
-    scalarplot!(p, real(allIL), imag(allIL); label = "calc", show = true, clear = false,
+    scalarplot!(vis, real(allIL), imag(allIL); label = "calc", show = true, clear = false,
                 color = :blue, linestyle = :solid)
 
     sum(allIL)
@@ -162,10 +179,13 @@ end
 using Test
 function runtests()
     testval = 57.92710286186797 + 23.163945443946027im
-    @test main(; unknown_storage = :sparse, assembly = :edgewise) ≈ testval &&
-          main(; unknown_storage = :dense, assembly = :edgewise) ≈ testval &&
-          main(; unknown_storage = :sparse, assembly = :cellwise) ≈ testval &&
-          main(; unknown_storage = :dense, assembly = :cellwise) ≈ testval
+    for unknown_storage in (:sparse, :dense)
+        for assembly in (:edgewise, :cellwise)
+            for time_embedding  in (:none, :builtin, :ordinarydiffeq)
+                @test main(; unknown_storage, assembly, time_embedding) ≈ testval
+            end
+        end
+    end
 end
 
 end
diff --git a/src/VoronoiFVM.jl b/src/VoronoiFVM.jl
@@ -72,6 +72,7 @@ abstract type AbstractSolutionArray{T,N} <: AbstractArray{T,N} end
 Base.getindex(a::AbstractSolutionArray, i::Int, j::Int)= getindex(a.u,i,j )
 Base.setindex!(a::AbstractSolutionArray,v, i::Int, j::Int) = setindex!(a.u,v,i,j)
 Base.size(a::AbstractSolutionArray)=size(a.u)
+solutionarray(a::AbstractSolutionArray)=a
 
 export AbstractSolutionArray
 

diff --git a/src/vfvm_densesolution.jl b/src/vfvm_densesolution.jl
@@ -21,6 +21,11 @@ end
 """
 DenseSolutionArray(u::Matrix{T}) where {T} =DenseSolutionArray{T,2}(u,nothing)
 
+"""
+    solutionarray(a::Matrix)
+"""
+solutionarray(u::Matrix{T}) where {T} =DenseSolutionArray{T,2}(u,nothing)
+
 """
  $(TYPEDSIGNATURES)
     

diff --git a/src/vfvm_diffeq_interface.jl b/src/vfvm_diffeq_interface.jl
@@ -9,7 +9,7 @@ function _eval_res_jac!(state, u, t)
     uhash = hash(u)
     if uhash != state.uhash
         ur = reshape(u, state.system)
-        eval_and_assemble(state.system, ur, ur, state.residual, state.matrix, state.dudp, value(t), Inf, 0.0, state.system.physics.data, zeros(0))
+        eval_and_assemble(state.system, ur, ur, state.residual, state.matrix, state.dudp, value(t), Inf, 0.0, state.system.physics.data, state.params)
         state.uhash = uhash
         state.history.nd += 1
     end
@@ -151,7 +151,9 @@ for more documentation.
 
 Defined in VoronoiFVM.jl.
 """
-function SciMLBase.ODEProblem(state::VoronoiFVM.SystemState, inival, tspan; callback = SciMLBase.CallbackSet())
+function SciMLBase.ODEProblem(state::VoronoiFVM.SystemState, inival, tspan;
+                              params=state.params, callback = SciMLBase.CallbackSet())
+    state.params.=params
     odefunction = SciMLBase.ODEFunction(state; jacval = dofs(inival), tjac = tspan[1])
     SciMLBase.ODEProblem(odefunction, dofs(inival), tspan, state, callback)
 end
@@ -174,7 +176,10 @@ by [solve()](https://diffeq.sciml.ai/stable/basics/common_solver_opts/).
 
 Defined in VoronoiFVM.jl.
 """
-SciMLBase.ODEProblem(sys::VoronoiFVM.System, inival, tspan;kwargs...)=SciMLBase.ODEProblem(SystemState(sys), inival, tspan; kwargs...)
+function SciMLBase.ODEProblem(sys::VoronoiFVM.System, inival, tspan;
+                              params=zeros(sys.num_parameters), kwargs...)
+    SciMLBase.ODEProblem(SystemState(sys), inival, tspan; params, kwargs...)
+end
 
 """
     reshape(ode_solution, system; times=nothing, state=nothing)

diff --git a/src/vfvm_geometryitems.jl b/src/vfvm_geometryitems.jl
@@ -55,7 +55,7 @@ Base.getindex(p::DParameters, i) = @inbounds p.val[p.offset + i]
 Return abstract vector of parameters passed via vector of unknonws. 
 This allows differentiation with respect to these parameters.
 """
-function parameters(u::AbstractNodeData{Tv}) where {Tv <: Number}
+function parameters(u::AbstractNodeData)
     DParameters(u.val, u.nspec)
 end
 
@@ -80,7 +80,7 @@ abstract type AbstractEdgeData{Tv <: Number} <: AbstractMatrix{Tv} end
 Base.size(u::AbstractEdgeData) = (u.n1, 2)
 Base.getindex(u::AbstractEdgeData, i, j) = @inbounds u.val[(j - 1) * u.n1 + i]
 
-function parameters(u::AbstractEdgeData{Tv}) where {Tv <: Number}
+function parameters(u::AbstractEdgeData)
     DParameters(u.val, u.n1 + u.n1)
 end
 
@@ -139,7 +139,7 @@ mutable struct Node{Tc, Tp, Ti} <: AbstractNode{Tc, Tp, Ti}
     embedparam::Float64
 
     """
-    parameters
+    parameters (deprecated)
     """
     params::Vector{Tp}
 
@@ -258,6 +258,9 @@ mutable struct BNode{Tv, Tc, Tp, Ti} <: AbstractNode{Tc, Tp, Ti}
     """
     embedparam::Float64
 
+    """
+    Parameters (deprecated)
+    """
     params::Vector{Tp}
 
     dirichlet_value::Vector{Tv}
@@ -356,6 +359,9 @@ mutable struct Edge{Tc, Tp, Ti} <: AbstractEdge{Tc, Tp, Ti}
     """
     embedparam::Float64
 
+    """
+    Parameters (deprecated)
+    """
     params::Vector{Tp}
 
     """
@@ -519,6 +525,9 @@ mutable struct BEdge{Tc, Tp, Ti} <: AbstractEdge{Tc, Tp, Ti}
     """
     embedparam::Float64
 
+    """
+    Parameters (deprecated)
+    """
     params::Vector{Tp}
 
     fac::Float64

diff --git a/src/vfvm_impedance.jl b/src/vfvm_impedance.jl
@@ -47,6 +47,7 @@ $(SIGNATURES)
 Construct impedance system from time domain system `sys` and steady state solution `U0`
 """
 function ImpedanceSystem(system::AbstractSystem{Tv, Tc, Ti}, U0::AbstractMatrix; λ0 = 0.0) where {Tv, Tc, Ti}
+    U0=solutionarray(U0)
     residual = copy(U0)
 
     if system.num_parameters > 0

diff --git a/src/vfvm_solver.jl b/src/vfvm_solver.jl
@@ -486,7 +486,6 @@ 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)
@@ -506,6 +505,8 @@ function CommonSolve.solve!(state::VoronoiFVM.SystemState;
     if !isnothing(data)
         state.data=data
     end
+    state.params.=params
+
 
     if haskey(kwargs, :times) && !isnothing(kwargs[:times])
         solve_transient!(state,

diff --git a/src/vfvm_sparsesolution.jl b/src/vfvm_sparsesolution.jl
@@ -28,6 +28,10 @@ end
 """
 SparseSolutionArray(a::SparseMatrixCSC{Tv,Ti}) where {Tv,Ti}=SparseSolutionArray{Tv,2,Ti}(a,nothing)
 
+"""
+    solutionarray(a::SparseMatrixCSC)
+"""
+solutionarray(a::SparseMatrixCSC{Tv,Ti}) where {Tv,Ti}=SparseSolutionArray{Tv,2,Ti}(a,nothing)
 
 """
 $(TYPEDEF)

diff --git a/src/vfvm_state.jl b/src/vfvm_state.jl
@@ -13,15 +13,15 @@ Type parameters:
 Type fields:
 $(TYPEDFIELDS)
 """
-mutable struct SystemState{Tv, TMatrix<:AbstractMatrix{Tv}, TSolArray<:AbstractMatrix{Tv}, TData}
+mutable struct SystemState{Tv, Tp, TMatrix<:AbstractMatrix{Tv}, TSolArray<:AbstractMatrix{Tv}, TData}
 
     """
     Related finite volume system
     """
     system::VoronoiFVM.System
 
     """
-    Parameter data 
+    User data 
     """
     data::TData
 
@@ -55,6 +55,11 @@ mutable struct SystemState{Tv, TMatrix<:AbstractMatrix{Tv}, TSolArray<:AbstractM
     """
     linear_cache::Union{Nothing, LinearSolve.LinearCache}
 
+    """
+    Parameter vector
+    """
+    params::Vector{Tp}
+
     """
     Hash value of latest unknowns vector the assembly was called with
     (used by differential equation interface)
@@ -85,9 +90,14 @@ Keyword arguments:
 """
 function SystemState(::Type{Tu}, system::AbstractSystem{Tv, Tc, Ti, Tm};
                      data=system.physics.data,
+                     params=zeros(system.num_parameters),
                      matrixtype=system.matrixtype) where {Tu,Tv,Tc, Ti, Tm}
     _complete!(system)
-
+
+    if (length(params)!=system.num_parameters)
+        error("length(params)!=system.num_parameters")
+    end
+
     nspec = size(system.node_dof, 1)
     n = num_dof(system)
 
@@ -123,7 +133,7 @@ function SystemState(::Type{Tu}, system::AbstractSystem{Tv, Tc, Ti, Tm};
     residual = unknowns(Tu, system)
     update = unknowns(Tu, system)
     dudp = [unknowns(Tu, system) for i = 1:(system.num_parameters)]
-    SystemState(system, data, solution, matrix, dudp, residual, update, nothing, zero(UInt64),  nothing)
+    SystemState(system, data, solution, matrix, dudp, residual, update, nothing, params, zero(UInt64),  nothing)
 end
 
 
@@ -154,7 +164,8 @@ function Base.similar(state::SystemState; data=state.data)
     residual=similar(state.residual)
     update=similar(state.update)
     linear_cache=nothing
+    params=similar(state.params)
     uhash=zero(UInt64)
     history=nothing
-    SystemState(system, data, solution, matrix, dudp, residual, update, linear_cache, uhash, history)
+    SystemState(system, data, solution, matrix, dudp, residual, update, linear_cache, params, uhash, history)
 end