Release Gridap 0.18

NEWS.md

@@ -5,7 +5,23 @@ All notable changes to this project will be documented in this file.
 The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/),
 and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 
-## [0.17.23] - 2024-01-28 
+## [0.18.0] - 2024-03-22
+
+### Breaking
+
+- ODE module extensive refactor. Breaking changes! See docs and PR for details. Since PR[965](https://github.com/gridap/Gridap.jl/pull/965).
+- Fixed name clash with `Statistics.mean`. Since PR[#988](https://github.com/gridap/Gridap.jl/pull/988).
+- Deprecated `SubVector` in favor of Julia's `view`. Since PR[#989](https://github.com/gridap/Gridap.jl/pull/989).
+
+### Added
+
+- Added some missing API methods to `Assemblers` and `MultiField`. Since PR[#985](https://github.com/gridap/Gridap.jl/pull/985).
+
+### Fixed
+
+- Fix when evaluating `\circ` operator with `CellState`. Since PR[#987](https://github.com/gridap/Gridap.jl/pull/987).
+
+## [0.17.23] - 2024-01-28
 
 ### Changed
 

Project.toml

@@ -1,7 +1,7 @@
 name = "Gridap"
 uuid = "56d4f2e9-7ea1-5844-9cf6-b9c51ca7ce8e"
 authors = ["Santiago Badia <[email protected]>", "Francesc Verdugo <[email protected]>", "Alberto F. Martin <[email protected]>"]
-version = "0.17.23"
+version = "0.18.0"
 
 [deps]
 AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"
@@ -26,6 +26,7 @@ Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 SparseMatricesCSR = "a0a7dd2c-ebf4-11e9-1f05-cf50bc540ca1"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
+Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 WriteVTK = "64499a7a-5c06-52f2-abe2-ccb03c286192"
 

docs/src/ODEs.md

@@ -91,7 +91,7 @@ A `TransientTrialFESpace` can be evaluated at any time derivative order, and the
 
 For example, the following creates a transient `FESpace` and evaluates its first two time derivatives.
 ```
-g(t::Real) = x -> x[1] + x[2] * t
+g(t) = x -> x[1] + x[2] * t
 V = FESpace(model, reffe, dirichlet_tags="boundary")
 U = TransientTrialFESpace (V, g)
 
@@ -155,6 +155,16 @@ TransientLinearFEOperator((stiffness, mass), res, U, V, constant_forms=(false, t
 ```
 If $\kappa$ is constant, the keyword `constant_forms` could be replaced by `(true, true)`.
 
+## The `TimeSpaceFunction` constructor
+Apply differential operators on a function that depends on time and space is somewhat cumbersome. Let `f` be a function of time and space, and `g(t) = x -> f(t, x)` (as in the prescription of the boundary conditions `g` above). Applying the operator $\partial_{t} - \Delta$  to `g` and evaluating at $(t, x)$ is written `∂t(g)(t)(x) - Δ(g(t))(x)`.
+
+The constructor `TimeSpaceFunction` allows for simpler notations: let `h = TimeSpaceFunction(g)`. The object `h` is a functor that supports the notations 
+* `op(h)`: a `TimeSpaceFunction` representing both `t -> x -> op(f)(t, x)` and `(t, x) -> op(f)(t, x)`,
+* `op(h)(t)`: a function of space representing `x -> op(f)(t, x)`
+* `op(h)(t, x)`: the quantity `op(f)(t, x)` (this notation is equivalent to `op(h)(t)(x)`),
+
+for all spatial and temporal differential operator, i.e. `op` in `(time_derivative, gradient, symmetric_gradient, divergence, curl, laplacian)` and their symbolic aliases (`∂t`, `∂tt`, `∇`, ...). The operator above applied to `h` and evaluated at `(t, x)` can be conveniently written `∂t(h)(t, x) - Δ(h)(t, x)`.
+
 ## Solver and solution
 The next step is to choose an ODE solver (see below for a full list) and specify the boundary conditions. The solution can then be iterated over until the final time is reached.
 
@@ -252,6 +262,8 @@ By looking at the behaviour of the stability function at infinity, we find that
 ## Generalised- $\alpha$ scheme for first-order ODEs
 This scheme relies on the state vector $\\{\boldsymbol{s}(t)\\} = \\{\boldsymbol{u}(t), \partial_{t} \boldsymbol{u}(t)\\}$. In particular, it needs a nontrivial starting procedure that evaluates $\partial_{t} \boldsymbol{u}(t_{0})$ by enforcing a zero residual at $t_{0}$. The finaliser can still return the first vector of the state vectors. For convenience, let $\partial_{t} \boldsymbol{u}\_{n}$ denote the approximation $\partial_{t} \boldsymbol{u}(t_{n})$.
 
+> Alternatively, the initial velocity can be provided manually: when calling `solve(odeslvr, tfeop, t0, tF, uhs0)`, set `uhs0 = (u0, v0, a0)` instead of `uhs0 = (u0, v0)`. This is useful when enforcing a zero initial residual would lead to a singular system.
+
 This method extends the $\theta$-method by considering the two-point quadrature rule $$\boldsymbol{u}(t_{n+1}) = \boldsymbol{u}\_{n} + \int_{t_{n}}^{t_{n+1}} \partial_{t} \boldsymbol{u}(t) \ \mathrm{d} t \approx \boldsymbol{u}\_{n} + h_{n} [(1 - \gamma) \partial_{t} \boldsymbol{u}(t_{n}) + \gamma \partial_{t} \boldsymbol{u}(t_{n+1})],$$
 where $0 \leq \gamma \leq 1$ is a free parameter. The question is now how to estimate $\partial_{t} \boldsymbol{u}(t_{n+1})$. This is achieved by enforcing a zero residual at $t_{n + \alpha_{F}} \doteq (1 - \alpha_{F}) t_{n} + \alpha_{F} t_{n+1}$, where $0 \leq \alpha_{F} \leq 1$ is another free parameter. The value of $\boldsymbol{u}$ at that time, $\boldsymbol{u}\_{n + \alpha_{F}}$, is obtained by the same linear combination of $\boldsymbol{u}$ at $t_{n}$ and $t_{n+1}$. Regarding $\partial_{t} \boldsymbol{u}$, it is taken as a linear combination weighted by another free parameter $0 < \alpha_{M} \leq 1$ of the time derivative at times $t_{n}$ and $t_{n+1}$. Note that $\alpha_{M}$ cannot be zero. Altogether, we have defined the discrete operators
 ```math
@@ -399,6 +411,8 @@ We note that the first column of the matrix and the first weight are all zero, s
 ## Generalised- $\alpha$ scheme for second-order ODEs
 This scheme relies on the state vector $\\{\boldsymbol{s}(t)\\} = \\{\boldsymbol{u}(t), \partial_{t} \boldsymbol{u}(t), \partial_{tt} \boldsymbol{u}(t)\\}$. It needs a nontrivial starting procedure that evaluates $\partial_{tt} \boldsymbol{u}(t_{0})$ by enforcing a zero residual at $t_{0}$. The finaliser can still return the first vector of the state vectors. For convenience, let $\partial_{tt} \boldsymbol{u}\_{n}$ denote the approximation $\partial_{tt} \boldsymbol{u}(t_{n})$.
 
+> The initial acceleration can alternatively be provided manually: when calling `solve(odeslvr, tfeop, t0, tF, uhs0)`, set `uhs0 = (u0, v0, a0)` instead of `uhs0 = (u0, v0)`. This is useful when enforcing a zero initial residual would lead to a singular system.
+
 This method is built out of the following update rule
 ```math
 \begin{align*}

src/Algebra/AlgebraInterfaces.jl

@@ -298,6 +298,13 @@ function axpy_entries!(
   B
 end
 
+function axpy_entries!(α::Number, A::T, B::T) where {T<:AbstractBlockMatrix}
+  map(blocks(A), blocks(B)) do a, b
+    axpy_entries!(α, a, b)
+  end
+  B
+end
+
 #
 # Some API associated with assembly routines
 #

src/Arrays/Arrays.jl

@@ -105,7 +105,6 @@ export find_local_index
 
 export IdentityVector
 
-export SubVector
 export pair_arrays
 export unpair_arrays
 
@@ -162,7 +161,7 @@ include("KeyToValMaps.jl")
 
 include("FilteredArrays.jl")
 
-include("SubVectors.jl")
+include("SubVectors.jl") # Deprecated
 
 include("ArrayPairs.jl")
 

src/Arrays/PrintOpTrees.jl

@@ -53,10 +53,6 @@ function get_children(n::TreeNode, a::LazyArray)
   (similar_tree_node(n,a.maps),map(i->similar_tree_node(n,i),a.args)...)
 end
 
-function get_children(n::TreeNode, a::SubVector)
-  (similar_tree_node(n,a.vector),)
-end
-
 function get_children(n::TreeNode, a::Reindex)
   (similar_tree_node(n,a.values),)
 end

src/Arrays/SubVectors.jl

@@ -5,6 +5,8 @@
       pini::Int
       pend::Int
     end
+
+  `SubVector` is deprecated, use `view` instead.
 """
 struct SubVector{T,A<:AbstractVector{T}} <: AbstractVector{T}
   vector::A

src/CellData/CellData.jl

@@ -34,6 +34,7 @@ import Gridap.Geometry: get_triangulation
 import Gridap.TensorValues: inner, outer, double_contraction, symmetric_part
 import LinearAlgebra: det, tr, cross, dot, ⋅, rmul!
 import Base: inv, abs, abs2, *, +, -, /, adjoint, transpose, real, imag, conj
+import Statistics: mean
 
 export gradient, ∇
 export ∇∇

src/FESpaces/Assemblers.jl

@@ -197,7 +197,7 @@ end
 
 """
 """
-function assemble_matrix_add!(mat,a::Assembler, matdata)
+function assemble_matrix_add!(mat,a::Assembler,matdata)
   @abstractmethod
 end
 
@@ -215,11 +215,11 @@ end
 
 """
 """
-function assemble_matrix_and_vector!(A,b,a::Assembler, data)
+function assemble_matrix_and_vector!(A,b,a::Assembler,data)
   @abstractmethod
 end
 
-function assemble_matrix_and_vector_add!(A,b,a::Assembler, data)
+function assemble_matrix_and_vector_add!(A,b,a::Assembler,data)
   @abstractmethod
 end
 
@@ -279,6 +279,17 @@ end
 # Some syntactic sugar for assembling from anonymous functions
 # and objects from which one can collect cell matrices/vectors
 
+function allocate_matrix(f::Function,a::Assembler,U::FESpace,V::FESpace)
+  v = get_fe_basis(V)
+  u = get_trial_fe_basis(U)
+  allocate_matrix(a,collect_cell_matrix(U,V,f(u,v)))
+end
+
+function allocate_vector(f::Function,a::Assembler,V::FESpace)
+  v = get_fe_basis(V)
+  allocate_vector(a,collect_cell_vector(V,f(v)))
+end
+
 function assemble_matrix(f::Function,a::Assembler,U::FESpace,V::FESpace)
   v = get_fe_basis(V)
   u = get_trial_fe_basis(U)
@@ -313,6 +324,14 @@ function assemble_matrix_and_vector!(f::Function,b::Function,M::AbstractMatrix,r
   assemble_matrix_and_vector!(M,r,a,collect_cell_matrix_and_vector(U,V,f(u,v),b(v)))
 end
 
+function assemble_matrix_and_vector!(
+  f::Function,b::Function,M::AbstractMatrix,r::AbstractVector,a::Assembler,U::FESpace,V::FESpace,uhd
+)
+  v = get_fe_basis(V)
+  u = get_trial_fe_basis(U)
+  assemble_matrix_and_vector!(M,r,a,collect_cell_matrix_and_vector(U,V,f(u,v),b(v),uhd))
+end
+
 function assemble_matrix(f,a::Assembler,U::FESpace,V::FESpace)
   assemble_matrix(a,collect_cell_matrix(U,V,f))
 end

src/MultiField/MultiFieldCellFields.jl

@@ -40,3 +40,8 @@ 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)
+
+function LinearAlgebra.dot(a::MultiFieldCellField,b::MultiFieldCellField)
+  @check num_fields(a) == num_fields(b)
+  return sum(map(dot,a.single_fields,b.single_fields))
+end

src/MultiField/MultiFieldFEFunctions.jl

@@ -68,3 +68,8 @@ 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)
+
+function LinearAlgebra.dot(a::MultiFieldFEFunction,b::MultiFieldFEFunction)
+  @check num_fields(a) == num_fields(b)
+  return sum(map(dot,a.single_fe_functions,b.single_fe_functions))
+end

src/MultiField/MultiFieldFESpaces.jl

@@ -137,6 +137,10 @@ function FESpaces.get_free_dof_ids(f::MultiFieldFESpace,::BlockMultiFieldStyle{N
   return BlockArrays.blockedrange(block_num_dofs)
 end
 
+function FESpaces.zero_dirichlet_values(f::MultiFieldFESpace)
+  map(zero_dirichlet_values,f.spaces)
+end
+
 FESpaces.get_dof_value_type(f::MultiFieldFESpace{MS,CS,V}) where {MS,CS,V} = eltype(V)
 
 FESpaces.get_vector_type(f::MultiFieldFESpace) = f.vector_type
@@ -262,6 +266,18 @@ function FESpaces.FEFunction(fe::MultiFieldFESpace, free_values)
   MultiFieldFEFunction(free_values,fe,blocks)
 end
 
+function FESpaces.FEFunction(
+  fe::MultiFieldFESpace, free_values::AbstractVector, dir_values::Vector{<:AbstractVector}
+)
+  @check length(dir_values) == num_fields(fe)
+  blocks = map(1:length(fe.spaces)) do i
+    free_values_i = restrict_to_field(fe,free_values,i)
+    dir_values_i  = dir_values[i]
+    FEFunction(fe.spaces[i],free_values_i,dir_values_i)
+  end
+  MultiFieldFEFunction(free_values,fe,blocks)
+end
+
 function FESpaces.EvaluationFunction(fe::MultiFieldFESpace, free_values)
   blocks = map(1:length(fe.spaces)) do i
     free_values_i = restrict_to_field(fe,free_values,i)
@@ -297,7 +313,7 @@ function  _restrict_to_field(f,
   offsets = _compute_field_offsets(U)
   pini = offsets[field] + 1
   pend = offsets[field] + num_free_dofs(U[field])
-  SubVector(free_values,pini,pend)
+  view(free_values,pini:pend)
 end
 
 function  _restrict_to_field(f,
@@ -316,7 +332,7 @@ function  _restrict_to_field(f,
   offsets = compute_field_offsets(f,mfs)
   pini = offsets[field] + 1
   pend = offsets[field] + num_free_dofs(U[field])
-  return SubVector(block_free_values,pini,pend)
+  return view(block_free_values,pini:pend)
 end
 
 """
@@ -596,7 +612,18 @@ function FESpaces.interpolate_everywhere(objects, fe::MultiFieldFESpace)
   for (field, (U,object)) in enumerate(zip(fe.spaces,objects))
     free_values_i = restrict_to_field(fe,free_values,field)
     dirichlet_values_i = zero_dirichlet_values(U)
-    uhi = interpolate_everywhere!(object, free_values_i,dirichlet_values_i,U)
+    uhi = interpolate_everywhere!(object,free_values_i,dirichlet_values_i,U)
+    push!(blocks,uhi)
+  end
+  MultiFieldFEFunction(free_values,fe,blocks)
+end
+
+function FESpaces.interpolate_everywhere!(objects,free_values::AbstractVector,dirichlet_values::Vector,fe::MultiFieldFESpace)
+  blocks = SingleFieldFEFunction[]
+  for (field, (U,object)) in enumerate(zip(fe.spaces,objects))
+    free_values_i = restrict_to_field(fe,free_values,field)
+    dirichlet_values_i = dirichlet_values[field]
+    uhi = interpolate_everywhere!(object,free_values_i,dirichlet_values_i,U)
     push!(blocks,uhi)
   end
   MultiFieldFEFunction(free_values,fe,blocks)

src/ODEs/ODESolvers/GeneralizedAlpha1.jl

@@ -43,18 +43,27 @@ function GeneralizedAlpha1(
   GeneralizedAlpha1(sysslvr, dt, αf, αm, γ)
 end
 
+# Default allocate_odecache without velocity
+function allocate_odecache(
+  odeslvr::GeneralizedAlpha1, odeop::ODEOperator,
+  t0::Real, us0::NTuple{1,AbstractVector}
+)
+  u0 = us0[1]
+  allocate_odecache(odeslvr, odeop, t0, (u0, u0))
+end
+
 ##################
 # Nonlinear case #
 ##################
 function allocate_odecache(
   odeslvr::GeneralizedAlpha1, odeop::ODEOperator,
-  t0::Real, us0::NTuple{1,AbstractVector}
+  t0::Real, us0::NTuple{2,AbstractVector}
 )
-  u0 = us0[1]
-  us0N = (u0, u0)
+  u0, v0 = us0[1], us0[2]
+  us0N = (u0, v0)
   odeopcache = allocate_odeopcache(odeop, t0, us0N)
 
-  uα, vα = copy(u0), copy(u0)
+  uα, vα = copy(u0), copy(v0)
 
   sysslvrcache = nothing
   odeslvrcache = (uα, vα, sysslvrcache)
@@ -104,6 +113,24 @@ function ode_start(
   (state0, odecache)
 end
 
+function ode_start(
+  odeslvr::GeneralizedAlpha1, odeop::ODEOperator,
+  t0::Real, us0::NTuple{2,AbstractVector},
+  odecache
+)
+  # Unpack inputs
+  u0, v0 = us0[1], us0[2]
+
+  # Allocate state
+  s0, s1 = copy(u0), copy(v0)
+
+  # Update state
+  state0 = (s0, s1)
+
+  # Pack outputs
+  (state0, odecache)
+end
+
 function ode_march!(
   stateF::NTuple{2,AbstractVector},
   odeslvr::GeneralizedAlpha1, odeop::ODEOperator,
@@ -163,13 +190,13 @@ end
 ###############
 function allocate_odecache(
   odeslvr::GeneralizedAlpha1, odeop::ODEOperator{<:AbstractLinearODE},
-  t0::Real, us0::NTuple{1,AbstractVector}
+  t0::Real, us0::NTuple{2,AbstractVector}
 )
-  u0 = us0[1]
-  us0N = (u0, u0)
+  u0, v0 = us0[1], us0[2]
+  us0N = (u0, v0)
   odeopcache = allocate_odeopcache(odeop, t0, us0N)
 
-  uα, vα = zero(u0), zero(u0)
+  uα, vα = zero(u0), zero(v0)
 
   constant_stiffness = is_form_constant(odeop, 0)
   constant_mass = is_form_constant(odeop, 1)
@@ -229,6 +256,24 @@ function ode_start(
   (state0, odecache)
 end
 
+function ode_start(
+  odeslvr::GeneralizedAlpha1, odeop::ODEOperator{<:AbstractLinearODE},
+  t0::Real, us0::NTuple{2,AbstractVector},
+  odecache
+)
+  # Unpack inputs
+  u0, v0 = us0[1], us0[2]
+
+  # Allocate state
+  s0, s1 = copy(u0), copy(v0)
+
+  # Update state
+  state0 = (s0, s1)
+
+  # Pack outputs
+  (state0, odecache)
+end
+
 function ode_march!(
   stateF::NTuple{2,AbstractVector},
   odeslvr::GeneralizedAlpha1, odeop::ODEOperator{<:AbstractLinearODE},