diff --git a/NEWS.md b/NEWS.md index 234a9a52b..7795d370f 100644 --- a/NEWS.md +++ b/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 diff --git a/Project.toml b/Project.toml index 065fa733c..8073aeab3 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.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" diff --git a/docs/src/ODEs.md b/docs/src/ODEs.md index 501611e39..2e6bbb4bc 100644 --- a/docs/src/ODEs.md +++ b/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*} diff --git a/src/Algebra/AlgebraInterfaces.jl b/src/Algebra/AlgebraInterfaces.jl index 41ed27f4f..7548cd9cb 100644 --- a/src/Algebra/AlgebraInterfaces.jl +++ b/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 # diff --git a/src/Arrays/Arrays.jl b/src/Arrays/Arrays.jl index 821304d96..ce7452453 100644 --- a/src/Arrays/Arrays.jl +++ b/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") diff --git a/src/Arrays/PrintOpTrees.jl b/src/Arrays/PrintOpTrees.jl index 718f38981..b4791f678 100644 --- a/src/Arrays/PrintOpTrees.jl +++ b/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 diff --git a/src/Arrays/SubVectors.jl b/src/Arrays/SubVectors.jl index e8c94f791..c742a888a 100644 --- a/src/Arrays/SubVectors.jl +++ b/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 diff --git a/src/CellData/CellData.jl b/src/CellData/CellData.jl index 9461e34d8..19624ec27 100644 --- a/src/CellData/CellData.jl +++ b/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 ∇∇ diff --git a/src/FESpaces/Assemblers.jl b/src/FESpaces/Assemblers.jl index 6d14c29ab..4aa7836d3 100644 --- a/src/FESpaces/Assemblers.jl +++ b/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 diff --git a/src/MultiField/MultiFieldCellFields.jl b/src/MultiField/MultiFieldCellFields.jl index ff2378fc0..c324dbf0e 100644 --- a/src/MultiField/MultiFieldCellFields.jl +++ b/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 diff --git a/src/MultiField/MultiFieldFEFunctions.jl b/src/MultiField/MultiFieldFEFunctions.jl index 6578e8fdc..12789bc8d 100644 --- a/src/MultiField/MultiFieldFEFunctions.jl +++ b/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 diff --git a/src/MultiField/MultiFieldFESpaces.jl b/src/MultiField/MultiFieldFESpaces.jl index 47cf02210..580c038c7 100644 --- a/src/MultiField/MultiFieldFESpaces.jl +++ b/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) diff --git a/src/ODEs/ODESolvers/GeneralizedAlpha1.jl b/src/ODEs/ODESolvers/GeneralizedAlpha1.jl index 96dd235ef..ec7ac7e0f 100644 --- a/src/ODEs/ODESolvers/GeneralizedAlpha1.jl +++ b/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}, diff --git a/src/ODEs/ODESolvers/GeneralizedAlpha2.jl b/src/ODEs/ODESolvers/GeneralizedAlpha2.jl index 6d50a01a2..07052cea0 100644 --- a/src/ODEs/ODESolvers/GeneralizedAlpha2.jl +++ b/src/ODEs/ODESolvers/GeneralizedAlpha2.jl @@ -68,18 +68,27 @@ function Newmark(sysslvr::NonlinearSolver, dt::Real, γ::Real, β::Real) GeneralizedAlpha2(sysslvr, dt, αf, αm, γ, β) end +# Default allocate_odecache without acceleration +function allocate_odecache( + odeslvr::GeneralizedAlpha2, odeop::ODEOperator, + t0::Real, us0::NTuple{2,AbstractVector} +) + u0, v0 = us0[1], us0[2] + allocate_odecache(odeslvr, odeop, t0, (u0, v0, v0)) +end + ################## # Nonlinear case # ################## function allocate_odecache( odeslvr::GeneralizedAlpha2, odeop::ODEOperator, - t0::Real, us0::NTuple{2,AbstractVector} + t0::Real, us0::NTuple{3,AbstractVector} ) - u0, v0 = us0[1], us0[2] - us0N = (u0, v0, v0) + u0, v0, a0 = us0[1], us0[2], us0[3] + us0N = (u0, v0, a0) odeopcache = allocate_odeopcache(odeop, t0, us0N) - uα, vα, aα = copy(u0), copy(v0), copy(v0) + uα, vα, aα = copy(u0), copy(v0), copy(a0) sysslvrcache = nothing odeslvrcache = (uα, vα, aα, sysslvrcache) @@ -129,6 +138,24 @@ function ode_start( (state0, odecache) end +function ode_start( + odeslvr::GeneralizedAlpha2, odeop::ODEOperator, + t0::Real, us0::NTuple{3,AbstractVector}, + odecache +) + # Unpack inputs + u0, v0, a0 = us0[1], us0[2], us0[3] + + # Allocate state + s0, s1, s2 = copy(u0), copy(v0), copy(a0) + + # Update state + state0 = (s0, s1, s2) + + # Pack outputs + (state0, odecache) +end + function ode_march!( stateF::NTuple{3,AbstractVector}, odeslvr::GeneralizedAlpha2, odeop::ODEOperator, @@ -191,13 +218,13 @@ end ############### function allocate_odecache( odeslvr::GeneralizedAlpha2, odeop::ODEOperator{<:AbstractLinearODE}, - t0::Real, us0::NTuple{2,AbstractVector} + t0::Real, us0::NTuple{3,AbstractVector} ) - u0, v0 = us0[1], us0[2] - us0N = (u0, v0, v0) + u0, v0, a0 = us0[1], us0[2], us0[3] + us0N = (u0, v0, a0) odeopcache = allocate_odeopcache(odeop, t0, us0N) - uα, vα, aα = zero(u0), zero(v0), zero(v0) + uα, vα, aα = zero(u0), zero(v0), zero(a0) constant_stiffness = is_form_constant(odeop, 0) constant_damping = is_form_constant(odeop, 1) @@ -258,6 +285,24 @@ function ode_start( (state0, odecache) end +function ode_start( + odeslvr::GeneralizedAlpha2, odeop::ODEOperator{<:AbstractLinearODE}, + t0::Real, us0::NTuple{3,AbstractVector}, + odecache +) + # Unpack inputs + u0, v0, a0 = us0[1], us0[2], us0[3] + + # Allocate state + s0, s1, s2 = copy(u0), copy(v0), copy(a0) + + # Update state + state0 = (s0, s1, s2) + + # Pack outputs + (state0, odecache) +end + function ode_march!( stateF::NTuple{3,AbstractVector}, odeslvr::GeneralizedAlpha2, odeop::ODEOperator{<:AbstractLinearODE}, diff --git a/src/ODEs/ODEs.jl b/src/ODEs/ODEs.jl index c172124d8..daec78a48 100644 --- a/src/ODEs/ODEs.jl +++ b/src/ODEs/ODEs.jl @@ -35,6 +35,7 @@ const ε = 100 * eps() include("TimeDerivatives.jl") +export TimeSpaceFunction export time_derivative export ∂t export ∂tt diff --git a/src/ODEs/TimeDerivatives.jl b/src/ODEs/TimeDerivatives.jl index df66902ac..d14885900 100644 --- a/src/ODEs/TimeDerivatives.jl +++ b/src/ODEs/TimeDerivatives.jl @@ -1,3 +1,44 @@ +##################### +# TimeSpaceFunction # +##################### +""" + struct TimeSpaceFunction{F} <: Function end + +`TimeSpaceFunction` allows for convenient ways to apply differential operators to +functions that depend on time and space. More precisely, if `f` is a function that, to +a given time, returns a function of space (i.e. `f` is evaluated at time `t` and position +`x` via `f(t)(x)`), then `F = TimeSpaceFunction(f)` supports the following syntax: +* `op(F)`: a `TimeSpaceFunction` representing both `t -> x -> op(f)(t)(x)` and `(t, x) -> op(f)(t)(x)`, +* `op(F)(t)`: a function of space representing `x -> op(f)(t)(x)` +* `op(F)(t, x)`: the quantity `op(f)(t)(x)` (this notation is equivalent to `op(F)(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`, `∇`, ...). +""" +struct TimeSpaceFunction{F} <: Function + f::F +end + +(ts::TimeSpaceFunction)(t) = ts.f(t) +(ts::TimeSpaceFunction)(t, x) = ts.f(t)(x) + +################################## +# Spatial differential operators # +################################## +# Using the rule spatial_op(ts)(t)(x) = spatial_op(ts(t))(x) +for op in (:(Fields.gradient), :(Fields.symmetric_gradient), :(Fields.divergence), + :(Fields.curl), :(Fields.laplacian)) + @eval begin + function ($op)(ts::TimeSpaceFunction) + function _op(t) + _op_t(x) = $op(ts(t))(x) + _op_t + end + TimeSpaceFunction(_op) + end + end +end + ############################# # time_derivative interface # ############################# @@ -64,28 +105,37 @@ end # Specialisation for `Function` # ################################# function time_derivative(f::Function) - function dfdt(x, t) - z = zero(return_type(f, x, t)) - _time_derivative(f, x, t, z) + function dfdt(t) + ft = f(t) + function dfdt_t(x) + T = return_type(ft, x) + _time_derivative(T, f, t, x) + end + dfdt_t end - # Extend definition to include restrictions - _dfdt(x, t) = dfdt(x, t) - _dfdt(x::VectorValue) = t -> dfdt(x, t) - _dfdt(t::Real) = x -> dfdt(x, t) - return _dfdt + dfdt +end + +function _time_derivative(T::Type{<:Real}, f, t, x) + partial(t) = f(t)(x) + ForwardDiff.derivative(partial, t) end -function _time_derivative(f, x, t, z) - ForwardDiff.derivative(t -> f(x, t), t) +function _time_derivative(T::Type{<:VectorValue}, f, t, x) + partial(t) = get_array(f(t)(x)) + VectorValue(ForwardDiff.derivative(partial, t)) end -function _time_derivative(f, x, t, z::VectorValue) - VectorValue(ForwardDiff.derivative(t -> get_array(f(x, t)), t)) - # VectorValue(ForwardDiff.derivative(t -> f(x, t), t)) +function _time_derivative(T::Type{<:TensorValue}, f, t, x) + partial(t) = get_array(f(t)(x)) + TensorValue(ForwardDiff.derivative(partial, t)) end -function _time_derivative(f, x, t, z::TensorValue) - TensorValue(ForwardDiff.derivative(t -> get_array(f(x, t)), t)) +########################################## +# Specialisation for `TimeSpaceFunction` # +########################################## +function time_derivative(ts::TimeSpaceFunction) + TimeSpaceFunction(time_derivative(ts.f)) end ############################### diff --git a/src/ODEs/TransientFEOperators.jl b/src/ODEs/TransientFEOperators.jl index 3244a42d1..cd6284aca 100644 --- a/src/ODEs/TransientFEOperators.jl +++ b/src/ODEs/TransientFEOperators.jl @@ -220,10 +220,10 @@ end # Constructor with manual jacobians function TransientFEOperator( res::Function, jacs::Tuple{Vararg{Function}}, - trial, test + trial, test; + assembler=SparseMatrixAssembler(trial, test) ) order = length(jacs) - 1 - assembler = SparseMatrixAssembler(trial, test) TransientFEOpFromWeakForm( res, jacs, assembler, trial, test, order @@ -234,33 +234,36 @@ end function TransientFEOperator( res::Function, jac::Function, - trial, test + trial, test; + assembler=SparseMatrixAssembler(trial, test) ) TransientFEOperator( - res, (jac,), - trial, test + res, (jac,), trial, test; + assembler ) end function TransientFEOperator( res::Function, jac::Function, jac_t::Function, - trial, test + trial, test; + assembler=SparseMatrixAssembler(trial, test) ) TransientFEOperator( - res, (jac, jac_t), - trial, test + res, (jac, jac_t), trial, test; + assembler ) end function TransientFEOperator( res::Function, jac::Function, jac_t::Function, jac_tt::Function, - trial, test + trial, test; + assembler=SparseMatrixAssembler(trial, test) ) TransientFEOperator( - res, (jac, jac_t, jac_tt), - trial, test + res, (jac, jac_t, jac_tt), trial, test; + assembler ) end @@ -268,7 +271,8 @@ end function TransientFEOperator( res::Function, trial, test; - order::Integer=1 + order::Integer=1, + assembler=SparseMatrixAssembler(trial, test) ) function jac_0(t, u, du, v) function res_0(y) @@ -291,7 +295,10 @@ function TransientFEOperator( jacs = (jacs..., jac_k) end - TransientFEOperator(res, jacs, trial, test) + TransientFEOperator( + res, jacs, trial, test; + assembler + ) end # TransientFEOperator interface @@ -338,15 +345,18 @@ end # Constructor with manual jacobians function TransientQuasilinearFEOperator( mass::Function, res::Function, jacs::Tuple{Vararg{Function}}, - trial, test + trial, test; + assembler=SparseMatrixAssembler(trial, test) ) order = length(jacs) - 1 if order == 0 @warn default_linear_msg - return TransientLinearFEOperator((mass,), res, jacs, trial, test) + return TransientLinearFEOperator( + (mass,), res, jacs, trial, test; + assembler + ) end - assembler = SparseMatrixAssembler(trial, test) TransientQuasilinearFEOpFromWeakForm( mass, res, jacs, assembler, trial, test, order @@ -357,31 +367,37 @@ end function TransientQuasilinearFEOperator( mass::Function, res::Function, jac::Function, - trial, test + trial, test; + assembler=SparseMatrixAssembler(trial, test) ) @warn default_linear_msg - TransientLinearFEOperator(mass, res, jac, trial, test) + TransientLinearFEOperator( + mass, res, jac, trial, test; + assembler + ) end function TransientQuasilinearFEOperator( mass::Function, res::Function, jac::Function, jac_t::Function, - trial, test + trial, test; + assembler=SparseMatrixAssembler(trial, test) ) TransientQuasilinearFEOperator( - mass, res, (jac, jac_t), - trial, test + mass, res, (jac, jac_t), trial, test; + assembler ) end function TransientQuasilinearFEOperator( mass::Function, res::Function, jac::Function, jac_t::Function, jac_tt::Function, - trial, test + trial, test; + assembler=SparseMatrixAssembler(trial, test) ) TransientQuasilinearFEOperator( - mass, res, (jac, jac_t, jac_tt), - trial, test + mass, res, (jac, jac_t, jac_tt), trial, test; + assembler ) end @@ -389,11 +405,15 @@ end function TransientQuasilinearFEOperator( mass::Function, res::Function, trial, test; - order::Integer=1 + order::Integer=1, + assembler=SparseMatrixAssembler(trial, test) ) if order == 0 @warn default_linear_msg - return TransientLinearFEOperator(mass, res, trial, test) + return TransientLinearFEOperator( + mass, res, trial, test; + assembler + ) end jacs = () @@ -428,7 +448,8 @@ function TransientQuasilinearFEOperator( jacs = (jacs..., jac_N) TransientQuasilinearFEOperator( - mass, res, jacs, trial, test + mass, res, jacs, trial, test; + assembler ) end @@ -481,16 +502,18 @@ function TransientSemilinearFEOperator( mass::Function, res::Function, jacs::Tuple{Vararg{Function}}, trial, test; constant_mass::Bool=false, + assembler=SparseMatrixAssembler(trial, test) ) order = length(jacs) - 1 if order == 0 @warn default_linear_msg + constant_forms = (constant_mass,) return TransientLinearFEOperator( - (mass,), res, jacs, trial, test; constant_mass + (mass,), res, jacs, trial, test; + constant_forms, assembler ) end - assembler = SparseMatrixAssembler(trial, test) TransientSemilinearFEOpFromWeakForm( mass, res, jacs, constant_mass, assembler, trial, test, order @@ -503,9 +526,14 @@ function TransientSemilinearFEOperator( jac::Function, trial, test; constant_mass::Bool=false, + assembler=SparseMatrixAssembler(trial, test) ) @warn default_linear_msg - TransientLinearFEOperator(mass, res, jac, trial, test; constant_mass) + constant_forms = (constant_mass,) + TransientLinearFEOperator( + mass, res, jac, trial, test; + constant_forms, assembler + ) end function TransientSemilinearFEOperator( @@ -513,11 +541,11 @@ function TransientSemilinearFEOperator( jac::Function, jac_t::Function, trial, test; constant_mass::Bool=false, + assembler=SparseMatrixAssembler(trial, test) ) TransientSemilinearFEOperator( - mass, res, (jac, jac_t), - trial, test; - constant_mass + mass, res, (jac, jac_t), trial, test; + constant_mass, assembler ) end @@ -526,11 +554,11 @@ function TransientSemilinearFEOperator( jac::Function, jac_t::Function, jac_tt::Function, trial, test; constant_mass::Bool=false, + assembler=SparseMatrixAssembler(trial, test) ) TransientSemilinearFEOperator( - mass, res, (jac, jac_t, jac_tt), - trial, test; - constant_mass + mass, res, (jac, jac_t, jac_tt), trial, test; + constant_mass, assembler ) end @@ -539,11 +567,16 @@ function TransientSemilinearFEOperator( mass::Function, res::Function, trial, test; order::Integer=1, - constant_mass::Bool=false + constant_mass::Bool=false, + assembler=SparseMatrixAssembler(trial, test) ) if order == 0 @warn default_linear_msg - return TransientLinearFEOperator(mass, res, trial, test; constant_mass) + constant_forms = (constant_mass,) + return TransientLinearFEOperator( + mass, res, trial, test; + constant_forms, assembler + ) end # When the operator is semilinear, the mass term can be omitted in the @@ -579,7 +612,7 @@ function TransientSemilinearFEOperator( TransientSemilinearFEOperator( mass, res, jacs, trial, test; - constant_mass + constant_mass, assembler ) end @@ -635,10 +668,10 @@ end function TransientLinearFEOperator( forms::Tuple{Vararg{Function}}, res::Function, jacs::Tuple{Vararg{Function}}, trial, test; - constant_forms::Tuple{Vararg{Bool}}=ntuple(_ -> false, length(forms)) + constant_forms::Tuple{Vararg{Bool}}=ntuple(_ -> false, length(forms)), + assembler=SparseMatrixAssembler(trial, test) ) order = length(jacs) - 1 - assembler = SparseMatrixAssembler(trial, test) TransientLinearFEOpFromWeakForm( forms, res, jacs, constant_forms, assembler, trial, test, order @@ -652,7 +685,8 @@ end function TransientLinearFEOperator( forms::Tuple{Vararg{Function}}, res::Function, trial, test; - constant_forms::Tuple{Vararg{Bool}}=ntuple(_ -> false, length(forms)) + constant_forms::Tuple{Vararg{Bool}}=ntuple(_ -> false, length(forms)), + assembler=SparseMatrixAssembler(trial, test) ) # When the operator is linear, the jacobians are the forms themselves order = length(forms) - 1 @@ -660,7 +694,7 @@ function TransientLinearFEOperator( TransientLinearFEOperator( forms, res, jacs, trial, test; - constant_forms + constant_forms, assembler ) end @@ -668,33 +702,36 @@ end function TransientLinearFEOperator( mass::Function, res::Function, trial, test; - constant_forms::NTuple{1,Bool}=(false,) + constant_forms::NTuple{1,Bool}=(false,), + assembler=SparseMatrixAssembler(trial, test) ) TransientLinearFEOperator( - (mass,), res, - trial, test; constant_forms + (mass,), res, trial, test; + constant_forms, assembler ) end function TransientLinearFEOperator( stiffness::Function, mass::Function, res::Function, trial, test; - constant_forms::NTuple{2,Bool}=(false, false) + constant_forms::NTuple{2,Bool}=(false, false), + assembler=SparseMatrixAssembler(trial, test) ) TransientLinearFEOperator( - (stiffness, mass), res, - trial, test; constant_forms + (stiffness, mass), res, trial, test; + constant_forms, assembler ) end function TransientLinearFEOperator( stiffness::Function, damping::Function, mass::Function, res::Function, trial, test; - constant_forms::NTuple{3,Bool}=(false, false, false) + constant_forms::NTuple{3,Bool}=(false, false, false), + assembler=SparseMatrixAssembler(trial, test) ) TransientLinearFEOpFromWeakForm( - (stiffness, damping, mass), res, - trial, test; constant_forms + (stiffness, damping, mass), res, trial, test; + constant_forms, assembler ) end diff --git a/test/ArraysTests/SubVectorsTests.jl b/test/ArraysTests/SubVectorsTests.jl deleted file mode 100644 index 3a9e34726..000000000 --- a/test/ArraysTests/SubVectorsTests.jl +++ /dev/null @@ -1,27 +0,0 @@ -module SubVectorsTests -using Gridap -using Gridap.Arrays - -a = collect(1:10) -pini = 1 -pend = 5 -b = SubVector(a,pini,pend) -test_array(b,collect(pini:pend)) - -pini = 7 -pend = 9 -b = SubVector(a,pini,pend) -test_array(b,collect(pini:pend)) - -inds = Int32[1,4,3] -a = rand(10) -sa = Gridap.Arrays.SubVector(a,1,5) -vsa = view(sa,inds) - -b = rand(10) -vb = view(b,inds) - -x = Base.broadcasted(+,vb,vsa) -Base.materialize!(vb,x) - -end # module diff --git a/test/ArraysTests/runtests.jl b/test/ArraysTests/runtests.jl index 1b4ba1a4c..076a3f863 100644 --- a/test/ArraysTests/runtests.jl +++ b/test/ArraysTests/runtests.jl @@ -26,8 +26,6 @@ using Test @testset "IdentityVectors" begin include("IdentityVectorsTests.jl") end -@testset "SubVectors" begin include("SubVectorsTests.jl") end - @testset "ArrayPairs" begin include("ArrayPairsTests.jl") end @testset "AppendedArrays" begin include("AppendedArraysTests.jl") end diff --git a/test/MultiFieldTests/MultiFieldCellFieldsTests.jl b/test/MultiFieldTests/MultiFieldCellFieldsTests.jl index 1d3292129..efbea1f0e 100644 --- a/test/MultiFieldTests/MultiFieldCellFieldsTests.jl +++ b/test/MultiFieldTests/MultiFieldCellFieldsTests.jl @@ -35,6 +35,8 @@ _cf1, _cf2 = cf @test cf1 === _cf1 @test cf2 === _cf2 +_cfdot = cf⋅cf + order = 2 domain = (0,1,0,1) diff --git a/test/MultiFieldTests/MultiFieldFEFunctionsTests.jl b/test/MultiFieldTests/MultiFieldFEFunctionsTests.jl index 8f3245efe..8b87ab611 100644 --- a/test/MultiFieldTests/MultiFieldFEFunctionsTests.jl +++ b/test/MultiFieldTests/MultiFieldFEFunctionsTests.jl @@ -36,4 +36,6 @@ uh, ph = xh cell_values = get_cell_dof_values(xh,trian) @test isa(cell_values[1],ArrayBlock) +xh2 = xh⋅xh + end # module diff --git a/test/MultiFieldTests/MultiFieldFESpacesTests.jl b/test/MultiFieldTests/MultiFieldFESpacesTests.jl index 786ed6785..21ee00846 100644 --- a/test/MultiFieldTests/MultiFieldFESpacesTests.jl +++ b/test/MultiFieldTests/MultiFieldFESpacesTests.jl @@ -60,6 +60,9 @@ uh, ph = xh @test isa(uh,FEFunction) @test isa(ph,FEFunction) +dir_values = zero_dirichlet_values(Y) +@test all(map((dv,Yi) -> dv == zero_dirichlet_values(Yi),dir_values,Y)) + cell_isconstr = get_cell_isconstrained(X,trian) @test cell_isconstr == Fill(false,num_cells(model)) diff --git a/test/ODEsTests/ODEProblemsTests/Order1ODETests.jl b/test/ODEsTests/ODEProblemsTests/Order1ODETests.jl index 544e1f33d..645116e9a 100644 --- a/test/ODEsTests/ODEProblemsTests/Order1ODETests.jl +++ b/test/ODEsTests/ODEProblemsTests/Order1ODETests.jl @@ -9,7 +9,7 @@ using Gridap.ODEs include("../ODEOperatorsMocks.jl") include("../ODESolversMocks.jl") -# M u̇ + M * Diag(λ) u = M * f(t), +# M u̇ + M * Diag(-λ) u = M * f(t), # where f(t) = exp(Diag(α) * t) t0 = 0.0 dt = 1.0e-3 @@ -87,28 +87,43 @@ maxiter = 100 sysslvr_l = LUSolver() sysslvr_nl = NonlinearSolverMock(rtol, atol, maxiter) -# odeslvrs = ( -# ForwardEuler(sysslvr_nl, dt), -# ThetaMethod(sysslvr_nl, dt, 0.2), -# MidPoint(sysslvr_nl, dt), -# ThetaMethod(sysslvr_nl, dt, 0.8), -# BackwardEuler(sysslvr_nl, dt), -# GeneralizedAlpha1(sysslvr_nl, dt, 0.0), -# GeneralizedAlpha1(sysslvr_nl, dt, 0.5), -# GeneralizedAlpha1(sysslvr_nl, dt, 1.0), -# ) -# for tableau in available_tableaus -# global odeslvrs -# odeslvr = RungeKutta(sysslvr_nl, sysslvr_l, dt, tableau) -# odeslvrs = (odeslvrs..., odeslvr) -# end - -# us0 = (u0,) -# for odeslvr in odeslvrs -# for odeop in odeops -# test_solver(odeslvr, odeop, us0, tol) -# end -# end +odeslvrs = ( + ForwardEuler(sysslvr_nl, dt), + ThetaMethod(sysslvr_nl, dt, 0.2), + MidPoint(sysslvr_nl, dt), + ThetaMethod(sysslvr_nl, dt, 0.8), + BackwardEuler(sysslvr_nl, dt), + GeneralizedAlpha1(sysslvr_nl, dt, 0.0), + GeneralizedAlpha1(sysslvr_nl, dt, 0.5), + GeneralizedAlpha1(sysslvr_nl, dt, 1.0), +) +for tableau in available_tableaus + global odeslvrs + odeslvr = RungeKutta(sysslvr_nl, sysslvr_l, dt, tableau) + odeslvrs = (odeslvrs..., odeslvr) +end + +us0 = (u0,) +for odeslvr in odeslvrs + for odeop in odeops + test_solver(odeslvr, odeop, us0, tol) + end +end + +# Tests with initial velocity +odeslvrs = ( + GeneralizedAlpha1(sysslvr_nl, dt, 0.0), + GeneralizedAlpha1(sysslvr_nl, dt, 0.5), + GeneralizedAlpha1(sysslvr_nl, dt, 1.0), +) + +v0 = exp.(α .* t0) - spdiagm(-λ) * u0 +us0 = (u0, v0) +for odeslvr in odeslvrs + for odeop in odeops + test_solver(odeslvr, odeop, us0, tol) + end +end # Solvers for `IMEXODEOperator`s odeops = ( diff --git a/test/ODEsTests/ODEProblemsTests/Order2ODETests.jl b/test/ODEsTests/ODEProblemsTests/Order2ODETests.jl index b663d11a8..239965a36 100644 --- a/test/ODEsTests/ODEProblemsTests/Order2ODETests.jl +++ b/test/ODEsTests/ODEProblemsTests/Order2ODETests.jl @@ -9,6 +9,8 @@ using Gridap.ODEs include("../ODEOperatorsMocks.jl") include("../ODESolversMocks.jl") +# M ü + M * Diag(-(λ .+ μ)) u̇ + M * Diag(λ .* μ) u = M * f(t), +# where f(t) = exp(Diag(α) * t) t0 = 0.0 dt = 1.0e-3 tF = t0 + 10 * dt @@ -106,4 +108,13 @@ for odeslvr in odeslvrs end end +# Tests with initial acceleration +a0 = exp.(α .* t0) - spdiagm(-(λ .+ μ)) * v0 - spdiagm(λ .* μ) * u0 +us0 = (u0, v0, a0) +for odeslvr in odeslvrs + for odeop in odeops + test_solver(odeslvr, odeop, us0, tol) + end +end + end # module Order2ODETests diff --git a/test/ODEsTests/TimeDerivativesTests.jl b/test/ODEsTests/TimeDerivativesTests.jl index 10a74b0e6..612970dbd 100644 --- a/test/ODEsTests/TimeDerivativesTests.jl +++ b/test/ODEsTests/TimeDerivativesTests.jl @@ -8,94 +8,100 @@ using Gridap using Gridap.ODEs # First time derivative, scalar-valued -f1(x, t) = 5 * x[1] * x[2] + x[2]^2 * t^3 -∂tf1(x, t) = 3 * x[2]^2 * t^2 +f1(t) = x -> 5 * x[1] * x[2] + x[2]^2 * t^3 +∂tf1(t) = x -> 3 * x[2]^2 * t^2 -f2(x, t) = t^2 -∂tf2(x, t) = 2 * t +f2(t) = x -> t^2 +∂tf2(t) = x -> 2 * t -f3(x, t) = x[1]^2 -∂tf3(x, t) = zero(x[1]) +f3(t) = x -> x[1]^2 +∂tf3(t) = x -> zero(x[1]) -f4(x, t) = x[1]^t^2 -∂tf4(x, t) = 2 * t * log(x[1]) * f4(x, t) +f4(t) = x -> x[1]^t^2 +∂tf4(t) = x -> 2 * t * log(x[1]) * f4(t)(x) for (f, ∂tf) in ((f1, ∂tf1), (f2, ∂tf2), (f3, ∂tf3), (f4, ∂tf4),) - dtf = (x, t) -> ForwardDiff.derivative(t -> f(x, t), t) + dtf(t) = x -> ForwardDiff.derivative(t -> f(t)(x), t) tv = rand(Float64) xv = Point(rand(Float64, 2)...) - @test ∂t(f)(xv, tv) ≈ ∂tf(xv, tv) - @test ∂t(f)(xv, tv) ≈ dtf(xv, tv) - @test ∂t(f)(xv, tv) ≈ ∂t(f)(xv)(tv) ≈ ∂t(f)(tv)(xv) + @test ∂t(f)(tv)(xv) ≈ ∂tf(tv)(xv) + @test ∂t(f)(tv)(xv) ≈ dtf(tv)(xv) + + F = TimeSpaceFunction(f) + @test F(tv)(xv) ≈ f(tv)(xv) + @test ∂t(F)(tv)(xv) ≈ ∂tf(tv)(xv) end # First time derivative, vector-valued -f1(x, t) = VectorValue(5 * x[1] * x[2], x[2]^2 * t^3) -∂tf1(x, t) = VectorValue(zero(x[1]), x[2]^2 * 3 * t^2) +f1(t) = x -> VectorValue(5 * x[1] * x[2], x[2]^2 * t^3) +∂tf1(t) = x -> VectorValue(zero(x[1]), x[2]^2 * 3 * t^2) -f2(x, t) = VectorValue(x[1]^2, zero(x[2])) -∂tf2(x, t) = VectorValue(zero(x[1]), zero(x[2])) +f2(t) = x -> VectorValue(x[1]^2, zero(x[2])) +∂tf2(t) = x -> VectorValue(zero(x[1]), zero(x[2])) -f3(x, t) = VectorValue(x[1]^2, t) -∂tf3(x, t) = VectorValue(zero(x[1]), one(t)) +f3(t) = x -> VectorValue(x[1]^2, t) +∂tf3(t) = x -> VectorValue(zero(x[1]), one(t)) for (f, ∂tf) in ((f1, ∂tf1), (f2, ∂tf2), (f3, ∂tf3),) - dtf = (x, t) -> VectorValue(ForwardDiff.derivative(t -> get_array(f(x, t)), t)) + dtf(t) = x -> VectorValue(ForwardDiff.derivative(t -> get_array(f(t)(x)), t)) tv = rand(Float64) xv = Point(rand(Float64, 2)...) - @test ∂t(f)(xv, tv) ≈ ∂tf(xv, tv) - @test ∂t(f)(xv, tv) ≈ dtf(xv, tv) - @test ∂t(f)(xv, tv) ≈ ∂t(f)(xv)(tv) ≈ ∂t(f)(tv)(xv) + @test ∂t(f)(tv)(xv) ≈ ∂tf(tv)(xv) + @test ∂t(f)(tv)(xv) ≈ dtf(tv)(xv) + + F = TimeSpaceFunction(f) + @test F(tv)(xv) ≈ f(tv)(xv) + @test ∂t(F)(tv)(xv) ≈ ∂tf(tv)(xv) end # First time derivative, tensor-valued -f1(x, t) = TensorValue(x[1] * t, x[2] * t, x[1] * x[2], x[1] * t^2) -∂tf1(x, t) = TensorValue(x[1], x[2], zero(x[1]), 2 * x[1] * t) +f1(t) = x -> TensorValue(x[1] * t, x[2] * t, x[1] * x[2], x[1] * t^2) +∂tf1(t) = x -> TensorValue(x[1], x[2], zero(x[1]), 2 * x[1] * t) for (f, ∂tf) in ((f1, ∂tf1),) - dtf = (x, t) -> TensorValue(ForwardDiff.derivative(t -> get_array(f(x, t)), t)) + dtf(t) = x -> TensorValue(ForwardDiff.derivative(t -> get_array(f(t)(x)), t)) tv = rand(Float64) xv = Point(rand(Float64, 2)...) - @test ∂t(f)(xv, tv) ≈ ∂tf(xv, tv) - @test ∂t(f)(xv, tv) ≈ dtf(xv, tv) - @test ∂t(f)(xv, tv) ≈ ∂t(f)(xv)(tv) ≈ ∂t(f)(tv)(xv) + @test ∂t(f)(tv)(xv) ≈ ∂tf(tv)(xv) + @test ∂t(f)(tv)(xv) ≈ dtf(tv)(xv) + + F = TimeSpaceFunction(f) + @test F(tv)(xv) ≈ f(tv)(xv) + @test ∂t(F)(tv)(xv) ≈ ∂tf(tv)(xv) end # Spatial derivatives -# f(x, t) = VectorValue(x[1]^2, t) -# ∇f(x, t) = ∇(y -> f(y, t))(x) +ft(t) = x -> x[1]^2 * t + x[2] +f = TimeSpaceFunction(ft) -# tv = rand(Float64) -# xv = Point(rand(Float64, 2)...) -# ∇f(xv, tv) - -# TODO -# @santiagobadia : Is there any way to make this transparent to the user -# I guess not unless we create a type for these analytical (space-only or -# space-time via a trait) functions -# Probably a try-catch? +tv = rand(Float64) +xv = Point(rand(Float64, 2)...) +@test ∇(f)(tv)(xv) ≈ Point(2 * xv[1] * tv, 1.0) # Second time derivative, scalar-valued -f1(x, t) = t^2 -∂tf1(x, t) = 2 * t -∂ttf1(x, t) = 2 * one(t) +f1(t) = x -> t^2 +∂tf1(t) = x -> 2 * t +∂ttf1(t) = x -> 2 * one(t) -f2(x, t) = x[1] * t^2 -∂tf2(x, t) = 2 * x[1] * t -∂ttf2(x, t) = 2 * x[1] +f2(t) = x -> x[1] * t^2 +∂tf2(t) = x -> 2 * x[1] * t +∂ttf2(t) = x -> 2 * x[1] for (f, ∂tf, ∂ttf) in ((f1, ∂tf1, ∂ttf1), (f2, ∂tf2, ∂ttf2),) - dtf = (x, t) -> ForwardDiff.derivative(t -> f(x, t), t) - dttf = (x, t) -> ForwardDiff.derivative(t -> dtf(x, t), t) + dtf(t) = x -> ForwardDiff.derivative(t -> f(t)(x), t) + dttf(t) = x -> ForwardDiff.derivative(t -> dtf(t)(x), t) tv = rand(Float64) xv = Point(rand(Float64, 2)...) - @test ∂tt(f)(xv, tv) ≈ ∂ttf(xv, tv) - @test ∂tt(f)(xv, tv) ≈ dttf(xv, tv) - @test ∂tt(f)(xv, tv) ≈ ∂tt(f)(xv)(tv) ≈ ∂tt(f)(tv)(xv) + @test ∂tt(f)(tv)(xv) ≈ ∂ttf(tv)(xv) + @test ∂tt(f)(tv)(xv) ≈ dttf(tv)(xv) + + F = TimeSpaceFunction(f) + @test F(tv)(xv) ≈ f(tv)(xv) + @test ∂tt(F)(tv)(xv) ≈ ∂ttf(tv)(xv) end end # module TimeDerivativesTests diff --git a/test/ODEsTests/TransientCellFieldsTests.jl b/test/ODEsTests/TransientCellFieldsTests.jl index bbe8a6a6b..e47be226f 100644 --- a/test/ODEsTests/TransientCellFieldsTests.jl +++ b/test/ODEsTests/TransientCellFieldsTests.jl @@ -12,8 +12,7 @@ using Gridap.MultiField using Gridap.ODEs using Gridap.ODEs: TransientMultiFieldCellField -f(x, t) = sum(x) -f(t::Real) = x -> f(x, t) +f(t) = x -> sum(x) domain = (0, 1, 0, 1) partition = (5, 5) diff --git a/test/ODEsTests/TransientFEOperatorsSolutionsTests.jl b/test/ODEsTests/TransientFEOperatorsSolutionsTests.jl index f6bc81ff2..3bb7024aa 100644 --- a/test/ODEsTests/TransientFEOperatorsSolutionsTests.jl +++ b/test/ODEsTests/TransientFEOperatorsSolutionsTests.jl @@ -15,12 +15,8 @@ using Gridap.ODEs include("ODESolversMocks.jl") # Analytical functions -u(x, t) = x[1] * (1 - x[2]) * (1 + t) -u(t::Real) = x -> u(x, t) -u(x) = t -> u(x, t) - -∂tu(x, t) = ∂t(u)(x, t) -∂tu(t::Real) = x -> ∂tu(x, t) +ut(t) = x -> x[1] * (1 - x[2]) * (1 + t) +u = TimeSpaceFunction(ut) # Geometry domain = (0, 1, 0, 1) @@ -108,7 +104,8 @@ end # First-order # ############### order = 1 -f(t) = x -> ∂t(u)(x, t) - Δ(u(t))(x) +ft(t) = x -> ∂t(u)(t, x) - Δ(u)(t, x) +f = TimeSpaceFunction(ft) mass(t, ∂ₜu, v) = ∫(∂ₜu ⋅ v) * dΩ mass(t, u, ∂ₜu, v) = mass(t, ∂ₜu, v) @@ -203,7 +200,8 @@ end # Second-order # ################ order = 2 -f(t) = x -> ∂tt(u)(x, t) + ∂t(u)(x, t) - Δ(u(t))(x) +ft(t) = x -> ∂tt(u)(t, x) + ∂t(u)(t, x) - Δ(u)(t, x) +f = TimeSpaceFunction(ft) mass(t, ∂ₜₜu, v) = ∫(∂ₜₜu ⋅ v) * dΩ mass(t, u, ∂ₜₜu, v) = mass(t, ∂ₜₜu, v) diff --git a/test/ODEsTests/TransientFEProblemsTests/FreeSurfacePotentialFlowTests.jl b/test/ODEsTests/TransientFEProblemsTests/FreeSurfacePotentialFlowTests.jl index 675d1be19..9da67f79d 100644 --- a/test/ODEsTests/TransientFEProblemsTests/FreeSurfacePotentialFlowTests.jl +++ b/test/ODEsTests/TransientFEProblemsTests/FreeSurfacePotentialFlowTests.jl @@ -23,10 +23,8 @@ tF = 10 * dt # 2 * π tol = 1.0e-2 # Exact solution -ϕₑ(x, t) = ω / k * ξ * (cosh(k * (x[2]))) / sinh(k * H) * sin(k * x[1] - ω * t) -ηₑ(x, t) = ξ * cos(k * x[1] - ω * t) -ϕₑ(t::Real) = x -> ϕₑ(x, t) -ηₑ(t::Real) = x -> ηₑ(x, t) +ϕₑ(t) = x -> ω / k * ξ * (cosh(k * (x[2]))) / sinh(k * H) * sin(k * x[1] - ω * t) +ηₑ(t) = x -> ξ * cos(k * x[1] - ω * t) # Domain domain = (0, L, 0, H) @@ -84,7 +82,7 @@ UΓ0 = U_Γ(t0) X0 = X(t0) uh0 = interpolate_everywhere(ϕₑ(t0), U0) uhΓ0 = interpolate_everywhere(ηₑ(t0), UΓ0) -xh0 = interpolate_everywhere([uh0, uhΓ0], X0); +xh0 = interpolate_everywhere([uh0, uhΓ0], X0) xhs0 = (xh0,) function test_flow_operator(op) diff --git a/test/ODEsTests/TransientFEProblemsTests/HeatEquationDGTests.jl b/test/ODEsTests/TransientFEProblemsTests/HeatEquationDGTests.jl index 9a650d8b4..1c0aa3167 100644 --- a/test/ODEsTests/TransientFEProblemsTests/HeatEquationDGTests.jl +++ b/test/ODEsTests/TransientFEProblemsTests/HeatEquationDGTests.jl @@ -10,12 +10,8 @@ using Gridap.FESpaces using Gridap.ODEs # Analytical functions -u(x, t) = x[1] * (1 - x[2]) * (1 + t) -u(t::Real) = x -> u(x, t) -u(x) = t -> u(x, t) - -∂tu(x, t) = ∂t(u)(x, t) -∂tu(t::Real) = x -> ∂tu(x, t) +ut(t) = x -> x[1] * (1 - x[2]) * (1 + t) +u = TimeSpaceFunction(ut) # Geometry domain = (0, 1, 0, 1) @@ -42,7 +38,8 @@ dΛ = Measure(Λ, degree) nΛ = get_normal_vector(Λ) # FE operator -f(t) = x -> ∂t(u)(x, t) - Δ(u(t))(x) +ft(t) = x -> ∂t(u)(t, x) - Δ(u)(t, x) +f = TimeSpaceFunction(ft) h = 1 / 5 γ = order * (order + 1) diff --git a/test/ODEsTests/TransientFEProblemsTests/HeatEquationMultiFieldTests.jl b/test/ODEsTests/TransientFEProblemsTests/HeatEquationMultiFieldTests.jl index 61bd123c6..234ceb412 100644 --- a/test/ODEsTests/TransientFEProblemsTests/HeatEquationMultiFieldTests.jl +++ b/test/ODEsTests/TransientFEProblemsTests/HeatEquationMultiFieldTests.jl @@ -10,12 +10,8 @@ using Gridap.FESpaces using Gridap.ODEs # Analytical functions -u(x, t) = x[1] * (1 - x[2]) * (1 + t) -u(t::Real) = x -> u(x, t) -u(x) = t -> u(x, t) - -∂tu(x, t) = ∂t(u)(x, t) -∂tu(t::Real) = x -> ∂tu(x, t) +ut(t) = x -> x[1] * (1 - x[2]) * (1 + t) +u = TimeSpaceFunction(ut) # Geometry domain = (0, 1, 0, 1) @@ -37,7 +33,8 @@ degree = 2 * order dΩ = Measure(Ω, degree) # FE operator -f(t) = x -> ∂t(u)(x, t) - Δ(u(t))(x) +ft(t) = x -> ∂t(u)(t, x) - Δ(u)(t, x) +f = TimeSpaceFunction(ft) _mass(t, ∂ₜu, v) = ∫(∂ₜu ⋅ v) * dΩ _mass(t, u, ∂ₜu, v) = _mass(t, ∂ₜu, v) _stiffness(t, u, v) = ∫(∇(u) ⊙ ∇(v)) * dΩ diff --git a/test/ODEsTests/TransientFEProblemsTests/HeatEquationNeumannTests.jl b/test/ODEsTests/TransientFEProblemsTests/HeatEquationNeumannTests.jl index f9a75c90f..1b6b68131 100644 --- a/test/ODEsTests/TransientFEProblemsTests/HeatEquationNeumannTests.jl +++ b/test/ODEsTests/TransientFEProblemsTests/HeatEquationNeumannTests.jl @@ -10,12 +10,8 @@ using Gridap.FESpaces using Gridap.ODEs # Analytical functions -u(x, t) = x[1] * (1 - x[2]) * (1 + t) -u(t::Real) = x -> u(x, t) -u(x) = t -> u(x, t) - -∂tu(x, t) = ∂t(u)(x, t) -∂tu(t::Real) = x -> ∂tu(x, t) +ut(t) = x -> x[1] * (1 - x[2]) * (1 + t) +u = TimeSpaceFunction(ut) # Geometry domain = (0, 1, 0, 1) @@ -40,7 +36,8 @@ dΓ = Measure(Γ, degree) nΓ = get_normal_vector(Γ) # FE operator -f(t) = x -> ∂t(u)(x, t) - Δ(u(t))(x) +ft(t) = x -> ∂t(u)(t, x) - Δ(u)(t, x) +f = TimeSpaceFunction(ft) mass(t, ∂ₜu, v) = ∫(∂ₜu ⋅ v) * dΩ mass(t, u, ∂ₜu, v) = mass(t, ∂ₜu, v) stiffness(t, u, v) = ∫(∇(u) ⊙ ∇(v)) * dΩ diff --git a/test/ODEsTests/TransientFEProblemsTests/HeatEquationScalarTests.jl b/test/ODEsTests/TransientFEProblemsTests/HeatEquationScalarTests.jl index 7038afa73..67ed3254c 100644 --- a/test/ODEsTests/TransientFEProblemsTests/HeatEquationScalarTests.jl +++ b/test/ODEsTests/TransientFEProblemsTests/HeatEquationScalarTests.jl @@ -10,11 +10,8 @@ using Gridap.FESpaces using Gridap.ODEs # Analytical functions -u(x, t) = x[1] * (1 - x[2]) * (1 + t) -∂tu(x, t) = ∂t(u)(x, t) - -u(t::Real) = x -> u(x, t) -∂tu(t::Real) = x -> ∂tu(x, t) +ut(t) = x -> x[1] * (1 - x[2]) * (1 + t) +u = TimeSpaceFunction(ut) # Geometry domain = (0, 1, 0, 1) @@ -33,8 +30,8 @@ degree = 2 * order dΩ = Measure(Ω, degree) # FE operator -f(t) = x -> ∂t(u)(x, t) - Δ(u(t))(x) - +ft(t) = x -> ∂t(u)(t, x) - Δ(u)(t, x) +f = TimeSpaceFunction(ft) mass(t, ∂ₜu, v) = ∫(∂ₜu ⋅ v) * dΩ mass(t, u, ∂ₜu, v) = mass(t, ∂ₜu, v) stiffness(t, u, v) = ∫(∇(u) ⊙ ∇(v)) * dΩ @@ -129,6 +126,21 @@ for odeslvr in odeslvrs end end +# Test GeneralizedAlpha1 with initial velocity +odeslvrs = ( + GeneralizedAlpha1(sysslvr_nl, dt, 0.0), + GeneralizedAlpha1(sysslvr_nl, dt, 0.5), + GeneralizedAlpha1(sysslvr_nl, dt, 1.0), +) + +vh0 = interpolate_everywhere(∂t(u)(t0), U0) +uhs0 = (uh0, vh0) +for odeslvr in odeslvrs + for tfeop in tfeops + test_transient_heat_scalar(odeslvr, tfeop, uhs0) + end +end + # Solvers for IMEX decompositions tfeops = ( tfeop_imex_man, diff --git a/test/ODEsTests/TransientFEProblemsTests/HeatEquationVectorTests.jl b/test/ODEsTests/TransientFEProblemsTests/HeatEquationVectorTests.jl index 22f88a6bf..5b9cd79e6 100644 --- a/test/ODEsTests/TransientFEProblemsTests/HeatEquationVectorTests.jl +++ b/test/ODEsTests/TransientFEProblemsTests/HeatEquationVectorTests.jl @@ -10,11 +10,8 @@ using Gridap.FESpaces using Gridap.ODEs # Analytical functions -u(x, t) = VectorValue(x[1] * (1 - x[2]), (1 - x[1]) * x[2]) * (1 + t) -∂tu(x, t) = ∂t(u)(x, t) - -u(t::Real) = x -> u(x, t) -∂tu(t::Real) = x -> ∂tu(x, t) +ut(t) = x -> VectorValue(x[1] * (1 - x[2]), (1 - x[1]) * x[2]) * (1 + t) +u = TimeSpaceFunction(ut) # Geometry domain = (0, 1, 0, 1) @@ -33,8 +30,8 @@ degree = 2 * order dΩ = Measure(Ω, degree) # FE operator -f(t) = x -> ∂t(u)(x, t) - Δ(u(t))(x) - +ft(t) = x -> ∂t(u)(t, x) - Δ(u)(t, x) +f = TimeSpaceFunction(ft) mass(t, ∂ₜu, v) = ∫(∂ₜu ⋅ v) * dΩ mass(t, u, ∂ₜu, v) = mass(t, ∂ₜu, v) stiffness(t, u, v) = ∫(∇(u) ⊙ ∇(v)) * dΩ diff --git a/test/ODEsTests/TransientFEProblemsTests/SecondOrderEquationTests.jl b/test/ODEsTests/TransientFEProblemsTests/SecondOrderEquationTests.jl index c5213610e..8902aa346 100644 --- a/test/ODEsTests/TransientFEProblemsTests/SecondOrderEquationTests.jl +++ b/test/ODEsTests/TransientFEProblemsTests/SecondOrderEquationTests.jl @@ -8,15 +8,9 @@ using Gridap using Gridap.Algebra using Gridap.FESpaces using Gridap.ODEs - # Analytical functions -u(x, t) = (1 - x[1]) * x[2] * (t^2 + 3.0) -∂tu(x, t) = ∂t(u)(x, t) -∂ttu(x, t) = ∂tt(u)(x, t) - -u(t::Real) = x -> u(x, t) -∂tu(t::Real) = x -> ∂tu(x, t) -∂ttu(t::Real) = x -> ∂ttu(x, t) +ut(t) = x -> (1 - x[1]) * x[2] * (t^2 + 3.0) +u = TimeSpaceFunction(ut) # Geometry domain = (0, 1, 0, 1) @@ -36,8 +30,8 @@ dΩ = Measure(Ω, degree) # FE operator order = 2 -f(t) = x -> ∂tt(u)(x, t) + ∂t(u)(x, t) - Δ(u(t))(x) - +ft(t) = x -> ∂tt(u)(t, x) + ∂t(u)(t, x) - Δ(u)(t, x) +f = TimeSpaceFunction(ft) mass(t, ∂ₜₜu, v) = ∫(∂ₜₜu ⋅ v) * dΩ mass(t, u, ∂ₜₜu, v) = mass(t, ∂ₜₜu, v) damping(t, ∂ₜu, v) = ∫(∂ₜu ⋅ v) * dΩ @@ -96,8 +90,7 @@ dt = 0.1 U0 = U(t0) uh0 = interpolate_everywhere(u(t0), U0) -∂tuh0 = interpolate_everywhere(∂tu(t0), U0) -∂ttuh0 = interpolate_everywhere(∂ttu(t0), U0) +∂tuh0 = interpolate_everywhere(∂t(u)(t0), U0) tol = 1.0e-6 sysslvr_l = LUSolver() @@ -125,4 +118,13 @@ for odeslvr in odeslvrs end end +# Test with initial acceleration +∂ttuh0 = interpolate_everywhere(∂tt(u)(t0), U0) +uhs0 = (uh0, ∂tuh0, ∂ttuh0) +for odeslvr in odeslvrs + for tfeop in tfeops + test_tfeop_order2(odeslvr, tfeop, uhs0) + end +end + end # module Order2FETests diff --git a/test/ODEsTests/TransientFEProblemsTests/StokesEquationTests.jl b/test/ODEsTests/TransientFEProblemsTests/StokesEquationTests.jl index c9f9702a4..c78227991 100644 --- a/test/ODEsTests/TransientFEProblemsTests/StokesEquationTests.jl +++ b/test/ODEsTests/TransientFEProblemsTests/StokesEquationTests.jl @@ -10,12 +10,11 @@ using Gridap.FESpaces using Gridap.ODEs # Analytical functions -u(x, t) = VectorValue(x[1], x[2]) * t -u(t::Real) = x -> u(x, t) +ut(t) = x -> VectorValue(x[1], x[2]) * t +u = TimeSpaceFunction(ut) -p(x, t) = (x[1] - x[2]) * t -p(t::Real) = x -> p(x, t) -q(x) = t -> p(x, t) +pt(t) = x -> (x[1] - x[2]) * t +p = TimeSpaceFunction(pt) # Geometry domain = (0, 1, 0, 1) @@ -41,8 +40,9 @@ degree = 2 * order dΩ = Measure(Ω, degree) # FE operator -f(t) = x -> ∂t(u)(t)(x) - Δ(u(t))(x) + ∇(p(t))(x) -g(t) = x -> (∇ ⋅ u(t))(x) +ft(t) = x -> ∂t(u)(t, x) - Δ(u)(t, x) + ∇(p)(t, x) +f = TimeSpaceFunction(ft) +g = ∇ ⋅ u mass(t, ∂ₜu, v) = ∫(∂ₜu ⋅ v) * dΩ stiffness(t, u, v) = ∫(∇(u) ⊙ ∇(v)) * dΩ forcing(t, (v, q)) = ∫(f(t) ⋅ v) * dΩ + ∫(g(t) * q) * dΩ diff --git a/test/ODEsTests/TransientFESpacesTests.jl b/test/ODEsTests/TransientFESpacesTests.jl index a2ee48d3b..6130d1c99 100644 --- a/test/ODEsTests/TransientFESpacesTests.jl +++ b/test/ODEsTests/TransientFESpacesTests.jl @@ -6,23 +6,8 @@ using Gridap using Gridap.Fields using Gridap.ODEs -u1(x, t) = (x[1] + x[2]) * t -u1(t::Real) = x -> u1(x, t) - -∂tu1(t) = x -> x[1] + x[2] -ODEs.∂t(::typeof(u1)) = ∂tu1 - -∂ttu1(t) = x -> zero(x[1]) -ODEs.∂tt(::typeof(u1)) = ∂ttu1 - -u2(x, t) = x[1] * t^2 + x[2] * t -u2(t::Real) = x -> u2(x, t) - -∂tu2(t) = x -> 2 * t * x[1] + x[2] -ODEs.∂t(::typeof(u2)) = ∂tu2 - -∂ttu2(t) = x -> 2 * x[1] -ODEs.∂tt(::typeof(u2)) = ∂ttu2 +u1(t) = x -> (x[1] + x[2]) * t +u2(t) = x -> x[1] * t^2 + x[2] * t domain = (0, 1, 0, 1) partition = (5, 5)