Update SurfaceIntegralOutput for UserVariable

flow360/component/simulation/outputs/outputs.py

@@ -119,7 +119,12 @@ class _OutputBase(Flow360BaseModel):
     def _validate_improper_surface_field_usage(cls, value: UniqueItemList):
         if any(
             output_type in cls.__name__
-            for output_type in ["SurfaceProbeOutput", "SurfaceOutput", "SurfaceSliceOutput"]
+            for output_type in [
+                "SurfaceProbeOutput",
+                "SurfaceOutput",
+                "SurfaceSliceOutput",
+                "SurfaceIntegralOutput",
+            ]
         ):
             return value
         for output_item in value.items:

flow360/component/simulation/translator/solver_translator.py

@@ -97,6 +97,8 @@
 )
 from flow360.component.simulation.unit_system import LengthType
 from flow360.component.simulation.user_code.core.types import Expression, UserVariable
+from flow360.component.simulation.user_code.functions import math
+from flow360.component.simulation.user_code.variables import solution
 from flow360.component.simulation.user_defined_dynamics.user_defined_dynamics import (
     UserDefinedDynamic,
 )
@@ -765,7 +767,7 @@ def translate_streamline_output(output_params: list):
 
 
 def translate_output(input_params: SimulationParams, translated: dict):
-    # pylint: disable=too-many-branches,too-many-statements
+    # pylint: disable=too-many-branches,too-many-statements,too-many-locals
     """Translate output settings."""
     outputs = input_params.outputs
 
@@ -832,6 +834,33 @@ def translate_output(input_params: SimulationParams, translated: dict):
             outputs, TimeAverageProbeOutput, inject_probe_info
         )
     if has_instance_in_list(outputs, SurfaceIntegralOutput):
+        for output in outputs:
+            if not isinstance(output, SurfaceIntegralOutput):
+                continue
+            output_fields_processed = []
+            for output_field in output.output_fields.items:
+                if isinstance(output_field, UserVariable):
+                    expression = Expression.model_validate(output_field.value)
+                    if expression.length == 0:
+                        expression_processed = expression * math.magnitude(
+                            solution.node_area_vector
+                        )
+                    else:
+                        expression_processed = [
+                            expression[i] * math.magnitude(solution.node_area_vector)
+                            for i in range(expression.length)
+                        ]
+
+                    output_fields_processed.append(
+                        UserVariable(
+                            name=output_field.name + "_integral",
+                            value=expression_processed,
+                        )
+                    )
+                else:
+                    output_fields_processed.append(output_field)
+            output.output_fields.items = output_fields_processed
+
         integral_output = translate_monitor_output(
             outputs, SurfaceIntegralOutput, inject_surface_list_info
         )

flow360/component/simulation/translator/user_expression_utils.py

@@ -89,6 +89,10 @@
     + "{___CfVec[i] = wallShearStress[i] / (0.5 * MachRef * MachRef);}",
     "solution.Cf": "double ___Cf;"
     + "___Cf = magnitude(wallShearStress) / (0.5 * MachRef * MachRef);",
+    "solution.node_unit_normal": "double ___node_unit_normal[3];"
+    + "double ___normalMag = magnitude(nodeNormals);"
+    + "for (int i = 0; i < 3; i++){___node_unit_normal[i] = "
+    + "nodeNormals[i] / ___normalMag;}",
     "solution.node_forces_per_unit_area": "double ___node_forces_per_unit_area[3];"
     + "double ___normalMag = magnitude(nodeNormals);"
     + "for (int i = 0; i < 3; i++){___node_forces_per_unit_area[i] = "

flow360/component/simulation/user_code/core/context.py

@@ -113,7 +113,8 @@ def _import_solution(_) -> Any:
             "CfVec",
             "Cf",
             "heat_flux",
-            "node_normals",
+            "node_area_vector",
+            "node_unit_normal",
             "node_forces_per_unit_area",
             "y_plus",
             "wall_shear_stress_magnitude",

flow360/component/simulation/user_code/variables/solution.py

@@ -194,12 +194,18 @@
     solver_name="heatFlux",
     variable_type="Surface",
 )
-node_normals = SolverVariable(
-    name="solution.node_normals",
-    value=[float("NaN"), float("NaN"), float("NaN")],
+node_area_vector = SolverVariable(
+    name="solution.node_area_vector",
+    value=[float("NaN"), float("NaN"), float("NaN")] * u.m**2,
     solver_name="nodeNormals",
     variable_type="Surface",
 )
+node_unit_normal = SolverVariable(
+    name="solution.node_unit_normal",
+    value=[float("NaN"), float("NaN"), float("NaN")],
+    solver_name="___node_unit_normal",
+    variable_type="Surface",
+)
 node_forces_per_unit_area = SolverVariable(
     name="solution.node_forces_per_unit_area",
     value=[float("NaN"), float("NaN"), float("NaN")] * u.Pa,

flow360/component/simulation/validation/validation_output.py

@@ -5,23 +5,36 @@
 from typing import List, Literal, Union, get_args, get_origin
 
 from flow360.component.simulation.models.volume_models import Fluid
-from flow360.component.simulation.outputs.outputs import AeroAcousticOutput
+from flow360.component.simulation.outputs.outputs import (
+    AeroAcousticOutput,
+    SurfaceIntegralOutput,
+)
 from flow360.component.simulation.time_stepping.time_stepping import Steady
 
 
 def _check_output_fields(params):
     """Check the specified output fields for each output item is valid."""
 
+    # pylint: disable=too-many-branches
     if params.outputs is None:
         return params
 
-    has_surface_integral_output = any(
-        output.output_type == "SurfaceIntegralOutput" for output in params.outputs
-    )
+    has_legacy_user_defined_field_in_surface_integral_output = False
+    for output in params.outputs:
+        if isinstance(output, SurfaceIntegralOutput):
+            for output_field in output.output_fields.items:
+                if isinstance(output_field, str):
+                    has_legacy_user_defined_field_in_surface_integral_output = True
+                    break
     has_user_defined_fields = len(params.user_defined_fields) > 0
 
-    if has_surface_integral_output and has_user_defined_fields is False:
-        raise ValueError("`SurfaceIntegralOutput` can only be used with `UserDefinedField`.")
+    if (
+        has_legacy_user_defined_field_in_surface_integral_output
+        and has_user_defined_fields is False
+    ):
+        raise ValueError(
+            "The legacy string output fields in `SurfaceIntegralOutput` must be used with `UserDefinedField`."
+        )
 
     def extract_literal_values(annotation):
         origin = get_origin(annotation)

tests/simulation/params/test_validators_params.py

@@ -83,6 +83,9 @@
 from flow360.component.simulation.simulation_params import SimulationParams
 from flow360.component.simulation.time_stepping.time_stepping import Steady, Unsteady
 from flow360.component.simulation.unit_system import SI_unit_system
+from flow360.component.simulation.user_code.core.types import UserVariable
+from flow360.component.simulation.user_code.functions import math
+from flow360.component.simulation.user_code.variables import solution
 from flow360.component.simulation.user_defined_dynamics.user_defined_dynamics import (
     UserDefinedDynamic,
 )
@@ -1152,7 +1155,7 @@ def test_output_fields_with_user_defined_fields():
                 ],
             )
 
-    msg = "`SurfaceIntegralOutput` can only be used with `UserDefinedField`."
+    msg = "The legacy string output fields in `SurfaceIntegralOutput` must be used with `UserDefinedField`."
     with pytest.raises(ValueError, match=re.escape(msg)):
         with SI_unit_system:
             _ = SimulationParams(
@@ -1165,6 +1168,24 @@ def test_output_fields_with_user_defined_fields():
                 ]
             )
 
+    with SI_unit_system:
+        params = SimulationParams(
+            outputs=[
+                SurfaceIntegralOutput(
+                    name="MassFluxIntegral",
+                    output_fields=[
+                        UserVariable(
+                            name="MassFluxProjected",
+                            value=-1
+                            * solution.density
+                            * math.dot(solution.velocity, solution.node_area_vector),
+                        )
+                    ],
+                    surfaces=[surface_1],
+                )
+            ]
+        )
+
     msg = (
         "In `outputs`[1] IsosurfaceOutput:, Cpp is not a valid iso field name. Allowed fields are "
     )

tests/simulation/translator/ref/Flow360_user_variable.json

@@ -53,6 +53,23 @@
         "outputFields": [],
         "outputFormat": "paraview"
     },
+    "monitorOutput": {
+        "monitors": {
+            "MassFluxIntegral": {
+                "animationFrequency": 1,
+                "animationFrequencyOffset": 0,
+                "computeTimeAverages": false,
+                "outputFields": [
+                    "MassFluxProjected_integral"
+                ],
+                "surfaces": [
+                    "VOLUME/LEFT"
+                ],
+                "type": "surfaceIntegral"
+            }
+        },
+        "outputFields": []
+    },
     "navierStokesSolver": {
         "CFLMultiplier": 1.0,
         "absoluteTolerance": 1e-10,
@@ -250,6 +267,10 @@
             "expression": "double ___velocity[3];___velocity[0] = primitiveVars[1] * velocityScale;___velocity[1] = primitiveVars[2] * velocityScale;___velocity[2] = primitiveVars[3] * velocityScale;iso_field_velocity_mag = (pow(((pow(___velocity[0], 2) + pow(___velocity[1], 2)) + pow(___velocity[2], 2)), 0.5) * 200.0);",
             "name": "iso_field_velocity_mag"
         },
+        {
+            "expression": "double ___node_unit_normal[3];double ___normalMag = magnitude(nodeNormals);for (int i = 0; i < 3; i++){___node_unit_normal[i] = nodeNormals[i] / ___normalMag;}double ___velocity[3];___velocity[0] = primitiveVars[1] * velocityScale;___velocity[1] = primitiveVars[2] * velocityScale;___velocity[2] = primitiveVars[3] * velocityScale;MassFluxProjected_integral = ((((-1 * primitiveVars[0]) * (((___velocity[0] * ___node_unit_normal[0]) + (___velocity[1] * ___node_unit_normal[1])) + (___velocity[2] * ___node_unit_normal[2]))) * pow(((pow(nodeNormals[0], 2) + pow(nodeNormals[1], 2)) + pow(nodeNormals[2], 2)), 0.5)) * 200000.0);",
+            "name": "MassFluxProjected_integral"
+        },
         {
             "expression": "double ___CfVec[3]; for (int i = 0; i < 3; i++){___CfVec[i] = wallShearStress[i] / (0.5 * MachRef * MachRef);}exp_res = (exp(___CfVec[0]) * 1);",
             "name": "exp_res"

tests/simulation/translator/test_output_translation.py

@@ -832,15 +832,15 @@ def surface_integral_output_config(vel_in_km_per_hr):
                     "animationFrequency": 1,
                     "animationFrequencyOffset": 0,
                     "computeTimeAverages": False,
-                    "outputFields": ["My_field_1", "velocity_in_km_per_hr"],
+                    "outputFields": ["My_field_1", "velocity_in_km_per_hr_integral"],
                     "surfaces": ["zoneName/surface1", "surface2"],
                     "type": "surfaceIntegral",
                 },
                 "prb 122": {
                     "animationFrequency": 1,
                     "animationFrequencyOffset": 0,
                     "computeTimeAverages": False,
-                    "outputFields": ["My_field_2", "velocity_in_km_per_hr"],
+                    "outputFields": ["My_field_2", "velocity_in_km_per_hr_integral"],
                     "surfaces": ["surface21", "surface22"],
                     "type": "surfaceIntegral",
                 },
@@ -1049,7 +1049,7 @@ def test_monitor_output(
                 "animationFrequency": 1,
                 "animationFrequencyOffset": 0,
                 "computeTimeAverages": False,
-                "outputFields": ["My_field_1", "velocity_in_km_per_hr"],
+                "outputFields": ["My_field_1", "velocity_in_km_per_hr_integral"],
                 "surfaces": ["zoneName/surface1", "surface2"],
                 "type": "surfaceIntegral",
             },
@@ -1067,7 +1067,7 @@ def test_monitor_output(
                 "animationFrequency": 1,
                 "animationFrequencyOffset": 0,
                 "computeTimeAverages": False,
-                "outputFields": ["My_field_2", "velocity_in_km_per_hr"],
+                "outputFields": ["My_field_2", "velocity_in_km_per_hr_integral"],
                 "surfaces": ["surface21", "surface22"],
                 "type": "surfaceIntegral",
             },

tests/simulation/translator/test_solver_translator.py

@@ -46,6 +46,7 @@
     Isosurface,
     IsosurfaceOutput,
     SliceOutput,
+    SurfaceIntegralOutput,
     SurfaceOutput,
     UserDefinedField,
     VolumeOutput,
@@ -123,6 +124,7 @@
 
 assertions = unittest.TestCase("__init__")
 
+import flow360.component.simulation.user_code.core.context as context
 from flow360.component.simulation.framework.updater_utils import compare_values
 from flow360.component.simulation.models.volume_models import (
     AngleExpression,
@@ -134,6 +136,11 @@
 from flow360.component.simulation.time_stepping.time_stepping import Unsteady
 
 
+@pytest.fixture(autouse=True)
+def reset_context():
+    clear_context()
+
+
 @pytest.fixture()
 def get_om6Wing_tutorial_param():
     my_wall = Surface(name="1")
@@ -660,14 +667,6 @@ def test_liquid_simulation_translation():
     translate_and_compare(param, mesh_unit=1 * u.m, ref_json_file="Flow360_liquid_rotation_dd.json")
 
 
-import flow360.component.simulation.user_code.core.context as context
-
-
-@pytest.fixture()
-def reset_context():
-    clear_context()
-
-
 def test_param_with_user_variables():
     some_dependent_variable_a = UserVariable(
         name="some_dependent_variable_a", value=[1.0 * u.m / u.s, 2.0 * u.m / u.s, 3.0 * u.m / u.s]
@@ -712,6 +711,10 @@ def test_param_with_user_variables():
     my_temperature = UserVariable(
         name="my_temperature", value=(solution.temperature + (-10 * u.K)) * 1.8
     )
+    surface_integral_variable = UserVariable(
+        name="MassFluxProjected",
+        value=-1 * solution.density * math.dot(solution.velocity, solution.node_unit_normal),
+    )
     iso_field_pressure = UserVariable(
         name="iso_field_pressure",
         value=0.5 * solution.Cp * solution.density * math.magnitude(solution.velocity) ** 2,
@@ -796,6 +799,11 @@ def test_param_with_user_variables():
                         ),
                     ],
                 ),
+                SurfaceIntegralOutput(
+                    name="MassFluxIntegral",
+                    output_fields=[surface_integral_variable],
+                    entities=Surface(name="VOLUME/LEFT"),
+                ),
                 SurfaceOutput(
                     name="surface_output",
                     entities=Surface(name="fluid/body"),