From 44d1c2a4229080f8be19da2b08c94af97eee21f9 Mon Sep 17 00:00:00 2001
From: Alexandre Magueresse <alexandre.magueresse@monash.edu>
Date: Fri, 12 Apr 2024 17:04:53 +1000
Subject: [PATCH 1/2] Changed sign for the residual in
 TransientLinearFEOperator

---
 docs/src/ODEs.md                                          | 8 ++++----
 src/ODEs/ODEOperators.jl                                  | 2 +-
 src/ODEs/ODEOpsFromTFEOps.jl                              | 4 +++-
 src/ODEs/TransientFEOperators.jl                          | 2 +-
 test/ODEsTests/ODEOperatorsMocks.jl                       | 8 ++++----
 test/ODEsTests/ODEOperatorsTests.jl                       | 4 ++--
 test/ODEsTests/ODEProblemsTests/Order1ODETests.jl         | 2 +-
 test/ODEsTests/ODEProblemsTests/Order2ODETests.jl         | 2 +-
 test/ODEsTests/ODESolversTests.jl                         | 2 +-
 test/ODEsTests/TransientFEOperatorsSolutionsTests.jl      | 4 ++--
 .../HeatEquationMultiFieldTests.jl                        | 2 +-
 .../TransientFEProblemsTests/HeatEquationNeumannTests.jl  | 2 +-
 .../TransientFEProblemsTests/HeatEquationScalarTests.jl   | 2 +-
 .../TransientFEProblemsTests/HeatEquationVectorTests.jl   | 2 +-
 .../TransientFEProblemsTests/SecondOrderEquationTests.jl  | 2 +-
 15 files changed, 25 insertions(+), 23 deletions(-)

diff --git a/docs/src/ODEs.md b/docs/src/ODEs.md
index 6ed8da9fd..3701457b9 100644
--- 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
index 70ebfe6a1..0727277ad 100644
--- 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
index 56bb775ec..dd56e2725 100644
--- 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
index 0824cb556..0d3455a87 100644
--- 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
index 88d0f6562..4ce6ebd85 100644
--- 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
index 4e448d6db..411a68c1e 100644
--- 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
index 645116e9a..53f9ae548 100644
--- 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
index 239965a36..62e872259 100644
--- 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
index 0273bc814..5055d45c6 100644
--- 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
index 3bb7024aa..e176928fb 100644
--- 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
index 234ceb412..ed7c0616b 100644
--- 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
index 1b6b68131..167b77d03 100644
--- 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
index 67ed3254c..d7d0ddd57 100644
--- 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
index 5b9cd79e6..86c75ab06 100644
--- 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
index 8902aa346..6cc04f36f 100644
--- 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Ω

From 5801cb786da853957c6a0f6852284ae38e9ab270 Mon Sep 17 00:00:00 2001
From: Alexandre Magueresse <alexandre.magueresse@monash.edu>
Date: Fri, 12 Apr 2024 17:40:35 +1000
Subject: [PATCH 2/2] Updated NEWS.md

---
 NEWS.md | 4 ++++
 1 file changed, 4 insertions(+)

diff --git a/NEWS.md b/NEWS.md
index 0409808b5..d928c1198 100644
--- 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).