Optimization in transient linear FE operator for zero Dirichlet value case

[oriolcg]

Hi @JordiManyer, I propose an optimization for the case of LinearODEOperator with zero Dirichlet values. See explanation below (AI-generated).

---

A linear transient problem in Gridap is written as

Σ_{k=0}^{N} form_k(t, ∂tᵏu, v) = res(t, v)

and built with TransientLinearFEOperator(forms, res, U, V; constant_forms). When a form is declared constant, allocate_odeopcache assembles its matrix A_k once and stores it in odeopcache.const_forms. For a linear problem with constant forms, every step still performs a full quadrature pass over all cells for every bilinear form, only to obtain a
vector that equals Σ_k A_k ∂tᵏu.

Since the system matrix is factorised once, this assembly becomes the dominant per-step cost (tested time of 3-4x slower than only doing algebraic operations).

Why Gridap cannot simply use A_k · ∂tᵏu

The identity

form_k(t, ∂tᵏuₕ, v) = A_k · us[k+1]     (on the free test dofs)

holds only if ∂tᵏuₕ is fully determined by its free values, i.e. if the space of the k-th time derivative has zero Dirichlet values. A_k only stores the free × free block; with non-zero Dirichlet values the residual also contains the lifting term A_k^{free,Dir} · ∂tᵏu_D(t), which is not stored. constant_forms only states that the bilinear form does not depend on time; it says nothing about the boundary data. So the stock code takes the generic, always-correct path.

The proposed change

In Algebra.residual!(r, ::ODEOpFromTFEOp{<:AbstractLinearODE}, …)
(src/ODEs/ODEOpsFromTFEOps.jl), each form is treated separately:

for k in 0:order
  const_form = odeopcache.const_forms[k+1]
  if !isnothing(const_form) && _has_zero_dirichlet_values(odeopcache.Us[k+1])
    mul!(r, const_form, us[k+1], true, true)      # stored matrix × free values
  else
    dc = dc + forms[k+1](t, ∂tkuh, v)             # re-assembled, as before
  end
  ...
end

with

_has_zero_dirichlet_values(U) = false             # unknown space: be safe
_has_zero_dirichlet_values(U::SingleFieldFESpace) =
  (num_dirichlet_dofs(U) == 0) || iszero(get_dirichlet_dof_values(U))
_has_zero_dirichlet_values(U::MultiFieldFESpace) = all(_has_zero_dirichlet_values, U)

The forcing res(t, v) and any form that does not satisfy the condition are assembled exactly as before. The assembly call is also skipped when nothing is left to assemble (num_domains(dc) == 0, as jacobian_add! already does).

Because the decision is taken per form, the common cases work as follows:

Dirichlet data of the trial space	stiffness (k=0)	damping (k=1)	mass (k=2)
none, or zero	matvec	matvec	matvec
constant, non-zero (e.g. u = 1 on ∂Ω)	assembled	matvec (∂t u_D = 0)	matvec
time-dependent, affine in t (e.g. u_D = g(x)(1 + t))	assembled	assembled	matvec (∂tt u_D = 0)
time-dependent, general (e.g. ramp(t)·sin(ωt))	assembled	assembled	assembled
non-transient trial space (∂t spaces are homogeneous)	depends on U	matvec	matvec

(In each case the table gives the path actually taken in practice, as reported by _has_zero_dirichlet_values on the evaluated spaces.)

Nothing changes for non-constant forms, for quasilinear / semilinear / nonlinear operators, or for the Jacobian.

[codecov]

Codecov Report

❌ Patch coverage is 92.85714% with 1 line in your changes missing coverage. Please review.
✅ Project coverage is 89.58%. Comparing base (853eaa2) to head (2878ef6).
⚠️ Report is 3 commits behind head on master.

Files with missing lines	Patch %	Lines
src/ODEs/ODEOpsFromTFEOps.jl	92.85%	1 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##           master    #1345      +/-   ##
==========================================
+ Coverage   89.52%   89.58%   +0.05%     
==========================================
  Files         243      243              
  Lines       32142    32138       -4     
==========================================
+ Hits        28776    28791      +15     
+ Misses       3366     3347      -19     

Flag	Coverage Δ
drivers	37.16% <0.00%> (+<0.01%)	⬆️
extensions	6.14% <0.00%> (+<0.01%)	⬆️
unit-adaptivity	37.75% <0.00%> (+<0.01%)	⬆️
unit-basics	14.21% <0.00%> (+<0.01%)	⬆️
unit-celldata	20.31% <0.00%> (+<0.01%)	⬆️
unit-fespaces-1	35.37% <0.00%> (+<0.01%)	⬆️
unit-fespaces-2	39.57% <0.00%> (+<0.01%)	⬆️
unit-fields	18.87% <0.00%> (+0.05%)	⬆️
unit-geometry	26.68% <0.00%> (+<0.01%)	⬆️
unit-multifield	28.92% <0.00%> (+<0.01%)	⬆️
unit-odes	26.69% <92.85%> (+0.03%)	⬆️
unit-referencefes	38.51% <0.00%> (+<0.01%)	⬆️
unit-visualization	10.99% <0.00%> (+<0.01%)	⬆️

Flags with carried forward coverage won't be shown. Click here to find out more.

☔ View full report in Codecov by Harness.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[JordiManyer]

Hi @oriolcg thanks for the contribution! From a mathematical point of view, I am ok with the change. I left a comment on some aesthetic changes, please have a look. Thank you!

[JordiManyer]

Maybe we can use this function instead of !isnothing(const_form)? Also the axpy_values? I believe we are doing the same here than what you want to do above.

Then, can we create a helper like _has_zero_dirichlet_values(odeopcache.Us[k+1]) that instead of acting on the space itself acts on the operator like is_form_constant(odeop, k)? Under the hood we can do the same thing, but lets try to keep things clean.

[oriolcg]

I modified according to the suggestion. We cannot use axpy_entries as this is a mat-vec product. We still need to check the underlying FEspace, so we cannot just rely on the operator alone. I added a wrapper function is_form_constant_and_zero_dirichlet. See latest commit.

[JordiManyer]

Alright, all good! Thanks for the contribution @oriolcg !