From 5c2cb4967ceb6ef9c1571edbe456267c2544d1ea Mon Sep 17 00:00:00 2001 From: Christian Merdon Date: Mon, 14 Aug 2023 17:13:21 +0200 Subject: [PATCH] some fixes, addded constructor for LinearOperatorFromMatrix, e.g. for more efficient custom time-dependent problems with fixed mass matrix --- docs/src/pdesolvers.md | 16 ++++++-- examples/Example103_BurgersEquation.jl | 15 +++---- src/common_operators/bilinear_operator.jl | 2 +- src/common_operators/linear_operator.jl | 48 +++++++++++++++++++++++ 4 files changed, 69 insertions(+), 12 deletions(-) diff --git a/docs/src/pdesolvers.md b/docs/src/pdesolvers.md index 4b364cca..b60fbed4 100644 --- a/docs/src/pdesolvers.md +++ b/docs/src/pdesolvers.md @@ -15,7 +15,18 @@ easily generated by FESpace{FEType}(grid::ExtendableGrid) where FEType denotes t A list of available FEType can be found in the documentation of ExtendableFEMBase. ## Solve (stationary) -If solve is applied to a PDEDescription and an array of FESpaces (that specify the ansatz spaces for the unknowns) an investigation of the PDEDescription is performed that decides if the problem is nonlinear (and has to be solved by a fixed-point algorithm) or if it can be solved directly in one step. +If solve is applied to a PDEDescription and an array of FESpaces (that specify the ansatz spaces for the unknowns) an investigation of the PDEDescription is performed that decides if the problem is nonlinear (and has to be solved by a fixed-point algorithm) or if it can be solved directly in one step. The solver attempts to bring the (nonlinear) residual to the prescribed target tolerance. In a time-dependent context this +is equivalent to seek a stationary solution and probably is not always what one wants. + +## Solve (time-dependent) + +For time-dependent (non-stationary) problems the user currently has these options: + +- add a pre-defined time derivative operator (i.e. BackwardDifference) and iterate the problem yourself +- add custom time derivatives to the problem (i.e. a mass matrix as a BilinearOperator and necessary LinearOperators for evaluating the previous time step(s), if more than one previous time step needs to be remembered, further unknowns must be registered) +- reframe the ProblemDescription as an ODE problem and evolve it via DifferentialEquations + +## Detailed headers ```@autodocs Modules = [ExtendableFEM] @@ -23,9 +34,6 @@ Pages = ["solvers.jl"] Order = [:type, :function] ``` -## Solve (time-dependent) - -work in progress ## Export and Postprocessing diff --git a/examples/Example103_BurgersEquation.jl b/examples/Example103_BurgersEquation.jl index 5c0cdcfc..efc7ebb4 100644 --- a/examples/Example103_BurgersEquation.jl +++ b/examples/Example103_BurgersEquation.jl @@ -65,11 +65,11 @@ function main(; p = GridVisualizer(; Plotter = Plotter, layout = (1,1), clear = true, resolution = (800,800)) scalarplot!(p[1,1], xgrid, nodevalues_view(sol[u])[1], flimits = (-0.75,2), levels = 0, title = "u_h (t = 0)") - if (use_diffeq) - ## generate mass matrix - M = FEMatrix(FES) - assemble!(M, BilinearOperator([id(1)])) + ## generate mass matrix + M = FEMatrix(FES) + assemble!(M, BilinearOperator([id(1)])) + if (use_diffeq) ## generate ODE problem prob = generate_ode(DifferentialEquations, SC, (0.0, T); mass_matrix = M.entries.cscmatrix) @@ -80,10 +80,11 @@ function main(; ## get final solution sol.entries .= de_sol[end] else - ## solve time-dependent problem (implicit Euler) - assign_operator!(PD, BilinearOperator([id(u)]; factor = 1/τ, store = true, kwargs...)) - assign_operator!(PD, LinearOperator([id(u)], [id(u)]; factor = 1/τ, kwargs...)) + ## add backward Euler time derivative + assign_operator!(PD, BilinearOperator(M, [u]; factor = 1/τ, kwargs...)) + assign_operator!(PD, LinearOperator(M, [u], [u]; factor = 1/τ, kwargs...)) + ## iterate tspan t = 0 for it = 1 : Int(floor(T/τ)) t += τ diff --git a/src/common_operators/bilinear_operator.jl b/src/common_operators/bilinear_operator.jl index 5db6bd94..60263786 100644 --- a/src/common_operators/bilinear_operator.jl +++ b/src/common_operators/bilinear_operator.jl @@ -851,7 +851,7 @@ function ExtendableFEM.assemble!(A, b, sol, O::BilinearOperatorFromMatrix{UT,MT} O.A = O.parameters[:callback!](O.A, [b[j] for j in ind_test], sol) end if MT <: FEMatrix - for (j,ij) in enumrate(ind_test), (k,ik) in enumerate(ind_ansatz) + for (j,ij) in enumerate(ind_test), (k,ik) in enumerate(ind_ansatz) addblock!(A[j,k], O.A[ij,ik]; factor = O.parameters[:factor]) end else diff --git a/src/common_operators/linear_operator.jl b/src/common_operators/linear_operator.jl index 9ea5bd1c..7634247a 100644 --- a/src/common_operators/linear_operator.jl +++ b/src/common_operators/linear_operator.jl @@ -4,6 +4,13 @@ mutable struct LinearOperatorFromVector{UT <: Union{Unknown, Integer}, bT} <: Ab parameters::Dict{Symbol,Any} end +mutable struct LinearOperatorFromMatrix{UT <: Union{Unknown, Integer}, MT} <: AbstractOperator + u_test::Array{UT,1} + u_args::Array{UT,1} + A::MT + parameters::Dict{Symbol,Any} +end + mutable struct LinearOperator{Tv <: Real, UT <: Union{Unknown, Integer}, KFT <: Function, ST} <: AbstractOperator u_test::Array{UT,1} ops_test::Array{DataType,1} @@ -144,6 +151,27 @@ function LinearOperator(b, u_test; kwargs...) return LinearOperatorFromVector{typeof(u_test[1]), typeof(b)}(u_test, b, parameters) end +""" +```` +function LinearOperator( + A, + u_test, + u_args; + kwargs...) +```` + +Generates a linear form from a user-provided matrix A, which can be an AbstractMatrix or a FEMatrix with +multiple blocks. The arguments u_args specify which coefficients of the current solution +should be multiplied with the matrix and u_test specifies where to put the +(blocks of the) resulting vector in the system right-hand side. + +""" +function LinearOperator(A, u_test, u_args; kwargs...) + parameters=Dict{Symbol,Any}( k => v[1] for (k,v) in default_linop_kwargs()) + _update_params!(parameters, kwargs) + return LinearOperatorFromMatrix{typeof(u_test[1]), typeof(A)}(u_test, u_args, A, parameters) +end + """ ```` @@ -586,4 +614,24 @@ function ExtendableFEM.assemble!(A, b, sol, O::LinearOperatorFromVector{UT,bT}, @assert length(ind_test) == 1 addblock!(b[ind_test[1]], O.b; factor = O.parameters[:factor]) end +end + + + +function ExtendableFEM.assemble!(A, b, sol, O::LinearOperatorFromMatrix{UT,MT}, SC::SolverConfiguration; kwargs...) where {UT,MT} + if UT <: Integer + ind_test = O.u_test + ind_args = O.u_args + elseif UT <: Unknown + ind_test = [get_unknown_id(SC, u) for u in O.u_test] + ind_args = [get_unknown_id(SC, u) for u in O.u_args] + end + if MT <: FEMatrix + for (j,ij) in enumerate(ind_test), k in ind_args + addblock_matmul!(b[j], O.A[ij, k], sol[k]; factor = O.parameters[:factor]) + end + else + @assert length(ind_test) == 1 && length(ind_args) == 1 + addblock_matmul!(b[ind_test[1]], O.A, sol[ind_args[1]]; factor = O.parameters[:factor]) + end end \ No newline at end of file