[ConvectionDiffusionApplication] Adding crosswind stabilization to Eulerian Convection–Diffusion element

applications/ConvectionDiffusionApplication/custom_elements/eulerian_conv_diff.cpp

@@ -103,6 +103,9 @@ namespace Kratos
         this-> GetNodalValues(Variables,rCurrentProcessInfo);
         double h = this->ComputeH(DN_DX);
 
+        array_1d<double,TDim> grad_phi_halfstep = prod(trans(DN_DX), 0.5*(Variables.phi+Variables.phi_old));
+        const double norm_grad = norm_2(grad_phi_halfstep);
+
         //Computing the divergence
         for (unsigned int i = 0; i < TNumNodes; i++)
         {
@@ -141,6 +144,22 @@ namespace Kratos
             //terms which multiply the gradient of phi
             noalias(aux2) += (1.0+tau*Variables.beta*Variables.div_v)*outer_prod(N, a_dot_grad);
             noalias(aux2) += tau*outer_prod(a_dot_grad, a_dot_grad);
+
+            //cross-wind term
+            if(norm_grad > 1e-3 && norm_vel > 1e-9)
+            {
+                const double C = rCurrentProcessInfo[CROSS_WIND_STABILIZATION_FACTOR];
+                const double time_derivative = Variables.dt_inv*(inner_prod(N,Variables.phi)-inner_prod(N,Variables.phi_old));
+                const double res = -time_derivative -inner_prod(vel_gauss, grad_phi_halfstep);
+
+                const double disc_capturing_coeff = 0.5*C*h*fabs(res/norm_grad);
+                BoundedMatrix<double,TDim,TDim> D = disc_capturing_coeff*( IdentityMatrix(TDim));
+                const double norm_vel_squared = norm_vel*norm_vel;
+                D += (std::max( disc_capturing_coeff - tau*norm_vel_squared , 0.0) - disc_capturing_coeff)/(norm_vel_squared) * outer_prod(vel_gauss,vel_gauss);
+
+                noalias(tmp) = prod(DN_DX,D);
+                noalias(aux2) += prod(tmp,trans(DN_DX));
+            }
         }
 
         //adding the second and third term in the formulation

applications/ConvectionDiffusionApplication/python_scripts/convection_diffusion_transient_solver.py

@@ -39,7 +39,8 @@ def GetDefaultParameters(cls):
             "time_integration_method" : "implicit",
             "transient_parameters" : {
                 "dynamic_tau": 1.0,
-                "theta"    : 0.5
+                "theta"    : 0.5,
+                "cross_wind_stabilization_factor" : 0.0
             }
         }""")
         this_defaults.AddMissingParameters(super().GetDefaultParameters())
@@ -50,6 +51,7 @@ def _CreateScheme(self):
         # Variable defining the temporal scheme (0: Forward Euler, 1: Backward Euler, 0.5: Crank-Nicolson)
         self.GetComputingModelPart().ProcessInfo[KratosMultiphysics.TIME_INTEGRATION_THETA] = self.settings["transient_parameters"]["theta"].GetDouble()
         self.GetComputingModelPart().ProcessInfo[KratosMultiphysics.DYNAMIC_TAU] = self.settings["transient_parameters"]["dynamic_tau"].GetDouble()
+        self.GetComputingModelPart().ProcessInfo[KratosMultiphysics.CROSS_WIND_STABILIZATION_FACTOR] = self.settings["transient_parameters"]["cross_wind_stabilization_factor"].GetDouble()
 
         # As the time integration is managed by the element, we set a "fake" scheme to perform the solution update
         if not self.main_model_part.IsDistributed():

applications/ConvectionDiffusionApplication/tests/SourceTermTest/SourceTermTestProjectParameters.json

@@ -45,7 +45,8 @@
         },
         "transient_parameters" : {
             "dynamic_tau": 1.0,
-            "theta"    : 1.0
+            "theta"    : 1.0,
+            "cross_wind_stabilization_factor": 0.0
         }
     },
     "processes"        : {

applications/ConvectionDiffusionApplication/tests/cpp_tests/test_eulerian_conv_diff.cpp

@@ -256,4 +256,60 @@ namespace Kratos::Testing
         KRATOS_EXPECT_MATRIX_NEAR(LHS, expected_LHS, 1.0e-4)
     }
 
+    KRATOS_TEST_CASE_IN_SUITE(EulerianConvDiff2D3NCrossWindStabilization, KratosConvectionDiffusionFastSuite)
+    {
+        // Create the test element
+        Model model;
+        auto &r_test_model_part = model.CreateModelPart("TestModelPart");
+        ConvectionDiffusionTestingUtilities::SetEntityUnitTestModelPart(r_test_model_part);
+
+        // Fill the process info container
+        auto &r_process_info = r_test_model_part.GetProcessInfo();
+        r_process_info.SetValue(CROSS_WIND_STABILIZATION_FACTOR, 0.7);
+        r_process_info.SetValue(TIME_INTEGRATION_THETA, 1.0);
+        r_process_info.SetValue(DELTA_TIME, 0.1);
+        r_process_info.SetValue(DYNAMIC_TAU, 0.001);
+
+        // Element creation
+        r_test_model_part.CreateNewNode(1, 0.0, 0.0, 0.0);
+        r_test_model_part.CreateNewNode(2, 1.0, 0.0, 0.0);
+        r_test_model_part.CreateNewNode(3, 0.0, 1.0, 0.0);
+        std::vector<ModelPart::IndexType> elem_nodes{1, 2, 3};
+        r_test_model_part.CreateNewElement("EulerianConvDiff2D", 1, elem_nodes, r_test_model_part.pGetProperties(0));
+
+        // Set the nodal values
+        for (auto &i_node : r_test_model_part.Nodes()) {
+            i_node.FastGetSolutionStepValue(DENSITY) = 1.0;
+            i_node.FastGetSolutionStepValue(CONDUCTIVITY) = 1.0;
+            i_node.FastGetSolutionStepValue(SPECIFIC_HEAT) = 1.0;
+            array_1d<double,3> aux_vel = ZeroVector(3);
+            aux_vel[0] = i_node.X();
+            aux_vel[1] = i_node.Y();
+            i_node.FastGetSolutionStepValue(VELOCITY) = aux_vel;
+        }
+
+        // Test element
+        auto p_element = r_test_model_part.pGetElement(1);
+        Vector RHS = ZeroVector(3);
+        Matrix LHS = ZeroMatrix(3,3);
+        const auto& rConstProcessInfoRef = r_test_model_part.GetProcessInfo();
+        p_element->CalculateLocalSystem(LHS, RHS, rConstProcessInfoRef);
+
+        std::vector<double> expected_RHS = {0.0, 0.0, 0.0};
+        Matrix expected_LHS(3, 3);
+        expected_LHS(0, 0) = 1.71340;
+        expected_LHS(0, 1) = -0.12313;
+        expected_LHS(0, 2) = -0.12313;
+        expected_LHS(1, 0) = -0.19003;
+        expected_LHS(1, 1) = 1.47086;
+        expected_LHS(1, 2) = 0.48560;
+        expected_LHS(2, 0) = -0.19003;
+        expected_LHS(2, 1) = 0.48560;
+        expected_LHS(2, 2) = 1.47086;
+
+        KRATOS_EXPECT_VECTOR_NEAR(RHS, expected_RHS, 1.0e-4)
+        KRATOS_EXPECT_MATRIX_NEAR(LHS, expected_LHS, 1.0e-4)
+
+    }
+
 } // namespace Kratos::Testing.