# Assemble\_matrix on frequency dependent form

**URL:** <https://fenicsproject.discourse.group/t/assemble-matrix-on-frequency-dependent-form/12078>\
**Category:** General\
**Created:** [August 30, 2023, 3:14pm UTC](https://fenicsproject.discourse.group/t/assemble-matrix-on-frequency-dependent-form/12078 "2023-08-30T15:14:12Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![bay\_swiss](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/bay_swiss/32/8778_2.png) [@bay\_swiss](https://fenicsproject.discourse.group/u/bay_swiss)\
**Post date:** [August 30, 2023, 3:14pm UTC](https://fenicsproject.discourse.group/t/assemble-matrix-on-frequency-dependent-form/12078/1 "2023-08-30T15:14:12Z")

</div>

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

```auto
#... 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:

```auto
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

```auto
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:

```auto
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.

---

<div class="post-metadata">

**Author:** ![nate](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/nate/32/17_2.png) [@nate](https://fenicsproject.discourse.group/u/nate)\
**Post date:** [August 30, 2023, 3:23pm UTC](https://fenicsproject.discourse.group/t/assemble-matrix-on-frequency-dependent-form/12078/2 "2023-08-30T15:23:14Z")

</div>

Regarding your use of `dolfinx.fem.form` see [Speeding up time-dependent diffusion equation - #7](https://fenicsproject.discourse.group/t/speeding-up-time-dependent-diffusion-equation/11856/7).

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

```python
# 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()

```

---

<div class="post-metadata">

**Author:** ![bay\_swiss](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/bay_swiss/32/8778_2.png) [@bay\_swiss](https://fenicsproject.discourse.group/u/bay_swiss)\
**Post date:** [August 30, 2023, 3:54pm UTC](https://fenicsproject.discourse.group/t/assemble-matrix-on-frequency-dependent-form/12078/3 "2023-08-30T15:54:25Z")

</div>

> [@nate](#):
>
> Regarding your use of `dolfinx.fem.form` see [Speeding up time-dependent diffusion equation - #7](https://fenicsproject.discourse.group/t/speeding-up-time-dependent-diffusion-equation/11856/7).

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:

> [@Speeding up time-dependent diffusion equation](https://fenicsproject.discourse.group/t/speeding-up-time-dependent-diffusion-equation/11856/6):
>
> Trying to find more optimisations will likely involve splitting up your stiffness matrix into linear combinations of those with multiples of updating constants.

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

```auto
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?

---

<div class="post-metadata">

**Author:** ![nate](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/nate/32/17_2.png) [@nate](https://fenicsproject.discourse.group/u/nate)\
**Post date:** [August 30, 2023, 4:11pm UTC](https://fenicsproject.discourse.group/t/assemble-matrix-on-frequency-dependent-form/12078/4 "2023-08-30T16:11:11Z")

</div>

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

```auto
    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](https://petsc.org/main/manualpages/Mat/MatAXPY/)

```auto
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()

```

---

<div class="post-metadata">

**Author:** ![bay\_swiss](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/bay_swiss/32/8778_2.png) [@bay\_swiss](https://fenicsproject.discourse.group/u/bay_swiss)\
**Post date:** [August 30, 2023, 4:24pm UTC](https://fenicsproject.discourse.group/t/assemble-matrix-on-frequency-dependent-form/12078/5 "2023-08-30T16:24:06Z")

</div>

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?