Change linear form after its computation

Question

Hello,

In a 2D problem I want to solve I have a right hand side of the form
$$L = \nabla \cdot \left( A \nabla f\right) + \omega^2 f$$
with $A = A(x,y)$ a $2\times 2$ matrix and $f = \exp{(i\omega x)}$.

When $A = Id(2)$ we clearly have $L=0$, but numerically computing $L$ even with third degree elements produces a solution of order $10^{-2}$ instead of being exactly 0.

There is however a whole region where $A$ is not the identity but assumes values depending on the position.

I tried feeding $L$ a modified function $f$ that is exactly 0 except for the region where $A \neq Id(2)$, but this produces a jump at the interface between the regions.

So I was wondering if there is a way to change the linear form $L$ after it's been computed to set it exactly to 0 in the region where it should in fact be null.

Thanks in advance

Answer 1

1e-2 or 0 are both <<1 ? Does the error not go down with mesh refinement? Why do you think this is the dominant source of errors in your formulation?

Comments

The error does indeed go down with mesh refinement, problem is my machine is only powerful enough to go up to a smaller and smaller value of mesh precision as the degree of the elements augments.

Plus, the computation with higher degrees requires time, and it seems meaningless to me to enforce such hard computation when I could just change the values of $L$.

I don't know if this trick will give me better results, but I hope it does.

Answer 2

Instead of using an analytic expression for $f$ and then having to use very accurate quadrature rules, why not computing directly $f_h$ as the solution of the following eigenvalue problem?

Find $f_h$ such that
$$ ( \nabla f_h, \nabla v_h) = \lambda (f_h, v_h) \forall v_h \in V_h$$

where you set $\lambda = \omega^2$.

You can then use the dolfin interface to SLEPC to solve the eigenproblem, and let $f_h$ be the eigenvector relative to $\lambda = \omega^2$.

Best,

Umbe

Comments

You should then define the Linear form $L: V_h \mapsto R $ as

$$ L(v_h) = - ( A \nabla f_h, \nabla v_h) + \omega^2 (f_h, v_h) $$

---

Hi Umberto,

Thank you for you answer.

If I got it correctly, you are suggesting to get $f_h$ directly on the grid via the solution of another eigenvalue problem, instead of defining it analitically and intepolating it on the grid, is that right?

I'll give it a try.

---

After all, the problem was not in the linear form. I have tried to implement your suggestion, but I think it was adding too much unnecessary complication.

If the problem had really been my RHS function maybe it'd have been more useful.

Thanks anyway!