How to assemble 'matrix' only once with dolfinx 0.7.3?

Hello everyone,

Dolfinx version is 0.7.3 and have been installed with conda.
Other version I need in the minimal code (at the bottom of the post) are :

  • numpy : 1.26.4
  • mpi4py : 3.1.6
  • ufl : 2.23.2.0
  • pyvista : 0.44.2

For the context, I’m solving an acoustic equation with Neumann or Robin boundary condition (but not dirichlet).

I solve it in the laplace domain according to time (like some people do in the harmonic domain).
Consequently, I must solve one problem for each Laplace parameter s.

A weak formulation gives something like :

∀s ,  a(u,v) + s²b(u,v)  = L(v) 

which in term of matrix gives :

∀s , (K + s²M).U = L

Then, for each s-value, the problem is solved using an iteratived domain decomposition method.
This implie to assemble the FE “matrix” too many times.

Consequently, i would like to “save/compute once” the Finite Element objects.
For one s-value (for one frequency), i have tried to follow the Solving a time-dependent problem — FEniCS 22 tutorial tutorial.

Q1) Problem of version to assemble the second member

Due to the version difference,

# Create the problem  A_form ~ (K+s²M)
problem = LinearProblem(A_form, L_form, bcs = [ ])

# Solve the problem 
uh = problem.solve()

# Create a function
func = fem.Function(V)

I can’t use the ‘.x.petsc_vec’ method for a solution of a problem (uh) or on a function (func).
I can assemble the bilinear form once as in the tutorial

  • How can I code something like in the tutorial for the second member ?

I do not plan to change version now since we are in code validation step.

Q2) My actual best choice to optimize my code

If i define the dolfinx form like :

# Bilinear form 
a_form = inner(nabla_grad(u),nabla_grad(v))*dx
b_form = inner(u,v)*dx

# Linear form
L_form = inner(f,v)*dx

# Or just  general A_form of the problem
A_form = inner(nabla_grad(u),nabla_grad(v))*dx + s**2 * inner(u,v)*dx

To optimize my code, I would like to :

  1. Assemble the two form once and for all cause these two objects do not depend on the frequency/s-value and i dont have dirichlet bc.
  • af = assemble(form(a_form)),
  • bf = assemble(form(b_form))
  1. Sum them at each frequency/s-value such that :

    af + s²bf
    
  • Is it possible to do something like that ?

Here is the minimal code I’m using to achieve my goal :

## Library

import numpy as np
from dolfinx import fem, la, mesh, plot
from mpi4py import MPI
import ufl
from dolfinx.fem.petsc import LinearProblem
import pyvista as pv

## Parameters

s = 3       # Laplace parameter

## Domain parameter

length, height = 10, 3
Nx, Ny = 80, 60
extent = [[0.0, 0.0], [length, height]]
domain = mesh.create_rectangle(MPI.COMM_WORLD, extent, [Nx, Ny], mesh.CellType.quadrilateral)

## Finite Element Space/functions

V = fem.FunctionSpace(domain,("P",2,(3,)))
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)

## Define a manufactured solution 

def Solution(x):
    X = np.zeros_like(x)
    X[0] = x[0]
    X[1] = np.cos(x[0]+x[2]**2)
    X[2] = 0.0
    return X

sol = fem.Function(V)
sol.interpolate(Solution)

## Solve the problem directly

# Linear form L(v)
L_form = (ufl.inner(ufl.nabla_grad(sol),ufl.nabla_grad(v)) + s**2 *ufl.inner(sol,v))* ufl.dx 

# Bilinear form A(u,v)
A_form = (ufl.inner(ufl.nabla_grad(u),ufl.nabla_grad(v)) + s**2 * ufl.inner(u,v))* ufl.dx 

# Create the problem
problem = LinearProblem(A_form, L_form, bcs = [ ])

# Solve the problem 
uh = problem.solve()


## Avoid multiple assembly

    ## Assemble the bilinear form

compiled_a = fem.form(A_form)
A = fem.petsc.assemble_matrix(compiled_a,bcs=[])
A.assemble()

    ## Assemble RHS : the linear form

compiled_L = fem.form(L_form)
L = fem.petsc.assemble_vector(compiled_L)


Thanks in advance.

I guess you are close. You already know how to use fem.petsc.assemble_matrix to assemble a matrix. Now you just need to call it twice, once on the form with af and another time on the form with bf. The resulting type is petsc4py.PETSc.Mat, and you can look into their documentation on how to sum them (the method will usually be called axpy). Finally, call the petsc4py.PETSc.KSP solver yourself rather than using LinearProblem (there are several demos in dolfinx that use that, search for KSP in the demos of your version)

Dear francesco-ballarin,

Thank you for your fast reply.

My phd supervisor recently started use FEniCSx (version 9.3) and has followed the time-dependant tutorial to save / assemble matrix once.

As a conclusion, saving Finite Element matrix may save 0.1-0.2 sec in my case which is not worth it.

Thanks again.