validation rule for sign of EM fields

Author: MasanariHosokawa

imas_validator/assets/rulesets/generic/equilibrium.py

@@ -1,50 +1,101 @@
 """Generic rules applying to the equilibrium IDS"""
 
-from imas_validator.common.cocos import compute_COCOS
+import numpy as np
 
-dd_cocos = {"3": 11, "4": 17}
 
-
-@validator("equilibrium")
-def validate_cocos(ids):
+@validator("equilibrium", version="<4.0.0")
+def validate_sign_DD3(ids):
     """
-    Validate that COCOS computed corresponds to the one in DD.
+    Validate that the equilibrium IDS has consistent sign conventions.
     """
 
-    ver = str(ids.ids_properties.version_put.data_dictionary.value)[0]
-    ref = dd_cocos.get(ver)
-    if ref is None:
-        raise ValueError(f"Unsupported DD version_put: {ver}")
+    # Loop through each time slice
+    for itime, time_slice in enumerate(ids.time_slice):
+
+        # Get the time for this time slice
+        if ids.ids_properties.homogeneous_time == 0:
+            time = time_slice.time
+        else:
+            time = ids.time[itime]
+
+        # Set sign of Ip
+        if time_slice.global_quantities.ip.has_value:
+            ipsign = np.sign(time_slice.global_quantities.ip)
+        else:
+            ipsign = -9
+
+        # Set sign of B0
+        if ids.vacuum_toroidal_field.b0.has_value:
+            b0sign = np.sign(
+                np.interp(time, ids.time, ids.vacuum_toroidal_field.b0)
+            )
+        else:
+            b0sign = -9
+
+        #
+        # Validate the sign conventions for EM fields
+        #
+
+        if time_slice.profiles_1d.f.has_value:
+            assert np.sign(time_slice.profiles_1d.f) == b0sign
+
+        if time_slice.profiles_1d.phi.has_value:
+            # phi at center found as zero in DINA, CORSICA and etc.
+            assert np.sign(time_slice.profiles_1d.phi[1:]) == b0sign
+
+        psi_axis = time_slice.global_quantities.psi_axis
+        psi_boundary = time_slice.global_quantities.psi_boundary
+        if psi_axis.has_value and psi_boundary.has_value:
+            assert np.sign(psi_boundary - psi_axis) == ipsign
+
+        if time_slice.global_quantities.q_95.has_value:
+            assert np.sign(time_slice.global_quantities.q_95) == ipsign * b0sign
 
-    # time_slice[:]
+
+@validator("equilibrium", version=">=4.0.0")
+def validate_sign(ids):
+    """
+    Validate that the equilibrium IDS has consistent sign conventions.
+    """
+
+    # Loop through each time slice
     for itime, time_slice in enumerate(ids.time_slice):
 
-        # time_slice[:].profiles_2d
-        for i1, profiles_2d in enumerate(time_slice.profiles_2d):
-
-            if not (
-                ids.vacuum_toroidal_field.b0.has_value
-                and time_slice.global_quantities.ip.has_value
-                and time_slice.global_quantities.magnetic_axis.z.has_value
-                and time_slice.profiles_1d.psi.has_value
-                and time_slice.profiles_1d.q.has_value
-                and time_slice.profiles_1d.dpressure_dpsi.has_value
-                and profiles_2d.b_field_z.has_value
-                and profiles_2d.grid_type.index == 1
-                and profiles_2d.psi.has_value
-                and profiles_2d.r.has_value
-            ):
-                continue
-
-            # Compute COCOS
-            try:
-                # ids._obj instead of ids since '&' operator not supported
-                # yet in IMAS-Validator
-                cocos = compute_COCOS(ids._obj, itime, i1)["COCOS"]
-            except Exception:
-                cocos = None
-
-            assert ref == cocos, (
-                f"COCOS mismatch for time_slice {itime}, profiles_2d {i1}, "
-                f"DD version_put/computed: {ref}/{cocos}"
+        # Get the time for this time slice
+        if ids.ids_properties.homogeneous_time == 0:
+            time = time_slice.time
+        else:
+            time = ids.time[itime]
+
+        # Set sign of Ip
+        if time_slice.global_quantities.ip.has_value:
+            ipsign = np.sign(time_slice.global_quantities.ip)
+        else:
+            ipsign = -9
+
+        # Set sign of B0
+        if ids.vacuum_toroidal_field.b0.has_value:
+            b0sign = np.sign(
+                np.interp(time, ids.time, ids.vacuum_toroidal_field.b0)
             )
+        else:
+            b0sign = -9
+
+        #
+        # Validate the sign conventions for EM fields
+        #
+
+        if time_slice.profiles_1d.f.has_value:
+            assert np.sign(time_slice.profiles_1d.f) == b0sign
+
+        if time_slice.profiles_1d.phi.has_value:
+            # phi at center found as zero in DINA, CORSICA and etc.
+            assert np.sign(time_slice.profiles_1d.phi[1:]) == b0sign
+
+        psi_axis = time_slice.global_quantities.psi_axis
+        psi_boundary = time_slice.global_quantities.psi_boundary
+        if psi_axis.has_value and psi_boundary.has_value:
+            assert np.sign(psi_boundary - psi_axis) == -ipsign
+
+        if time_slice.global_quantities.q_95.has_value:
+            assert np.sign(time_slice.global_quantities.q_95) == ipsign * b0sign