Separate prepending code to declaration and computation parts

flow360/component/simulation/translator/solver_translator.py

@@ -77,6 +77,9 @@
 from flow360.component.simulation.primitives import Box, SurfacePair
 from flow360.component.simulation.simulation_params import SimulationParams
 from flow360.component.simulation.time_stepping.time_stepping import Steady, Unsteady
+from flow360.component.simulation.translator.user_expression_utils import (
+    udf_prepending_code,
+)
 from flow360.component.simulation.translator.utils import (
     _get_key_name,
     convert_tuples_to_lists,
@@ -97,95 +100,6 @@
 )
 from flow360.exceptions import Flow360TranslationError
 
-udf_prepending_code = {
-    "solution.Cp": "double Cp = (primitiveVars[4] - pressureFreestream) / (0.5 * MachRef * MachRef);",
-    "solution.Cpt": "double MachUser = sqrt(primitiveVars[1] * primitiveVars[1] + "
-    + "primitiveVars[2] * primitiveVars[2] + primitiveVars[3] * primitiveVars[3])"
-    + "/sqrt(1.4 * primitiveVars[4] / primitiveVars[0]);"
-    + "double Cpt = (1.4 * primitiveVars[4] * pow(1.0 + (1.4 - 1.0) / 2. * MachUser * MachUser,"
-    + "1.4 / (1.4 - 1.0)) - pow(1.0 + (1.4 - 1.0) / 2. * MachRef * MachRef,"
-    + "1.4 / (1.4 - 1.0))) / (0.5 * 1.4 * MachRef * MachRef);",
-    "solution.grad_density": "double gradDensity[3] = {gradPrimitive[0][0], "
-    + "gradPrimitive[0][1], gradPrimitive[0][2]};",
-    "solution.grad_velocity_x": "double gradVelocityX[3] = {gradPrimitive[1][0], "
-    + " gradPrimitive[1][1], gradPrimitive[1][2]};",
-    "solution.grad_velocity_y": "double gradVelocityY[3] = {gradPrimitive[2][0],"
-    + "gradPrimitive[2][1], gradPrimitive[2][2]};",
-    "solution.grad_velocity_z": "double gradVelocityZ[3] = {gradPrimitive[3][0], "
-    + "gradPrimitive[3][1], gradPrimitive[3][2]};",
-    "solution.grad_pressure": "double gradPressure[3] = {gradPrimitive[4][0], "
-    + "gradPrimitive[4][1], gradPrimitive[4][2]};",
-    "solution.Mach": "double Mach = usingLiquidAsMaterial ? 0 : "
-    + "sqrt(primitiveVars[1] * primitiveVars[1] + "
-    + "primitiveVars[2] * primitiveVars[2] + "
-    + "primitiveVars[3] * primitiveVars[3]) / "
-    + "sqrt(1.4 * primitiveVars[4] / primitiveVars[0]);",
-    "solution.mut_ratio": "double mutRatio;mutRatio = mut / mu;",
-    "solution.velocity": "double velocity[3];"
-    + "velocity[0] = primitiveVars[1] * velocityScale;"
-    + "velocity[1] = primitiveVars[2] * velocityScale;"
-    + "velocity[2] = primitiveVars[3] * velocityScale;",
-    "solution.qcriterion": "double qcriterion;"
-    + "double ux = gradPrimitive[1][0];"
-    + "double uy = gradPrimitive[1][1];"
-    + "double uz = gradPrimitive[1][2];"
-    + "double vx = gradPrimitive[2][0];"
-    + "double vy = gradPrimitive[2][1];"
-    + "double vz = gradPrimitive[2][2];"
-    + "double wx = gradPrimitive[3][0];"
-    + "double wy = gradPrimitive[3][1];"
-    + "double wz = gradPrimitive[3][2];"
-    + "double str11 = ux;"
-    + "double str22 = vy;"
-    + "double str33 = wz;"
-    + "double str12 = 0.5 * (uy + vx);"
-    + "double str13 = 0.5 * (uz + wx);"
-    + "double str23 = 0.5 * (vz + wy);"
-    + "double str_norm = str11 * str11 + str22 * str22 + str33 * str33 + "
-    + "2 * (str12 * str12) + 2 * (str13 * str13) + 2 * (str23 * str23);"
-    + "double omg12 = 0.5 * (uy - vx);"
-    + "double omg13 = 0.5 * (uz - wx);"
-    + "double omg23 = 0.5 * (vz - wy);"
-    + "double omg_norm = 2 * (omg12 * omg12) + 2 * (omg13 * omg13) + 2 * (omg23 * omg23);"
-    + "qcriterion = 0.5 * (omg_norm - str_norm) * (velocityScale * velocityScale);",
-    "solution.entropy": "double entropy;entropy = log(primitiveVars[4] / (1.0 / 1.4) / pow(primitiveVars[0], 1.4));",
-    "solution.temperature": f"double epsilon = {np.finfo(np.float64).eps};"
-    "double temperature = (primitiveVars[0] < epsilon && HeatEquation_solution != nullptr) ? "
-    "HeatEquation_solution[0] : primitiveVars[4] / (primitiveVars[0] * (1.0 / 1.4));",
-    "solution.vorticity": "double vorticity[3];"
-    + "vorticity[0] = (gradPrimitive[3][1] - gradPrimitive[2][2]) * velocityScale;"
-    + "vorticity[1] = (gradPrimitive[1][2] - gradPrimitive[3][0]) * velocityScale;"
-    + "vorticity[2] = (gradPrimitive[2][0] - gradPrimitive[1][1]) * velocityScale;",
-    "solution.CfVec": "double CfVec[3];"
-    + "for (int i = 0; i < 3; i++)"
-    + "{CfVec[i] = wallShearStress[i] / (0.5 * MachRef * MachRef);}",
-    "solution.Cf": "double Cf;Cf = magnitude(wallShearStress) / (0.5 * MachRef * MachRef);",
-    "solution.node_forces_per_unit_area": "double nodeForcesPerUnitArea[3];"
-    + "double normalMag = magnitude(nodeNormals);"
-    + "for (int i = 0; i < 3; i++){nodeForcesPerUnitArea[i] = "
-    + "((primitiveVars[4] - pressureFreestream) * nodeNormals[i] / normalMag + wallViscousStress[i])"
-    + " * (velocityScale * velocityScale);}",
-    "solution.heat_transfer_coefficient_static_temperature": "double temperature = "
-    + "primitiveVars[4] / (primitiveVars[0] * 1.0 / 1.4);"
-    + f"double temperatureSafeDivide; double epsilon = {np.finfo(np.float64).eps};"
-    + "temperatureSafeDivide = (temperature - 1.0 < 0) ? "
-    + "temperature - 1.0 - epsilon : "
-    + "temperature - 1.0 + epsilon;"
-    + "double heatTransferCoefficientStaticTemperature = "
-    + "abs(temperature - 1.0) > epsilon ? "
-    + "- heatFlux / temperatureSafeDivide :  1.0 / epsilon;",
-    "solution.heat_transfer_coefficient_total_temperature": "double temperature = "
-    + "primitiveVars[4] / (primitiveVars[0] * 1.0 / 1.4);"
-    + "double temperatureTotal = 1.0 + (1.4 - 1.0) / 2.0 * MachRef * MachRef;"
-    + f"double temperatureSafeDivide; double epsilon = {np.finfo(np.float64).eps};"
-    + "temperatureSafeDivide = (temperature - temperatureTotal < 0) ? "
-    + "temperature - temperatureTotal - epsilon : "
-    + "temperature - temperatureTotal + epsilon;"
-    + "double heatTransferCoefficientTotalTemperature = "
-    + "abs(temperature - temperatureTotal) > epsilon ? "
-    + "temperatureTotal = - heatFlux / temperatureSafeDivide :  1.0 / epsilon;",
-}
-
 
 def dump_dict(input_params):
     """Dumping param/model to dictionary."""
@@ -660,6 +574,22 @@ def _compute_coefficient_and_offset(source_unit: u.Unit, target_unit: u.Unit):
 
         return coefficient, offset
 
+    def _prepare_prepending_code(expression: Expression):
+        prepending_code = []
+        for name in expression.solver_variable_names():
+            if not udf_prepending_code.get(name):
+                continue
+            if name.split(".")[-1] == variable.name:
+                # Avoid duplicate declaration if the intermediate variable name is
+                # the same as the solver_name.
+                prepending_code.append(udf_prepending_code[name]["computation"])
+                continue
+            prepending_code.append(
+                udf_prepending_code[name]["declaration"] + udf_prepending_code[name]["computation"]
+            )
+        prepending_code = "".join(prepending_code)
+        return prepending_code
+
     expression: Expression = variable.value
 
     requested_unit: Union[u.Unit, None] = expression.get_output_units(input_params=input_params)
@@ -677,12 +607,7 @@ def _compute_coefficient_and_offset(source_unit: u.Unit, target_unit: u.Unit):
         )
 
     expression_length = expression.length
-    prepending_code = [
-        udf_prepending_code[name]
-        for name in expression.solver_variable_names()
-        if udf_prepending_code.get(name)
-    ]
-    prepending_code = "".join(prepending_code)
+    prepending_code = _prepare_prepending_code(expression=expression)
 
     if expression_length == 1:
         expression = expression.evaluate(raise_on_non_evaluable=False, force_evaluate=False)

flow360/component/simulation/translator/user_expression_utils.py

@@ -0,0 +1,180 @@
+"""Utilities for user expression translation."""
+
+import numpy as np
+
+udf_prepending_code = {
+    "solution.Cp": {
+        "declaration": "double Cp;",
+        "computation": "Cp = (primitiveVars[4] - pressureFreestream) / (0.5 * MachRef * MachRef);",
+    },
+    "solution.Cpt": {
+        "declaration": "double Cpt;",
+        "computation": "double MachTmp = sqrt(primitiveVars[1] * primitiveVars[1] + "
+        + "primitiveVars[2] * primitiveVars[2] + primitiveVars[3] * primitiveVars[3]) / "
+        + "sqrt(1.4 * primitiveVars[4] / primitiveVars[0]);"
+        + "Cpt = (1.4 * primitiveVars[4] * pow(1.0 + (1.4 - 1.0) / 2. * MachTmp * MachTmp,"
+        + "1.4 / (1.4 - 1.0)) - pow(1.0 + (1.4 - 1.0) / 2. * MachRef * MachRef,"
+        + "1.4 / (1.4 - 1.0))) / (0.5 * 1.4 * MachRef * MachRef);",
+    },
+    "solution.grad_density": {
+        "declaration": "double gradDensity[3];",
+        "computation": " gradDensity[0] = gradPrimitive[0][0];"
+        + "gradDensity[1] = gradPrimitive[0][1];"
+        + "gradDensity[2] = gradPrimitive[0][2];",
+    },
+    "solution.grad_velocity_x": {
+        "declaration": "double gradVelocityX[3];",
+        "computation": "gradVelocityX[0] = gradPrimitive[1][0] * velocityScale;"
+        + "gradVelocityX[1] = gradPrimitive[1][1] * velocityScale;"
+        + "gradVelocityX[2] = gradPrimitive[1][2] * velocityScale;",
+    },
+    "solution.grad_velocity_y": {
+        "declaration": "double gradVelocityY[3];",
+        "computation": "gradVelocityY[0] = gradPrimitive[2][0] * velocityScale;"
+        + "gradVelocityY[1] = gradPrimitive[2][1] * velocityScale;"
+        + "gradVelocityY[2] = gradPrimitive[2][2] * velocityScale;",
+    },
+    "solution.grad_velocity_z": {
+        "declaration": "double gradVelocityZ[3];",
+        "computation": "gradVelocityZ[0] = gradPrimitive[3][0] * velocityScale;"
+        + "gradVelocityZ[1] = gradPrimitive[3][1] * velocityScale;"
+        + "gradVelocityZ[2] = gradPrimitive[3][2] * velocityScale;",
+    },
+    "solution.grad_pressure": {
+        "declaration": "double gradPressure[3];",
+        "computation": "gradPressure[0] = gradPrimitive[4][0]; "
+        + "gradPressure[1] = gradPrimitive[4][1]; "
+        + "gradPressure[2] = gradPrimitive[4][2];",
+    },
+    "solution.Mach": {
+        "declaration": "double Mach;",
+        "computation": "Mach = usingLiquidAsMaterial ? 0 : "
+        + "sqrt(primitiveVars[1] * primitiveVars[1] + "
+        + "primitiveVars[2] * primitiveVars[2] + "
+        + "primitiveVars[3] * primitiveVars[3]) / "
+        + "sqrt(1.4 * primitiveVars[4] / primitiveVars[0]);",
+    },
+    "solution.mut_ratio": {
+        "declaration": "double mutRatio;",
+        "computation": "mutRatio = mut / mu;",
+    },
+    "solution.nu_hat": {
+        "declaration": "double nuHat;",
+        "computation": "nuHat = SpalartAllmaras_solution * velocityScale;",
+    },
+    "solution.turbulence_kinetic_energy": {
+        "declaration": "double turbulenceKineticEnergy;",
+        "computation": "turbulenceKineticEnergy = kOmegaSST_solution[0] * pow(velocityScale, 2);",
+    },
+    "solution.specific_rate_of_dissipation": {
+        "declaration": "double specificRateOfDissipation;",
+        "computation": "specificRateOfDissipation = kOmegaSST_solution[1] * velocityScale;",
+    },
+    "solution.velocity": {
+        "declaration": "double velocity[3];",
+        "computation": "velocity[0] = primitiveVars[1] * velocityScale;"
+        + "velocity[1] = primitiveVars[2] * velocityScale;"
+        + "velocity[2] = primitiveVars[3] * velocityScale;",
+    },
+    "solution.velocity_magnitude": {
+        "declaration": "double velocityMagnitude;",
+        "computation": "double velocityTmp[3];velocityTmp[0] = primitiveVars[1] * velocityScale;"
+        + "velocityTmp[1] = primitiveVars[2] * velocityScale;"
+        + "velocityTmp[2] = primitiveVars[3] * velocityScale;"
+        + "velocityMagnitude = magnitude(velocityTmp);",
+    },
+    "solution.qcriterion": {
+        "declaration": "double qcriterion;",
+        "computation": "double ux = gradPrimitive[1][0];"
+        + "double uy = gradPrimitive[1][1];"
+        + "double uz = gradPrimitive[1][2];"
+        + "double vx = gradPrimitive[2][0];"
+        + "double vy = gradPrimitive[2][1];"
+        + "double vz = gradPrimitive[2][2];"
+        + "double wx = gradPrimitive[3][0];"
+        + "double wy = gradPrimitive[3][1];"
+        + "double wz = gradPrimitive[3][2];"
+        + "double str11 = ux;"
+        + "double str22 = vy;"
+        + "double str33 = wz;"
+        + "double str12 = 0.5 * (uy + vx);"
+        + "double str13 = 0.5 * (uz + wx);"
+        + "double str23 = 0.5 * (vz + wy);"
+        + "double str_norm = str11 * str11 + str22 * str22 + str33 * str33 + "
+        + "2 * (str12 * str12) + 2 * (str13 * str13) + 2 * (str23 * str23);"
+        + "double omg12 = 0.5 * (uy - vx);"
+        + "double omg13 = 0.5 * (uz - wx);"
+        + "double omg23 = 0.5 * (vz - wy);"
+        + "double omg_norm = 2 * (omg12 * omg12) + 2 * (omg13 * omg13) + 2 * (omg23 * omg23);"
+        + "qcriterion = 0.5 * (omg_norm - str_norm) * (velocityScale * velocityScale);",
+    },
+    "solution.entropy": {
+        "declaration": "double entropy;",
+        "computation": "entropy = log(primitiveVars[4] / (1.0 / 1.4) / pow(primitiveVars[0], 1.4));",
+    },
+    "solution.temperature": {
+        "declaration": "double temperature;",
+        "computation": f"double epsilon = {np.finfo(np.float64).eps};"
+        "temperature = (primitiveVars[0] < epsilon && HeatEquation_solution != nullptr) ? "
+        "HeatEquation_solution[0] : primitiveVars[4] / (primitiveVars[0] * (1.0 / 1.4));",
+    },
+    "solution.vorticity": {
+        "declaration": "double vorticity[3];",
+        "computation": "vorticity[0] = (gradPrimitive[3][1] - gradPrimitive[2][2]) * velocityScale;"
+        + "vorticity[1] = (gradPrimitive[1][2] - gradPrimitive[3][0]) * velocityScale;"
+        + "vorticity[2] = (gradPrimitive[2][0] - gradPrimitive[1][1]) * velocityScale;",
+    },
+    "solution.vorticity_magnitude": {
+        "declaration": "double vorticityMagnitude;",
+        "computation": "double vorticityTmp[3];"
+        + "vorticityTmp[0] = (gradPrimitive[3][1] - gradPrimitive[2][2]) * velocityScale;"
+        + "vorticityTmp[1] = (gradPrimitive[1][2] - gradPrimitive[3][0]) * velocityScale;"
+        + "vorticityTmp[2] = (gradPrimitive[2][0] - gradPrimitive[1][1]) * velocityScale;"
+        + "vorticityMagnitude = magnitude(vorticityTmp);",
+    },
+    "solution.CfVec": {
+        "declaration": "double CfVec[3];",
+        "computation": "for (int i = 0; i < 3; i++)"
+        + "{CfVec[i] = wallShearStress[i] / (0.5 * MachRef * MachRef);}",
+    },
+    "solution.Cf": {
+        "declaration": "double Cf;",
+        "computation": "Cf = magnitude(wallShearStress) / (0.5 * MachRef * MachRef);",
+    },
+    "solution.node_forces_per_unit_area": {
+        "declaration": "double nodeForcesPerUnitArea[3];",
+        "computation": "double normalMag = magnitude(nodeNormals);"
+        + "for (int i = 0; i < 3; i++){nodeForcesPerUnitArea[i] = "
+        + "((primitiveVars[4] - pressureFreestream) * nodeNormals[i] / normalMag + wallViscousStress[i])"
+        + " * (velocityScale * velocityScale);}",
+    },
+    "solution.heat_transfer_coefficient_static_temperature": {
+        "declaration": "double heatTransferCoefficientStaticTemperature;",
+        "computation": "double temperatureTmp = "
+        + "primitiveVars[4] / (primitiveVars[0] * 1.0 / 1.4);"
+        + f"double epsilon = {np.finfo(np.float64).eps};"
+        + "double temperatureSafeDivide = (temperatureTmp - 1.0 < 0) ? "
+        + "temperatureTmp - 1.0 - epsilon : "
+        + "temperatureTmp - 1.0 + epsilon;"
+        + "heatTransferCoefficientStaticTemperature = "
+        + "abs(temperatureTmp - 1.0) > epsilon ? "
+        + "- wallHeatFlux / temperatureSafeDivide :  1.0 / epsilon;",
+    },
+    "solution.heat_transfer_coefficient_total_temperature": {
+        "declaration": "double heatTransferCoefficientTotalTemperature;",
+        "computation": "double temperatureTmp = "
+        + "primitiveVars[4] / (primitiveVars[0] * 1.0 / 1.4);"
+        + "double temperatureTotal = 1.0 + (1.4 - 1.0) / 2.0 * MachRef * MachRef;"
+        + f"double epsilon = {np.finfo(np.float64).eps};"
+        + "double temperatureSafeDivide = (temperatureTmp - temperatureTotal < 0) ? "
+        + "temperatureTmp - temperatureTotal - epsilon : "
+        + "temperatureTmp - temperatureTotal + epsilon;"
+        + "double heatTransferCoefficientTotalTemperature = "
+        + "abs(temperatureTmp - temperatureTotal) > epsilon ? "
+        + "temperatureTotal = - wallHeatFlux / temperatureSafeDivide :  1.0 / epsilon;",
+    },
+    "solution.wall_shear_stress_magnitude": {
+        "declaration": "double wallShearStressMagnitude;",
+        "computation": "wallShearStressMagnitude = magnitude(wallShearStress);",
+    },
+}

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

@@ -94,7 +94,6 @@ def _import_solution(_) -> Any:
             "grad_pressure",
             "Mach",
             "mut",
-            "mu",
             "mut_ratio",
             "nu_hat",
             "turbulence_kinetic_energy",
@@ -115,7 +114,7 @@ def _import_solution(_) -> Any:
             "node_normals",
             "node_forces_per_unit_area",
             "y_plus",
-            "wall_shear_stress",
+            "wall_shear_stress_magnitude",
             "heat_transfer_coefficient_static_temperature",
             "heat_transfer_coefficient_total_temperature",
         ],

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

@@ -647,11 +647,11 @@ def get_output_units(self, input_params=None):
         def get_unit_from_unit_system(expression: Expression, unit_system_name: str):
             """Derive the default output unit based on the value's dimensionality and current unit system"""
             numerical_value = expression.evaluate(raise_on_non_evaluable=False, force_evaluate=True)
-            if not isinstance(numerical_value, (u.unyt_array, u.unyt_quantity, list)):
-                # Pure dimensionless constant
-                return None
             if isinstance(numerical_value, list):
                 numerical_value = numerical_value[0]
+            if not isinstance(numerical_value, (u.unyt_array, u.unyt_quantity)):
+                # Pure dimensionless constant
+                return None
 
             if unit_system_name in ("SI", "SI_unit_system"):
                 return numerical_value.in_base("mks").units