Renamed TimeSlicing into TimeSpaceFunction and updated constructor

Author: AlexandreMagueresse

docs/src/ODEs.md

@@ -91,8 +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.
 ```
-gtx(t, x) = x[1] + x[2] * t
-g = time_slicing(gtx)
+g(t) = x -> x[1] + x[2] * t
 V = FESpace(model, reffe, dirichlet_tags="boundary")
 U = TransientTrialFESpace (V, g)
 
@@ -106,13 +105,6 @@ U0 = U(t0)
 ∂ttU0 = ∂ttU(t0)
 ```
 
-## The `time_slicing` constructor
-Note that the constructor `time_slicing` was used above to prescribe the boundary condition, instead of a more natural function `f(t, x) = gtx(t, x)` or functor `F(t) = x -> gtx(t, x)`. It would be possible to provide `f` instead of `g` in this example, but it would not be possible to do so for time-dependent coefficients inside an integrand, simply because `Gridap` expects an integrand to be a function of space only. This `time_slicing` constructor unifies these two scenarios and is motivated by performance and convenience reasons:
-* With the functor approach, `h(t)` is an anonymous function of `x`. This leads to performance drops when evaluating `h(t)(x)`, even when the function `h(t)` is instantiated only once, i.e. creating `ht = h(t)` and calling `ht(x)`. An even more significant loss in performance happens when computing (space or time) derivatives with automatic differentiation. With the `time_slicing` constructor, `g(t)` is a named function of `x`, and brings the same performance as if calling `gtx`.
-* With the functor approach, applying spatial differential operators is somewhat cumbersome, for example one would have to write `∂t(h)(t)(x) - Δ(h(t))(x)`. With the `time_slicing` constructor, one can simply write `∂t(g)(t, x) - Δ(g)(t, x)`.
-
-As a summary, `g = time_slicing(gtx)` allows the syntax `op(g)`, `op(g)(t)` and `op(g)(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`, `∇`, ...).
-
 ## Cell fields
 The time-dependent equivalent of `CellField` is `TransientCellField`. It stores the cell field itself together with its derivatives up to the order of the ODE.
 
@@ -163,7 +155,15 @@ TransientLinearFEOperator((stiffness, mass), res, U, V, constant_forms=(false, t
 ```
 If $\kappa$ is constant, the keyword `constant_forms` could be replaced by `(true, true)`.
 
-**Important** Note that in these examples, for optimal performance, `κ` should be a `time_slicing`.
+## 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.

src/Exports.jl

@@ -182,7 +182,6 @@ using Gridap.CellData: ∫; export ∫
 @publish Visualization createpvd
 @publish Visualization savepvd
 
-@publish ODEs time_slicing
 @publish ODEs ∂t
 @publish ODEs ∂tt
 @publish ODEs ForwardEuler

src/ODEs/ODEs.jl

@@ -35,9 +35,7 @@ const ε = 100 * eps()
 
 include("TimeDerivatives.jl")
 
-export TimeSlicing
-export time_slicing
-
+export TimeSpaceFunction
 export time_derivative
 export ∂t
 export ∂tt

src/ODEs/TimeDerivatives.jl

@@ -1,47 +1,40 @@
-################
-# Time slicing #
-################
+#####################
+# TimeSpaceFunction #
+#####################
 """
-    struct TimeSlicing{F} <: Function end
-
-`TimeSlicing` is an operator that can be applied to a function of time and space, and
-that, to a given time, produces the function of space at that time, allowing for the
-following syntax: if `f(t, x)` is defined, then `TimeSlicing(f)(t)(x) = f(t, x)`.
-
-This is needed prevent performance drops with automatic differentiation (e.g. when taking
-the time derivative of Dirichlet boundary conditions), and with integration (i.e. when
-an integrand contains a time-dependent coefficient).
+    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 TimeSlicing{F} <: Function
+struct TimeSpaceFunction{F} <: Function
   f::F
 end
 
-function (ts::TimeSlicing)(t)
-  _ftx(x) = ts.f(t, x)
-  _ftx
-end
-
-function (ts::TimeSlicing)(t, x)
-  ts.f(t, x)
-end
-
-"""
-    const time_slicing = TimeSlicing
-
-Alias for `TimeSlicing`.
-"""
-const time_slicing = TimeSlicing
+(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)
+# 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::TimeSlicing)
-      _op(t, x) = $op(ts(t))(x)
-      TimeSlicing(_op)
+    function ($op)(ts::TimeSpaceFunction)
+      function _op(t)
+        _op_t(x) = $op(ts(t))(x)
+        _op_t
+      end
+      TimeSpaceFunction(_op)
     end
   end
 end
@@ -112,33 +105,37 @@ end
 # Specialisation for `Function` #
 #################################
 function time_derivative(f::Function)
-  function dfdt(t, x)
-    T = return_type(f, t, x)
-    _time_derivative(T, f, t, x)
+  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
   dfdt
 end
 
 function _time_derivative(T::Type{<:Real}, f, t, x)
-  partial(t) = f(t, x)
+  partial(t) = f(t)(x)
   ForwardDiff.derivative(partial, t)
 end
 
 function _time_derivative(T::Type{<:VectorValue}, f, t, x)
-  partial(t) = get_array(f(t, x))
+  partial(t) = get_array(f(t)(x))
   VectorValue(ForwardDiff.derivative(partial, t))
 end
 
 function _time_derivative(T::Type{<:TensorValue}, f, t, x)
-  partial(t) = get_array(f(t, x))
+  partial(t) = get_array(f(t)(x))
   TensorValue(ForwardDiff.derivative(partial, t))
 end
 
-####################################
-# Specialisation for `TimeSlicing` #
-####################################
-function time_derivative(ts::TimeSlicing)
-  TimeSlicing(time_derivative(ts.f))
+##########################################
+# Specialisation for `TimeSpaceFunction` #
+##########################################
+function time_derivative(ts::TimeSpaceFunction)
+  TimeSpaceFunction(time_derivative(ts.f))
 end
 
 ###############################

test/ODEsTests/TimeDerivativesTests.jl

@@ -8,100 +8,100 @@ using Gridap
 using Gridap.ODEs
 
 # First time derivative, scalar-valued
-f1(t, x) = 5 * x[1] * x[2] + x[2]^2 * t^3
-∂tf1(t, x) = 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(t, x) = t^2
-∂tf2(t, x) = 2 * t
+f2(t) = x -> t^2
+∂tf2(t) = x -> 2 * t
 
-f3(t, x) = x[1]^2
-∂tf3(t, x) = zero(x[1])
+f3(t) = x -> x[1]^2
+∂tf3(t) = x -> zero(x[1])
 
-f4(t, x) = x[1]^t^2
-∂tf4(t, x) = 2 * t * log(x[1]) * f4(t, x)
+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(t, x) = ForwardDiff.derivative(t -> f(t, x), t)
+  dtf(t) = x -> ForwardDiff.derivative(t -> f(t)(x), t)
 
   tv = rand(Float64)
   xv = Point(rand(Float64, 2)...)
-  @test ∂t(f)(tv, xv) ≈ ∂tf(tv, xv)
-  @test ∂t(f)(tv, xv) ≈ dtf(tv, xv)
+  @test ∂t(f)(tv)(xv) ≈ ∂tf(tv)(xv)
+  @test ∂t(f)(tv)(xv) ≈ dtf(tv)(xv)
 
-  slicing = time_slicing(f)
-  @test slicing(tv)(xv) ≈ f(tv, xv)
-  @test ∂t(slicing)(tv)(xv) ≈ ∂tf(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(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)
+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(t, x) = VectorValue(x[1]^2, zero(x[2]))
-∂tf2(t, x) = 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(t, x) = VectorValue(x[1]^2, t)
-∂tf3(t, x) = 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(t, x) = VectorValue(ForwardDiff.derivative(t -> get_array(f(t, x)), 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)(tv, xv) ≈ ∂tf(tv, xv)
-  @test ∂t(f)(tv, xv) ≈ dtf(tv, xv)
+  @test ∂t(f)(tv)(xv) ≈ ∂tf(tv)(xv)
+  @test ∂t(f)(tv)(xv) ≈ dtf(tv)(xv)
 
-  slicing = time_slicing(f)
-  @test slicing(tv)(xv) ≈ f(tv, xv)
-  @test ∂t(slicing)(tv)(xv) ≈ ∂tf(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(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)
+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(t, x) = TensorValue(ForwardDiff.derivative(t -> get_array(f(t, x)), 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)(tv, xv) ≈ ∂tf(tv, xv)
-  @test ∂t(f)(tv, xv) ≈ dtf(tv, xv)
+  @test ∂t(f)(tv)(xv) ≈ ∂tf(tv)(xv)
+  @test ∂t(f)(tv)(xv) ≈ dtf(tv)(xv)
 
-  slicing = time_slicing(f)
-  @test slicing(tv)(xv) ≈ f(tv, xv)
-  @test ∂t(slicing)(tv)(xv) ≈ ∂tf(tv, xv)
+  F = TimeSpaceFunction(f)
+  @test F(tv)(xv) ≈ f(tv)(xv)
+  @test ∂t(F)(tv)(xv) ≈ ∂tf(tv)(xv)
 end
 
 # Spatial derivatives
-ftx(t, x) = x[1]^2 * t + x[2]
-f = time_slicing(ftx)
+ft(t) = x -> x[1]^2 * t + x[2]
+f = TimeSpaceFunction(ft)
 
 tv = rand(Float64)
 xv = Point(rand(Float64, 2)...)
-@test ∇(f)(tv, xv) ≈ Point(2 * xv[1] * tv, 1.0)
+@test ∇(f)(tv)(xv) ≈ Point(2 * xv[1] * tv, 1.0)
 
 # Second time derivative, scalar-valued
-f1(t, x) = t^2
-∂tf1(t, x) = 2 * t
-∂ttf1(t, x) = 2 * one(t)
+f1(t) = x -> t^2
+∂tf1(t) = x -> 2 * t
+∂ttf1(t) = x -> 2 * one(t)
 
-f2(t, x) = x[1] * t^2
-∂tf2(t, x) = 2 * x[1] * t
-∂ttf2(t, x) = 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(t, x) = ForwardDiff.derivative(t -> f(t, x), t)
-  dttf(t, x) = ForwardDiff.derivative(t -> dtf(t, x), 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)(tv, xv) ≈ ∂ttf(tv, xv)
-  @test ∂tt(f)(tv, xv) ≈ dttf(tv, xv)
+  @test ∂tt(f)(tv)(xv) ≈ ∂ttf(tv)(xv)
+  @test ∂tt(f)(tv)(xv) ≈ dttf(tv)(xv)
 
-  slicing = time_slicing(f)
-  @test slicing(tv)(xv) ≈ f(tv, xv)
-  @test ∂tt(slicing)(tv)(xv) ≈ ∂ttf(tv, xv)
+  F = TimeSpaceFunction(f)
+  @test F(tv)(xv) ≈ f(tv)(xv)
+  @test ∂tt(F)(tv)(xv) ≈ ∂ttf(tv)(xv)
 end
 
 end # module TimeDerivativesTests

test/ODEsTests/TransientCellFieldsTests.jl

@@ -12,8 +12,7 @@ using Gridap.MultiField
 using Gridap.ODEs
 using Gridap.ODEs: TransientMultiFieldCellField
 
-ftx(t, x) = sum(x)
-f = time_slicing(ftx)
+f(t) = x -> sum(x)
 
 domain = (0, 1, 0, 1)
 partition = (5, 5)

test/ODEsTests/TransientFEOperatorsSolutionsTests.jl

@@ -15,8 +15,8 @@ using Gridap.ODEs
 include("ODESolversMocks.jl")
 
 # Analytical functions
-utx(t, x) = x[1] * (1 - x[2]) * (1 + t)
-u = time_slicing(utx)
+ut(t) = x -> x[1] * (1 - x[2]) * (1 + t)
+u = TimeSpaceFunction(ut)
 
 # Geometry
 domain = (0, 1, 0, 1)
@@ -104,8 +104,8 @@ end
 # First-order #
 ###############
 order = 1
-ftx(t, x) = ∂t(u)(t)(x) - Δ(u(t))(x)
-f = time_slicing(ftx)
+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)
@@ -200,8 +200,8 @@ end
 # Second-order #
 ################
 order = 2
-ftx(t, x) = ∂tt(u)(t)(x) + ∂t(u)(t)(x) - Δ(u(t))(x)
-f = time_slicing(ftx)
+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)

test/ODEsTests/TransientFEProblemsTests/FreeSurfacePotentialFlowTests.jl

@@ -23,10 +23,8 @@ tF = 10 * dt # 2 * π
 tol = 1.0e-2
 
 # Exact solution
-ϕₑtx(t, x) = ω / k * ξ * (cosh(k * (x[2]))) / sinh(k * H) * sin(k * x[1] - ω * t)
-ϕₑ = time_slicing(ϕₑtx)
-ηₑtx(t, x) = ξ * cos(k * x[1] - ω * t)
-ηₑ = time_slicing(ηₑtx)
+ϕₑ(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)