Adding new hydraulic utilities #12280

uxuech · 2024-04-15T14:23:03Z

📝 Description
In this PR, new hydraulic utilities have been added, each thoroughly tested with corresponding tests. Below are the main changes:
🆕 Changelog

Corner detection to prevent velocity peaks.
Enabling the option to turn off gravity in air elements.
Hydraulic inlet under subcritical conditions (possibility of having a submerged inlet).
Assigning artificial viscosity to uncut elements (this has been moved from the hydraulic solver in order to optimize simulation times)."

.../FluidDynamicsHydraulicsApplication/custom_utilities/hydraulic_fluid_auxiliary_utilities.cpp

rubenzorrilla · 2024-04-16T05:44:29Z

.../FluidDynamicsHydraulicsApplication/custom_utilities/hydraulic_fluid_auxiliary_utilities.cpp

+    });
+}
+
+void HydraulicFluidAuxiliaryUtilities::TurnOffGravityOnAirElements(ModelPart &rModelPart){


How do you turn it on when an element changes back from air to water?

rubenzorrilla · 2024-04-16T05:45:05Z

.../FluidDynamicsHydraulicsApplication/custom_utilities/hydraulic_fluid_auxiliary_utilities.cpp

+            }
+        if (FluidAuxiliaryUtilities::IsPositive(distances)){
+            for (unsigned int i_nodes = 0; i_nodes < r_geom.PointsNumber(); i_nodes++)
+            {


Please be consistent in the formatting, you either put the { in the same line or in a new one.

.../FluidDynamicsHydraulicsApplication/custom_utilities/hydraulic_fluid_auxiliary_utilities.cpp

rubenzorrilla · 2024-04-16T05:47:16Z

.../FluidDynamicsHydraulicsApplication/custom_utilities/hydraulic_fluid_auxiliary_utilities.cpp

+void HydraulicFluidAuxiliaryUtilities::FixCornerNodeVelocity(ModelPart &rModelPart, double MaximumAngle)
+{
+    // Obtain for each condition each neighbor condition
+    block_for_each(rModelPart.Nodes(), [&](NodeType &rNode){ rNode.SetValue(AUX_INDEX, 0); });


I think there's a method in the variable utils for doing this.

.../FluidDynamicsHydraulicsApplication/custom_utilities/hydraulic_fluid_auxiliary_utilities.cpp

rubenzorrilla · 2024-04-16T05:49:01Z

.../FluidDynamicsHydraulicsApplication/custom_utilities/hydraulic_fluid_auxiliary_utilities.cpp

+    double acceptable_cos = cos(angle_corner_rad);
+    double max_height = block_for_each<MaxReduction<double>>(rModelPart.Nodes(), [&](NodeType &rNode){return rNode.Z();});


I think both can be const.

.../FluidDynamicsHydraulicsApplication/custom_utilities/hydraulic_fluid_auxiliary_utilities.cpp

rubenzorrilla · 2024-04-16T05:52:56Z

.../FluidDynamicsHydraulicsApplication/custom_utilities/hydraulic_fluid_auxiliary_utilities.cpp

+            array_1d<double, 2> ids_neig;
+            array_1d<double, 2> ids_cond;


Suggested change

array_1d<double, 2> ids_neig;

array_1d<double, 2> ids_cond;

array_1d<double, 2> ids_neig;

array_1d<double, 2> ids_cond;

Allocate these two outside the loop.

rubenzorrilla · 2024-04-16T05:54:03Z

.../FluidDynamicsHydraulicsApplication/custom_utilities/hydraulic_fluid_auxiliary_utilities.cpp

+                                double repeated = edgelist[edge][i].GetValue(AUX_INDEX);
+                                repeated += 1;
+                                edgelist[edge][i].SetValue(AUX_INDEX, repeated);


You can avoid the get/set by directly using the get with += 1.

rubenzorrilla · 2024-04-16T05:54:32Z

.../FluidDynamicsHydraulicsApplication/custom_utilities/hydraulic_fluid_auxiliary_utilities.cpp

        }
+    }
+    block_for_each(rModelPart.Nodes(),[&](NodeType &rNode){
+        if (rNode.GetValue(AUX_INDEX)*0.5>=3 && rNode.Is(MARKER)){


Suggested change

if (rNode.GetValue(AUX_INDEX)*0.5>=3 && rNode.Is(MARKER)){

if (rNode.GetValue(AUX_INDEX)*0.5>=3 && rNode.Is(MARKER)){

Why this?

rubenzorrilla · 2024-04-16T05:54:47Z

.../FluidDynamicsHydraulicsApplication/custom_utilities/hydraulic_fluid_auxiliary_utilities.cpp

+    }
+    block_for_each(rModelPart.Nodes(),[&](NodeType &rNode){
+        if (rNode.GetValue(AUX_INDEX)*0.5>=3 && rNode.Is(MARKER)){
+            Vector veloctiy = ZeroVector(3);


Allocate the vector outside the loop.

rubenzorrilla · 2024-04-16T05:55:13Z

.../FluidDynamicsHydraulicsApplication/custom_utilities/hydraulic_fluid_auxiliary_utilities.cpp

+            double value= rNode.GetValue(AUX_INDEX);
+            rNode.SetValue(AUX_INDEX,value*0.5);


Same comment, you can avoid the get/set.

rubenzorrilla · 2024-04-16T05:55:49Z

.../FluidDynamicsHydraulicsApplication/custom_utilities/hydraulic_fluid_auxiliary_utilities.cpp

+
+bool HydraulicFluidAuxiliaryUtilities::MaximumWaterDepthChange(const ModelPart &rModelPart){
+    // make sure the Ids of the Nodes start with one
+    const double &min_Z = block_for_each<MinReduction<double>>(rModelPart.Nodes(), [&](const auto &rNode){ return rNode.Z(); });


Why taking a reference in here?

rubenzorrilla · 2024-04-16T05:57:13Z

.../FluidDynamicsHydraulicsApplication/custom_utilities/hydraulic_fluid_auxiliary_utilities.cpp

+        if (std::abs(rNode.Z() - min_Z) < 1.0e-10 && !maximum_water_depth_change){
+            const double inlet_phi = rNode.GetValue(AUX_DISTANCE);
+            const double domain_phi = rNode.FastGetSolutionStepValue(DISTANCE);
+            if (std::abs(domain_phi) > std::abs(inlet_phi)){


Are you sure you want to check this? Shouldn't you check that the absolute value of the difference is below a threshold?

rubenzorrilla · 2024-04-16T05:57:31Z

.../FluidDynamicsHydraulicsApplication/custom_utilities/hydraulic_fluid_auxiliary_utilities.cpp

+    return maximum_water_depth_change;
+}
+
+void HydraulicFluidAuxiliaryUtilities::CalculateArtificialViscosity(ModelPart &rModelPart, double LimiterCoefficient)


Please follow the style guide for the arguments.

Besides, I just realized that MaximumArtificialViscosity would be a much convenient variable name for this magnitude. As far as I understand, this is not a limiting coefficient but a viscosity.

rubenzorrilla · 2024-04-16T05:58:44Z

.../FluidDynamicsHydraulicsApplication/custom_utilities/hydraulic_fluid_auxiliary_utilities.cpp

+void HydraulicFluidAuxiliaryUtilities::CalculateArtificialViscosity(ModelPart &rModelPart, double LimiterCoefficient)
+{
+    const auto &r_process_info = rModelPart.GetProcessInfo();
+    block_for_each(rModelPart.Elements(), [&](Element &rElement){


Can you please respect indentation in this loop? It would ease quite a lot going through the code inside the block_for_each...

.../FluidDynamicsHydraulicsApplication/custom_utilities/hydraulic_fluid_auxiliary_utilities.cpp

rubenzorrilla · 2024-04-16T06:04:32Z

...ns/FluidDynamicsHydraulicsApplication/custom_utilities/hydraulic_fluid_auxiliary_utilities.h

+    struct NodeConditionDataStruct
+    {
+        NodeType Node;                 // Condition normal
+        double CornerNode; // Gauss point shape functions values
+    };


I think this is used nowhere...

...ns/FluidDynamicsHydraulicsApplication/custom_utilities/hydraulic_fluid_auxiliary_utilities.h

rubenzorrilla · 2024-04-16T06:07:08Z

...ications/FluidDynamicsHydraulicsApplication/custom_python/add_custom_utilities_to_python.cpp

@@ -26,12 +26,26 @@ void AddCustomUtilitiesToPython(pybind11::module& m)
 {


Please, avoid all the non-necessary changes in the format as these complicate very much differentiating the new additions.

...cations/FluidDynamicsHydraulicsApplication/tests/test_hydraulic_fluid_auxiliary_utilities.py

rubenzorrilla · 2024-04-16T06:10:44Z

...cations/FluidDynamicsHydraulicsApplication/tests/test_hydraulic_fluid_auxiliary_utilities.py

+        Kratos.VariableUtils().AddDof(Kratos.VELOCITY_Z,self.model_part)
+        KratosFluidHydraulics.HydraulicFluidAuxiliaryUtilities.FixCornerNodeVelocity(self.model_part , corner_angle)
+        tested_node = self.model_part.GetNode(33)
+        self.assertAlmostEqual(tested_node.GetSolutionStepValue(Kratos.BODY_FORCE_Z),0.0,12)


What are you checking in here? Shouldn't you check the fixity status?

rubenzorrilla

Besides some things that I consider tiny errors, there is still room for improvement.

…es/hydraulic_fluid_auxiliary_utilities.cpp Co-authored-by: Rubén Zorrilla <rubenzorrillamartinez@hotmail.com>

…es/hydraulic_fluid_auxiliary_utilities.h Co-authored-by: Rubén Zorrilla <rubenzorrillamartinez@hotmail.com>

…raulic_fluid_auxiliary_utilities.py Co-authored-by: Rubén Zorrilla <rubenzorrillamartinez@hotmail.com>

…es/hydraulic_fluid_auxiliary_utilities.cpp Co-authored-by: Rubén Zorrilla <rubenzorrillamartinez@hotmail.com>

…es/hydraulic_fluid_auxiliary_utilities.h Co-authored-by: Rubén Zorrilla <rubenzorrillamartinez@hotmail.com>

…es/hydraulic_fluid_auxiliary_utilities.cpp Co-authored-by: Rubén Zorrilla <rubenzorrillamartinez@hotmail.com>

rubenzorrilla · 2024-04-24T09:02:44Z

@uxuech any advance in this regard?

adding new hydraulic utilities

a075d84

uxuech requested a review from rubenzorrilla April 15, 2024 14:23

uxuech enabled auto-merge April 15, 2024 14:23

adding new utilities

6007772