diff --git a/NEWS.md b/NEWS.md index 5c344a2a2..f5a5d1c75 100644 --- a/NEWS.md +++ b/NEWS.md @@ -11,13 +11,30 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 - Implemented real/imag for VectorValues +### Fixed + +- `BlockMultiFieldStyle` available for `TransientMultiFieldFESpaces` since PR [#946](/~https://github.com/gridap/Gridap.jl/pull/946). + +## [0.17.20] - 2023-10-01 + +### Added + +- Block assembly now generalised to work with `AbstractBlockArrays`, to include changes in GridapDistributed. Since PR [939](/~https://github.com/gridap/Gridap.jl/pull/939). +- Implici-Explicit Runge-Kutta ODE solvers. Since PR [#919](/~https://github.com/gridap/Gridap.jl/pull/919). + +### Fixed + +- Using Broadcasting(\circ) instead of \circ in one of the lazy_maps used to transform a coarse field into a fine field. Since PR [#938](/~https://github.com/gridap/Gridap.jl/pull/938). +- Better infinite norm computation in `Algebra._check_convergence`. Now works for any `AbstractArray` type, including `PVector`. Since PR [#940](/~https://github.com/gridap/Gridap.jl/pull/940). +- Updated Runge-Kutta solver. Since PR [#919](/~https://github.com/gridap/Gridap.jl/pull/919). + ## [0.17.19] - 2023-08-23 ### Fixed - Reimplemented `DomainStyle` for `CellQuadrature` to fix breaking low-level Poisson tutorial. Since PR [#937](/~https://github.com/gridap/Gridap.jl/pull/937). -## [0.17.18] - 2023-08-15 +## [0.17.18] - 2023-08-15 ### Added diff --git a/Project.toml b/Project.toml index 9e89d5dba..53306bba5 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "Gridap" uuid = "56d4f2e9-7ea1-5844-9cf6-b9c51ca7ce8e" authors = ["Santiago Badia ", "Francesc Verdugo ", "Alberto F. Martin "] -version = "0.17.19" +version = "0.17.20" [deps] AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c" diff --git a/src/Adaptivity/AdaptedTriangulations.jl b/src/Adaptivity/AdaptedTriangulations.jl index 733c8760c..03ce692a1 100644 --- a/src/Adaptivity/AdaptedTriangulations.jl +++ b/src/Adaptivity/AdaptedTriangulations.jl @@ -278,7 +278,7 @@ function change_domain_o2n(f_coarse,ctrian::Triangulation{Dc},ftrian::AdaptedTri ref_coord_map = get_n2o_reference_coordinate_map(glue) # f_fine(x_fine) = f_f2c ∘ Φ^(-1)(x_fine) = f_coarse(x_coarse) - fine_mface_to_field = lazy_map(∘,f_f2c,ref_coord_map) + fine_mface_to_field = lazy_map(Broadcasting(∘),f_f2c,ref_coord_map) ### New Model -> New Triangulation fine_tface_to_field = lazy_map(Reindex(fine_mface_to_field),fglue.tface_to_mface) diff --git a/src/Algebra/NLSolvers.jl b/src/Algebra/NLSolvers.jl index c83a044ea..3990e5543 100644 --- a/src/Algebra/NLSolvers.jl +++ b/src/Algebra/NLSolvers.jl @@ -77,21 +77,13 @@ function _solve_nr!(x,A,b,dx,ns,nls,op) end function _check_convergence(nls,b) - m0 = _inf_norm(b) - (false, m0) + m0 = maximum(abs,b) + return (false, m0) end function _check_convergence(nls,b,m0) - m = _inf_norm(b) - m < nls.tol * m0 -end - -function _inf_norm(b) - m = 0 - for bi in b - m = max(m,abs(bi)) - end - m + m = maximum(abs,b) + return m < nls.tol * m0 end # Default concrete implementation diff --git a/src/MultiField/BlockSparseMatrixAssemblers.jl b/src/MultiField/BlockSparseMatrixAssemblers.jl index 8430bf54d..b4a336d74 100644 --- a/src/MultiField/BlockSparseMatrixAssemblers.jl +++ b/src/MultiField/BlockSparseMatrixAssemblers.jl @@ -80,8 +80,9 @@ function BlockSparseMatrixAssembler(trial::MultiFieldFESpace, @notimplemented msg end -function BlockSparseMatrixAssembler(trial::MultiFieldFESpace{<:BlockMultiFieldStyle{NB,SB,P}}, - test::MultiFieldFESpace{<:BlockMultiFieldStyle{NB,SB,P}}, +function BlockSparseMatrixAssembler(::BlockMultiFieldStyle{NB,SB,P}, + trial, + test, matrix_builder, vector_builder, strategy=FESpaces.DefaultAssemblyStrategy()) where {NB,SB,P} @@ -103,12 +104,13 @@ function BlockSparseMatrixAssembler(trial::MultiFieldFESpace{<:BlockMultiFieldSt return BlockSparseMatrixAssembler{NB,NV,SB,P}(block_assemblers) end -function FESpaces.SparseMatrixAssembler(mat, - vec, - trial::MultiFieldFESpace{<:BlockMultiFieldStyle}, - test::MultiFieldFESpace{<:BlockMultiFieldStyle}, - strategy::AssemblyStrategy=DefaultAssemblyStrategy()) - return BlockSparseMatrixAssembler(trial,test,SparseMatrixBuilder(mat),ArrayBuilder(vec),strategy) +function FESpaces.SparseMatrixAssembler(mat,vec, + trial::MultiFieldFESpace{MS}, + test ::MultiFieldFESpace{MS}, + strategy::AssemblyStrategy=DefaultAssemblyStrategy() + ) where MS <: BlockMultiFieldStyle + mfs = MultiFieldStyle(test) + return BlockSparseMatrixAssembler(mfs,trial,test,SparseMatrixBuilder(mat),ArrayBuilder(vec),strategy) end # BlockArrays extensions @@ -224,7 +226,7 @@ for T in (:AddEntriesMap,:TouchEntriesMap) cache end - function Fields.evaluate!(cache, k::$T,A::$MT,v::MatrixBlock,I::VectorBlock,J::VectorBlock) + function Fields.evaluate!(cache,k::$T,A::$MT,v::MatrixBlock,I::VectorBlock,J::VectorBlock) ni,nj = size(v.touched) for j in 1:nj for i in 1:ni @@ -272,22 +274,25 @@ end # In place assembly modifications (for dispatching) # We convert from BlockArray to ArrayBlock to be able to expland the blocks. # After assembly we convert back to BlockArray automatically. -function FESpaces.assemble_vector_add!(b::BlockVector,a::BlockSparseMatrixAssembler,vecdata) - b1 = ArrayBlock(b.blocks,fill(true,size(b.blocks))) +function FESpaces.assemble_vector_add!(b::AbstractBlockVector,a::BlockSparseMatrixAssembler,vecdata) + b1 = ArrayBlock(blocks(b),fill(true,blocksize(b))) b2 = expand_blocks(a,b1) FESpaces.assemble_vector_add!(b2,a,vecdata) end -function FESpaces.assemble_matrix_add!(mat::BlockMatrix,a::BlockSparseMatrixAssembler,matdata) - m1 = ArrayBlock(mat.blocks,fill(true,size(mat.blocks))) +function FESpaces.assemble_matrix_add!(mat::AbstractBlockMatrix,a::BlockSparseMatrixAssembler,matdata) + m1 = ArrayBlock(blocks(mat),fill(true,blocksize(mat))) m2 = expand_blocks(a,m1) FESpaces.assemble_matrix_add!(m2,a,matdata) end -function FESpaces.assemble_matrix_and_vector_add!(A::BlockMatrix,b::BlockVector,a::BlockSparseMatrixAssembler,data) - m1 = ArrayBlock(A.blocks,fill(true,size(A.blocks))) +function FESpaces.assemble_matrix_and_vector_add!(A::AbstractBlockMatrix, + b::AbstractBlockVector, + a::BlockSparseMatrixAssembler, + data) + m1 = ArrayBlock(blocks(A),fill(true,blocksize(A))) m2 = expand_blocks(a,m1) - b1 = ArrayBlock(b.blocks,fill(true,size(b.blocks))) + b1 = ArrayBlock(blocks(b),fill(true,blocksize(b))) b2 = expand_blocks(a,b1) FESpaces.assemble_matrix_and_vector_add!(m2,b2,a,data) end \ No newline at end of file diff --git a/src/MultiField/MultiFieldFESpaces.jl b/src/MultiField/MultiFieldFESpaces.jl index 436e28979..0df1528a7 100644 --- a/src/MultiField/MultiFieldFESpaces.jl +++ b/src/MultiField/MultiFieldFESpaces.jl @@ -34,12 +34,12 @@ function BlockMultiFieldStyle(NB::Integer) return BlockMultiFieldStyle(NB,SB) end -function BlockMultiFieldStyle(::BlockMultiFieldStyle{NB,SB,P},spaces::Vector{<:SingleFieldFESpace}) where {NB,SB,P} +function BlockMultiFieldStyle(::BlockMultiFieldStyle{NB,SB,P},spaces) where {NB,SB,P} @check length(spaces) == sum(SB) return BlockMultiFieldStyle(NB,SB,P) end -function BlockMultiFieldStyle(::BlockMultiFieldStyle{0,0,0},spaces::Vector{<:SingleFieldFESpace}) +function BlockMultiFieldStyle(::BlockMultiFieldStyle{0,0,0},spaces) NB = length(spaces) return BlockMultiFieldStyle(NB) end diff --git a/src/ODEs/Exports.jl b/src/ODEs/Exports.jl index fcb94c335..60090ccb8 100644 --- a/src/ODEs/Exports.jl +++ b/src/ODEs/Exports.jl @@ -10,6 +10,7 @@ end @publish_gridapodes ODETools MidPoint @publish_gridapodes ODETools ThetaMethod @publish_gridapodes ODETools RungeKutta +@publish_gridapodes ODETools IMEXRungeKutta @publish_gridapodes ODETools Newmark @publish_gridapodes ODETools GeneralizedAlpha @publish_gridapodes ODETools ∂t @@ -23,3 +24,5 @@ end @publish_gridapodes TransientFETools TransientAffineFEOperator @publish_gridapodes TransientFETools TransientConstantFEOperator @publish_gridapodes TransientFETools TransientConstantMatrixFEOperator +@publish_gridapodes TransientFETools TransientRungeKuttaFEOperator +@publish_gridapodes TransientFETools TransientIMEXRungeKuttaFEOperator diff --git a/src/ODEs/ODETools/IMEXRungeKutta.jl b/src/ODEs/ODETools/IMEXRungeKutta.jl new file mode 100644 index 000000000..a6211ce41 --- /dev/null +++ b/src/ODEs/ODETools/IMEXRungeKutta.jl @@ -0,0 +1,267 @@ +""" +Implicit-Explicit Runge-Kutta ODE solver. + +This struct defines an ODE solver for the system of ODEs + + M(u,t)du/dt = f(u,t) + g(u,t) + +where `f` is a nonlinear function of `u` and `t` that will treated implicitly and + `g` is a nonlinear function of `u` and `t` that will be treated explicitly. + The ODE is solved using an implicit-explicit Runge-Kutta method. +""" +struct IMEXRungeKutta <: ODESolver + nls_stage::NonlinearSolver + nls_update::NonlinearSolver + dt::Float64 + tableau::IMEXButcherTableau + function IMEXRungeKutta(nls_stage::NonlinearSolver, nls_update::NonlinearSolver, dt, type::Symbol) + bt = IMEXButcherTableau(type) + new(nls_stage, nls_update, dt, bt) + end +end + +""" +solve_step!(uf,odesol,op,u0,t0,cache) + +Solve one step of the ODE problem defined by `op` using the ODE solver `odesol` + with initial solution `u0` at time `t0`. The solution is stored in `uf` and + the final time in `tf`. The cache is used to store the solution of the + nonlinear system of equations and auxiliar variables. +""" +function solve_step!(uf::AbstractVector, + solver::IMEXRungeKutta, + op::ODEOperator, + u0::AbstractVector, + t0::Real, + cache) + + # Unpack variables + dt = solver.dt + s = solver.tableau.s + aᵢ = solver.tableau.aᵢ + bᵢ = solver.tableau.bᵢ + aₑ = solver.tableau.aₑ + bₑ = solver.tableau.bₑ + c = solver.tableau.c + d = solver.tableau.d + + # Create cache if not there + if cache === nothing + ode_cache = allocate_cache(op) + vi = similar(u0) + fi = Vector{typeof(u0)}(undef,0) + gi = Vector{typeof(u0)}(undef,0) + for i in 1:s + push!(fi,similar(u0)) + push!(gi,similar(u0)) + end + nls_stage_cache = nothing + nls_update_cache = nothing + else + ode_cache, vi, fi, gi, nls_stage_cache, nls_update_cache = cache + end + + # Create RKNL stage operator + nlop_stage = IMEXRungeKuttaStageNonlinearOperator(op,t0,dt,u0,ode_cache,vi,fi,gi,0,aᵢ,aₑ) + + # Compute intermediate stages + for i in 1:s + + # Update time + ti = t0 + c[i]*dt + ode_cache = update_cache!(ode_cache,op,ti) + update!(nlop_stage,ti,fi,gi,i) + + if(aᵢ[i,i]==0) + # Skip stage solve if a_ii=0 => u_i=u_0, f_i = f_0, gi = g_0 + @. uf = u0 + else + # solve at stage i + nls_stage_cache = solve!(uf,solver.nls_stage,nlop_stage,nls_stage_cache) + end + + # Update RHS at stage i using solution at u_i + rhs!(nlop_stage, uf) + explicit_rhs!(nlop_stage, uf) + + end + + # Update final time + tf = t0+dt + + # Skip final update if not necessary + if !(c[s]==1.0 && aᵢ[s,:] == bᵢ && aₑ[s,:] == bₑ) + + # Create RKNL final update operator + ode_cache = update_cache!(ode_cache,op,tf) + nlop_update = IMEXRungeKuttaUpdateNonlinearOperator(op,tf,dt,u0,ode_cache,vi,fi,gi,s,bᵢ,bₑ) + + # solve at final update + nls_update_cache = solve!(uf,solver.nls_update,nlop_update,nls_update_cache) + + end + + # Update final cache + cache = (ode_cache, vi, fi, gi, nls_stage_cache, nls_update_cache) + + return (uf, tf, cache) + +end + +""" +IMEXRungeKuttaStageNonlinearOperator <: NonlinearOperator + +Nonlinear operator for the implicit-explicit Runge-Kutta stage. + At a given stage `i` it represents the nonlinear operator A(t,u_i,(u_i-u_n)/dt) such that +```math +A(t,u_i,(u_i-u_n)/dt) = M(u_i,t)(u_i-u_n)/Δt - ∑aᵢ[i,j] * f(u_j,t_j) - ∑aₑ[i,j] * g(u_j,t_j) = 0 +``` +""" +mutable struct IMEXRungeKuttaStageNonlinearOperator <: RungeKuttaNonlinearOperator + odeop::ODEOperator + ti::Float64 + dt::Float64 + u0::AbstractVector + ode_cache + vi::AbstractVector + fi::Vector{AbstractVector} + gi::Vector{AbstractVector} + i::Int + aᵢ::Matrix{Float64} + aₑ::Matrix{Float64} +end + +""" +IMEXRungeKuttaUpdateNonlinearOperator <: NonlinearOperator + +Nonlinear operator for the implicit-explicit Runge-Kutta final update. + At the final update it represents the nonlinear operator A(t,u_t,(u_t-u_n)/dt) such that + ```math + A(t,u_f,(u_f-u_n)/dt) = M(u_f,t)(u_f-u_n)/Δt - ∑aᵢ[i,j] * f(u_j,t_j) - ∑aₑ[i,j] * g(u_j,t_j) = 0 + ``` +""" +mutable struct IMEXRungeKuttaUpdateNonlinearOperator <: RungeKuttaNonlinearOperator + odeop::ODEOperator + ti::Float64 + dt::Float64 + u0::AbstractVector + ode_cache + vi::AbstractVector + fi::Vector{AbstractVector} + gi::Vector{AbstractVector} + s::Int + bᵢ::Vector{Float64} + bₑ::Vector{Float64} +end + +""" +residual!(b,op::IMEXRungeKuttaStageNonlinearOperator,x) + +Compute the residual of the IMEXR Runge-Kutta nonlinear operator `op` at `x` and +store it in `b` for a given stage `i`. +```math +b = A(t,x,(x-x₀)/dt) = ∂ui/∂t - ∑aᵢ[i,j] * f(xj,tj) +``` + +Uses the vector b as auxiliar variable to store the residual of the left-hand side of +the i-th stage ODE operator, then adds the corresponding contribution from right-hand side +at all earlier stages. +```math +b = M(ui,ti)∂u/∂t +b - ∑_{j<=i} aᵢ_ij * f(uj,tj) - ∑_{j u_i=u_0, f_i = f_0 + # Update time + ti = t0 + c[i]*dt + ode_cache = update_cache!(ode_cache,op,ti) + update!(nlop_stage,ti,fi,i) + if(a[i,i]==0) - @assert c[i] == 0 - ti = t0 - update!(nlop,ti,fi,i) - fi[i] = get_fi(u0,nlop,nl_cache) + # Skip stage solve if a_ii=0 => u_i=u_0, f_i = f_0 + @. uf = u0 else # solve at stage i - ti = t0 + c[i]*dt - ode_cache = update_cache!(ode_cache,op,ti) - update!(nlop,ti,fi,i) - nl_cache = solve!(uf,solver.nls,nlop,nl_cache) - fi[i] = get_fi(uf,nlop,nl_cache) + nls_stage_cache = solve!(uf,solver.nls_stage,nlop_stage,nls_stage_cache) end + # Update RHS at stage i using solution at u_i + rhs!(nlop_stage, uf) + end - # update - uf = u0 - for i in 1:s - uf = uf + dt*b[i]*fi[i] + # Update final time + tf = t0+dt + + # Skip final update if not necessary + if !(c[s]==1.0 && a[s,:] == b) + + # Create RKNL final update operator + ode_cache = update_cache!(ode_cache,op,tf) + nlop_update = RungeKuttaUpdateNonlinearOperator(op,tf,dt,u0,ode_cache,vi,fi,s,b) + + # solve at final update + nls_update_cache = solve!(uf,solver.nls_update,nlop_update,nls_update_cache) + end - cache = (ode_cache, vi, fi, nl_cache) + # Update final cache + cache = (ode_cache, vi, fi, nls_stage_cache, nls_update_cache) - tf = t0 + dt return (uf,tf,cache) end +abstract type RungeKuttaNonlinearOperator <: NonlinearOperator end + """ Nonlinear operator that represents the Runge-Kutta nonlinear operator at a given time step and stage, i.e., A(t,u_i,(u_i-u_n)/dt) """ -mutable struct RungeKuttaNonlinearOperator <: NonlinearOperator +mutable struct RungeKuttaStageNonlinearOperator <: RungeKuttaNonlinearOperator odeop::ODEOperator ti::Float64 dt::Float64 u0::AbstractVector ode_cache vi::AbstractVector - fi::AbstractVector + fi::Vector{AbstractVector} i::Int a::Matrix end -function residual!(b::AbstractVector,op::RungeKuttaNonlinearOperator,x::AbstractVector) - # A(t,ui,∂ui/∂t) = ∂ui/∂t - a_ii * f(ui,ti) - ∑_{j 0.0) - ui = x - vi = op.vi - vi = (x-op.u0)/(op.a[op.i,op.i]*op.dt) - residual!(b,op.odeop,op.ti,(ui,vi),op.ode_cache) - for j in 1:op.i-1 - b .= b - op.a[op.i,j]/op.a[op.i,op.i] * op.fi[j] +""" +Nonlinear operator that represents the Runge-Kutta nonlinear operator at the +final updated of a given time step, i.e., A(t,u_f,(u_f-u_n)/dt) +""" +mutable struct RungeKuttaUpdateNonlinearOperator <: RungeKuttaNonlinearOperator + odeop::ODEOperator + ti::Float64 + dt::Float64 + u0::AbstractVector + ode_cache + vi::AbstractVector + fi::Vector{AbstractVector} + s::Int + b::Vector{Number} +end + + + +""" +Compute the residual of the Runge-Kutta nonlinear operator at stage i. +```math +A(t,ui,∂ui/∂t) = ∂ui/∂t - a_ii * f(ui,ti) - ∑_{j 0.0) +function jacobian!(A::AbstractMatrix,op::RungeKuttaStageNonlinearOperator,x::AbstractVector) + # @assert (abs(op.a[op.i,op.i]) > 0.0) ui = x vi = op.vi - vi = (x-op.u0)/(op.a[op.i,op.i]*op.dt) + @. vi = (x-op.u0)/(op.dt) z = zero(eltype(A)) fillstored!(A,z) - jacobians!(A,op.odeop,op.ti,(ui,vi),(1.0,1.0/(op.a[op.i,op.i]*op.dt)),op.ode_cache) + jacobians!(A,op.odeop,op.ti,(ui,vi),(op.a[op.i,op.i],1.0/op.dt),op.ode_cache) +end + +function jacobian!(A::AbstractMatrix,op::RungeKuttaUpdateNonlinearOperator,x::AbstractVector) + uf = x + vf = op.vi + @. vf = (x-op.u0)/(op.dt) + z = zero(eltype(A)) + fillstored!(A,z) + jacobian!(A,op.odeop,op.ti,(uf,vf),2,1.0/(op.dt),op.ode_cache) +end + +function jacobian!(A::AbstractMatrix,op::RungeKuttaStageNonlinearOperator,x::AbstractVector, + i::Integer,γᵢ::Real) + ui = x + vi = op.vi + @. vi = (x-op.u0)/(op.dt) + z = zero(eltype(A)) + fillstored!(A,z) + jacobian!(A,op.odeop,op.ti,(ui,vi),i,γᵢ,op.ode_cache) end function allocate_residual(op::RungeKuttaNonlinearOperator,x::AbstractVector) @@ -193,28 +213,19 @@ function zero_initial_guess(op::RungeKuttaNonlinearOperator) x0 end -function get_fi(x::AbstractVector, op::RungeKuttaNonlinearOperator, cache::Nothing) - ui = x - vi = op.vi - if(op.a[op.i,op.i]==0.0) - vi=zero(x) - else - vi = (x-op.u0)/(op.a[op.i,op.i]*op.dt) - end - b=similar(x) - residual!(b,op.odeop,op.ti,(ui,vi),op.ode_cache) - (vi-b) # store fi for future stages +function rhs!(op::RungeKuttaNonlinearOperator, x::AbstractVector) + u = x + v = op.vi + @. v = (x-op.u0)/(op.dt) + f = op.fi + rhs!(f[op.i],op.odeop,op.ti,(u,v),op.ode_cache) end -function get_fi(x::AbstractVector, op::RungeKuttaNonlinearOperator, cache) - ui = x - vi = op.vi - if(op.a[op.i,op.i]==0.0) - vi=zero(x) - else - vi = (x-op.u0)/(op.a[op.i,op.i]*op.dt) - end - residual!(cache.b,op.odeop,op.ti,(ui,vi),op.ode_cache) - (vi-cache.b) # store fi for future stages + +function lhs!(b::AbstractVector, op::RungeKuttaNonlinearOperator, x::AbstractVector) + u = x + v = op.vi + @. v = (x-op.u0)/(op.dt) + lhs!(b,op.odeop,op.ti,(u,v),op.ode_cache) end function update!(op::RungeKuttaNonlinearOperator,ti::Float64,fi::AbstractVector,i::Int) @@ -222,3 +233,9 @@ function update!(op::RungeKuttaNonlinearOperator,ti::Float64,fi::AbstractVector, op.fi = fi op.i = i end + +function _mass_matrix!(A,odeop,t,ode_cache,u) + z = zero(eltype(A)) + fillstored!(A,z) + jacobian!(A,odeop,t,(u,u),2,1.0,ode_cache) +end diff --git a/src/ODEs/ODETools/Tableaus.jl b/src/ODEs/ODETools/Tableaus.jl new file mode 100644 index 000000000..1fd6c2dc6 --- /dev/null +++ b/src/ODEs/ODETools/Tableaus.jl @@ -0,0 +1,192 @@ +abstract type ButcherTableauType end + +struct BE_1_0_1 <: ButcherTableauType end +struct CN_2_0_2 <: ButcherTableauType end +struct SDIRK_2_0_2 <: ButcherTableauType end +struct SDIRK_2_0_3 <: ButcherTableauType end +struct ESDIRK_3_1_2 <: ButcherTableauType end +struct TRBDF2_3_2_3 <: ButcherTableauType end + +""" +Butcher tableau +""" +struct ButcherTableau{T <: ButcherTableauType} + s::Int # stages + p::Int # embedded order + q::Int # order + a::Matrix # A_ij + b::Vector # b_j + c::Vector # c_i + d::Vector # d_j (embedded) +end + +# Butcher Tableaus constructors +""" +Backward-Euler + +number of stages: 1 +embedded method: no +order: 1 +""" +function ButcherTableau(::BE_1_0_1) + s = 1 + p = 0 + q = 1 + a = reshape([1.0],1,1) + b = [1.0] + c = [1.0] + d = [0.0] + ButcherTableau{BE_1_0_1}(s,p,q,a,b,c,d) +end + +""" +Crank-Nicolson (equivalent to trapezoidal rule) + +number of stages: 2 +embedded method: no +order: 2 +""" +function ButcherTableau(type::CN_2_0_2) + s = 2 + p = 0 + q = 2 + a = [0.0 0.0; 0.5 0.5] + b = [0.5, 0.5] + c = [0.0, 1.0] + d = [0.0, 0.0] + ButcherTableau{CN_2_0_2}(s,p,q,a,b,c,d) +end + +""" +Qin and Zhang's SDIRK + +number of stages: 2 +embedded method: no +order: 2 +""" +function ButcherTableau(type::SDIRK_2_0_2) + s = 2 + p = 0 + q = 2 + a = [0.25 0.0; 0.5 0.25] + b = [0.5, 0.5] + c = [0.25, 0.75] + d = [0.0, 0.0] + ButcherTableau{SDIRK_2_0_2}(s,p,q,a,b,c,d) +end + +""" +3rd order SDIRK + +number of stages: 2 +embedded method: no +order: 3 +""" +function ButcherTableau(type::SDIRK_2_0_3) + s = 2 + p = 0 + q = 3 + γ = (3-√(3))/6 + a = [γ 0.0; 1-2γ γ] + b = [0.5, 0.5] + c = [γ, 1-γ] + d = [0.0, 0.0] + ButcherTableau{SDIRK_2_0_3}(s,p,q,a,b,c,d) +end + +function ButcherTableau(type::ESDIRK_3_1_2) +s = 3 +p = 1 +q = 2 +γ = (2-√(2))/2 +b₂ = (1 − 2γ)/(4γ) +b̂₂ = γ*(−2 + 7γ − 5(γ^2) + 4(γ^3)) / (2(2γ − 1)) +b̂₃ = −2*(γ^2)*(1 − γ + γ^2) / (2γ − 1) +a = [0.0 0.0 0.0; γ γ 0.0; (1 − b₂ − γ) b₂ γ] +b = [(1 − b₂ − γ), b₂, γ] +c = [0.0, 2γ, 1.0] +d = [(1 − b̂₂ − b̂₃), b̂₂, b̂₃] +ButcherTableau{ESDIRK_3_1_2}(s,p,q,a,b,c,d) +end + +function ButcherTableau(type::TRBDF2_3_2_3) + s = 3 + p = 2 + q = 3 + aux = 2.0-√2.0 + a = [0.0 0.0 0.0; aux/2 aux/2 0.0; √2/4 √2/4 aux/2] + b = [√2/4, √2/4, aux/2] + c = [0.0, aux, 1.0] + d = [(1.0-(√2/4))/3, ((3*√2)/4+1.0)/3, aux/6] + ButcherTableau{TRBDF2_3_2_3}(s,p,q,a,b,c,d) +end + +function ButcherTableau(type::Symbol) + eval(:(ButcherTableau($type()))) +end + +abstract type IMEXButcherTableauType end + +struct IMEX_FE_BE_2_0_1 <: IMEXButcherTableauType end +struct IMEX_Midpoint_2_0_2 <: IMEXButcherTableauType end + +""" +Implicit-Explicit Butcher tableaus +""" +struct IMEXButcherTableau{T <: IMEXButcherTableauType} + s::Int # stages + p::Int # embedded order + q::Int # order + aᵢ::Matrix # A_ij implicit + aₑ::Matrix # A_ij explicit + bᵢ::Vector # b_j implicit + bₑ::Vector # b_j explicit + c::Vector # c_i + d::Vector # d_j (embedded) +end + +# IMEX Butcher Tableaus constructors +""" +IMEX Forward-Backward-Euler + +number of stages: 2 +embedded method: no +order: 1 +""" +function IMEXButcherTableau(::IMEX_FE_BE_2_0_1) + s = 2 + p = 0 + q = 1 + aᵢ = [0.0 0.0; 0.0 1.0] + aₑ = [0.0 0.0; 1.0 0.0] + bᵢ = [0.0, 1.0] + bₑ = [0.0, 1.0] + c = [0.0, 1.0] + d = [0.0, 0.0] + IMEXButcherTableau{IMEX_FE_BE_2_0_1}(s,p,q,aᵢ,aₑ,bᵢ,bₑ,c,d) +end + +""" +IMEX Midpoint + +number of stages: 2 +embedded method: no +order: 2 +""" +function IMEXButcherTableau(::IMEX_Midpoint_2_0_2) + s = 2 + p = 0 + q = 2 + aᵢ = [0.0 0.0; 0.0 0.5] + aₑ = [0.0 0.0; 0.5 0.0] + bᵢ = [0.0, 1.0] + bₑ = [0.0, 1.0] + c = [0.0, 0.5] + d = [0.0, 0.0] + IMEXButcherTableau{IMEX_Midpoint_2_0_2}(s,p,q,aᵢ,aₑ,bᵢ,bₑ,c,d) +end + + +function IMEXButcherTableau(type::Symbol) + eval(:(IMEXButcherTableau($type()))) +end diff --git a/src/ODEs/TransientFETools/ODEOperatorInterfaces.jl b/src/ODEs/TransientFETools/ODEOperatorInterfaces.jl index 939d53b77..b31dc9cae 100644 --- a/src/ODEs/TransientFETools/ODEOperatorInterfaces.jl +++ b/src/ODEs/TransientFETools/ODEOperatorInterfaces.jl @@ -121,3 +121,57 @@ function jacobians!( xh=TransientCellField(EvaluationFunction(Xh[1],xhF[1]),dxh) jacobians!(J,op.feop,t,xh,γ,ode_cache) end + +""" +It provides the Left hand side, RHS, of LHS(t,uh,∂tuh) = RHS(t,uh) for a given (t,uh,∂tuh,...,∂t^Nuh) +""" +function lhs!( + lhs::AbstractVector, + op::ODEOpFromFEOp, + t::Real, + xhF::Tuple{Vararg{AbstractVector}}, + ode_cache) + Xh, = ode_cache + dxh = () + for i in 2:get_order(op)+1 + dxh = (dxh...,EvaluationFunction(Xh[i],xhF[i])) + end + xh=TransientCellField(EvaluationFunction(Xh[1],xhF[1]),dxh) + lhs!(lhs,op.feop,t,xh,ode_cache) +end + +""" +It provides the Right hand side, RHS, of LHS(t,uh,∂tuh) = RHS(t,uh) for a given (t,uh,∂tuh,...,∂t^Nuh) +""" +function rhs!( + rhs::AbstractVector, + op::ODEOpFromFEOp, + t::Real, + xhF::Tuple{Vararg{AbstractVector}}, + ode_cache) + Xh, = ode_cache + dxh = () + for i in 2:get_order(op)+1 + dxh = (dxh...,EvaluationFunction(Xh[i],xhF[i])) + end + xh=TransientCellField(EvaluationFunction(Xh[1],xhF[1]),dxh) + rhs!(rhs,op.feop,t,xh,ode_cache) +end + +""" +It provides the explicit right hand side, E_RHS, of LHS(t,uh,∂tuh) = I_RHS(t,uh) + E_RHS(t,uh) for a given (t,uh,∂tuh,...,∂t^Nuh) +""" +function explicit_rhs!( + explicit_rhs::AbstractVector, + op::ODEOpFromFEOp, + t::Real, + xhF::Tuple{Vararg{AbstractVector}}, + ode_cache) + Xh, = ode_cache + dxh = () + for i in 2:get_order(op)+1 + dxh = (dxh...,EvaluationFunction(Xh[i],xhF[i])) + end + xh=TransientCellField(EvaluationFunction(Xh[1],xhF[1]),dxh) + explicit_rhs!(explicit_rhs,op.feop,t,xh,ode_cache) +end diff --git a/src/ODEs/TransientFETools/TransientFEOperators.jl b/src/ODEs/TransientFETools/TransientFEOperators.jl index 34047967c..557fb0f85 100644 --- a/src/ODEs/TransientFETools/TransientFEOperators.jl +++ b/src/ODEs/TransientFETools/TransientFEOperators.jl @@ -103,6 +103,7 @@ Transient FE operator that is defined by a transient Weak form """ struct TransientFEOperatorFromWeakForm{C} <: TransientFEOperator{C} res::Function + rhs::Function jacs::Tuple{Vararg{Function}} assem_t::Assembler trials::Tuple{Vararg{Any}} @@ -113,40 +114,44 @@ end function TransientConstantFEOperator(m::Function,a::Function,b::Function, trial,test) res(t,u,v) = -1.0 * b(v) + rhs(t,u,v) = b(v) jac(t,u,du,v) = a(du,v) jac_t(t,u,dut,v) = m(dut,v) assem_t = SparseMatrixAssembler(trial,test) - TransientFEOperatorFromWeakForm{Constant}(res,(jac,jac_t),assem_t,(trial,∂t(trial)),test,1) + TransientFEOperatorFromWeakForm{Constant}(res,rhs,(jac,jac_t),assem_t,(trial,∂t(trial)),test,1) end function TransientConstantMatrixFEOperator(m::Function,a::Function,b::Function, trial,test) res(t,u,v) = m(∂t(u),v) + a(u,v) - b(t,v) + rhs(t,u,v) = b(t,v) - a(u,v) jac(t,u,du,v) = a(du,v) jac_t(t,u,dut,v) = m(dut,v) assem_t = SparseMatrixAssembler(trial,test) - TransientFEOperatorFromWeakForm{ConstantMatrix}(res,(jac,jac_t),assem_t,(trial,∂t(trial)),test,1) + TransientFEOperatorFromWeakForm{ConstantMatrix}(res,rhs,(jac,jac_t),assem_t,(trial,∂t(trial)),test,1) end function TransientAffineFEOperator(m::Function,a::Function,b::Function, trial,test) res(t,u,v) = m(t,∂t(u),v) + a(t,u,v) - b(t,v) + rhs(t,u,v) = b(t,v) - a(t,u,v) jac(t,u,du,v) = a(t,du,v) jac_t(t,u,dut,v) = m(t,dut,v) assem_t = SparseMatrixAssembler(trial,test) - TransientFEOperatorFromWeakForm{Affine}(res,(jac,jac_t),assem_t,(trial,∂t(trial)),test,1) + TransientFEOperatorFromWeakForm{Affine}(res,rhs,(jac,jac_t),assem_t,(trial,∂t(trial)),test,1) end function TransientFEOperator(res::Function,jac::Function,jac_t::Function, trial,test) assem_t = SparseMatrixAssembler(trial,test) - TransientFEOperatorFromWeakForm{Nonlinear}(res,(jac,jac_t),assem_t,(trial,∂t(trial)),test,1) + TransientFEOperatorFromWeakForm{Nonlinear}(res,rhs_error,(jac,jac_t),assem_t,(trial,∂t(trial)),test,1) end function TransientConstantFEOperator(m::Function,c::Function,a::Function,b::Function, trial,test) res(t,u,v) = -1.0 * b(v) + rhs(t,u,v) = b(v) jac(t,u,du,v) = a(du,v) jac_t(t,u,dut,v) = c(dut,v) jac_tt(t,u,dutt,v) = m(dutt,v) @@ -154,12 +159,13 @@ function TransientConstantFEOperator(m::Function,c::Function,a::Function,b::Func trial_t = ∂t(trial) trial_tt = ∂t(trial_t) TransientFEOperatorFromWeakForm{Constant}( - res,(jac,jac_t,jac_tt),assem_t,(trial,trial_t,trial_tt),test,2) + res,rhs,(jac,jac_t,jac_tt),assem_t,(trial,trial_t,trial_tt),test,2) end function TransientConstantMatrixFEOperator(m::Function,c::Function,a::Function,b::Function, trial,test) res(t,u,v) = m(∂tt(u),v) + c(∂t(u),v) + a(u,v) - b(t,v) + rhs(t,u,v) = b(t,v) - c(∂t(u),v) - a(u,v) jac(t,u,du,v) = a(du,v) jac_t(t,u,dut,v) = c(dut,v) jac_tt(t,u,dutt,v) = m(dutt,v) @@ -167,12 +173,13 @@ function TransientConstantMatrixFEOperator(m::Function,c::Function,a::Function,b trial_t = ∂t(trial) trial_tt = ∂t(trial_t) TransientFEOperatorFromWeakForm{ConstantMatrix}( - res,(jac,jac_t,jac_tt),assem_t,(trial,trial_t,trial_tt),test,2) + res,rhs,(jac,jac_t,jac_tt),assem_t,(trial,trial_t,trial_tt),test,2) end function TransientAffineFEOperator(m::Function,c::Function,a::Function,b::Function, trial,test) res(t,u,v) = m(t,∂tt(u),v) + c(t,∂t(u),v) + a(t,u,v) - b(t,v) + rhs(t,u,v) = b(t,v) - c(t,∂t(u),v) - a(t,u,v) jac(t,u,du,v) = a(t,du,v) jac_t(t,u,dut,v) = c(t,dut,v) jac_tt(t,u,dutt,v) = m(t,dutt,v) @@ -180,7 +187,7 @@ function TransientAffineFEOperator(m::Function,c::Function,a::Function,b::Functi trial_t = ∂t(trial) trial_tt = ∂t(trial_t) TransientFEOperatorFromWeakForm{Affine}( - res,(jac,jac_t,jac_tt),assem_t,(trial,trial_t,trial_tt),test,2) + res,rhs,(jac,jac_t,jac_tt),assem_t,(trial,trial_t,trial_tt),test,2) end function TransientFEOperator(res::Function,jac::Function,jac_t::Function, @@ -189,7 +196,7 @@ function TransientFEOperator(res::Function,jac::Function,jac_t::Function, trial_t = ∂t(trial) trial_tt = ∂t(trial_t) TransientFEOperatorFromWeakForm{Nonlinear}( - res,(jac,jac_t,jac_tt),assem_t,(trial,trial_t,trial_tt),test,2) + res,rhs_error,(jac,jac_t,jac_tt),assem_t,(trial,trial_t,trial_tt),test,2) end function TransientFEOperator(res::Function,trial,test;order::Integer=1) @@ -215,22 +222,79 @@ function TransientFEOperator(res::Function,trial,test;order::Integer=1) TransientFEOperator(res,jacs...,trial,test) end -function SparseMatrixAssembler( - trial::Union{TransientTrialFESpace,TransientMultiFieldTrialFESpace}, - test::FESpace) - SparseMatrixAssembler(evaluate(trial,nothing),test) +function allocate_residual( + op::TransientFEOperatorFromWeakForm, + t0::Real, + uh::T, + cache) where T + V = get_test(op) + v = get_fe_basis(V) + dxh = () + for i in 1:get_order(op) + dxh = (dxh...,uh) + end + xh = TransientCellField(uh,dxh) + vecdata = collect_cell_vector(V,op.res(t0,xh,v)) + allocate_vector(op.assem_t,vecdata) end -get_assembler(op::TransientFEOperatorFromWeakForm) = op.assem_t +function residual!( + b::AbstractVector, + op::TransientFEOperatorFromWeakForm, + t::Real, + xh::T, + cache) where T + V = get_test(op) + v = get_fe_basis(V) + vecdata = collect_cell_vector(V,op.res(t,xh,v)) + assemble_vector!(b,op.assem_t,vecdata) + b +end -get_test(op::TransientFEOperatorFromWeakForm) = op.test +""" +Transient FE operator that is defined by a transient Weak form with the +form: LHS(t,u,∂u/∂t,...) ∂u/∂t = RHS(t,u,∂u/∂t,...). Used in Runge-Kutta schemes +""" +struct TransientRKFEOperatorFromWeakForm{C} <: TransientFEOperator{C} + lhs::Function + rhs::Function + jacs::Tuple{Vararg{Function}} + assem_t::Assembler + trials::Tuple{Vararg{Any}} + test::FESpace + order::Integer +end -get_trial(op::TransientFEOperatorFromWeakForm) = op.trials[1] +function TransientRungeKuttaFEOperator(lhs::Function,rhs::Function,jac::Function, + jac_t::Function,trial,test) + assem_t = SparseMatrixAssembler(trial,test) + TransientRKFEOperatorFromWeakForm{Nonlinear}(lhs,rhs,(jac,jac_t),assem_t,(trial,∂t(trial)),test,1) +end -get_order(op::TransientFEOperatorFromWeakForm) = op.order +function TransientRungeKuttaFEOperator(lhs::Function,rhs::Function,trial,test) + res(t,u,v) = lhs(t,u,v) - rhs(t,u,v) + function jac_0(t,x,dx0,dv) + function res_0(y) + x0 = TransientCellField(y,x.derivatives) + res(t,x0,dv) + end + jacobian(res_0,x.cellfield) + end + jacs = (jac_0,) + function jac_t(t,x,dxt,dv) + function res_t(y) + derivatives = (y,x.derivatives[2:end]...) + xt = TransientCellField(x.cellfield,derivatives) + res(t,xt,dv) + end + jacobian(res_t,x.derivatives[1]) + end + jacs = (jac_0,jac_t) + TransientRungeKuttaFEOperator(lhs,rhs,jacs...,trial,test) +end function allocate_residual( - op::TransientFEOperatorFromWeakForm, + op::TransientRKFEOperatorFromWeakForm, t0::Real, uh::T, cache) where T @@ -241,25 +305,159 @@ function allocate_residual( dxh = (dxh...,uh) end xh = TransientCellField(uh,dxh) - vecdata = collect_cell_vector(V,op.res(t0,xh,v)) + vecdata = collect_cell_vector(V,op.lhs(t0,xh,v)) allocate_vector(op.assem_t,vecdata) end -function residual!( +function lhs!( b::AbstractVector, - op::TransientFEOperatorFromWeakForm, + op::TransientRKFEOperatorFromWeakForm, t::Real, xh::T, cache) where T V = get_test(op) v = get_fe_basis(V) - vecdata = collect_cell_vector(V,op.res(t,xh,v)) + vecdata = collect_cell_vector(V,op.lhs(t,xh,v)) assemble_vector!(b,op.assem_t,vecdata) b end +function rhs!( + rhs::AbstractVector, + op::TransientRKFEOperatorFromWeakForm, + t::Real, + xh::T, + cache) where T + V = get_test(op) + v = get_fe_basis(V) + vecdata = collect_cell_vector(V,op.rhs(t,xh,v)) + assemble_vector!(rhs,op.assem_t,vecdata) + rhs +end + +# IMEX-RK Transient FE operators +""" +Transient FE operator that is defined by a transient Weak form with the +form: LHS(t,u,∂u/∂t,...) ∂u/∂t = I_RHS(t,u,∂u/∂t,...) + E_RHS(t,u,∂u/∂t,...). +Used in Implicit-Explicit Runge-Kutta schemes +""" +struct TransientIMEXRKFEOperatorFromWeakForm{C} <: TransientFEOperator{C} + lhs::Function + rhs::Function + explicit_rhs::Function + jacs::Tuple{Vararg{Function}} + assem_t::Assembler + trials::Tuple{Vararg{Any}} + test::FESpace + order::Integer +end + +function TransientIMEXRungeKuttaFEOperator(lhs::Function,rhs::Function, + explicit_rhs::Function,jac::Function,jac_t::Function,trial,test) + assem_t = SparseMatrixAssembler(trial,test) + TransientIMEXRKFEOperatorFromWeakForm{Nonlinear}(lhs,rhs,explicit_rhs,(jac,jac_t),assem_t,(trial,∂t(trial)),test,1) +end + +function TransientIMEXRungeKuttaFEOperator(lhs::Function,rhs::Function, + explicit_rhs::Function,trial,test) + res(t,u,v) = lhs(t,u,v) - rhs(t,u,v) + function jac_0(t,x,dx0,dv) + function res_0(y) + x0 = TransientCellField(y,x.derivatives) + res(t,x0,dv) + end + jacobian(res_0,x.cellfield) + end + jacs = (jac_0,) + function jac_t(t,x,dxt,dv) + function res_t(y) + derivatives = (y,x.derivatives[2:end]...) + xt = TransientCellField(x.cellfield,derivatives) + res(t,xt,dv) + end + jacobian(res_t,x.derivatives[1]) + end + jacs = (jac_0,jac_t) + TransientIMEXRungeKuttaFEOperator(lhs,rhs,explicit_rhs,jacs...,trial,test) +end + +function allocate_residual( + op::TransientIMEXRKFEOperatorFromWeakForm, + t0::Real, + uh::T, + cache) where T + V = get_test(op) + v = get_fe_basis(V) + dxh = () + for i in 1:get_order(op) + dxh = (dxh...,uh) + end + xh = TransientCellField(uh,dxh) + vecdata = collect_cell_vector(V,op.lhs(t0,xh,v)) + allocate_vector(op.assem_t,vecdata) +end + +function lhs!( + b::AbstractVector, + op::TransientIMEXRKFEOperatorFromWeakForm, + t::Real, + xh::T, + cache) where T + V = get_test(op) + v = get_fe_basis(V) + vecdata = collect_cell_vector(V,op.lhs(t,xh,v)) + assemble_vector!(b,op.assem_t,vecdata) + b +end + +function rhs!( + rhs::AbstractVector, + op::TransientIMEXRKFEOperatorFromWeakForm, + t::Real, + xh::T, + cache) where T + V = get_test(op) + v = get_fe_basis(V) + vecdata = collect_cell_vector(V,op.rhs(t,xh,v)) + assemble_vector!(rhs,op.assem_t,vecdata) + rhs +end + +function explicit_rhs!( + explicit_rhs::AbstractVector, + op::TransientIMEXRKFEOperatorFromWeakForm, + t::Real, + xh::T, + cache) where T + V = get_test(op) + v = get_fe_basis(V) + vecdata = collect_cell_vector(V,op.explicit_rhs(t,xh,v)) + assemble_vector!(explicit_rhs,op.assem_t,vecdata) + explicit_rhs +end + +# Common functions + +TransientFEOperatorsFromWeakForm = Union{TransientFEOperatorFromWeakForm, +TransientRKFEOperatorFromWeakForm, TransientIMEXRKFEOperatorFromWeakForm} + +function SparseMatrixAssembler( + trial::Union{TransientTrialFESpace,TransientMultiFieldTrialFESpace}, + test::FESpace) + SparseMatrixAssembler(evaluate(trial,nothing),test) +end + +get_assembler(op::TransientFEOperatorsFromWeakForm) = op.assem_t + +get_test(op::TransientFEOperatorsFromWeakForm) = op.test + +get_trial(op::TransientFEOperatorsFromWeakForm) = op.trials[1] + +get_order(op::TransientFEOperatorsFromWeakForm) = op.order + + function allocate_jacobian( - op::TransientFEOperatorFromWeakForm, + op::TransientFEOperatorsFromWeakForm, t0::Real, uh::CellField, cache) @@ -270,7 +468,7 @@ end function jacobian!( A::AbstractMatrix, - op::TransientFEOperatorFromWeakForm, + op::TransientFEOperatorsFromWeakForm, t::Real, xh::T, i::Integer, @@ -283,7 +481,7 @@ end function jacobians!( A::AbstractMatrix, - op::TransientFEOperatorFromWeakForm, + op::TransientFEOperatorsFromWeakForm, t::Real, xh::TransientCellField, γ::Tuple{Vararg{Real}}, @@ -294,7 +492,7 @@ function jacobians!( A end -function fill_initial_jacobians(op::TransientFEOperatorFromWeakForm,t0::Real,uh) +function fill_initial_jacobians(op::TransientFEOperatorsFromWeakForm,t0::Real,uh) dxh = () for i in 1:get_order(op) dxh = (dxh...,uh) @@ -308,7 +506,7 @@ function fill_initial_jacobians(op::TransientFEOperatorFromWeakForm,t0::Real,uh) end function fill_jacobians( - op::TransientFEOperatorFromWeakForm, + op::TransientFEOperatorsFromWeakForm, t::Real, xh::T, γ::Tuple{Vararg{Real}}) where T @@ -339,7 +537,7 @@ function _vcat_matdata(_matdata) end function _matdata_jacobian( - op::TransientFEOperatorFromWeakForm, + op::TransientFEOperatorsFromWeakForm, t::Real, xh::T, i::Integer, @@ -351,6 +549,12 @@ function _matdata_jacobian( matdata = collect_cell_matrix(Uh,V,γᵢ*op.jacs[i](t,xh,du,v)) end +function rhs_error(t::Real,xh,v) + error("The \"rhs\" function is not defined for this TransientFEOperator. + Please, try to use another type of TransientFEOperator that supports this + functionality.") +end + # Tester function test_transient_fe_operator(op::TransientFEOperator,uh) diff --git a/src/ODEs/TransientFETools/TransientFESpaces.jl b/src/ODEs/TransientFETools/TransientFESpaces.jl index b13ebec12..b6312b96b 100644 --- a/src/ODEs/TransientFETools/TransientFESpaces.jl +++ b/src/ODEs/TransientFETools/TransientFESpaces.jl @@ -34,7 +34,7 @@ end Allocate the space to be used as first argument in evaluate! """ function allocate_trial_space(U::TransientTrialFESpace) - HomogeneousTrialFESpace(U.space) + HomogeneousTrialFESpace(U.space) end """ @@ -43,14 +43,14 @@ Time evaluation allocating Dirichlet vals function evaluate(U::TransientTrialFESpace,t::Real) Ut = allocate_trial_space(U) evaluate!(Ut,U,t) - Ut + return Ut end """ We can evaluate at `nothing` when we do not care about the Dirichlet vals """ function evaluate(U::TransientTrialFESpace,t::Nothing) - U.Ud0 + return U.Ud0 end evaluate(U::TrialFESpace,t::Nothing) = U @@ -80,6 +80,11 @@ Time 2nd derivative of the Dirichlet functions ∂tt(U::MultiFieldFESpace) = MultiFieldFESpace(∂tt.(U.spaces)) ∂tt(t::T) where T<:Number = zero(T) +zero_free_values(f::TransientTrialFESpace) = zero_free_values(f.space) +has_constraints(f::TransientTrialFESpace) = has_constraints(f.space) +get_dof_value_type(f::TransientTrialFESpace) = get_dof_value_type(f.space) +get_vector_type(f::TransientTrialFESpace) = get_vector_type(f.space) + # Testing the interface function test_transient_trial_fe_space(Uh) @@ -100,68 +105,122 @@ end # Define the TransientTrialFESpace interface for stationary spaces -function evaluate!(Ut::FESpace,U::FESpace,t::Real) - U -end +evaluate!(Ut::FESpace,U::FESpace,t::Real) = U +allocate_trial_space(U::FESpace) = U +evaluate(U::FESpace,t::Real) = U +evaluate(U::FESpace,t::Nothing) = U -function allocate_trial_space(U::FESpace) - U +@static if VERSION >= v"1.3" + (U::FESpace)(t) = U end -function evaluate(U::FESpace,t::Real) - U +# Define the interface for MultiField + +struct TransientMultiFieldTrialFESpace{MS<:MultiFieldStyle,CS<:ConstraintStyle,V} + vector_type::Type{V} + spaces::Vector + multi_field_style::MS + constraint_style::CS + function TransientMultiFieldTrialFESpace( + ::Type{V}, + spaces::Vector, + multi_field_style::MultiFieldStyle) where V + @assert length(spaces) > 0 + + MS = typeof(multi_field_style) + if any( map(has_constraints,spaces) ) + constraint_style = Constrained() + else + constraint_style = UnConstrained() + end + CS = typeof(constraint_style) + new{MS,CS,V}(V,spaces,multi_field_style,constraint_style) + end end -function evaluate(U::FESpace,t::Nothing) - U +# Default constructors +function TransientMultiFieldFESpace(spaces::Vector; + style = ConsecutiveMultiFieldStyle()) + Ts = map(get_dof_value_type,spaces) + T = typeof(*(map(zero,Ts)...)) + if isa(style,BlockMultiFieldStyle) + style = BlockMultiFieldStyle(style,spaces) + VT = typeof(mortar(map(zero_free_values,spaces))) + else + VT = Vector{T} + end + TransientMultiFieldTrialFESpace(VT,spaces,style) end -@static if VERSION >= v"1.3" - (U::FESpace)(t) = U +function TransientMultiFieldFESpace(::Type{V},spaces::Vector) where V + TransientMultiFieldTrialFESpace(V,spaces,ConsecutiveMultiFieldStyle()) end -# Define the interface for MultiField +function TransientMultiFieldFESpace(spaces::Vector{<:SingleFieldFESpace}; + style = ConsecutiveMultiFieldStyle()) + MultiFieldFESpace(spaces,style=style) +end -struct TransientMultiFieldTrialFESpace - spaces::Vector +function TransientMultiFieldFESpace(::Type{V},spaces::Vector{<:SingleFieldFESpace}) where V + MultiFieldFESpace(V,spaces,ConsecutiveMultiFieldStyle()) end + Base.iterate(m::TransientMultiFieldTrialFESpace) = iterate(m.spaces) Base.iterate(m::TransientMultiFieldTrialFESpace,state) = iterate(m.spaces,state) Base.getindex(m::TransientMultiFieldTrialFESpace,field_id::Integer) = m.spaces[field_id] Base.length(m::TransientMultiFieldTrialFESpace) = length(m.spaces) - -function TransientMultiFieldFESpace(spaces::Vector) - TransientMultiFieldTrialFESpace(spaces) -end - -function TransientMultiFieldFESpace(spaces::Vector{<:SingleFieldFESpace}) - MultiFieldFESpace(spaces) -end - function evaluate!(Ut::T,U::TransientMultiFieldTrialFESpace,t::Real) where T spaces_at_t = [evaluate!(Uti,Ui,t) for (Uti,Ui) in zip(Ut,U)] - MultiFieldFESpace(spaces_at_t) + mfs = MultiFieldStyle(U) + return MultiFieldFESpace(spaces_at_t;style=mfs) end function allocate_trial_space(U::TransientMultiFieldTrialFESpace) spaces = allocate_trial_space.(U.spaces) - MultiFieldFESpace(spaces) + mfs = MultiFieldStyle(U) + return MultiFieldFESpace(spaces;style=mfs) end function evaluate(U::TransientMultiFieldTrialFESpace,t::Real) Ut = allocate_trial_space(U) evaluate!(Ut,U,t) - Ut + return Ut end function evaluate(U::TransientMultiFieldTrialFESpace,t::Nothing) - MultiFieldFESpace([evaluate(fesp,nothing) for fesp in U.spaces]) + spaces = [evaluate(fesp,nothing) for fesp in U.spaces] + mfs = MultiFieldStyle(U) + MultiFieldFESpace(spaces;style=mfs) end (U::TransientMultiFieldTrialFESpace)(t) = evaluate(U,t) function ∂t(U::TransientMultiFieldTrialFESpace) spaces = ∂t.(U.spaces) - TransientMultiFieldFESpace(spaces) + mfs = MultiFieldStyle(U) + TransientMultiFieldFESpace(spaces;style=mfs) +end + +function zero_free_values(f::TransientMultiFieldTrialFESpace{<:BlockMultiFieldStyle{NB,SB,P}}) where {NB,SB,P} + block_ranges = get_block_ranges(NB,SB,P) + block_num_dofs = map(range->sum(map(num_free_dofs,f.spaces[range])),block_ranges) + block_vtypes = map(range->get_vector_type(first(f.spaces[range])),block_ranges) + return mortar(map(allocate_vector,block_vtypes,block_num_dofs)) +end + +get_dof_value_type(f::TransientMultiFieldTrialFESpace{MS,CS,V}) where {MS,CS,V} = eltype(V) +get_vector_type(f::TransientMultiFieldTrialFESpace) = f.vector_type +ConstraintStyle(::Type{TransientMultiFieldTrialFESpace{S,B,V}}) where {S,B,V} = B() +ConstraintStyle(::TransientMultiFieldTrialFESpace) = ConstraintStyle(typeof(f)) +MultiFieldStyle(::Type{TransientMultiFieldTrialFESpace{S,B,V}}) where {S,B,V} = S() +MultiFieldStyle(f::TransientMultiFieldTrialFESpace) = MultiFieldStyle(typeof(f)) + +function SparseMatrixAssembler(mat,vec, + trial::TransientMultiFieldTrialFESpace{MS}, + test ::TransientMultiFieldTrialFESpace{MS}, + strategy::AssemblyStrategy=DefaultAssemblyStrategy() + ) where MS <: BlockMultiFieldStyle + mfs = MultiFieldStyle(test) + return BlockSparseMatrixAssembler(mfs,trial,test,SparseMatrixBuilder(mat),ArrayBuilder(vec),strategy) end diff --git a/src/ODEs/TransientFETools/TransientFETools.jl b/src/ODEs/TransientFETools/TransientFETools.jl index 4cc6f801d..08f862407 100644 --- a/src/ODEs/TransientFETools/TransientFETools.jl +++ b/src/ODEs/TransientFETools/TransientFETools.jl @@ -39,6 +39,8 @@ export TransientFEOperator export TransientAffineFEOperator export TransientConstantFEOperator export TransientConstantMatrixFEOperator +export TransientRungeKuttaFEOperator +export TransientIMEXRungeKuttaFEOperator using Gridap.FESpaces: Assembler using Gridap.FESpaces: SparseMatrixAssembler import Gridap.ODEs.ODETools: allocate_cache @@ -52,6 +54,9 @@ import Gridap.ODEs.ODETools: allocate_jacobian import Gridap.ODEs.ODETools: residual! import Gridap.ODEs.ODETools: jacobian! import Gridap.ODEs.ODETools: jacobians! +import Gridap.ODEs.ODETools: lhs! +import Gridap.ODEs.ODETools: rhs! +import Gridap.ODEs.ODETools: explicit_rhs! import Gridap.ODEs.ODETools: OperatorType using Gridap.ODEs.ODETools: Nonlinear using Gridap.ODEs.ODETools: Affine @@ -107,6 +112,13 @@ import Gridap.CellData: gradient import Gridap.CellData: ∇∇ import Gridap.CellData: change_domain import Gridap.FESpaces: BasisStyle +using Gridap.FESpaces: Constrained, UnConstrained, AssemblyStrategy +using Gridap.MultiField: ConsecutiveMultiFieldStyle, BlockSparseMatrixAssembler +import Gridap.MultiField: ConstraintStyle, MultiFieldStyle, BlockMultiFieldStyle +import Gridap.FESpaces: zero_free_values, has_constraints, SparseMatrixAssembler +import Gridap.FESpaces: get_dof_value_type, get_vector_type + +using BlockArrays include("TransientFESpaces.jl") diff --git a/test/AdaptivityTests/FineToCoarseFieldsTests.jl b/test/AdaptivityTests/FineToCoarseFieldsTests.jl index 3d37210ed..e09b54f8d 100644 --- a/test/AdaptivityTests/FineToCoarseFieldsTests.jl +++ b/test/AdaptivityTests/FineToCoarseFieldsTests.jl @@ -91,4 +91,13 @@ eh2 = u_c - u_fc eh3 = u_c - u_fc2 @test sum(∫(eh3⋅eh3)*dΩ_c) < 1.e-12 +modelH=CartesianDiscreteModel((0,1,0,1),(1,1)) +modelh=refine(modelH,2) +reffe=LagrangianRefFE(Float64,QUAD,1) +XH = TestFESpace(modelH,reffe) +xH = get_fe_basis(XH) +xHh = change_domain(xH,get_triangulation(modelh),ReferenceDomain()) +evaluate(Gridap.CellData.get_data(xHh)[1], + [Point(0.0,0.0),Point(0.5,0.5)]) + end \ No newline at end of file diff --git a/test/ODEsTests/ODEsTests/ODEOperatorMocks.jl b/test/ODEsTests/ODEsTests/ODEOperatorMocks.jl index 0c2cdee35..cfec64c2b 100644 --- a/test/ODEsTests/ODEsTests/ODEOperatorMocks.jl +++ b/test/ODEsTests/ODEsTests/ODEOperatorMocks.jl @@ -16,6 +16,9 @@ import Gridap.ODEs.ODETools: jacobian! import Gridap.ODEs.ODETools: jacobians! import Gridap.ODEs.ODETools: allocate_jacobian import Gridap.ODEs.ODETools: residual! +import Gridap.ODEs.ODETools: rhs! +import Gridap.ODEs.ODETools: explicit_rhs! +import Gridap.ODEs.ODETools: lhs! using SparseArrays: spzeros struct ODEOperatorMock{T<:Real,C} <: ODEOperator{C} @@ -35,6 +38,28 @@ function residual!(r::AbstractVector,op::ODEOperatorMock,t::Real,x::NTuple{2,Abs r end +function rhs!(r::AbstractVector,op::ODEOperatorMock,t::Real,x::NTuple{2,AbstractVector},ode_cache) + u,u_t = x + r .= 0 + r[1] = op.a * u[1] + r[2] = op.b * u[1] + op.c * u[2] + r +end + +function explicit_rhs!(r::AbstractVector,op::ODEOperatorMock,t::Real,x::NTuple{2,AbstractVector},ode_cache) + u,u_t = x + r .= 0 + r +end + +function lhs!(r::AbstractVector,op::ODEOperatorMock,t::Real,x::NTuple{2,AbstractVector},ode_cache) + u,u_t = x + r .= 0 + r[1] = u_t[1] + r[2] = u_t[2] + r +end + function residual!(r::AbstractVector,op::ODEOperatorMock,t::Real,x::NTuple{3,AbstractVector},ode_cache) u,u_t,u_tt = x r .= 0 @@ -57,9 +82,9 @@ function jacobian!(J::AbstractMatrix, @assert get_order(op) == 1 @assert 0 < i <= get_order(op)+1 if i==1 - J[1,1] += -op.a - J[2,1] += -op.b - J[2,2] += -op.c + J[1,1] += -op.a*γᵢ + J[2,1] += -op.b*γᵢ + J[2,2] += -op.c*γᵢ elseif i==2 J[1,1] += 1.0*γᵢ J[2,2] += 1.0*γᵢ @@ -77,9 +102,9 @@ function jacobian!(J::AbstractMatrix, @assert get_order(op) == 2 @assert 0 < i <= get_order(op)+1 if i==1 - J[1,1] += -op.a - J[2,1] += -op.b - J[2,2] += -op.c + J[1,1] += -op.a*γᵢ + J[2,1] += -op.b*γᵢ + J[2,2] += -op.c*γᵢ elseif i==2 J[1,1] += op.b*γᵢ J[2,2] += op.a*γᵢ diff --git a/test/ODEsTests/ODEsTests/ODESolversTests.jl b/test/ODEsTests/ODEsTests/ODESolversTests.jl index 20aafb256..3794a1168 100644 --- a/test/ODEsTests/ODEsTests/ODESolversTests.jl +++ b/test/ODEsTests/ODEsTests/ODESolversTests.jl @@ -95,67 +95,84 @@ _J = jacobian(sop,x) @test all(_r .== [ -11.0 -11.0]) @test all(_J .== [ 9.0 0.0; 0.0 9.0]) -# BackwardEuler tests - ls = LUSolver() + +# BackwardEuler tests odesol = BackwardEuler(ls,dt) uf = copy(u0) uf.=1.0 cache = nothing -# Juno.@enter solve_step!(uf,odesol,op,u0,t0,cache) uf, tf, cache = solve_step!(uf,odesol,op,u0,t0,cache) @test tf==t0+dt @test all(uf.≈1+11/9) - @test test_ode_solver(odesol,op,u0,t0,tf) -test_ode_solver(odesol,op,u0,t0,tf) # Affine and nonlinear solvers - op = ODEOperatorMock{Float64,Nonlinear}(1.0,0.0,1.0,1) cache = nothing uf, tf, cache = solve_step!(uf,odesol,op,u0,t0,cache) @test tf==t0+dt @test all(uf.≈1+11/9) +@test test_ode_solver(odesol,op,u0,t0,tf) op = ODEOperatorMock{Float64,Affine}(1.0,0.0,1.0,1) cache = nothing uf, tf, cache = solve_step!(uf,odesol,op,u0,t0,cache) @test tf==t0+dt @test all(uf.≈1+11/9) +@test test_ode_solver(odesol,op,u0,t0,tf) -# RK tests +# ThetaMethod tests +odesolθ = ThetaMethod(ls,dt,0.5) +ufθ = copy(u0) +ufθ.=1.0 +ufθ, tf, cache = solve_step!(ufθ,odesolθ,op,u0,t0,nothing) +@test all(ufθ.≈(dt/(1-0.5dt) + 1)*u0) +@test test_ode_solver(odesol,op,u0,t0,tf) -ls = LUSolver() -# BE equivalent -odesol = RungeKutta(ls,dt,:BE_1_0_1) -uf = copy(u0) -uf.=1.0 +# RK tests +# RK: BE equivalent +odesol = RungeKutta(ls,ls,dt,:BE_1_0_1) cache = nothing uf, tf, cache = solve_step!(uf,odesol,op,u0,t0,cache) @test tf==t0+dt @test all(uf.≈1+11/9) -# SDIRK 2nd order -odesol = RungeKutta(ls,dt,:SDIRK_2_1_2) -uf = copy(u0) -uf.=1.0 +@test test_ode_solver(odesol,op,u0,t0,tf) + +# RK: CN 2nd order +odesol = RungeKutta(ls,ls,dt,:CN_2_0_2) cache = nothing uf, tf, cache = solve_step!(uf,odesol,op,u0,t0,cache) @test tf==t0+dt -@test all(uf.≈u0*(1.0+dt/(2*(1-dt))+dt*(1-2*dt)/(2*(1-dt)^2))) -# TRBDF (2nd order with some 0 on the diagonal) -odesol = RungeKutta(ls,dt,:TRBDF2_3_3_2) -uf.=1.0 +@test all(uf.≈ufθ) +@test test_ode_solver(odesol,op,u0,t0,tf) + +# RK: SDIRK 2nd order +odesol = RungeKutta(ls,ls,dt,:SDIRK_2_0_2) +cache = nothing +uf, tf, cache = solve_step!(uf,odesol,op,u0,t0,cache) +@test tf==t0+dt +@test all(uf.≈u0*1.1051939513477975) +@test test_ode_solver(odesol,op,u0,t0,tf) + +# RK: TRBDF (2nd order with some 0 on the diagonal) +odesol = RungeKutta(ls,ls,dt,:TRBDF2_3_2_3) cache = nothing uf, tf, cache = solve_step!(uf,odesol,op,u0,t0,cache) @test tf==t0+dt @test all(uf.≈u0*1.105215241) +@test test_ode_solver(odesol,op,u0,t0,tf) +# IMEX RK tests (explicit part = 0) +odesol = IMEXRungeKutta(ls,ls,dt,:IMEX_FE_BE_2_0_1) +cache = nothing +uf, tf, cache = solve_step!(uf,odesol,op,u0,t0,cache) +@test tf==t0+dt +@test all(uf.≈1+11/9) @test test_ode_solver(odesol,op,u0,t0,tf) -test_ode_solver(odesol,op,u0,t0,tf) -# Newmark test +# Newmark test op_const = ODEOperatorMock{Float64,Constant}(1.0,0.0,0.0,2) op_const_mat = ODEOperatorMock{Float64,ConstantMatrix}(1.0,0.0,0.0,2) op_affine = ODEOperatorMock{Float64,Affine}(1.0,0.0,0.0,2) @@ -182,31 +199,23 @@ for op in ops end # GeneralizedAlpha test - op = ODEOperatorMock{Float64,Nonlinear}(1.0,0.0,1.0,1) -ls = LUSolver() ρ∞ = 1.0 # Equivalent to θ-method with θ=0.5 αf = 1.0/(1.0 + ρ∞) αm = 0.5 * (3-ρ∞) / (1+ρ∞) γ = 0.5 + αm - αf odesolα = GeneralizedAlpha(ls,dt,ρ∞) -odesolθ = ThetaMethod(ls,dt,0.5) ufα = copy(u0) -ufθ = copy(u0) v0 = 0.0*ones(2) vf = copy(v0) ufα.=1.0 -ufθ.=1.0 (ufα, vf), tf, cache = solve_step!((ufα,vf),odesolα,op,(u0,v0),t0,nothing) -ufθ, tf, cache = solve_step!(ufθ,odesolθ,op,u0,t0,nothing) @test tf==t0+dt @test all(ufα.≈ufθ) @test all(vf.≈ 1/(γ*dt) * (ufα-u0) + (1-1/γ)*v0) - # GeneralizedAlpha ∂tt test op = ODEOperatorMock{Float64,Nonlinear}(0.0,0.0,0.0,2) -ls = LUSolver() γ = 0.5 β = 0.25 ρ∞ = 1.0 # Equivalent to Newmark(0.5, 0.25) diff --git a/test/ODEsTests/TransientFEsTests/TransientBlockMultiFieldStyleTests.jl b/test/ODEsTests/TransientFEsTests/TransientBlockMultiFieldStyleTests.jl new file mode 100644 index 000000000..6c6e4dbe3 --- /dev/null +++ b/test/ODEsTests/TransientFEsTests/TransientBlockMultiFieldStyleTests.jl @@ -0,0 +1,103 @@ +module TransientBlockMultiFieldStyleTests +using Test, BlockArrays, SparseArrays, LinearAlgebra + +using Gridap +using Gridap.FESpaces, Gridap.ReferenceFEs, Gridap.MultiField +using Gridap.ODEs.TransientFETools + +function main(n_spaces,mfs,weakform,Ω,dΩ,U,V) + mass, biform, liform = weakform + res(t,x,y) = mass(t,∂t(x),y) + biform(t,x,y) - liform(t,y) + jac(t,x,dx,y) = biform(t,dx,y) + jac_t(t,xt,dxt,y) = mass(t,dxt,y) + + ############################################################################################ + # Normal assembly + + Y = MultiFieldFESpace(fill(V,n_spaces)) + X = TransientMultiFieldFESpace(fill(U,n_spaces)) + + u = get_trial_fe_basis(X(0.0)) + v = get_fe_basis(Y) + uₜ = TransientCellField(u,(u,)) + + matdata_jac = collect_cell_matrix(X(0),Y,jac(0,uₜ,u,v)) + matdata_jac_t = collect_cell_matrix(X(0),Y,jac_t(0,uₜ,u,v)) + matdata_jacs = (matdata_jac,matdata_jac_t) + matdata = TransientFETools._vcat_matdata(matdata_jacs) + vecdata = collect_cell_vector(Y,liform(0,v)) + + assem = SparseMatrixAssembler(X(0),Y) + A1 = assemble_matrix(assem,matdata) + b1 = assemble_vector(assem,vecdata) + + ############################################################################################ + # Block MultiFieldStyle + + Yb = MultiFieldFESpace(fill(V,n_spaces);style=mfs) + Xb = TransientMultiFieldFESpace(fill(U,n_spaces);style=mfs) + test_fe_space(Yb) + test_fe_space(Xb(0)) + + ub = get_trial_fe_basis(Xb(0)) + vb = get_fe_basis(Yb) + ubₜ = TransientCellField(ub,(ub,)) + + bmatdata_jac = collect_cell_matrix(Xb(0),Yb,jac(0,ubₜ,ub,vb)) + bmatdata_jac_t = collect_cell_matrix(Xb(0),Yb,jac_t(0,ubₜ,ub,vb)) + bmatdata_jacs = (bmatdata_jac,bmatdata_jac_t) + bmatdata = TransientFETools._vcat_matdata(bmatdata_jacs) + bvecdata = collect_cell_vector(Yb,liform(0,vb)) + + ############################################################################################ + # Block Assembly + + assem_blocks = SparseMatrixAssembler(Xb,Yb) + + A1_blocks = assemble_matrix(assem_blocks,bmatdata) + b1_blocks = assemble_vector(assem_blocks,bvecdata) + @test A1 ≈ A1_blocks + @test b1 ≈ b1_blocks + + y1_blocks = similar(b1_blocks) + mul!(y1_blocks,A1_blocks,b1_blocks) + y1 = similar(b1) + mul!(y1,A1,b1) + @test y1_blocks ≈ y1 + + A3_blocks = allocate_matrix(assem_blocks,bmatdata) + b3_blocks = allocate_vector(assem_blocks,bvecdata) + assemble_matrix!(A3_blocks,assem_blocks,bmatdata) + assemble_vector!(b3_blocks,assem_blocks,bvecdata) + @test A3_blocks ≈ A1 + @test b3_blocks ≈ b1_blocks + +end + +############################################################################################ + +sol(x,t) = sum(x) +sol(t::Real) = x->sol(x,t) + +model = CartesianDiscreteModel((0.0,1.0,0.0,1.0),(5,5)) +Ω = Triangulation(model) + +reffe = LagrangianRefFE(Float64,QUAD,1) +V = FESpace(Ω, reffe; dirichlet_tags="boundary") +U = TransientTrialFESpace(V,sol) + +dΩ = Measure(Ω, 2) +mass2(t,(u1t,u2t),(v1,v2)) = ∫(u1t⋅v1)*dΩ +biform2(t,(u1,u2),(v1,v2)) = ∫(∇(u1)⋅∇(v1) + u2⋅v2 - u1⋅v2)*dΩ +liform2(t,(v1,v2)) = ∫(v1 - v2)*dΩ +mass3(t,(u1t,u2t,u3t),(v1,v2,v3)) = ∫(u1t⋅v1)*dΩ +biform3(t,(u1,u2,u3),(v1,v2,v3)) = ∫(∇(u1)⋅∇(v1) + u2⋅v2 - u1⋅v2 - u3⋅v2 - u2⋅v3)*dΩ +liform3(t,(v1,v2,v3)) = ∫(v1 - v2 + 2.0*v3)*dΩ + +for (n_spaces,weakform) in zip([2,3],[(mass2,biform2,liform2),(mass3,biform3,liform3)]) + for mfs in [BlockMultiFieldStyle(),BlockMultiFieldStyle(2,(1,n_spaces-1))] + main(n_spaces,mfs,weakform,Ω,dΩ,U,V) + end +end + +end # module diff --git a/test/ODEsTests/TransientFEsTests/TransientFEOperatorsTests.jl b/test/ODEsTests/TransientFEsTests/TransientFEOperatorsTests.jl index 1574ceac3..de08cddf8 100644 --- a/test/ODEsTests/TransientFEsTests/TransientFEOperatorsTests.jl +++ b/test/ODEsTests/TransientFEsTests/TransientFEOperatorsTests.jl @@ -14,67 +14,88 @@ using Gridap.FESpaces: get_algebraic_operator u(x,t) = (1.0-x[1])*x[1]*(1.0-x[2])*x[2]*(t+3.0) u(t::Real) = x -> u(x,t) v(x) = t -> u(x,t) - f(t) = x -> ∂t(u)(x,t)-Δ(u(t))(x) +∂tu(x,t) = ∂t(u)(x,t) +∂tu(t::Real) = x -> ∂tu(x,t) +# Domain and triangulations domain = (0,1,0,1) partition = (2,2) model = CartesianDiscreteModel(domain,partition) - order = 2 - reffe = ReferenceFE(lagrangian,Float64,order) V0 = FESpace( model, reffe, conformity=:H1, dirichlet_tags="boundary") - U = TransientTrialFESpace(V0,u) - Ω = Triangulation(model) degree = 2*order dΩ = Measure(Ω,degree) +# Affine FE operator a(u,v) = ∫(∇(v)⊙∇(u))dΩ m(u,v) = ∫(v*u)dΩ b(v,t) = ∫(v*f(t))dΩ - res(t,u,v) = a(u,v) + m(∂t(u),v) - b(v,t) +lhs(t,u,v) = m(∂t(u),v) +rhs(t,u,v) = b(v,t) - a(u,v) +irhs(t,u,v) = b(v,t) - a(u,v)#∫( -1.0*(∇(v)⊙∇(u)))dΩ +erhs(t,u,v) = ∫( 0.0*(∇(v)⊙∇(u)))dΩ#b(v,t) jac(t,u,du,v) = a(du,v) jac_t(t,u,dut,v) = m(dut,v) - op = TransientFEOperator(res,jac,jac_t,U,V0) +opRK = TransientRungeKuttaFEOperator(lhs,rhs,jac,jac_t,U,V0) +opIMEXRK = TransientIMEXRungeKuttaFEOperator(lhs,irhs,erhs,jac,jac_t,U,V0) +# Time stepping t0 = 0.0 tF = 1.0 dt = 0.1 - +# Initial solution U0 = U(0.0) uh0 = interpolate_everywhere(u(0.0),U0) +∂tuh0 = interpolate_everywhere(∂tu(0.0),U0) -ls = LUSolver() -ode_solver = ThetaMethod(ls,dt,θ) +function test_ode_solver(ode_solver,op,xh0) + sol_t = solve(ode_solver,op,xh0,t0,tF) -sol_t = solve(ode_solver,op,uh0,t0,tF) + l2(w) = w*w + + tol = 1.0e-6 + _t_n = t0 -l2(w) = w*w + for (uh_tn, tn) in sol_t + # global _t_n + _t_n += dt + @test tn≈_t_n + e = u(tn) - uh_tn + el2 = sqrt(sum( ∫(l2(e))dΩ )) + @test el2 < tol + end -tol = 1.0e-6 -_t_n = t0 + @test length( [uht for uht in sol_t] ) == ceil((tF - t0)/dt) -for (uh_tn, tn) in sol_t - global _t_n - _t_n += dt - @test tn≈_t_n - e = u(tn) - uh_tn - el2 = sqrt(sum( ∫(l2(e))dΩ )) - @test el2 < tol end -@test length( [uht for uht in sol_t] ) == ceil((tF - t0)/dt) +# Linear solver +ls = LUSolver() +# ODE solvers +ode_solvers = [] +push!(ode_solvers,(ThetaMethod(ls,dt,θ),op,uh0)) +push!(ode_solvers,(BackwardEuler(ls,dt),op,uh0)) +push!(ode_solvers,(MidPoint(ls,dt),op,uh0)) +push!(ode_solvers,(GeneralizedAlpha(ls,dt,1.0),op,(uh0,∂tuh0))) +push!(ode_solvers,(RungeKutta(ls,ls,dt,:BE_1_0_1),opRK,uh0)) +push!(ode_solvers,(RungeKutta(ls,ls,dt,:CN_2_0_2),opRK,uh0)) +push!(ode_solvers,(RungeKutta(ls,ls,dt,:SDIRK_2_0_2),opRK,uh0)) +push!(ode_solvers,(IMEXRungeKutta(ls,ls,dt,:IMEX_FE_BE_2_0_1),opIMEXRK,uh0)) +for ode_solver in ode_solvers + test_ode_solver(ode_solver...) +end # u0 = get_free_dof_values(uh0) @@ -88,6 +109,8 @@ nl_cache = nothing # tf = t0+dt +# Nonlinear ThetaMethod +ode_solver = ThetaMethod(ls,dt,θ) ode_solver.θ == 0.0 ? dtθ = dt : dtθ = dt*ode_solver.θ tθ = t0+dtθ ode_cache = update_cache!(ode_cache,odeop,tθ) @@ -129,14 +152,14 @@ function extract_matrix_vector(a,fst) return Ast, bst end -A,rhs = extract_matrix_vector(a,fst) +A,vec = extract_matrix_vector(a,fst) gst(x) = u(tf)(x) m(u,v) = ∫(u*v)dΩ M,_ = extract_matrix_vector(m,gst) -@test rhs ≈ h +@test vec ≈ h @test A+M/(θ*dt) ≈ K rhs