Feature request: Assembly of high-order polynomial nonlinear terms (e.g., trilinear forms) for reduced-order modelling

[tiagomrns]

The feature request is related to the following problem
My team is working on reduced-order modelling (ROM) techniques for nonlinear problems, similar to what is done in the MORFEproject/MORFE2.0. In such methods, one often needs to precompute tensors that represent high-order polynomial nonlinearities. For example, given a weak form that contains a bilinear term like ∫ (u·∇)u · v dΩ, discretisation yields a trilinear form C(U, U, V) which can be represented as a third-order tensor H such that H(U, U, V) = H_{ijk} U_j U_k V_i. More generally, one may need terms of the form G(U1, U2) (bilinear) or H(U1, U2, U3) (trilinear) where each argument comes from a different finite element space (or the same).

Currently, Gridap provides excellent tools for assembling matrices from bilinear forms via assemble_matrix (which handles two FESpaces). However, for forms involving more than two FESpaces, there is no direct way to assemble the corresponding multi-dimensional arrays (tensors). This would be extremely valuable for ROM workflows, where one precomputes these tensors to enable efficient online evaluations.

Describing a possible solution
We propose adding functions similar to assemble_matrix but capable of handling an arbitrary number of FESpace arguments. For instance:

• assemble_tensor(a, trial_fe1, trial_fe2, test_fe) for a trilinear form a(u1, u2, v) returning a 3‑D array (or a sparse representation thereof).
• assemble_tensor(a, trial_fe1, trial_fe2, trial_fe3, test_fe) for quadrilinear forms, etc.

These functions would iterate over the degrees of freedom of each space and evaluate the form for combinations of basis functions, storing the result in a multi-dimensional container. For efficiency, sparsity patterns could be derived from the underlying meshes and finite element spaces, similar to how assemble_matrix uses a SparsityPattern.

Optionally, the user could also specify a test_fe and multiple trial_fe spaces in a flexible way (e.g., using a tuple of spaces). The output could be a dense or sparse tensor, depending on the problem size and sparsity.

Describing alternatives we have considered
One could manually loop over degrees of freedom and call evaluate on the weak form for each combination, but this is cumbersome, error‑prone, and does not leverage Gridap’s efficient cell‑wise assembly and sparsity handling. Another alternative is to use the existing assemble_matrix repeatedly by freezing some arguments, but that would be extremely inefficient for high‑order tensors.

Additional context
In reduced‑order modelling, after projecting the full‑order model onto a reduced basis, the nonlinear terms become low‑dimensional polynomial expressions. For a quadratic nonlinearity, one precomputes a third‑order tensor H (size N_r × N × N, where N_r is the number of test basis functions and N the number of trial basis functions). This tensor is then used online to compute the reduced nonlinear term as H_{ijk} q_j q_k. Similarly, for a cubic nonlinearity, a fourth‑order tensor is needed.

The current infrastructure in Gridap (e.g., FESpace, SparseMatrixAssemblers.jl) already contains many of the building blocks: cell iterators, local assembly loops, and sparsity pattern construction. Extending this to higher‑order tensors would likely involve:

1. Creating a new assembler type that can handle multiple trial spaces.
2. Implementing a routine to compute a sparsity pattern for the tensor (based on the support of basis functions from each space).
3. Looping over cells and local degrees of freedom to fill the global tensor entries.

We believe this feature would greatly enhance Gridap’s capabilities for nonlinear and reduced‑order modelling, and we’d be happy to help with design discussions or even contribute if given some guidance.

Thank you for considering this request!

[JordiManyer]

Hi @tiagomrns , thank you for using the library!

Let me first say I know very little about ROMS, so take all I say with a grain of salt. I would say that things are a bit more complicated than you think. Under the hood, I am quite confident we are assuming (at least in some places) that we have a Nquad x Ntest x Ntrial array being computed in each cell (for a bilinear form). So this would need to be expanded somehow, which would (as a minimum) entail recasting the bases into the appropriate shapes.

In any case, we develop things as we need them (we are researchers after all), which means we might never develop this ourselves. We are happy to accept contributions, though, so feel free to drop a PR if you manage to make it work. I would suggest looking at GridapROMs as well as the publication attached to see if something there might help you. Maybe @nichomueller can be of more help.

Cheers!

[nichomueller]

Hi @tiagomrns

As @JordiManyer just mentioned, the GridapROMs library implements quite a few ROM examples you might find useful, including for steady and unsteady Navier-Stokes equations. If you have questions in regards to the ROM side of things, I can help you; though I'm afraid I don't have an answer for your specific doubt. At a glance, I would say you could implement a higher-dimensional structure similar to an AbstractParamArray, which I implemented in GridapROMs. I refer to this paper for more details, but it's basically an interface for arrays related to parameterised PDEs (a parametric dimension is added to the vectors/matrices one would obtain by integrating a non-parameterised problem in Gridap).

Cheers :)

[tiagomrns]

Thank you for the quick responses, @JordiManyer and @nichomueller :)

The idea is to extend the existing assembly machinery to handle an arbitrary number of FESpace arguments.
We propose the following design:

1. New assembler type: SparseTensorAssembler{N} storing a tensor builder, test space, and N trial spaces.
2. Tensor storage: A sparse representation of an (N+1)-dimensional array (e.g., a dictionary of keys, or a custom SparseTensor type with coordinates and values).
3. Local evaluation: The multilinear form is evaluated cell‑wise on all combinations of test and trial basis functions, producing an array of size (n_test, n_trial1, …, n_trialN).
4. Assembly loops: Two‑pass symbolic/numeric loops similar to assemble_matrix, but with multiple index sets.

Code sketch

struct SparseTensorAssembler{N,T} <: Assembler
  tensor_builder   # builds a sparse tensor of order N+1
  test_space::FESpace
  trial_spaces::NTuple{N,FESpace}
  strategy::AssemblyStrategy
end

function SparseTensorAssembler(trial_spaces::FESpace..., test_space::FESpace)
  T = get_dof_value_type(test_space)
  builder = SparseTensorBuilder(T, test_space, trial_spaces...)  # to be defined
  SparseTensorAssembler(builder, test_space, trial_spaces, DefaultAssemblyStrategy())
end

function assemble_tensor(a::SparseTensorAssembler, tensordata)
  t1 = nz_counter(a.tensor_builder, (get_test_dofs(a), get_trial_dofs(a)...))
  symbolic_loop_tensor!(t1, a, tensordata)
  t2 = nz_allocation(t1)
  numeric_loop_tensor!(t2, a, tensordata)
  t3 = create_from_nz(t2)
  t3
end

function symbolic_loop_tensor!(T, a::SparseTensorAssembler, tensordata)
  strategy = get_assembly_strategy(a)
  for (cell_tensor, cell_test_dofs, cell_trial_dofs...) in zip(tensordata...)
    test_dofs = map_cell_rows(strategy, cell_test_dofs)
    trial_dofs_per_space = map(c -> map_cell_cols(strategy, c), cell_trial_dofs)
    # Use caching and touch! to mark non‑zero entries
  end
  T
end

function numeric_loop_tensor!(T, a::SparseTensorAssembler, tensordata)
  # Similar, using AddEntriesMap for tensor
end

The SparseTensorBuilder would implement nz_counter, nz_allocation, create_from_nz, and support entry‑wise operations like TouchEntriesMap and AddEntriesMap. The weak‑form API would need to allow multilinear forms with multiple trial spaces.

This extension would enable assembling tensors for terms like H(U1, U2, V) directly.