Fix PML instability in Gmsh CPML MPI example#1274
Merged
danielpeter merged 2 commits intoMar 31, 2026
Merged
Conversation
Three bugs caused numerical blowup with PML boundary conditions when using external Gmsh meshes: 1. meshio2spec2d.py: pml_transition_layer default was True, which removed the innermost PML row via get_cpml_cells_except_damping(). For a 3-element PML this left only 2 layers - insufficient for absorption. Changed default to False. 2. meshio2spec2d.py: When CPML is detected (PML_X physical group exists), all outer PML boundaries now auto-set their absorbing flags to True. This ensures Dirichlet BCs are applied at outer PML edges as required by determine_boundary_abs_points_PML() in pml_init.F90. 3. create_geom_with_pygmsh.py: When top_pml=False, the top domain rectangle was tagged 'Top' which created a spurious '_Top' inner boundary. Introduced 'TopFree' tag that adds the outer edge to the Top physical group without creating an inner boundary. Updated create_mesh_topo.ipynb and regenerated MESH files accordingly.
Previously, the offset was hardcoded to the minimum cell ID of M1, which causes wrong material assignments when another material has a smaller minimum cell ID. Now the offset is the global minimum cell ID across all material keys (M1, M2, ...).
Member
|
thanks for this update :) |
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
2 participants
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Summary
Fixes numerical blowup (simulation instability) when using PML boundary conditions with external Gmsh meshes in the
Gmsh_example_CPML_MPIexample.Problem
The SPECFEM2D simulation blows up when using CPML with Gmsh-generated external meshes, particularly at the top boundary. The internal mesher works correctly for the same configuration.
Root Cause
Three bugs were identified in the Gmsh mesh converter (
meshio2spec2d.py) and geometry builder (create_geom_with_pygmsh.py):Bug 1:
pml_transition_layerdefault removes PML elementsmeshio2spec2d.py: Thewrite()method defaultedpml_transition_layer=True, which calledget_cpml_cells_except_damping()to remove the innermost PML row. For a 3-element PML (NELEM_PML_THICKNESS=3), this reduced the effective PML to 2 layers — insufficient for stable absorption, causing late-time blowup.Fix: Changed default to
pml_transition_layer=False.Bug 2: Missing absorbing boundary flags for CPML
meshio2spec2d.py: When CPML is used, all outer PML boundaries need absorbing edges for Dirichlet boundary conditions (seepml_init.F90: determine_boundary_abs_points_PML()). Thetop_absflag defaulted toFalse, meaning the top absorbing boundary was missing when CPML was detected.Fix: Auto-set all
*_absflags toTruewhen CPML physical groups are detected.Bug 3: Spurious inner boundary when
top_pml=Falsecreate_geom_with_pygmsh.py: Whentop_pml=False, the top domain rectangle was tagged"Top", which created a spurious_Topinner boundary (the bottom edge of the top domain rectangle). For multi-rectangle models, this incorrectly excluded side PML elements near the domain-top interface.Fix: Introduced
"TopFree"tag that adds the outer edge to theTopphysical group without creating an inner boundary.Files Changed
meshio2spec2d.py: Defaultpml_transition_layerchanged toFalse; auto-detect absorbing flags when CPML is usedcreate_geom_with_pygmsh.py: New"TopFree"boundary tag for non-PML top boundarycreate_mesh_topo.ipynb: Updatedwrite()call to usepml_transition_layer=FalseMESH/*: Regenerated mesh files with fixed codeDATA/STATIONS: Updated station fileTesting
The simulation was run with 4 MPI ranks for 10,000 timesteps (
DT=0.002, total time 20s) and completed stably. Final wavefield amplitudes remained bounded (elastic max norm: 4.22e-04, acoustic max: 188.0).