some fixes, addded constructor for LinearOperatorFromMatrix, e.g. for…

… more efficient

custom time-dependent problems with fixed mass matrix

docs/src/pdesolvers.md

@@ -15,17 +15,25 @@ 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]
 Pages = ["solvers.jl"]
 Order   = [:type, :function]
 ```
 
-## Solve (time-dependent)
-
-work in progress
 
 ## Export and Postprocessing
 

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 += τ   

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

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