Merge pull request #996 from gridap/ode-residual-sign

Change sign of residual in `TransientLinearFEOperator`
gridap · Apr 12, 2024 · 532e8de · 532e8de · JordiManyer · Apr 14, 2024
2 parents 80fabe9 + 5801cb7
commit 532e8de
Show file tree

Hide file tree

Showing 16 changed files with 29 additions and 23 deletions.
diff --git a/NEWS.md b/NEWS.md
@@ -7,6 +7,10 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ## [0.18.1] - 2024-04-12
 
+### Changed
+
+- Changed the sign of the residual in `TransientLinearFEOperator` to align with the conventions of `AffineFEOperator`. Since PR[#996](https://github.com/gridap/Gridap.jl/pull/996).
+
 ### Fixed
 
 - Bugfix in `restrict_to_field` for `BlockMultiFieldStyle`. Since PR[#993](https://github.com/gridap/Gridap.jl/pull/993).

diff --git a/docs/src/ODEs.md b/docs/src/ODEs.md
@@ -62,9 +62,9 @@ We call the matrix ``\boldsymbol{M}: \mathbb{R} \to \mathbb{R}^{d \times d}`` th
 In particular, a semilinear ODE is a quasilinear ODE.
 * **Linear**. The residual is linear with respect to all time derivatives, i.e.
 ```math
-\boldsymbol{r}(t, \partial_{t}^{0} \boldsymbol{u}, \ldots, \partial_{t}^{n} \boldsymbol{u}) = \sum_{0 \leq k \leq n} \boldsymbol{A}_{k}(t) \partial_{t}^{k} \boldsymbol{u} + \boldsymbol{f}(t).
+\boldsymbol{r}(t, \partial_{t}^{0} \boldsymbol{u}, \ldots, \partial_{t}^{n} \boldsymbol{u}) = \sum_{0 \leq k \leq n} \boldsymbol{A}_{k}(t) \partial_{t}^{k} \boldsymbol{u} - \boldsymbol{f}(t).
 ```
-We refer to the matrix ``\boldsymbol{A}_{k}: \mathbb{R} \to \mathbb{R}^{d \times d}`` as the ``k``-th linear form of the residual. We may still define the mass matrix ``\boldsymbol{M} = \boldsymbol{A}_{n}``. In particular, a linear ODE is a semilinear ODE.
+We refer to the matrix ``\boldsymbol{A}_{k}: \mathbb{R} \to \mathbb{R}^{d \times d}`` as the ``k``-th linear form of the residual. We may still define the mass matrix ``\boldsymbol{M} = \boldsymbol{A}_{n}``. Note that the term independent of $u$, i.e. the forcing term, is subtracted from the residual. This aligns with standard conventions, and in particular with those of `AffineFEOperator` (see example in _Finite element operators_ below, in the construction of a `TransientLinearFEOperator`). In particular, a linear ODE is a semilinear ODE.
 
 > Note that for residuals of order zero (i.e. "standard" systems of equations), the definitions of quasilinear, semilinear, and linear coincide.
 
@@ -159,7 +159,7 @@ The time-dependent analog of `FEOperator` is `TransientFEOperator`. It has the f
 * `TransientSemilinearFEOperator(mass, res, jacs, trial, test; constant_mass)` and `TransientSemilinearFEOperator(mass, res, trial, test; order, constant_mass)` for the version with automatic jacobians. (The jacobian with respect to ``\partial_{t}^{n} \boldsymbol{u}`` is simply the mass term). The mass and residual are expected to have the signatures `mass(t, dtNu, v)` and `residual(t, u, v)`, where here again `dtNu` is the highest-order derivative. In particular, the mass is specified as a linear form of `dtNu`.
 * `TransientLinearFEOperator(forms, res, jacs, trial, test; constant_forms)` and `TransientLinearFEOperator(forms, res, trial, test; constant_forms)` for the version with automatic jacobians. (In fact, the jacobians are simply the forms themselves). The forms and residual are expected to have the signatures `form_k(t, dtku, v)` and `residual(t, v)`, i.e. `form_k` is a linear form of the ``k``-th order derivative, and the residual does not depend on `u`.
 
-It is important to note that all the terms are gathered in the residual, including the forcing term. In the common case where the forcing term is on the right-hand side, it will need to be made negative in this description.
+It is important to note that all the terms are gathered in the residual, including the forcing term. In the common case where the ODE is linear, the residual is only the forcing term, and it is subtracted from the bilinear forms (see example below).
 
 Here, in the signature of the residual, `t` is the time at which the residual is evaluated, `u` is a function in the trial space, and `v` is a test function. Time derivatives of `u` can be included in the residual via the `∂t` operator, applied as many times as needed, or using the shortcut `∂t(u, N)`.
 
@@ -194,7 +194,7 @@ TransientSemilinearFEOperator(mass, res, U, V, constant_mass=true)
 ```
 stiffness(t, u, v) = ∫( ∇(v) ⋅ (κ(t) ⋅ ∇(u)) ) dΩ
 mass(t, dtu, v) = ∫( v ⋅ dtu ) dΩ
-res(t, u, v) = -(∫( v ⋅ f(t) ) dΩ)
+res(t, u, v) = ∫( v ⋅ f(t) ) dΩ
 TransientLinearFEOperator((stiffness, mass), res, U, V, constant_forms=(false, true))
 ```
 If ``\kappa`` is constant, the keyword `constant_forms` could be replaced by `(true, true)`.

diff --git a/src/ODEs/ODEOperators.jl b/src/ODEs/ODEOperators.jl
@@ -44,7 +44,7 @@ struct SemilinearODE <: AbstractSemilinearODE end
 
 ODE operator whose residual is linear with respect to all time derivatives, i.e.
 ```
-residual(t, ∂t^0[u], ..., ∂t^N[u]) = ∑_{0 ≤ k ≤ N} A_k(t) ∂t^k[u] + res(t),
+residual(t, ∂t^0[u], ..., ∂t^N[u]) = ∑_{0 ≤ k ≤ N} A_k(t) ∂t^k[u] - res(t),
 ```
 where `N` is the order of the ODE operator, and `∂t^k[u]` is the `k`-th-order
 time derivative of `u`.

diff --git a/src/ODEs/ODEOpsFromTFEOps.jl b/src/ODEs/ODEOpsFromTFEOps.jl
@@ -260,7 +260,9 @@ function Algebra.residual!(
 
   # Residual
   res = get_res(odeop.tfeop)
-  dc = res(t, uh, v)
+  # Need a negative sign here:
+  # residual(t, u, v) = ∑_{0 ≤ k ≤ N} form_k(t, ∂t^k[u], v) - res(t, v)
+  dc = (-1) * res(t, uh, v)
 
   # Forms
   order = get_order(odeop)

diff --git a/src/ODEs/TransientFEOperators.jl b/src/ODEs/TransientFEOperators.jl
@@ -643,7 +643,7 @@ get_assembler(tfeop::TransientSemilinearFEOpFromWeakForm) = tfeop.assembler
 
 Transient `FEOperator` defined by a transient weak form
 ```
-residual(t, u, v) = ∑_{0 ≤ k ≤ N} form_k(t, ∂t^k[u], v) + res(t, v) = 0,
+residual(t, u, v) = ∑_{0 ≤ k ≤ N} form_k(t, ∂t^k[u], v) - res(t, v) = 0,
 ```
 where `N` is the order of the operator, `form_k` is linear in `∂t^k[u]` and
 does not depend on the other time derivatives of `u`, and the `form_k` and

diff --git a/test/ODEsTests/ODEOperatorsMocks.jl b/test/ODEsTests/ODEOperatorsMocks.jl
@@ -14,7 +14,7 @@ using Gridap.ODEs
 
 Mock linear ODE of arbitrary order
 ```
-∑_{0 ≤ k ≤ N} form_k(t) ∂t^k u + forcing(t) = 0.
+∑_{0 ≤ k ≤ N} form_k(t) ∂t^k u - forcing(t) = 0.
 ```
 """
 struct ODEOperatorMock{T} <: ODEOperator{T}
@@ -43,7 +43,7 @@ function Algebra.residual!(
 )
   order = get_order(odeop)
   !add && fill!(r, zero(eltype(r)))
-  axpy!(1, odeop.forcing(t), r)
+  axpy!(-1, odeop.forcing(t), r)
   for k in 0:order
     mat = odeop.forms[k+1](t)
     axpy!(1, mat * us[k+1], r)
@@ -58,7 +58,7 @@ function Algebra.residual!(
 )
   order = get_order(odeop)
   !add && fill!(r, zero(eltype(r)))
-  axpy!(1, odeop.forcing(t), r)
+  axpy!(-1, odeop.forcing(t), r)
   for k in 0:order-1
     mat = odeop.forms[k+1](t)
     axpy!(1, mat * us[k+1], r)
@@ -76,7 +76,7 @@ function Algebra.residual!(
 )
   order = get_order(odeop)
   !add && fill!(r, zero(eltype(r)))
-  axpy!(1, odeop.forcing(t), r)
+  axpy!(-1, odeop.forcing(t), r)
   for k in 0:order
     mat = odeop.forms[k+1](t)
     axpy!(1, mat * us[k+1], r)

diff --git a/test/ODEsTests/ODEOperatorsTests.jl b/test/ODEsTests/ODEOperatorsTests.jl
@@ -60,13 +60,13 @@ for N in 0:order_max
       for T_ex in Ts
         ex_odeop = ODEOperatorMock{T_ex}(ex_forms, ex_forcing)
         imex_odeop = IMEXODEOperator(im_odeop, ex_odeop)
-        # odeops = (odeops..., imex_odeop)
+        odeops = (odeops..., imex_odeop)
       end
     end
 
     # Compute expected residual
     f = forcing(t)
-    copy!(exp_r, f)
+    copy!(exp_r, -f)
     for (ui, formi) in zip(us, forms)
       form = formi(t)
       exp_r .+= form * ui

diff --git a/test/ODEsTests/ODEProblemsTests/Order1ODETests.jl b/test/ODEsTests/ODEProblemsTests/Order1ODETests.jl
@@ -29,7 +29,7 @@ nonzeros(mat0) .= 0
 form_zero(t) = mat0
 
 α = randn(num_eqs)
-forcing(t) = -M * exp.(α .* t)
+forcing(t) = M * exp.(α .* t)
 forcing_zero(t) = zeros(typeof(t), num_eqs)
 
 u0 = randn(num_eqs)

diff --git a/test/ODEsTests/ODEProblemsTests/Order2ODETests.jl b/test/ODEsTests/ODEProblemsTests/Order2ODETests.jl
@@ -32,7 +32,7 @@ nonzeros(mat0) .= 0
 form_zero(t) = mat0
 
 α = randn(num_eqs)
-forcing(t) = -M * exp.(α .* t)
+forcing(t) = M * exp.(α .* t)
 forcing_zero(t) = zeros(typeof(t), num_eqs)
 
 u0 = randn(num_eqs)

diff --git a/test/ODEsTests/ODESolversTests.jl b/test/ODEsTests/ODESolversTests.jl
@@ -51,7 +51,7 @@ for N in 1:order_max
 
   f = forcing(tx)
   fill!(exp_r, zero(eltype(exp_r)))
-  exp_r .+= f
+  exp_r .= -f
 
   m = last(forms)(tx)
   fillstored!(exp_J, zero(eltype(exp_J)))

diff --git a/test/ODEsTests/TransientFEOperatorsSolutionsTests.jl b/test/ODEsTests/TransientFEOperatorsSolutionsTests.jl
@@ -117,7 +117,7 @@ jac(t, u, du, v) = stiffness(t, du, v)
 jac_t(t, u, dut, v) = mass(t, u, dut, v)
 
 res_ql(t, u, v) = stiffness(t, u, v) - forcing(t, v)
-res_l(t, v) = (-1) * forcing(t, v)
+res_l(t, v) = forcing(t, v)
 
 # TODO could think of a simple and optimised way to create a zero residual or
 # jacobian without assembling the vector / matrix
@@ -215,7 +215,7 @@ jac_t(t, u, dut, v) = damping(t, dut, v)
 jac_tt(t, u, dutt, v) = mass(t, dutt, v)
 
 res_ql(t, u, v) = damping(t, ∂t(u), v) + stiffness(t, u, v) - forcing(t, v)
-res_l(t, v) = (-1) * forcing(t, v)
+res_l(t, v) = forcing(t, v)
 
 im_res(t, u, v) = ∫(0 * u * v) * dΩ
 im_jac(t, u, du, v) = ∫(0 * du * v) * dΩ

diff --git a/test/ODEsTests/TransientFEProblemsTests/HeatEquationMultiFieldTests.jl b/test/ODEsTests/TransientFEProblemsTests/HeatEquationMultiFieldTests.jl
@@ -42,7 +42,7 @@ _forcing(t, v) = ∫(f(t) ⋅ v) * dΩ
 
 _res(t, u, v) = _mass(t, ∂t(u), v) + _stiffness(t, u, v) - _forcing(t, v)
 _res_ql(t, u, v) = _stiffness(t, u, v) - _forcing(t, v)
-_res_l(t, v) = (-1) * _forcing(t, v)
+_res_l(t, v) = _forcing(t, v)
 
 mass(t, (∂ₜu1, ∂ₜu2), (v1, v2)) = _mass(t, ∂ₜu1, v1) + _mass(t, ∂ₜu2, v2)
 mass(t, (u1, u2), (∂ₜu1, ∂ₜu2), (v1, v2)) = _mass(t, u1, ∂ₜu1, v1) + _mass(t, u2, ∂ₜu2, v2)

diff --git a/test/ODEsTests/TransientFEProblemsTests/HeatEquationNeumannTests.jl b/test/ODEsTests/TransientFEProblemsTests/HeatEquationNeumannTests.jl
@@ -50,7 +50,7 @@ jac(t, u, du, v) = stiffness(t, du, v)
 jac_t(t, u, dut, v) = mass(t, dut, v)
 
 res_ql(t, u, v) = stiffness(t, u, v) - forcing(t, v)
-res_l(t, v) = (-1) * forcing(t, v)
+res_l(t, v) = forcing(t, v)
 
 args_man = ((jac, jac_t), U, V)
 tfeop_nl_man = TransientFEOperator(res, args_man...)

diff --git a/test/ODEsTests/TransientFEProblemsTests/HeatEquationScalarTests.jl b/test/ODEsTests/TransientFEProblemsTests/HeatEquationScalarTests.jl
@@ -42,7 +42,7 @@ jac(t, u, du, v) = stiffness(t, du, v)
 jac_t(t, u, dut, v) = mass(t, u, dut, v)
 
 res_ql(t, u, v) = stiffness(t, u, v) - forcing(t, v)
-res_l(t, v) = (-1) * forcing(t, v)
+res_l(t, v) = forcing(t, v)
 
 res0(t, u, v) = ∫(0 * u * v) * dΩ
 jac0(t, u, du, v) = ∫(0 * du * v) * dΩ

diff --git a/test/ODEsTests/TransientFEProblemsTests/HeatEquationVectorTests.jl b/test/ODEsTests/TransientFEProblemsTests/HeatEquationVectorTests.jl
@@ -42,7 +42,7 @@ jac(t, u, du, v) = stiffness(t, du, v)
 jac_t(t, u, dut, v) = mass(t, u, dut, v)
 
 res_ql(t, u, v) = stiffness(t, u, v) - forcing(t, v)
-res_l(t, v) = (-1) * forcing(t, v)
+res_l(t, v) = forcing(t, v)
 
 args_man = ((jac, jac_t), U, V)
 tfeop_nl_man = TransientFEOperator(res, args_man...)

diff --git a/test/ODEsTests/TransientFEProblemsTests/SecondOrderEquationTests.jl b/test/ODEsTests/TransientFEProblemsTests/SecondOrderEquationTests.jl
@@ -44,7 +44,7 @@ jac_t(t, u, dut, v) = damping(t, dut, v)
 jac_tt(t, u, dutt, v) = mass(t, dutt, v)
 
 res_ql(t, u, v) = damping(t, ∂t(u), v) + stiffness(t, u, v) - forcing(t, v)
-res_l(t, v) = (-1) * forcing(t, v)
+res_l(t, v) = forcing(t, v)
 
 res0(t, u, v) = ∫(0 * u * v) * dΩ
 jac0(t, u, du, v) = ∫(0 * du * v) * dΩ