Update conversion logic to fill equilibrium contour_tree

- Handle multiple x-points in boundary_[secondary_]separatrix
- Determine if O-point at magnetic axis is a local minimum or maximum of the psi map

Author: maarten-ic

imas/ids_convert.py

@@ -1098,26 +1098,36 @@ def _equilibrium_boundary_3to4(eq3: IDSToplevel, eq4: IDSToplevel, deepcopy: boo
                 and ts3.boundary_secondary_separatrix.psi.has_value
             ):
                 n_nodes = 3
-        ts4.contour_tree.node.resize(n_nodes)
-        # Magnetic axis (primary O-point)
         node = ts4.contour_tree.node
-        node[0].critical_type = 0  # minimum (?)
+        node.resize(n_nodes)
+        # Magnetic axis (primary O-point)
+        axis_is_psi_minimum = (
+            # Note the sign flip for psi due to the COCOS change between DD3 and DD4!
+            -ts3.global_quantities.psi_axis < -ts3.global_quantities.psi_boundary
+        )
+        node[0].critical_type = 0 if axis_is_psi_minimum else 2
         node[0].r = ts3.global_quantities.magnetic_axis.r
         node[0].z = ts3.global_quantities.magnetic_axis.z
         node[0].psi = -ts3.global_quantities.psi_axis  # COCOS change
 
         # X-points
         if n_nodes >= 2:
             if ts3.boundary_separatrix.type == 0:  # limiter plasma
-                node[1].critical_type = 2  # maximum (?)
+                node[1].critical_type = 2 if axis_is_psi_minimum else 0
                 node[1].r = ts3.boundary_separatrix.active_limiter_point.r
                 node[1].z = ts3.boundary_separatrix.active_limiter_point.z
             else:
                 node[1].critical_type = 1  # saddle-point (x-point)
                 if len(ts3.boundary_separatrix.x_point):
                     node[1].r = ts3.boundary_separatrix.x_point[0].r
                     node[1].z = ts3.boundary_separatrix.x_point[0].z
-                    # TODO: what if there are multiple x-points?
+                # Additional x-points. N.B. levelset is only stored on the first node
+                for i in range(1, len(ts3.boundary_separatrix.x_point)):
+                    node.resize(len(node) + 1, keep=True)
+                    node[-1].critical_type = 1
+                    node[-1].r = ts3.boundary_separatrix.x_point[i].r
+                    node[-1].z = ts3.boundary_separatrix.x_point[i].z
+                    node[-1].psi = -ts3.boundary_separatrix.psi
             node[1].psi = -ts3.boundary_separatrix.psi  # COCOS change
             node[1].levelset.r = copy(ts3.boundary_separatrix.outline.r)
             node[1].levelset.z = copy(ts3.boundary_separatrix.outline.z)
@@ -1127,7 +1137,13 @@ def _equilibrium_boundary_3to4(eq3: IDSToplevel, eq4: IDSToplevel, deepcopy: boo
             if len(ts3.boundary_secondary_separatrix.x_point):
                 node[2].r = ts3.boundary_secondary_separatrix.x_point[0].r
                 node[2].z = ts3.boundary_secondary_separatrix.x_point[0].z
-                # TODO: what if there are multiple x-points?
+                # Additional x-points. N.B. levelset is only stored on the first node
+                for i in range(1, len(ts3.boundary_secondary_separatrix.x_point)):
+                    node.resize(len(node) + 1, keep=True)
+                    node[-1].critical_type = 1
+                    node[-1].r = ts3.boundary_secondary_separatrix.x_point[i].r
+                    node[-1].z = ts3.boundary_secondary_separatrix.x_point[i].z
+                    node[-1].psi = -ts3.boundary_secondary_separatrix.psi
             node[2].psi = -ts3.boundary_secondary_separatrix.psi  # COCOS change
             node[2].levelset.r = copy(ts3.boundary_secondary_separatrix.outline.r)
             node[2].levelset.z = copy(ts3.boundary_secondary_separatrix.outline.z)

imas/test/test_ids_convert.py

@@ -569,7 +569,7 @@ def test_3to4_equilibrium_boundary():
             ts.boundary_separatrix.gap.resize(1)
         if i == 3:
             # Fill second_separatrix
-            ts.boundary_secondary_separatrix.psi = -1.0
+            ts.boundary_secondary_separatrix.psi = -1.1
             # Use limiter for time_slice[1], otherwise divertor:
             ts.boundary_secondary_separatrix.outline.r = [0.9, 3.1, 2.1, 0.9]
             ts.boundary_secondary_separatrix.outline.z = [0.9, 2.1, 3.1, 0.9]
@@ -588,10 +588,11 @@ def test_3to4_equilibrium_boundary():
     assert len(eq4.time_slice) == 5
     for i, ts in enumerate(eq4.time_slice):
         node = ts.contour_tree.node
-        assert len(node) == [1, 2, 2, 3, 2][i]
+        assert len(node) == [1, 2, 2, 3, 3][i]
         # Test magnetic axis
         assert node[0].critical_type == 0
         assert node[0].r == node[0].z == 2.0
+        assert node[0].psi == -1.0
         assert len(node[0].levelset.r) == len(node[0].levelset.z) == 0
         # boundary_separatrix
         if i == 1:  # node[1] is boundary for limiter plasma
@@ -602,18 +603,33 @@ def test_3to4_equilibrium_boundary():
             assert node[1].critical_type == 1
             assert node[1].r == node[1].z == 1.0
         if i > 0:
+            assert node[1].psi == 1.0
             assert numpy.array_equal(node[1].levelset.r, [1.0, 3.0, 2.0, 1.0])
             assert numpy.array_equal(node[1].levelset.z, [1.0, 2.0, 3.0, 1.0])
         # boundary_secondary_separatrix
         if i == 3:
             assert node[2].critical_type == 1
             assert node[2].r == 2.1
             assert node[2].z == 3.1
+            assert node[2].psi == 1.1
             assert numpy.array_equal(node[2].levelset.r, [0.9, 3.1, 2.1, 0.9])
             assert numpy.array_equal(node[2].levelset.z, [0.9, 2.1, 3.1, 0.9])
+        # Second x-point from boundary_separatrix
+        if i == 4:
+            assert node[2].critical_type == 1
+            assert node[2].r == 2.0
+            assert node[2].z == 3.0
+            assert node[2].psi == node[1].psi == 1.0
+            # Levelset is only filled for the main x-point (node[1])
+            assert not node[2].levelset.r.has_value
+            assert not node[2].levelset.z.has_value
 
     # not deepcopied, should share numpy arrays
-    assert (
-        eq342.time_slice[1].boundary_separatrix.outline.r.value
-        is eq4.time_slice[1].contour_tree.node[1].levelset.r.value
+    slice1_outline_r = eq342.time_slice[1].boundary_separatrix.outline.r.value
+    assert slice1_outline_r is eq4.time_slice[1].contour_tree.node[1].levelset.r.value
+
+    # deepcopy should create a copy of the numpy arrays
+    eq4_cp = convert_ids(eq342, "4.0.0", deepcopy=True)
+    assert not numpy.may_share_memory(
+        slice1_outline_r, eq4_cp.time_slice[1].contour_tree.node[1].levelset.r.value
     )