From d4724d4d17afecacb850e39b8af4211b8f7fe9b5 Mon Sep 17 00:00:00 2001 From: Alexandre Magueresse Date: Tue, 28 Nov 2023 23:27:49 +1100 Subject: [PATCH] Rewrote ForwardEuler and ThetaMethod in terms of the slope and rely on an AffineOperator whenever possible, removed BackwardEuler --- src/ODEs/ODEOperatorsFromFEOperators.jl | 41 +++-- src/ODEs/ODESolvers.jl | 49 +++-- src/ODEs/ODESolvers/BackwardEuler.jl | 100 ----------- src/ODEs/ODESolvers/ForwardEuler.jl | 109 ++++++------ src/ODEs/ODESolvers/ThetaMethod.jl | 167 +++++++++++++----- src/ODEs/ODEs.jl | 2 +- .../ThetaMethodSolversFamilyTests.jl | 24 ++- test/ODEsTests/runtests.jl | 4 +- 8 files changed, 256 insertions(+), 240 deletions(-) delete mode 100644 src/ODEs/ODESolvers/BackwardEuler.jl diff --git a/src/ODEs/ODEOperatorsFromFEOperators.jl b/src/ODEs/ODEOperatorsFromFEOperators.jl index aa6fdd361..9c9ded652 100644 --- a/src/ODEs/ODEOperatorsFromFEOperators.jl +++ b/src/ODEs/ODEOperatorsFromFEOperators.jl @@ -131,7 +131,7 @@ end function Algebra.residual!( r::AbstractVector, op::ODEOpFromFEOp{MassLinearODE}, t::Real, us::Tuple{Vararg{AbstractVector}}, - cache + cache; include_highest::Bool=true ) mass_stored = !isnothing(cache.jacs[end]) forcing_stored = !isnothing(cache.forcing) @@ -142,13 +142,16 @@ function Algebra.residual!( end fill!(r, zero(eltype(r))) - if mass_stored - mul!(cache.jacvec, cache.jacs[end], us[end]) - axpy!(1, cache.jacvec, r) - else - mass = get_mass(op.fe_op) - vecdata = collect_cell_vector(V, mass(t, uh, v)) - assemble_vector_add!(r, get_assembler(op.fe_op), vecdata) + + if include_highest + if mass_stored + mul!(cache.jacvec, cache.jacs[end], us[end]) + axpy!(1, cache.jacvec, r) + else + mass = get_mass(op.fe_op) + vecdata = collect_cell_vector(V, mass(t, uh, v)) + assemble_vector_add!(r, get_assembler(op.fe_op), vecdata) + end end if forcing_stored @@ -165,7 +168,7 @@ end function Algebra.residual!( r::AbstractVector, op::ODEOpFromFEOp{LinearODE}, t::Real, us::Tuple{Vararg{AbstractVector}}, - cache + cache; include_highest::Bool=true ) forms_stored = !any(isnothing, cache.jacs) forcing_stored = !isnothing(cache.forcing) @@ -175,9 +178,11 @@ function Algebra.residual!( uh = _make_uh_from_us(op, us, cache.Us) end - order = get_order(op) fill!(r, zero(eltype(r))) - for k in 0:order + + order = get_order(op) + order_max = include_highest ? order : order - 1 + for k in 0:order_max if !isnothing(cache.jacs[k+1]) mul!(cache.jacvec, cache.jacs[k+1], us[k+1]) axpy!(1, cache.jacvec, r) @@ -225,9 +230,7 @@ function Algebra.jacobian!( cache ) if !isnothing(cache.jacs[k+1]) - # TODO optimise the sum of sparse matrices - # This is surprisingly better than axpy!(γ, cache.jacs[k+1], J) - @. J += γ * cache.jacs[k+1] + axpy_sparse!(γ, cache.jacs[k+1], J) else uh = _make_uh_from_us(op, us, cache.Us) matdata = _matdata_jacobian(op.fe_op, t, uh, k, γ) @@ -253,9 +256,7 @@ function jacobians!( γ = γs[k+1] if !iszero(γ) if !isnothing(cache.jacs[k+1]) - # TODO optimise the sum of sparse matrices - # This is surprisingly better than axpy!(γ, cache.jacs[k+1], J) - @. J += γ * cache.jacs[k+1] + axpy_sparse!(γ, cache.jacs[k+1], J) else _matdata = (_matdata..., _matdata_jacobian(op.fe_op, t, uh, k, γ)) end @@ -329,3 +330,9 @@ function _make_uh_from_us(op, us, Us) end TransientCellFieldType(u, dus) end + +function axpy_sparse!(α, A, B) + # TODO optimise the sum of sparse matrices + # This is surprisingly better than axpy!(α, A, B) + @. B += α * A +end diff --git a/src/ODEs/ODESolvers.jl b/src/ODEs/ODESolvers.jl index 717485eb3..6038a0009 100644 --- a/src/ODEs/ODESolvers.jl +++ b/src/ODEs/ODESolvers.jl @@ -85,28 +85,47 @@ end # Import solvers # ################## """ - _discrete_time_derivative!( - u̇::AbstractVector, - u0::AbstractVector, u::AbstractVector, - dt::Real + _v_from_u( + v::AbstractVector, + u::AbstractVector, u0::AbstractVector, dt::Real ) -> AbstractVector -Compute the discrete time derivative `u̇ = (u - u0) / dt`. +Safely write `(u - u0) / dt` into `v`. """ -function _discrete_time_derivative!( - u̇::AbstractVector, - u0::AbstractVector, u::AbstractVector, - dt::Real +function _v_from_u( + v::AbstractVector, + u::AbstractVector, u0::AbstractVector, dt::Real ) - copy!(u̇, u) - axpy!(-1, u0, u̇) - rdiv!(u̇, dt) - u̇ + if u !== v + copy!(v, u) + end + axpy!(-1, u0, v) + rdiv!(v, dt) + v end -include("ODESolvers/ForwardEuler.jl") +""" + _u_from_v!( + u::AbstractVector, + u0::AbstractVector, dt::Real, v::AbstractVector + ) -> AbstractVector + +Safely write `u0 + dt * v` into `u`. +""" +function _u_from_v!( + u::AbstractVector, + u0::AbstractVector, dt::Real, v::AbstractVector +) + if u === v + axpby!(1, u0, dt, u) + else + copy!(u, u0) + axpy!(dt, v, u) + end + u +end -include("ODESolvers/BackwardEuler.jl") +include("ODESolvers/ForwardEuler.jl") include("ODESolvers/ThetaMethod.jl") diff --git a/src/ODEs/ODESolvers/BackwardEuler.jl b/src/ODEs/ODESolvers/BackwardEuler.jl deleted file mode 100644 index 9d5c1874f..000000000 --- a/src/ODEs/ODESolvers/BackwardEuler.jl +++ /dev/null @@ -1,100 +0,0 @@ -""" - struct BackwardEuler <: ODESolver end - -BackwardEuler ODE solver. -""" -struct BackwardEuler <: ODESolver - nls::NonlinearSolver - dt::Float64 -end - -function solve_step!( - uF::AbstractVector, - solver::BackwardEuler, op::ODEOperator, - u0::AbstractVector, t0::Real, - cache -) - if isnothing(cache) - ode_cache = allocate_cache(op) - u̇F = similar(u0) - nls_cache = nothing - else - ode_cache, u̇F, nls_cache = cache - end - - dt = solver.dt - tF = t0 + dt - - # Update Dirichlet boundary conditions - ode_cache = update_cache!(ode_cache, op, tF) - - # Create and solve discrete ODE operator - nl_op = BackwardEulerSolverOperator(op, ode_cache, tF, dt, u0, u̇F) - nls_cache = solve!(uF, solver.nls, nl_op, nls_cache) - - # Update cache - cache = (ode_cache, u̇F, nls_cache) - - (uF, tF, cache) -end - -""" - struct BackwardEulerSolverOperator <: NonlinearOperator end - -Backward Euler operator at a given time step, i.e. -```math -residual(t_n+1, u_n+1, (u_n+1 - u_n) / dt) = 0. -``` -""" -struct BackwardEulerSolverOperator <: NonlinearOperator - ode_op::ODEOperator - ode_cache - tF::Float64 - dt::Float64 - u0::AbstractVector - u̇F::AbstractVector -end - -function Algebra.zero_initial_guess(op::BackwardEulerSolverOperator) - u = similar(op.u0) - fill!(u, zero(eltype(u))) - u -end - -function Algebra.allocate_residual( - op::BackwardEulerSolverOperator, - u::AbstractVector -) - allocate_residual(op.ode_op, op.tF, (u, u), op.ode_cache) -end - -function Algebra.residual!( - r::AbstractVector, - op::BackwardEulerSolverOperator, - u::AbstractVector -) - uF = u - u0, u̇F, dt = op.u0, op.u̇F, op.dt - _discrete_time_derivative!(u̇F, u0, uF, dt) - residual!(r, op.ode_op, op.tF, (uF, u̇F), op.ode_cache) -end - -function Algebra.allocate_jacobian( - op::BackwardEulerSolverOperator, - u::AbstractVector -) - allocate_jacobian(op.ode_op, op.tF, (u, u), op.ode_cache) -end - -function Algebra.jacobian!( - J::AbstractMatrix, - op::BackwardEulerSolverOperator, - u::AbstractVector -) - fillstored!(J, zero(eltype(J))) - - uF = u - u0, u̇F, dt = op.u0, op.u̇F, op.dt - _discrete_time_derivative!(u̇F, u0, uF, dt) - jacobians!(J, op.ode_op, op.tF, (uF, u̇F), (1, inv(op.dt)), op.ode_cache) -end diff --git a/src/ODEs/ODESolvers/ForwardEuler.jl b/src/ODEs/ODESolvers/ForwardEuler.jl index 4781cbe40..97155866a 100644 --- a/src/ODEs/ODESolvers/ForwardEuler.jl +++ b/src/ODEs/ODESolvers/ForwardEuler.jl @@ -15,12 +15,11 @@ function solve_step!( cache ) if isnothing(cache) - ode_cache = allocate_cache(ode_op) - u̇F = similar(u0) + ode_cache = allocate_cache(ode_op, t0, (u0, u0)) nls_cache = nothing - solver_cache = allocate_solver_cache(ode_op, ode_cache, t0, u0) + solver_cache = allocate_solver_cache(solver, ode_op, ode_cache, t0, u0) else - ode_cache, u̇F, nls_cache, solver_cache = cache + ode_cache, nls_cache, solver_cache = cache end dt = solver.dt @@ -30,134 +29,138 @@ function solve_step!( ode_cache = update_cache!(ode_cache, ode_op, t0) # Create and solve discrete ODE operator - nl_op = ForwardEulerSolverOperator( + nl_op = ForwardEulerOperator( ode_op, ode_cache, solver_cache, - t0, dt, u0, u̇F + t0, u0 ) - nls_cache = solve!(uF, solver.nls, nl_op, nls_cache) + v = uF + nls_cache = solve!(v, solver.nls, nl_op, nls_cache) + uF = _u_from_v!(uF, u0, dt, v) # Update cache - cache = (ode_cache, u̇F, nls_cache, solver_cache) + cache = (ode_cache, nls_cache, solver_cache) (uF, tF, cache) end function allocate_solver_cache( - ode_op::ODEOperator, ode_cache, + solver::ForwardEuler, ode_op::ODEOperator, ode_cache, t::Real, u::AbstractVector ) nothing end function allocate_solver_cache( - ode_op::ODEOperator{<:AbstractMassLinearODE}, ode_cache, + solver::ForwardEuler, ode_op::ODEOperator{<:AbstractMassLinearODE}, ode_cache, t::Real, u::AbstractVector ) J = allocate_jacobian(ode_op, t, (u, u), ode_cache) r = allocate_residual(ode_op, t, (u, u), ode_cache) - J, r + (J, r) end """ - ForwardEulerSolverOperator( + ForwardEulerOperator( ode_op::ODEOperator, ode_cache, solver_cache, - t0::Real, dt::Real, u0::AbstractVector, u̇F::AbstractVector - ) -> Union{ForwardEulerSolverNonlinearOperator, ForwardEulerSolverLinearOperator} + t0::Real, u0::AbstractVector + ) -> Union{ForwardEulerNonlinearOperator, ForwardEulerLinearOperator} Return the linear or nonlinear Forward Euler operator, depending on the type of `ODEOperator`. """ -function ForwardEulerSolverOperator( +function ForwardEulerOperator( ode_op::ODEOperator, ode_cache, solver_cache, - t0::Real, dt::Real, u0::AbstractVector, u̇F::AbstractVector + t0::Real, u0::AbstractVector ) - ForwardEulerSolverNonlinearOperator( + ForwardEulerNonlinearOperator( ode_op, ode_cache, - t0, dt, u0, u̇F + t0, u0 ) end -function ForwardEulerSolverOperator( +function ForwardEulerOperator( ode_op::ODEOperator{<:AbstractMassLinearODE}, ode_cache, solver_cache, - t0::Real, dt::Real, u0::AbstractVector, u̇F::AbstractVector + t0::Real, u0::AbstractVector ) - ForwardEulerSolverLinearOperator( + ForwardEulerLinearOperator( ode_op, ode_cache, solver_cache, - t0, dt, u0, u̇F + t0, u0 ) end """ - struct ForwardEulerSolverNonlinearOperator <: NonlinearOperator + struct ForwardEulerNonlinearOperator <: NonlinearOperator Forward Euler nonlinear operator at a given time step, i.e. ```math residual(t_n, u_n, (u_n+1 - u_n) / dt) = 0. ``` """ -struct ForwardEulerSolverNonlinearOperator <: NonlinearOperator +struct ForwardEulerNonlinearOperator <: NonlinearOperator ode_op::ODEOperator ode_cache t0::Float64 - dt::Float64 u0::AbstractVector - u̇F::AbstractVector end -function Algebra.zero_initial_guess(op::ForwardEulerSolverNonlinearOperator) - u = similar(op.u0) - fill!(u, zero(eltype(u))) - u +function Algebra.zero_initial_guess(op::ForwardEulerNonlinearOperator) + v = similar(op.u0) + fill!(v, zero(eltype(v))) + v end function Algebra.allocate_residual( - op::ForwardEulerSolverNonlinearOperator, - u::AbstractVector + op::ForwardEulerNonlinearOperator, + v::AbstractVector ) - allocate_residual(op.ode_op, op.t0, (u, u), op.ode_cache) + allocate_residual(op.ode_op, op.t0, (op.u0, v), op.ode_cache) end function Algebra.residual!( r::AbstractVector, - op::ForwardEulerSolverNonlinearOperator, - u::AbstractVector + op::ForwardEulerNonlinearOperator, + v::AbstractVector ) - u0, u̇F, dt = op.u0, op.u̇F, op.dt - _discrete_time_derivative!(u̇F, u0, u, dt) - residual!(r, op.ode_op, op.t0, (u0, u̇F), op.ode_cache) + residual!(r, op.ode_op, op.t0, (op.u0, v), op.ode_cache) end function Algebra.allocate_jacobian( - op::ForwardEulerSolverNonlinearOperator, - u::AbstractVector + op::ForwardEulerNonlinearOperator, + v::AbstractVector ) - allocate_jacobian(op.ode_op, op.t0, (u, u), op.ode_cache) + allocate_jacobian(op.ode_op, op.t0, (op.u0, v), op.ode_cache) end function Algebra.jacobian!( J::AbstractMatrix, - op::ForwardEulerSolverNonlinearOperator, - u::AbstractVector + op::ForwardEulerNonlinearOperator, + v::AbstractVector ) fillstored!(J, zero(eltype(J))) - - u0, u̇F, dt = op.u0, op.u̇F, op.dt - _discrete_time_derivative!(u̇F, u0, u, dt) - jacobian!(J, op.ode_op, op.t0, (u0, u̇F), 1, inv(dt), op.ode_cache) + jacobian!(J, op.ode_op, op.t0, (op.u0, v), 1, 1, op.ode_cache) end """ - struct ForwardEulerSolverLinearOperator <: AffineOperator + ForwardEulerLinearOperator( + ode_op, ode_cache, solver_cache, + t0, u0 + ) -> AffineOperator Forward Euler linear operator at a given time step, i.e. -`mass(t_n, u_n) * (u_n+1 - u_n) / dt + res(t_n, u_n) = 0`. For simplicity, -we solve for `k_n = (u_n+1 - u_n) / dt` and we recover `u_n+1` from `k_n`, -`u_n` and `dt`. +```math +mass(t_n, u_n) * (u_n+1 - u_n) / dt + res(t_n, u_n) = 0. +``` """ -function ForwardEulerSolverLinearOperator( +function ForwardEulerLinearOperator( ode_op, ode_cache, solver_cache, - t0, dt, u0, u̇F + t0, u0 ) J, r = solver_cache + v = u0 + fillstored!(J, zero(eltype(J))) + jacobian!(J, ode_op, t0, (u0, v), 1, 1, ode_cache) + residual!(r, ode_op, t0, (u0, v), ode_cache, include_highest=false) + rmul!(r, -1) + AffineOperator(J, r) end diff --git a/src/ODEs/ODESolvers/ThetaMethod.jl b/src/ODEs/ODESolvers/ThetaMethod.jl index bd45b3a1e..bba64a4db 100644 --- a/src/ODEs/ODESolvers/ThetaMethod.jl +++ b/src/ODEs/ODESolvers/ThetaMethod.jl @@ -2,9 +2,6 @@ struct ThetaMethod <: ODESolver end θ-method ODE solver. -* θ = 0 Forward Euler -* θ = 1/2 Crank-Nicolson / MidPoint -* θ = 1 Backward Euler """ struct ThetaMethod <: ODESolver nls::NonlinearSolver @@ -12,17 +9,25 @@ struct ThetaMethod <: ODESolver θ::Float64 function ThetaMethod(nls, dt, θ) - if θ <= 0 + θ01 = clamp(θ, 0, 1) + if θ01 != θ + msg = """ + The parameter θ of ThetaMethod must be between zero and one. + Setting θ to $(θ01). + """ + println(msg) + end + + if iszero(θ01) ForwardEuler(nls, dt) - elseif θ >= 1 - BackwardEulerEuler(nls, dt) else - new(nls, dt, θ) + new(nls, dt, θ01) end end end MidPoint(nls, dt) = ThetaMethod(nls, dt, 0.5) +BackwardEuler(nls, dt) = ThetaMethod(nls, dt, 1) function solve_step!( uF::AbstractVector, @@ -31,11 +36,11 @@ function solve_step!( cache ) if isnothing(cache) - ode_cache = allocate_cache(ode_op) - u̇θ = similar(u0) + ode_cache = allocate_cache(ode_op, t0, (u0, u0)) nls_cache = nothing + solver_cache = allocate_solver_cache(solver, ode_op, ode_cache, t0, u0) else - ode_cache, u̇θ, nls_cache = cache + ode_cache, nls_cache, solver_cache = cache end dt, θ = solver.dt, solver.θ @@ -47,77 +52,149 @@ function solve_step!( ode_cache = update_cache!(ode_cache, ode_op, tθ) # Create and solve discrete ODE operator - nl_op = ThetaMethodSolverOperator(ode_op, ode_cache, tθ, dtθ, u0, u̇θ) - nls_cache = solve!(uF, solver.nls, nl_op, nls_cache) - - # The operator above solves for uθ - # Express uF in terms of uθ - axpy!(θ - 1, u0, uF) - rdiv!(uF, θ) + nl_op = ThetaMethodOperator( + ode_op, ode_cache, solver_cache, + tθ, u0, dtθ + ) + v = uF + nls_cache = solve!(v, solver.nls, nl_op, nls_cache) + uF = _u_from_v!(uF, u0, dt, v) # Update cache - cache = (ode_cache, u̇θ, nls_cache) + cache = (ode_cache, nls_cache, solver_cache) (uF, tF, cache) end +function allocate_solver_cache( + solver::ThetaMethod, ode_op::ODEOperator, ode_cache, + t::Real, u::AbstractVector +) + (similar(u),) +end + +function allocate_solver_cache( + solver::ThetaMethod, ode_op::ODEOperator{LinearODE}, ode_cache, + t::Real, u::AbstractVector +) + J = allocate_jacobian(ode_op, t, (u, u), ode_cache) + r = allocate_residual(ode_op, t, (u, u), ode_cache) + (J, r) +end + +""" + ThetaMethodOperator( + ode_op::ODEOperator, ode_cache, solver_cache, + tθ::Real, u0::AbstractVector, dtθ::Real + ) -> Union{ThetaMethodNonlinearOperator, ThetaMethodLinearOperator} + +Return the linear or nonlinear Theta Method operator, depending on the +type of `ODEOperator`. +""" +function ThetaMethodOperator( + ode_op::ODEOperator, ode_cache, solver_cache, + tθ::Real, u0::AbstractVector, dtθ::Real +) + u, = solver_cache + ThetaMethodNonlinearOperator( + ode_op, ode_cache, + tθ, u0, dtθ, u + ) +end + +function ThetaMethodOperator( + ode_op::ODEOperator{LinearODE}, ode_cache, solver_cache, + tθ::Real, u0::AbstractVector, dtθ::Real +) + ThetaMethodLinearOperator( + ode_op, ode_cache, solver_cache, + tθ, u0, dtθ + ) +end + """ - struct ThetaMethodSolverOperator <: NonlinearOperator end + struct ThetaMethodNonlinearOperator <: NonlinearOperator -Theta method operator at a given time step, i.e. +Theta method nonlinear operator at a given time step, i.e. ```math residual(t_n+θ, u_n+θ, (u_n+θ - u_n) / θ / dt) = 0. ``` """ -struct ThetaMethodSolverOperator <: NonlinearOperator +struct ThetaMethodNonlinearOperator <: NonlinearOperator ode_op::ODEOperator ode_cache tθ::Float64 - dtθ::Float64 u0::AbstractVector - u̇θ::AbstractVector + dtθ::Float64 + u::AbstractVector end -function Algebra.zero_initial_guess(op::ThetaMethodSolverOperator) - u = similar(op.u0) - fill!(u, zero(eltype(u))) - u +function Algebra.zero_initial_guess(op::ThetaMethodNonlinearOperator) + v = similar(op.u0) + fill!(v, zero(eltype(v))) + v end function Algebra.allocate_residual( - op::ThetaMethodSolverOperator, - u::AbstractVector + op::ThetaMethodNonlinearOperator, + v::AbstractVector ) - allocate_residual(op.ode_op, op.tθ, (u, u), op.ode_cache) + u = op.u + u = _u_from_v!(u, op.u0, op.dtθ, v) + allocate_residual(op.ode_op, op.tθ, (u, v), op.ode_cache) end function Algebra.residual!( r::AbstractVector, - op::ThetaMethodSolverOperator, - u::AbstractVector + op::ThetaMethodNonlinearOperator, + v::AbstractVector ) - uθ = u - u0, u̇θ, dtθ = op.u0, op.u̇θ, op.dtθ - _discrete_time_derivative!(u̇θ, u0, uθ, dtθ) - residual!(r, op.ode_op, op.tθ, (uθ, u̇θ), op.ode_cache) + u = op.u + u = _u_from_v!(u, op.u0, op.dtθ, v) + residual!(r, op.ode_op, op.tθ, (u, v), op.ode_cache) end function Algebra.allocate_jacobian( - op::ThetaMethodSolverOperator, - u::AbstractVector + op::ThetaMethodNonlinearOperator, + v::AbstractVector ) - allocate_jacobian(op.ode_op, op.tθ, (u, u), op.ode_cache) + u = op.u + u = _u_from_v!(u, op.u0, op.dtθ, v) + allocate_jacobian(op.ode_op, op.tθ, (u, v), op.ode_cache) end function Algebra.jacobian!( J::AbstractMatrix, - op::ThetaMethodSolverOperator, - u::AbstractVector + op::ThetaMethodNonlinearOperator, + v::AbstractVector ) fillstored!(J, zero(eltype(J))) + u = op.u + u = _u_from_v!(u, op.u0, op.dtθ, v) + jacobians!(J, op.ode_op, op.tθ, (u, v), (op.dtθ, 1), op.ode_cache) +end + +""" + ThetaMethodLinearOperator( + ode_op, ode_cache, solver_cache, + tθ, u0, dtθ + ) -> AffineOperator + +Theta method linear operator at a given time step, i.e. +```math +mass(t_n+θ) * (u_n+θ - u_n) / θ / dt + stiffness(t_n+θ) * u_n+θ + res(t_n+θ) = 0. +``` +""" +function ThetaMethodLinearOperator( + ode_op, ode_cache, solver_cache, + tθ, u0, dtθ +) + J, r = solver_cache + + fillstored!(J, zero(eltype(J))) + jacobians!(J, ode_op, tθ, (u0, u0), (dtθ, 1), ode_cache) + residual!(r, ode_op, tθ, (u0, u0), ode_cache, include_highest=false) + rmul!(r, -1) - uθ = u - u0, u̇θ, dtθ = op.u0, op.u̇θ, op.dtθ - _discrete_time_derivative!(u̇θ, u0, uθ, dtθ) - jacobians!(J, op.ode_op, op.tθ, (uθ, u̇θ), (1, inv(dtθ)), op.ode_cache) + AffineOperator(J, r) end diff --git a/src/ODEs/ODEs.jl b/src/ODEs/ODEs.jl index c21bc48b6..f5aafe35e 100644 --- a/src/ODEs/ODEs.jl +++ b/src/ODEs/ODEs.jl @@ -69,9 +69,9 @@ export solve_step! export solve export ForwardEuler -export BackwardEuler export ThetaMethod export MidPoint +export BackwardEuler export TableauType export ExplicitTableau diff --git a/test/ODEsTests/ODESolversAllTests/ThetaMethodSolversFamilyTests.jl b/test/ODEsTests/ODESolversAllTests/ThetaMethodSolversFamilyTests.jl index fcef60015..e0d738120 100644 --- a/test/ODEsTests/ODESolversAllTests/ThetaMethodSolversFamilyTests.jl +++ b/test/ODEsTests/ODESolversAllTests/ThetaMethodSolversFamilyTests.jl @@ -44,21 +44,31 @@ a(t, u, v) = ∫(∇(u) ⊙ ∇(v)) * dΩ b(t, v) = ∫(f(t) ⋅ v) * dΩ mass(t, u, v) = m(t, ∂t(u), v) -part_res(t, u, v) = a(t, u, v) - b(t, v) -full_res(t, u, v) = mass(t, u, v) + part_res(t, u, v) +stiffness(t, u, v) = a(t, u, v) +res(t, u, v) = mass(t, u, v) + stiffness(t, u, v) - b(t, v) +res_masslinear(t, u, v) = stiffness(t, u, v) - b(t, v) +res_linear(t, v) = (-1) * b(t, v) jac(t, u, du, v) = a(t, du, v) jac_t(t, u, dut, v) = m(t, dut, v) -op_nonlinear = TransientFEOperator(full_res, jac, jac_t, U, V) -op_masslinear = TransientMassLinearFEOperator(mass, part_res, jac, jac_t, U, V) -op_constmasslinear = TransientConstantMassLinearFEOperator( - mass, part_res, jac, jac_t, U, V +jacs_constant = (true, true) +op_nonlinear = TransientFEOperator( + res, jac, jac_t, U, V; + jacs_constant +) +op_masslinear = TransientMassLinearFEOperator( + mass, res_masslinear, jac, jac_t, U, V; + jacs_constant +) +op_linear = TransientLinearFEOperator( + mass, stiffness, res_linear, jac, jac_t, U, V; + jacs_constant ) ops = [ op_nonlinear, op_masslinear, - op_constmasslinear + op_linear ] # ODE solver diff --git a/test/ODEsTests/runtests.jl b/test/ODEsTests/runtests.jl index 2a13afa92..92c97debe 100644 --- a/test/ODEsTests/runtests.jl +++ b/test/ODEsTests/runtests.jl @@ -16,9 +16,9 @@ using Test # @time @testset "TransientFEOperators" begin include("TransientFEOperatorsTests.jl") end -@time @testset "TransientFESolutions" begin include("TransientFESolutionsTests.jl") end +# @time @testset "TransientFESolutions" begin include("TransientFESolutionsTests.jl") end -# @time @testset "ODESolversAll" begin include("ODESolversAllTests/runtests.jl") end +@time @testset "ODESolversAll" begin include("ODESolversAllTests/runtests.jl") end # @time @testset "TransientProblems" begin include("TransientProblemsTests/runtests.jl") end