Skip to content

Add warning for non-positive semidefinite sigma tensor for Lorentzian susceptibility#3230

Open
oskooi wants to merge 2 commits into
NanoComp:masterfrom
oskooi:lorentzian_sus_stability_check
Open

Add warning for non-positive semidefinite sigma tensor for Lorentzian susceptibility#3230
oskooi wants to merge 2 commits into
NanoComp:masterfrom
oskooi:lorentzian_sus_stability_check

Conversation

@oskooi

@oskooi oskooi commented Jun 14, 2026

Copy link
Copy Markdown
Collaborator

Related to #3229.

This PR adds a warning if the user-supplied anisotropic Lorentzian susceptibility tensor (sigma) is non positive semidefinite. Additionally, the guard introduced in #666 to skip update of the polarization fields at material interfaces (to ensure stability) is reorganized to enable loop vectorization.

Comment thread src/structure.cpp
if (!sbb || !sab) continue;
LOOP_OVER_VOL(gv, ca, idx) {
const realnum a = saa[idx], b = sbb[idx], o = sab[idx];
if (a > 0 && b > 0 && o * o > a * b * (1 + 1e-6)) return true;

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why is it checking "2x2" tensors, and not the full 3x3 matrix?

Comment thread src/susceptibility.cpp
// boundaries (see PR #666). sigma is static in time, so the
// predicate is the same every step; apply it branchlessly so the
// loop still vectorizes. Where s[i] == 0, p[i] and pp[i] are left
// unchanged, exactly as the original guarded update did.

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It would be nice to know how much of a performance difference this makes in two cases:

  1. the chunk is completely filled with σ ≠ 0 (so mask == 1 everywhere).
  2. the chunk is mostly σ = 0 (so that you end up doing a bunch of extra computation that would have otherwise been removed by the if).

@stevengj

stevengj commented Jun 18, 2026

Copy link
Copy Markdown
Collaborator

My other concern about this PR is that it is looking at σ on the Yee grid, for which different components of σ are stored at different Yee grid points. If σ is not constant in space, then it is problematic to simply treat the different components as part of the same σ tensor … you would need to define some effective σ tensor by some kind of spatial average of σ components.

One alternative is simply to check positive-definiteness at a higher level, e.g. in the list of materials/objects in meepgeom.cpp.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants