[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++)
         {
@@ -116,6 +119,11 @@ namespace Kratos
         BoundedMatrix<double,TNumNodes, TNumNodes> aux1 = ZeroMatrix(TNumNodes, TNumNodes); //terms multiplying dphi/dt
         BoundedMatrix<double,TNumNodes, TNumNodes> aux2 = ZeroMatrix(TNumNodes, TNumNodes); //terms multiplying phi
         bounded_matrix<double,TNumNodes, TDim> tmp;
+        // CrossWind variables
+        // Projected velocity
+        array_1d<double,TDim> u_proj;
+        // Diffusion contribution
+        BoundedMatrix<double,TDim,TDim> Dcw;
 
         // Gauss points and Number of nodes coincides in this case.
         for(unsigned int igauss=0; igauss<TNumNodes; igauss++)
@@ -141,6 +149,55 @@ 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
+            const double norm_vel2 = norm_vel * norm_vel;
+            if(Variables.C > 0.0 && norm_grad > 1e-3 && norm_vel2 > 1e-9)
+            {
+                // Temporal derivative
+                const double phi_gauss = inner_prod(N, Variables.phi);
+                const double phi_old_gauss = inner_prod(N, Variables.phi_old);
+                const double dphi_dt = Variables.dt_inv * (phi_gauss - phi_old_gauss);
+
+                // Convective term
+                const double conv = inner_prod(vel_gauss, grad_phi_halfstep);
+
+                // Residual
+                const double res = dphi_dt + conv;
+
+                // Projected velocity
+                u_proj = (conv/norm_grad) * (grad_phi_halfstep/norm_grad);
+
+                // Projected Peclet
+                const double Pe_proj = norm_2(u_proj) * h / (2.0 * Variables.conductivity + 1e-12);
+
+                // Limiter
+                const double alpha_c = std::max( 0.0, Variables.C - 1.0/(Pe_proj + 1e-12));
+
+                // Discontinuity capturing coefficient
+                bool force_diffusion = true; // to avoid k_c = 0.0 when alpha_c = 0.0 using Pe_proj, add to ProjectParameters?
+                double k_c;
+                if (alpha_c == 0.0 && force_diffusion)
+                {
+                    k_c = Variables.C * h * std::abs(res / norm_grad);
+                } else
+                {
+                    k_c = 0.5 * alpha_c * h * std::abs(res / norm_grad);
+                }
+
+                // Crosswind diffusion tensor
+                Dcw = k_c * IdentityMatrix(TDim);
+
+                // Removing diffusion in streamline direcction
+                const double k_streamline = tau * norm_vel2;
+                const double correction = std::max(k_c - k_streamline, 0.0) - k_c;
+
+                Dcw += (correction / norm_vel2) * outer_prod(vel_gauss, vel_gauss);
+
+                // Add diffusion contribution
+                noalias(tmp) = prod(DN_DX, Dcw);
+                noalias(aux2) += prod(tmp, trans(DN_DX));
+            }
         }
 
         //adding the second and third term in the formulation
@@ -174,6 +231,7 @@ namespace Kratos
     {
         KRATOS_TRY
 
+        rVariables.C = rCurrentProcessInfo[CROSS_WIND_STABILIZATION_FACTOR];
         rVariables.theta = rCurrentProcessInfo[TIME_INTEGRATION_THETA]; //Variable defining the temporal scheme (0: Forward Euler, 1: Backward Euler, 0.5: Crank-Nicolson)
         rVariables.dyn_st_beta = rCurrentProcessInfo[DYNAMIC_TAU];
         const double delta_t = rCurrentProcessInfo[DELTA_TIME];

applications/ConvectionDiffusionApplication/custom_elements/eulerian_conv_diff.h

@@ -119,6 +119,7 @@ class KRATOS_API(CONVECTION_DIFFUSION_APPLICATION) EulerianConvectionDiffusionEl
         double density;
         double beta;
         double div_v;
+        double C;
 
         array_1d<double,TNumNodes> phi;
         array_1d<double,TNumNodes> phi_old;

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.