How to avoid a null solution for the Helmholtz PDE with Dirichlet BC?

General Discussion
1

Hi all! I am a beginner to freefem, please excuse any gaps in my knowledge. So as part of my research I am doing numerical simulations of the (vector) Helmholz equation in the time domain:

-\nabla^2 \mathbf{E} - k^2 \mathbf{E} = 0

The weak form that I obtained was:

\int_\Omega \nabla_J \mathbf{E} :\nabla_J \Phi + \int_\Omega k^2 \Phi \cdot \mathbf{E} - \int_{\partial \Omega} \Phi \nabla_J \mathbf{E} = 0

Where \nabla_J is the Jacobian. I wanted to verify that my code for my numerical simulations worked correctly by verifying against a simple case with a known analytical solution, that of the unit square. For simplicity I am using the scalar variant of the Helmholz equation -\nabla^2 E - k^2 E = 0, where for the vector-valued solution I would have \mathbf{E}_x = E, \mathbf{E}_y = E.

It is known that the Helmholz equation solved on a unit square with Dirichlet BC E \big|_{\partial \Omega} = 0 has the analytical solution:

E(x, y) = \sin(n \pi x)\sin(m \pi y)

Where k^2 = \pi^2(m^2 + n^2). I believe there is a uniqueness theorem that says that for Dirichlet boundary conditions E \big|_{\partial \Omega} = 0 over all the boundaries of the domain, if an eigenvalue is specified and satisfies the Helmholtz equation, then the solution is guaranteed to be both unique and non-trivial (nonzero). So I substituted this in, started the simulation, and got…nothing but zeroes. Here is my (simplified) code for reference:

int squaren = 80;
mesh Th = square(squaren, squaren);
plot(Th, value=1, wait=1);

/*
For values of k^2 = (n^2 + m^2) pi^2
then the homogenous Dirichlet problem
yields a solution even with
E = 0 at the boundaries.

For all other values of k, the Dirichlet
problem (E = 0 at boundaries) yields a
null solution, different boundary conditions
are necessary instead.
*/

int m = 1;
int n = 1;
real ksquared = (m^2 + n^2) * pi^2;
func AnalyticalMag = sin(pi*x) * sin(pi*y);

fespace Vh(Th, P1);
Vh TestPhi;
Vh Emag;
Vh ExactEmag=AnalyticalMag;

macro grad(u) [dx(u), dy(u)] //

solve ElectricField(Emag, TestPhi, solver=sparsesolver)
    = int2d(Th)(
        grad(Emag)' * grad(TestPhi)
    )
    + int2d(Th)(
        ksquared * Emag * TestPhi
    )
    + on(1, 2, 3, 4, Emag=0);

plot(Emag, value=true, fill=true, dim=3, cmm="Electric field (spatial components)", wait=1);
plot(ExactEmag, value=true, fill=true, dim=3, cmm="Electric field (analytical solution)", wait=1);

The best I could do was to set E \big|_{\partial \Omega} = \epsilon where \epsilon \approx 0 which FreeFEM returns a nontrivial but still incorrect solution. Even more confusingly, after using FreeFEM’s EigenValue() in a separate validation script, no matter what I tried, I couldn’t get the analytical values of the eigenvalue k^2.

At this point I am thinking that perhaps the elliptic nature of the PDE is to blame. Maybe I should try some sort of Newton method? Some advice on this problem would be greatly appreciated.

2

Hi,
You first say you consider the Helmholtz equation in the time domain (with a time derivative), but then write an eigenvalue equation. If you want eigenvalues/eigenvectors you should indeed call the EigenValue() function. In your script you solve a problem without a source that has a solution Emag=0.
Please refer to the tutorial example for acoustics: Acoustics. It is for a Neumann boundary condition, but just add your on() condition to varf op().

1 Like
3

Ah, that makes sense. And sorry, “time domain” was inaccurate, I was only interested in the spatial domain, I simply meant the Helmholtz equation resulting from separation of variables from the time-domain wave equation. That was my accidental mistaken explanation.

For future reference to other people here, this is the result:

  • Solving the Helmholz equation in general on Dirichlet BC E \big |_{\partial \Omega} will always result in a null solution (E = 0) unless with a value of k^2 that is an eigenvalue of the Laplacian operator
    • If you want to solve for these solutions (eigenfunctions), use EigenValue(), the docs are here.
    • This will give you your desired number of eigenfunctions on Dirichlet boundaries as well as their associated eigenvalues
    • Don’t just plug in a value of k because due to the homogeneity and linearity of the PDE, an eigenvalue solution is non-unique; any valid solution u + u_e where u_e is an eigenfunction will result in another valid solution
  • Otherwise, you need to solve with non-Dirichlet boundary conditions, i.e. you can’t have E \big|_{\partial \Omega} = 0 as your boundary condition.
4

Hello, in order to avoid the nul solution, you have to inject energy inside your 2D cavity.
For example, you can add a non nul source term (a non nul RHS in the weak form) exciting one arbitrary place inside your cavity (tip: with a region label different to zero).

Be careful: the Laplacian gives minus*gradient’gradient not plusgradient’*gradient
See below your code re-written.
Regards

int squaren = 80;

mesh Th = square(squaren, squaren);

Th = change( Th, fregion=((x-1./3)^2+(y-1./4)^2)<0.001 ); // to define the exciting place

plot(Th, value=1, wait=0);

/*

For values of k^2 = (n^2 + m^2) pi^2

then the homogenous Dirichlet problem

yields a solution even with

E = 0 at the boundaries.

For all other values of k, the Dirichlet

problem (E = 0 at boundaries) yields a

null solution, different boundary conditions

are necessary instead.

*/

int m = 1;

int n = 1;

real ksquared = (m^2 + n^2) * pi^2;

func AnalyticalMag = sin(pix) * sin(piy);

fespace Vh(Th, P1);

Vh TestPhi;

Vh Emag;

Vh ExactEmag=AnalyticalMag;

macro grad(u) [dx(u), dy(u)] //

solve ElectricField(Emag, TestPhi, solver=sparsesolver)

= -int2d(Th)( grad(Emag) '* grad(TestPhi) )

  • int2d(Th)( ksquared * Emag * TestPhi )

  • int2d(Th)( ksquared * region * TestPhi ) // non null RHS (source term where region!=0)

  • on(1, 2, 3, 4, Emag=0);

Emag = Emag/Emag.max;

plot(Emag, value=true, fill=true, dim=3, cmm=“Electric field (spatial components)”, wait=1);

plot(ExactEmag, value=true, fill=true, dim=3, cmm=“Electric field (analytical solution)”, wait=0);

5

Thank you for taking the time to write this. That being said, I have a few questions. First, does the sign convention matter? Second, I am specifically looking for free-space solutions ie. no source term. Indeed, adding a source term would fix the issues but that’s not what I am looking for.

6

Hi, in your case the problem becomes a eigenvalue(vector) one.
Please see the code below.
About the minus sign, it comes naturally from the transformation of the strong form to the weak (or variational) form with the integration by parts.
In the first version of your code, it was omited, so your solution was a vanish wave instead of a periodic sinusoidal one.
regards

int squaren = 80;
mesh Th = square(squaren, squaren);
plot(Th, value=1, wait=0);

/*
For values of k^2 = (n^2 + m^2) pi^2
then the homogenous Dirichlet problem
yields a solution even with
E = 0 at the boundaries.
For all other values of k, the Dirichlet
problem (E = 0 at boundaries) yields a
null solution, different boundary conditions
are necessary instead.
*/

int m = 1;
int n = 1;
real ksquared = (m^2 + n^2) * pi^2;
func AnalyticalMag = sin(pix) * sin(piy);

fespace Vh(Th, P1);
Vh TestPhi;
Vh Emag;
Vh ExactEmag=AnalyticalMag;

macro grad(u) [dx(u), dy(u)] //

// solve ElectricField(Emag, TestPhi, solver=sparsesolver) =
// -int2d(Th)( grad(Emag) '* grad(TestPhi) )
// + int2d(Th)( ksquared * Emag * TestPhi )
// + on(1, 2, 3, 4, Emag=0);

plot(ExactEmag, value=true, fill=true, dim=3, cmm=“Electric field (analytical solution)”, wait=0);

// Operator = A - sigma B ; // the shifted matrix (A=stiffness matrix, B=eigenvector)
varf op(Emag,TestPhi) = -int2d(Th)( grad(Emag)'grad(TestPhi) + ksquaredEmagTestPhi ) + on(1, 2, 3, 4, Emag=0);
varf b(Emag,TestPhi) = int2d(Th)(EmagTestPhi); // no Boundary condition

matrix OP = op(Vh,Vh,solver=Crout,factorize=1);
matrix B = b(Vh,Vh,solver=CG,eps=1e-20);

int nev=16; // number of requested eigenvalues near ksquared
real[int] tev(nev); // to store the nev eigenvalue
Vh[int] teVH(nev); // to store the nev eigenvector

int k=EigenValue(OP,B,sym=true,sigma=ksquared,value=tev,vector=teVH,tol=1e-10,maxit=0,ncv=0);

for(int i=0; i<min(nev,k); i++) {
teVH[i] = teVH[i]/teVH[i].max;
plot(teVH[i], value=true, fill=true, dim=3, cmm=“Electric field (spatial components), eigenvalue=”+tev[i], wait=0);
}