Fix PML instability in Gmsh CPML MPI example

EXAMPLES/applications/meshing/Gmsh_example_CPML_MPI/create_geom_with_pygmsh.py

@@ -90,7 +90,7 @@ def add_pml_layers(self, geom, w_pml, lc_pml, top_pml=False):
 
         # add top tag if top_pml is False
         if not top_pml:
-            self.list_rects[irect_top].bound_tag = "Top"
+            self.list_rects[irect_top].bound_tag = "TopFree"
 
         # side of the main rectangles
         tmp_list_rects = self.list_rects.copy()
@@ -313,6 +313,11 @@ def build_points_edges(self, geom, debug=False):
             elif rect.bound_tag == "Top":
                 list_l_top.append(rect.list_lines[2])
                 list_l_top_inner.append(rect.list_lines[0])
+            elif rect.bound_tag == "TopFree":
+                # Top boundary of the domain when top_pml is False.
+                # Only add outer edge to Top group; do NOT create _Top inner
+                # boundary since there is no top PML layer.
+                list_l_top.append(rect.list_lines[2])
             elif rect.bound_tag == "Bottom":
                 list_l_bottom.append(rect.list_lines[0])
                 list_l_bottom_inner.append(rect.list_lines[2])

EXAMPLES/applications/meshing/Gmsh_example_CPML_MPI/create_mesh_topo.ipynb

@@ -422,7 +422,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 3,
+   "execution_count": null,
    "metadata": {},
    "outputs": [
     {
@@ -447,7 +447,7 @@
     "from meshio2spec2d import *\n",
     "\n",
     "mio2spec = Meshio2Specfem2D(mesh)\n",
-    "mio2spec.write(\"test\", pml_transition_layer=True)"
+    "mio2spec.write(\"test\", pml_transition_layer=False)"
    ]
   },
   {

EXAMPLES/applications/meshing/Gmsh_example_CPML_MPI/meshio2spec2d.py

@@ -208,6 +208,13 @@ def __init__(self, mesh, top_abs=False, bot_abs=True, left_abs=True, right_abs=T
         # check if PML_X is included in the physical groups
         if "PML_X" in self.mesh.cell_sets_dict:
             self.use_cpml = True
+            # When CPML is used, all outer PML boundaries need absorbing edges
+            # for Dirichlet boundary conditions (see pml_init.F90:
+            # determine_boundary_abs_points_PML). Auto-set all abs flags to True.
+            self.top_abs = True
+            self.bot_abs = True
+            self.left_abs = True
+            self.right_abs = True
 
         # check if second order elements are included in the mesh
         if "quad9" in self.mesh.cells_dict:
@@ -255,7 +262,14 @@ def write_material(self):
         arr_mflag = np.ones(self.n_cells, dtype=int) * -1
 
         # id offset for quad (subtract the number of lines)
-        self.cell_id_offset = int(np.min(self.mesh.cell_sets_dict["M1"][self.key_quad]))
+        # find the minimum cell id across all materials (not just M1,
+        # since M1 does not necessarily have the smallest cell id)
+        min_cell_id = int(np.min(self.mesh.cell_sets_dict[M_keys[0]][self.key_quad]))
+        for key in M_keys[1:]:
+            c = int(np.min(self.mesh.cell_sets_dict[key][self.key_quad]))
+            if c < min_cell_id:
+                min_cell_id = c
+        self.cell_id_offset = min_cell_id
 
         print("cell_id_offset: ", self.cell_id_offset)
 
@@ -407,7 +421,7 @@ def write_cpml(self):
         np.savetxt(self.fname_CPML, str_lines, fmt="%s")
 
 
-    def write(self, filename_out="TEST", pml_transition_layer=True):
+    def write(self, filename_out="TEST", pml_transition_layer=False):
 
         # measure time
         start_time = time.time()