Hi @Hugo,
Thanks for sharing your code and the detailed explanation of your approach. I see you're trying to increase the accuracy of your mode solver by reducing the simulation space through symmetry, which is a smart move. Let's go over your comments and figure out the issue.
1. Boundary Condition Setup: Dirichlet and Neumann on borders
You mentioned that you're applying Dirichlet boundary conditions on the top and right sides, and Neumann on the bottom and left sides. This is a valid approach for reducing the domain using symmetry.
However, I noticed something in your code that could be causing issues. It looks like you’re adding `1/d^2` to the diagonal terms for both the right and bottom edges, which isn’t how * Neumann boundary conditions* are typically applied in finite difference methods. Instead of modifying the diagonal, we need to modify the off-diagonal terms (the neighboring values) to mirror the values across the boundary.
For example, at the bottom edge (y = 0), a Neumann condition means that the value at the bottom should be the same as the value at the point just above it. So, we don’t touch the diagonal; we adjust the terms that represent neighbors.
2. Neumann Boundary Condition Misapplication
You wrote that you're adjusting the diagonal at the boundaries by adding `1/d^2`. The issue here is that Neumann conditions (i.e., zero derivative at the boundary) should be applied by doubling the off-diagonal terms that correspond to the boundary's neighbors, not modifying the diagonal term itself.
Here’s how the correction should look:
Instead of this:
if(ii>=l*(m-1)) % Right side
Phi(ii,ii) = Phi(ii,ii) + 1/d^2;
end
if(mod(ii,m)==0) % Bottom side
Phi(ii,ii) = Phi(ii,ii) + 1/d^2;
end
You should be doing this:
if(mod(ii, m) == 0) % Bottom edge
if ii + m <= l * m
Phi(ii + m, ii) = 2 / d^2; % mirror in y-direction
end
end
if(ii >= l * (m - 1)) % Right edge
if ii - 1 >= 1
Phi(ii - 1, ii) = 2 / d^2; % mirror in x-direction
end
end
This will correctly implement the Neumann boundary conditions and give you the right behavior for your mode solver near the edges.
3. Why the Full Field Results Differ
You mentioned that your reduced domain simulation doesn't match the full-field solution, and it's not just an *accuracy issue*. That’s because applying the Neumann condition incorrectly near the boundaries — by adjusting the diagonal instead of the off-diagonal neighbors — introduces errors in the way the field is being calculated. Once you correct this, you should get results that match those from the full-field simulation more closely.
4. Corner Treatment
For the bottom-right corner, where Neumann conditions from both axes meet, you’ll need to mirror the neighbors in both x and y directions. This means adjusting the off-diagonal elements in both directions, without touching the diagonal.
This is how we ensure that the boundary condition is correctly enforced in the corner, maintaining the zero derivative across both boundaries.
5. Additional Thoughts on the Eigenvalue Formulation
@Torsten also raised an insightful theoretical point: when you're solving an eigenvalue problem like this, you need to make sure the boundary conditions are applied correctly to fully define the problem. If the reduced domain doesn’t have sufficient boundary conditions, the eigenvalue solver could behave unpredictably, as there might be too many unknowns or conflicting constraints.
Once the boundary conditions are correctly implemented, you should see much more stable results for the eigenvalue problem.
So, once you update your boundary conditions as I’ve outlined (modifying the off-diagonal terms instead of the diagonal for Neumann BCs), your reduced-domain simulation should align better with the full-field solution. This change will also resolve the physics issue introduced by incorrectly applied Neumann BCs.
Let me know if you'd like further clarification from us or need help with testing the change.
