Merge pull request #965 from gridap/clean-rk-solvers

Refactoring of GridapODEs

src/Algebra/Algebra.jl

@@ -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

src/Algebra/AlgebraInterfaces.jl

@@ -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
 #

src/Arrays/Interface.jl

@@ -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)
@@ -272,5 +282,3 @@ function test_array(
   end
   true
 end
-
-

src/Exports.jl

@@ -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
 

src/MultiField/MultiFieldCellFields.jl

@@ -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)

src/MultiField/MultiFieldFEFunctions.jl

@@ -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)