Limiter update

buildMirge.sh

@@ -2,8 +2,6 @@
 
 # default branch for building mirgecom for this driver
 mirge_branch="production"
-#mirge_branch="limit-pocl"
-#mirge_branch="m-to-n-restart-update"
 # conda environment name
 conda_env="mirgeDriver.Y3prediction"
 

data/convert_nastran_to_gmsh.py

@@ -1,4 +1,5 @@
 import numpy as np
+import time
 
 def read_nastran_mesh(file_path):
     nodes = {}
@@ -80,7 +81,8 @@ def read_line(line, field_length=8):
             elif reading_elements:
                 if line.startswith("CQUAD") or\
                    line.startswith("CHEX") or\
-                   line.startswith("CROD"):
+                   line.startswith("CROD") or\
+                   line.startswith("CTRIA3"):
                     parts = read_line(line)
                     element_id = int(parts[1])
                     element_nodes = [int(parts[i]) for i in range(3, len(parts))]
@@ -158,6 +160,8 @@ def write_gmsh_mesh(file_path, nodes, elements, mesh_format,
         #for i, (element_nodes, element_tag) in enumerate(zip(elements, element_tags)):
             if len(element_nodes) == 2:
                    element_type = 1
+            elif len(element_nodes) == 3:
+                   element_type = 2
             elif len(element_nodes) == 4:
                    element_type = 3
             elif len(element_nodes) == 8:
@@ -172,13 +176,16 @@ def write_gmsh_mesh(file_path, nodes, elements, mesh_format,
         file.write("$EndElements\n")
 
 def convert_mesh(input_mesh, output_mesh):
-    # Example usage:
+    # Read the nastran file
+    begin_read_time = time.perf_counter()
     nodes, elements, element_tags, physical_names, physical_dim = read_nastran_mesh(input_mesh)
+    end_read_time = time.perf_counter()
+    print(f"Read mesh in {(end_read_time - begin_read_time):.1f} seconds.")
 
     # This is the mesh format meshmode expects
     mesh_format = "2.2 0 8"
 
-    print("Nodes:")
+    begin_minmax_time = time.perf_counter()
     first_node = 1e9
     last_node = -1
     for node_id, coordinates in nodes.items():
@@ -187,56 +194,72 @@ def convert_mesh(input_mesh, output_mesh):
         #print(f"Node {node_id}: {coordinates}")
     print(f"Read {len(nodes)} nodes, first node id {first_node}, last node id {last_node}")
 
-    print("\nElements:")
     first_element = 1e9
     last_element = -1
     for element_id, element_nodes in elements.items():
         first_element = min(first_element, element_id)
         last_element = max(last_element, element_id)
         #print(f"Element Nodes for Element {element_id}: {element_nodes}")
     print(f"Read {len(elements)} elements, first element id {first_element}, last element id {last_element}")
+    end_minmax_time = time.perf_counter()
+    print(f"Spent {(end_read_time - begin_read_time):.1f} seconds finding min/max node and element numbers.")
 
-    # renumber to start at 1
+    # renumber nodes to start at 1
+    begin_renumber_nodes_time = time.perf_counter()
     renumber_nodes = {}
     for new_id, old_id in enumerate(nodes.keys(), start=1):
         renumber_nodes[new_id] = nodes[old_id]
 
+    # update the element node numbers
     renumber_elements = {}
+    #renumber_elements = elements
     element_mapping = {}
-    for new_id, (element_id, node_ids) in enumerate(elements.items(), start=1):
-        renumber_elements[new_id] = [list(nodes.keys()).index(node_id) + 1 for node_id in node_ids]
-        element_mapping[element_id] = new_id
+
+    # Create a mapping for old node IDs to new IDs
+    old_to_new_node_mapping = {old_id: new_id for new_id, old_id in enumerate(nodes.keys(), start=1)}
+
+    # Renumbering nodes
+    for new_id, (old_id, coords) in enumerate(nodes.items(), start=1):
+        renumber_nodes[new_id] = coords  # Copy coordinates directly
+
+    # Renumbering elements
+    for new_id, (element_id, node_ids) in enumerate(list(elements.items()), start=1):
+        renumbered_node_ids = [old_to_new_node_mapping[node_id] for node_id in node_ids if node_id in old_to_new_node_mapping]
+        renumber_elements[new_id] = renumbered_node_ids  # Store the new node IDs
+        element_mapping[element_id] = new_id  # Map old element ID to new ID
 
     renumber_element_tags = {new_id: element_tags[element_id] for element_id, new_id in element_mapping.items()}
 
+    #renumber_element_tags = element_tags
+    end_renumber_nodes_time = time.perf_counter()
+
+    print(f"Spent {(end_renumber_nodes_time - begin_renumber_nodes_time):.1f} seconds renumbering nodes.")
+
     print("Renumber Nodes:")
     first_node = 1e9
     last_node = -1
     for node_id, coordinates in renumber_nodes.items():
         first_node = min(first_node, node_id)
         last_node = max(last_node, node_id)
         #print(f"Node {node_id}: {coordinates}")
-    print(f"Writing {len(nodes)} nodes, first node id {first_node}, last node id {last_node}")
+    print(f"Renumbered {len(nodes)} nodes, first node id {first_node}, last node id {last_node}")
 
-    print("\nRenumber Elements:")
     first_element = 1e9
     last_element = -1
     for element_id, element_nodes in renumber_elements.items():
         first_element = min(first_element, element_id)
         last_element = max(last_element, element_id)
         #print(f"Element Nodes for Element {element_id}: {element_nodes}")
-    print(f"Writing {len(elements)} elements, first element id {first_element}, last element id {last_element}")
-
-
-    #print("\nElement Tags:")
-    #for element_id, tags in element_tags.items():
-        #print(f"Element Tags for Element {element_id}: {tags}")
+    print(f"Renumbered {len(renumber_elements)} elements, first element id {first_element}, last element id {last_element}")
 
     print("\nPhysical Names:")
     for name_id, name in physical_names.items():
         print(f"Physical name {name_id}: {name}")
 
+    begin_write_time = time.perf_counter()
     write_gmsh_mesh(output_mesh, renumber_nodes, renumber_elements, mesh_format, physical_names, renumber_element_tags, physical_dim)
+    end_write_time = time.perf_counter()
+    print(f"Wrote new mesh in {(end_write_time - begin_write_time):.1f} seconds.")
 
 import sys
 def main():
@@ -247,7 +270,11 @@ def main():
     input_file = sys.argv[1]
     output_file = sys.argv[2]
 
+    total_time_begin = time.perf_counter()
     convert_mesh(input_mesh=input_file, output_mesh=output_file)
+    total_time_end = time.perf_counter()
+
+    print(f"Done converting mesh in {(total_time_end - total_time_begin):.1f} seconds.")
 
 if __name__ == "__main__":
     main()

smoke_test_ks/run_lazy.sh

@@ -1,4 +1,4 @@
 #!/bin/bash
 #mpirun -n 2 python -u -O -m mpi4py driver.py -i run_params.yaml --log --lazy
 # turn off optimizations for CI (enables asserts)
-mpirun -n 2 python -u -m mpi4py driver.py -i run_params.yaml --log --lazy
+mpirun -n 2 python -u -m mpi4py driver.py -i run_params.yaml --log --lazy --overintegration

smoke_test_ks/run_params.yaml

@@ -28,7 +28,7 @@ wall_material: 2
 # model control
 ##############
 nlimit: 1
-use_species_limiter: 1
+use_species_limiter: 2
 smooth_char_length_alpha: 0.025
 smooth_char_length: 5
 #smooth_char_length: 0

test/test_limiter.py

@@ -247,9 +247,10 @@ def test_positivity_preserving_limiter(actx_factory, order, dim,
     #
     from grudge.dof_desc import DD_VOLUME_ALL
     smin = 11
+    entropy_min = actx.np.zeros_like(fluid_cv.mass) + smin
     limited_cv = limit_fluid_state_lv(
         dcoll=dcoll, cv=fluid_cv, temperature_seed=tseed,
-        gas_model=gas_model, dd=DD_VOLUME_ALL, limiter_smin=smin)
+        gas_model=gas_model, dd=DD_VOLUME_ALL, entropy_min=entropy_min)
     limited_mass = limited_cv.mass
     print(f"{limited_mass=}")
     assert actx.to_numpy(actx.np.min(limited_mass)) >= 0.0
@@ -378,9 +379,10 @@ def test_positivity_preserving_limiter_multi(actx_factory, order, dim, nspecies,
     #
     from grudge.dof_desc import DD_VOLUME_ALL
     smin = 11
+    entropy_min = actx.np.zeros_like(fluid_cv.mass) + smin
     limited_cv = limit_fluid_state_lv(
         dcoll=dcoll, cv=fluid_cv, temperature_seed=tseed,
-        gas_model=gas_model, dd=DD_VOLUME_ALL, limiter_smin=smin)
+        gas_model=gas_model, dd=DD_VOLUME_ALL, entropy_min=entropy_min)
     limited_mass = limited_cv.mass
     print(f"{actx.to_numpy(limited_mass)=}")
     assert actx.to_numpy(actx.np.min(limited_mass)) >= 0.0

y3prediction/forward_step.py

@@ -0,0 +1,169 @@
+"""mirgecom driver initializer for 1d shock."""
+from functools import partial
+
+
+def get_mesh(dim, size, bl_ratio, interface_ratio,
+             transfinite=False, use_wall=True, use_quads=False,
+             use_gmsh=True):
+    """Generate a grid using `gmsh`."""
+
+    if use_gmsh:
+        from meshmode.mesh.io import (
+            generate_gmsh,
+            ScriptSource
+        )
+
+        # for 2D, the line segments/surfaces need to be specified clockwise to
+        # get the correct facing (right-handed) surface normals
+        my_string = (f"""
+            Point(1) = {{ 0., 0., 0., {size}}};
+            Point(2) = {{ 0.6, 0.0, 0., {size}}};
+            Point(3) = {{ 0.6, 0.2, 0., {size}}};
+            Point(4) = {{ 3.0, 0.2, 0., {size}}};
+            Point(5) = {{ 3.0, 1.0, 0., {size}}};
+            Point(6) = {{ 0.6, 1.0, 0., {size}}};
+            Point(7) = {{ 0.0, 1.0, 0., {size}}};
+            Point(8) = {{ 0.0, 0.2, 0., {size}}};
+
+            Point(9) = {{ 0.45, 0.0, 0., {size}}};
+            Point(10) = {{ 0.45, 0.2, 0., {size}}};
+            Point(11) = {{ 0.45, 1.0, 0., {size}}};
+            Point(12) = {{ 0.0, 0.35, 0., {size}}};
+            Point(13) = {{ 0.6, 0.35, 0., {size}}};
+            Point(14) = {{ 3.0, 0.35, 0., {size}}};
+            Point(15) = {{ 0.45, 0.35, 0., {size}}};
+        """)
+
+        my_string += ("""
+            Line(1) = {1, 9};
+            Line(2) = {9, 10};
+            Line(3) = {10, 8};
+            Line(4) = {8, 1};
+            Line(6) = {9, 2};
+            Line(7) = {2, 3};
+            Line(8) = {3, 10};
+            Line(9) = {10, 15};
+            Line(10) = {15, 12};
+            Line(11) = {12, 8};
+            Line(12) = {3, 13};
+            Line(13) = {13, 15};
+            Line(14) = {15, 11};
+            Line(15) = {11, 7};
+            Line(16) = {7, 12};
+            Line(17) = {13, 6};
+            Line(18) = {6, 11};
+            Line(19) = {3, 4};
+            Line(20) = {4, 14};
+            Line(21) = {14, 13};
+            Line(22) = {14, 5};
+            Line(23) = {5, 6};
+        """)
+
+        my_string += ("""
+            Line Loop(1) = {1, 2, 3, 4};
+            Line Loop(2) = {6, 7, 8, -2};
+            Line Loop(3) = {3, -11, -10, -9};
+            Line Loop(4) = {8, 9, -13, -12};
+            Line Loop(5) = {10, -16, -15, -14};
+            Line Loop(6) = {13, 14, -18, -17};
+            Line Loop(7) = {19, 20, 21, -12};
+            Line Loop(8) = {21, 17, -23, -22};
+        """)
+
+        my_string += ("""
+            Surface(1) = {1};
+            Surface(2) = {2};
+            Surface(3) = {3};
+            Surface(4) = {4};
+            Surface(5) = {5};
+            Surface(6) = {6};
+            Surface(7) = {7};
+            Surface(8) = {8};
+
+        """)
+
+        my_string += ("""
+            Physical Surface('fluid') = {1:8};
+            Physical Curve('inflow') = {4, 11, 16};
+            Physical Curve('outflow') = {20, 22};
+            Physical Curve('flow') = {4, 11, 16, 20, 22};
+            Physical Curve('slip_wall') = {1, 6, 15, 18, 23};
+            Physical Curve('noslip_wall') = {7, 19};
+        """)
+
+        if transfinite:
+            my_string += (f"""
+                Transfinite Curve {{1, 3, 10, 15}} = {0.45} / {size} + 1;
+                Transfinite Curve {{16, 14, 17, 22}} = {0.65} / {size} + 1;
+                Transfinite Curve {{11, 9, 12, 20}} = {0.15} / {size} + 1;
+                Transfinite Curve {{2, 4, 7}} = {0.2} / {size} + 1;
+                Transfinite Curve {{6, 8, 13, 18}} = {0.15} / {size} + 1;
+                Transfinite Curve {{19, 21, 23}} = {2.4} / {size} + 1;
+            """)
+            my_string += (
+                "Transfinite Surface {1:8} Right;"
+            )
+
+        else:
+            my_string += (f"""
+                // Create distance field from curves, excludes cavity
+                Field[1] = Distance;
+                Field[1].CurvesList = {{7, 19}};
+                Field[1].NumPointsPerCurve = 100000;
+
+                //Create threshold field that varies element size near boundaries
+                Field[2] = Threshold;
+                Field[2].InField = 1;
+                Field[2].SizeMin = {size} / {bl_ratio};
+                Field[2].SizeMax = {size};
+                Field[2].DistMin = 0.0002;
+                Field[2].DistMax = 0.005;
+                Field[2].StopAtDistMax = 1;
+
+                //  background mesh size
+                Field[3] = Box;
+                Field[3].XMin = 0.;
+                Field[3].XMax = 3.0;
+                Field[3].YMin = -1.0;
+                Field[3].YMax = 1.0;
+                Field[3].VIn = {size};
+
+                // Create distance field from curves, excludes cavity
+                Field[4] = Distance;
+                Field[4].CurvesList = {{2}};
+                Field[4].NumPointsPerCurve = 100000;
+
+                // take the minimum of all defined meshing fields
+                Field[100] = Min;
+                Field[100].FieldsList = {{2, 3, 5}};
+                Background Field = 100;
+            """)
+
+        my_string += ("""
+            Mesh.MeshSizeExtendFromBoundary = 0;
+            Mesh.MeshSizeFromPoints = 0;
+            Mesh.MeshSizeFromCurvature = 0;
+
+            Mesh.Algorithm = 5;
+            Mesh.OptimizeNetgen = 1;
+            Mesh.Smoothing = 100;
+        """)
+
+        if use_quads:
+            my_string += ("""
+            // Convert the triangles back to quads
+            Mesh.Algorithm = 6;
+            Mesh.Algorithm3D = 1;
+            Mesh.RecombinationAlgorithm = 2;
+            Mesh.RecombineAll = 1;
+            Recombine Surface {1, 2, 3};
+            """)
+
+        print(my_string)
+        mesh_construction_kwargs = {
+            "force_positive_orientation":  True,
+            "skip_element_orientation_test":  True}
+        return partial(generate_gmsh, ScriptSource(my_string, "geo"),
+                       force_ambient_dim=dim, dimensions=dim, target_unit="M",
+                       mesh_construction_kwargs=mesh_construction_kwargs,
+                       return_tag_to_elements_map=True)