Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Refactoring of GridapODEs #965

Merged
merged 42 commits into from
Mar 18, 2024
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
42 commits
Select commit Hold shift + click to select a range
f6ff614
Generalised ButcherTableaus and RungeKutta (not IMEX for now), furthe…
AlexandreMagueresse Nov 20, 2023
30801ea
Refactored ODETools and their tests (only interfaces for now)
AlexandreMagueresse Nov 21, 2023
2ceded7
Improved tests for ODETools, imported TransientFETools from current code
AlexandreMagueresse Nov 21, 2023
ed48d0e
Cleaned TransientFETools. MassLinear case and tests remaining
AlexandreMagueresse Nov 22, 2023
70d98e7
Restructured Gridap.ODEs and started to rewrite tests for the parts r…
AlexandreMagueresse Nov 22, 2023
489b5b5
Finished tests for the ODEs interfaces
AlexandreMagueresse Nov 22, 2023
37905b9
Cleaned and tested the Theta-method family solvers (added BackwardEuler)
AlexandreMagueresse Nov 22, 2023
3949ac9
Cleaned and tested Tableaus
AlexandreMagueresse Nov 22, 2023
3f3f9e3
Started the implementation of MassLinearFEOperator
AlexandreMagueresse Nov 22, 2023
bbc562b
First version of TransientMassLinearFEOperator, made documentation an…
AlexandreMagueresse Nov 27, 2023
8b8e018
Finished implementing and testing MassLinearODE and LinearODE operato…
AlexandreMagueresse Nov 28, 2023
d4724d4
Rewrote ForwardEuler and ThetaMethod in terms of the slope and rely o…
AlexandreMagueresse Nov 28, 2023
ebb4c5a
Introduced a DiscreteODEOperator interface for all solvers and wrote …
AlexandreMagueresse Nov 29, 2023
c2b9596
Implemented Diagonally-Implicit Runge-Kutta, defaulting to a LinearDi…
AlexandreMagueresse Nov 29, 2023
6c40b7c
Rewrote GeneralizedAlpha for first-order ODEs
AlexandreMagueresse Nov 29, 2023
dc50dcb
Extensive renaming, consistent types for first- and higher-order ODEs
AlexandreMagueresse Dec 2, 2023
b24f074
Rewrote Generalised alpha in terms of the highest derivative, wrote t…
AlexandreMagueresse Dec 2, 2023
1e808b7
Implemented and tested IMEXODEOperator and TransientIMEXFEOperator
AlexandreMagueresse Dec 3, 2023
21827cf
Replaced similar by zero, added fine-grained filters for the residual…
AlexandreMagueresse Dec 4, 2023
844a1e7
Finished implementing and testing IMEX RK, changed MassLinear to
AlexandreMagueresse Dec 4, 2023
2d15f87
Implemented IMEX in the general semilinear case, added and merged tes…
AlexandreMagueresse Dec 5, 2023
cc4e0e4
Passed all previously working tests (AD is still not fully working)
AlexandreMagueresse Dec 5, 2023
c5595fd
Redefined the constructors of TransientFEOperators and fixed some min…
AlexandreMagueresse Dec 6, 2023
57136d3
Removed filters for residuals (contribution of dirichlet values), arb…
AlexandreMagueresse Dec 7, 2023
3880c4a
Fixed behaviour of constant bilinear forms and extended tests
AlexandreMagueresse Dec 8, 2023
5113ca4
Pull and merged with master
AlexandreMagueresse Dec 8, 2023
fa8a3f1
Major fix: changed behaviour of linear solvers on s (in general the j…
AlexandreMagueresse Dec 8, 2023
9088da8
Merge branch 'master' of github.com:gridap/Gridap.jl into clean-rk-so…
AlexandreMagueresse Dec 8, 2023
a337e1f
Refactoring of ODESolvers to follow the framework of general linear m…
AlexandreMagueresse Dec 17, 2023
25f4748
Wrote some documentation and reactivated all the tests of the ODE module
AlexandreMagueresse Dec 18, 2023
5f502af
Fixed LaTeX syntax for github
AlexandreMagueresse Dec 18, 2023
91b5b5f
Aggregate residuals and jacobians into DomainContributions instead of…
AlexandreMagueresse Dec 19, 2023
9ab1e7b
Bug fix: Automatic Differentiation with MultiField
AlexandreMagueresse Dec 19, 2023
39e0599
Replaced TransientMultiFieldFESpace with MultiFieldFESpace
AlexandreMagueresse Dec 20, 2023
c26edc6
Fixed constructors of TransientFEOperator and only assemble jacobian …
AlexandreMagueresse Dec 20, 2023
b8ea043
Removed unnecessary specialisations for time derivatives
AlexandreMagueresse Dec 20, 2023
adc18df
Activating all tests back
AlexandreMagueresse Dec 20, 2023
288ef01
Compute constant jacobians within matrices with the same sparsity pat…
AlexandreMagueresse Dec 21, 2023
2cb91fd
Turned some tests back on
AlexandreMagueresse Dec 22, 2023
219c208
Merge branch 'master' of github.com:gridap/Gridap.jl into clean-rk-so…
JordiManyer Jan 13, 2024
e29dfc4
Merge branch 'master' of github.com:gridap/Gridap.jl into clean-rk-so…
amartinhuertas Jan 29, 2024
f578f90
Update ODEs.md
AlexandreMagueresse Mar 18, 2024
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
57 changes: 0 additions & 57 deletions docs/ODEs.md

This file was deleted.

452 changes: 452 additions & 0 deletions docs/src/ODEs.md

Large diffs are not rendered by default.

1 change: 1 addition & 0 deletions src/Algebra/Algebra.jl
Original file line number Diff line number Diff line change
Expand Up @@ -28,6 +28,7 @@ export allocate_in_domain
export allocate_in_range
export add_entries!
export muladd!
export axpy_entries!
export nz_counter
export nz_allocation
export create_from_nz
Expand Down
54 changes: 54 additions & 0 deletions src/Algebra/AlgebraInterfaces.jl
Original file line number Diff line number Diff line change
Expand Up @@ -244,6 +244,60 @@ else
end
end

"""
axpy_entries!(α::Number, A::T, B::T) where {T<: AbstractMatrix} -> T

Efficient implementation of axpy! for sparse matrices.
"""
function axpy_entries!(α::Number, A::T, B::T) where {T<:AbstractMatrix}
iszero(α) && return B

axpy!(α, A, B)
B
end

# For sparse matrices, it is surprisingly quicker to call `@. B += α * A` than
# `axpy!(α, A, B)`.` Calling axpy! on the nonzero values of A and B is the most
# efficient approach but this is only possible when A and B have the same
# sparsity pattern. The checks add some non-negligible overhead so we make them
# optional by adding a keyword.
const cannot_axpy_entries_msg = """
It is only possible to efficiently add two sparse matrices that have the same
sparsity pattern.
"""

function axpy_entries!(
α::Number, A::T, B::T;
check::Bool=true
) where {T<:SparseMatrixCSC}
iszero(α) && return B

if check
msg = cannot_axpy_entries_msg
@check rowvals(A) == rowvals(B) msg
@check all(nzrange(A, j) == nzrange(B, j) for j in axes(A, 2)) msg
end

axpy!(α, nonzeros(A), nonzeros(B))
B
end

function axpy_entries!(
α::Number, A::T, B::T;
check::Bool=true
) where {T<:Union{SparseMatrixCSR,SymSparseMatrixCSR}}
iszero(α) && return B

if check
msg = cannot_axpy_entries_msg
@check colvals(A) == colvals(B) msg
@check all(nzrange(A, j) == nzrange(B, j) for j in axes(A, 1)) msg
end

axpy!(α, nonzeros(A), nonzeros(B))
B
end

#
# Some API associated with assembly routines
#
Expand Down
12 changes: 10 additions & 2 deletions src/Arrays/Interface.jl
Original file line number Diff line number Diff line change
Expand Up @@ -168,6 +168,16 @@ function testvalue(::Type{T}) where T<:AbstractArray{E,N} where {E,N}
similar(T,tfill(0,Val(N))...)
end

# When the jacobian of a residual is obtained through automatic differentiation,
# the return type is BlockArray{<:SubArray} and the behaviour of testvalue
# does not allow broadcasting operations between BlockArray{<:AbstractMatrix}
# and BlockArray{<:SubArray}. This function returns a matrix of size a
# P-dimensional array where each dimension has length 0, i.e., (0, ..., 0).
function testvalue(::Type{<:SubArray{T,P,AT}}) where {T,P,AT}
a = testvalue(AT)
return SubArray(a, ntuple(_ -> 0:-1, P))
end

function testvalue(::Type{T}) where T<:Transpose{E,A} where {E,A}
a = testvalue(A)
Transpose(a)
Expand Down Expand Up @@ -272,5 +282,3 @@ function test_array(
end
true
end


20 changes: 19 additions & 1 deletion src/Exports.jl
Original file line number Diff line number Diff line change
Expand Up @@ -182,7 +182,25 @@ using Gridap.CellData: ∫; export ∫
@publish Visualization createpvd
@publish Visualization savepvd

include("ODEs/Exports.jl")
@publish ODEs ∂t
@publish ODEs ∂tt
@publish ODEs ForwardEuler
@publish ODEs ThetaMethod
@publish ODEs MidPoint
@publish ODEs BackwardEuler
@publish ODEs GeneralizedAlpha1
@publish ODEs ButcherTableau
@publish ODEs available_tableaus
@publish ODEs RungeKutta
# @publish ODEs GeneralizedAlpha2
# @publish ODEs Newmark
@publish ODEs TransientTrialFESpace
@publish ODEs TransientMultiFieldFESpace
@publish ODEs TransientFEOperator
@publish ODEs TransientIMEXFEOperator
@publish ODEs TransientSemilinearFEOperator
@publish ODEs TransientQuasilinearFEOperator
@publish ODEs TransientLinearFEOperator

# Deprecated / removed

Expand Down
1 change: 1 addition & 0 deletions src/MultiField/MultiFieldCellFields.jl
Original file line number Diff line number Diff line change
Expand Up @@ -39,3 +39,4 @@ num_fields(a::MultiFieldCellField) = length(a.single_fields)
Base.getindex(a::MultiFieldCellField,i::Integer) = a.single_fields[i]
Base.iterate(a::MultiFieldCellField) = iterate(a.single_fields)
Base.iterate(a::MultiFieldCellField,state) = iterate(a.single_fields,state)
Base.length(a::MultiFieldCellField) = num_fields(a)
1 change: 1 addition & 0 deletions src/MultiField/MultiFieldFEFunctions.jl
Original file line number Diff line number Diff line change
Expand Up @@ -67,3 +67,4 @@ num_fields(m::MultiFieldFEFunction) = length(m.single_fe_functions)
Base.iterate(m::MultiFieldFEFunction) = iterate(m.single_fe_functions)
Base.iterate(m::MultiFieldFEFunction,state) = iterate(m.single_fe_functions,state)
Base.getindex(m::MultiFieldFEFunction,field_id::Integer) = m.single_fe_functions[field_id]
Base.length(m::MultiFieldFEFunction) = num_fields(m)
30 changes: 0 additions & 30 deletions src/ODEs/Exports.jl

This file was deleted.

Loading
Loading