Assemble_matrix on frequency dependent form

General
1

Hello guys,
I am currently solving frequency dependent Helmholtz equation in this way:

  1. After setting up mesh, defining frequency and wavenumber as Constant(msh, PETSc.ScalarType()), I write the weak form and define the linear problem as follows
#... initial settings ...
# ...
omega = Constant(msh, PETSc.ScalarType(1))
k0 = Constant(msh, PETSc.ScalarType(1))
#...

V = FunctionSpace(msh, ("CG", 2)

# ...other settings

f  = Function(V)
# other stuff about f
a  = inner(grad(u), grad(v)) * dx - k0**2 * inner(u, v) * dx 
L  = inner(f, v) * dx

uh      = Function(V)
uh.name = "u"
problem = LinearProblem(a, L, u=uh, petsc_options={"ksp_type": "preonly", "pc_type": "lu","pc_factor_mat_solver_type": "mumps"}
  1. I update the value of k0 and omega inside a for loop and solve the problem for every frequency:
for freq in frequency_range:

    # Compute solution
    omega.value = 2*np.pi*freq
    k0.value = 2*np.pi*freq*c0
    problem.solve()

Now I need the LHS to change multiple times for some reason, so I need to assemble the system without using LinearProblem().
Question: Does the k0.value = something have effect if the matrix has been assembled with

A = assemble_matrix(form(a), [])
A.assemble()

?

If not, do I need to bring out of the matrix all the frequency dependent terms and build something similar to:

[{K_{a}}] + \omega^{2}[M_a]

where [{K_{a}}] and [M_a] are assembled like this:

K_a = assemble_matrix(form(inner(grad(u), grad(v)) * dx), [])
M_a = assemble_matrix(form(1/c0**2 * inner(u, v) * dx ), [])

K_a.assemble()
M_a.assemble()

thank you very much, as always.

2

Regarding your use of dolfinx.fem.form see Speeding up time-dependent diffusion equation - #7.

If you make a change to underlying Functions or Constants in the bilinear form, you’d need to update the matrices K_a and M_a by reassembling them. Something along the lines of:

# Change some constants and functions
...
# Reassemble
K_a = assemble_matrix(K_a_form, [])
K_a.assemble()
M_a = assemble_matrix(M_a_form, [])
M_a.assemble()
3

Yes, I understand the issue. This means that the optimal approach is to make the matrices (and vectors) not frequency dependent and assemble them outside of the for loop. This will also move dolfinx.fem.form outside of the loop.

then I would build the matrix A as linear combination of the assembled matrix as you say here:

that in my case means building the matrix A in this way:

K_a = assemble_matrix(K_a_form), [])
K_a.assemble()
M_a = assemble_matrix(M_a_form), [])
M_a.assemble()

# creating KSP solver
...
# frequency loop
for freq in frequency_range:

    omega.value = 2*np.pi*freq # this is the only frequency dependent value
    
    A = K_a - omega**2 * M_a

    # solving the problem  A * x = b

is it a correct approach?

4

Looks reasonable. After making sure your approach works, consider further optimisation by predeclaring the linear operator A as I believe the line

    A = K_a - omega**2 * M_a

will create a copy every loop.

I’m writing pseudocode here, something like the following which exploits MatAXPY — PETSc v3.19.4-1061-gd1b98e1ad7 documentation

A = K_a.copy()  # Needs to be the right shape and sparsity pattern
#
...
#
ksp.setOperator(A)
#
...
#
for freq in frequency_range:
  A.zeroEntries()
  A.axpy(1.0, K_a)
  A.axpy(-omega**2, M_a)
  A.assemble()
5

Thank you, I din’t know about the function axpy(), it’s very useful.

Another (hopefully the last) question:
Do I steel need to call the assemble() function, even though A is a linear combination of assembled matrices?