Merge pull request #11646 from KratosMultiphysics/fluid/weakly-compre…

…ssible-permeability

[Fluid] Weakly-compressible automatic differentiation enhancements

applications/FluidDynamicsApplication/custom_python/fluid_dynamics_python_application.cpp

@@ -62,6 +62,9 @@ PYBIND11_MODULE(KratosFluidDynamicsApplication,m)
     KRATOS_REGISTER_IN_PYTHON_VARIABLE(m,GAPS);
     KRATOS_REGISTER_IN_PYTHON_VARIABLE(m,DIVERGENCE);
 
+    // Darcy's flow variables
+    KRATOS_REGISTER_IN_PYTHON_VARIABLE(m, RESISTANCE);
+
     // Wall modelling
     KRATOS_REGISTER_IN_PYTHON_VARIABLE(m, SLIP_TANGENTIAL_CORRECTION_SWITCH)
 

applications/FluidDynamicsApplication/custom_utilities/weakly_compressible_navier_stokes_data.h

@@ -12,11 +12,18 @@
 //
 
 
-#if !defined(KRATOS_WEAKLY_COMPRESSIBLE_NAVIER_STOKES_DATA_H)
-#define KRATOS_WEAKLY_COMPRESSIBLE_NAVIER_STOKES_DATA_H
+#pragma once
 
+// System includes
+
+
+// External includes
+
+
+// Project includes
 #include "includes/constitutive_law.h"
 
+// Application includes
 #include "fluid_dynamics_application_variables.h"
 #include "custom_utilities/fluid_element_data.h"
 #include "utilities/element_size_calculator.h"
@@ -61,6 +68,7 @@ NodalScalarData SoundVelocity;
 double DynamicViscosity;
 double DeltaTime;      // Time increment
 double DynamicTau;     // Dynamic tau considered in ASGS stabilization coefficients
+double Resistance;     // Darcy's law resistance (permeability inverse)
 
 double bdf0;
 double bdf1;
@@ -102,10 +110,11 @@ void Initialize(const Element& rElement, const ProcessInfo& rProcessInfo) overri
     bdf1 = BDFVector[1];
     bdf2 = BDFVector[2];
 
-    ElementSize = ElementSizeCalculator<TDim,TNumNodes>::MinimumElementSize(r_geometry);
+    // Get data from element database
+    this->FillFromElementData(Resistance, RESISTANCE, rElement);
 
-    noalias(lhs) = ZeroMatrix(TNumNodes*(TDim+1),TNumNodes*(TDim+1));
-    noalias(rhs) = ZeroVector(TNumNodes*(TDim+1));
+    // Calculate element characteristic size
+    ElementSize = ElementSizeCalculator<TDim,TNumNodes>::MinimumElementSize(r_geometry);
 }
 
 void UpdateGeometryValues(
@@ -141,5 +150,3 @@ static int Check(const Element& rElement, const ProcessInfo& rProcessInfo)
 ///@}
 
 }
-
-#endif

applications/FluidDynamicsApplication/fluid_dynamics_application_variables.cpp

@@ -38,6 +38,9 @@ KRATOS_CREATE_3D_VARIABLE_WITH_COMPONENTS(CONVECTION_VELOCITY)
 KRATOS_CREATE_3D_VARIABLE_WITH_COMPONENTS(CONVECTION_SCALAR_GRADIENT)
 KRATOS_CREATE_3D_VARIABLE_WITH_COMPONENTS(AUXILIAR_VECTOR_VELOCITY)
 
+// Darcy's flow variables
+KRATOS_CREATE_VARIABLE(double, RESISTANCE)
+
 // Wall modelling
 KRATOS_CREATE_VARIABLE(bool, SLIP_TANGENTIAL_CORRECTION_SWITCH)
 

applications/FluidDynamicsApplication/fluid_dynamics_application_variables.h

@@ -48,6 +48,9 @@ KRATOS_DEFINE_APPLICATION_VARIABLE( FLUID_DYNAMICS_APPLICATION, Vector, GAPS )
 KRATOS_DEFINE_APPLICATION_VARIABLE( FLUID_DYNAMICS_APPLICATION, double, DIVERGENCE)
 KRATOS_DEFINE_APPLICATION_VARIABLE( FLUID_DYNAMICS_APPLICATION, double, FS_PRESSURE_GRADIENT_RELAXATION_FACTOR)
 
+// Darcy's flow variables
+KRATOS_DEFINE_APPLICATION_VARIABLE(FLUID_DYNAMICS_APPLICATION, double, RESISTANCE)
+
 // Wall modelling
 KRATOS_DEFINE_APPLICATION_VARIABLE(FLUID_DYNAMICS_APPLICATION, bool, SLIP_TANGENTIAL_CORRECTION_SWITCH)
 
@@ -118,6 +121,7 @@ KRATOS_DEFINE_3D_APPLICATION_VARIABLE_WITH_COMPONENTS( FLUID_DYNAMICS_APPLICATIO
 KRATOS_DEFINE_3D_APPLICATION_VARIABLE_WITH_COMPONENTS( FLUID_DYNAMICS_APPLICATION, DRAG_FORCE_CENTER )
 KRATOS_DEFINE_APPLICATION_VARIABLE( FLUID_DYNAMICS_APPLICATION, double, SMOOTHING_COEFFICIENT )
 KRATOS_DEFINE_APPLICATION_VARIABLE(FLUID_DYNAMICS_APPLICATION,double,SPECTRAL_RADIUS_LIMIT)
+
 // Two-phase flow with surface tension
 KRATOS_DEFINE_APPLICATION_VARIABLE( FLUID_DYNAMICS_APPLICATION, double, SURFACE_TENSION_COEFFICIENT )
 KRATOS_DEFINE_APPLICATION_VARIABLE( FLUID_DYNAMICS_APPLICATION, bool, SURFACE_TENSION )

applications/FluidDynamicsApplication/symbolic_generation/navier_stokes/README.md → applications/FluidDynamicsApplication/python_scripts/symbolic_generation/weakly_compressible_navier_stokes/README.md

@@ -1,18 +1,19 @@
-# Navier-Stokes and Stokes elements automatic differentiation
+# Weakly-compressible Navier-Stokes and Stokes elements automatic differentiation
 
 ## ELEMENT DESCRIPTION:
-Current directory contains the required files for the Automatic Differenctiation (AD) of both the _"symbolic_stokes"_ and _"symbolic_navier_stokes"_ elements of the FluidDynamicsApplication. These elements implement a **Stokes** and a **weakly-compressible Navier-Stokes** formulation for both **2D** and **3D** cases.
+Current directory contains the required files for the Automatic Differenctiation (AD) of both the _"symbolic_stokes"_ and _"weakly_compressible_navier_stokes"_ elements of the FluidDynamicsApplication. These elements implement a **Stokes** and a **weakly-compressible Navier-Stokes** formulation for both **2D** and **3D** cases.
 
 ## SYMBOLIC GENERATOR SETTINGS:
 *  Dimension to compute: This symbolic generator is valid for both **2D** and **3D** cases. The element has been implemented using a template argument for the problem dimension, so it is advised to set the _dim_to_compute_ flag as "_Both_". In this case the generated .cpp file will contain both **2D** and **3D** implementations.
 *  Formulation: As mentioned, the symbolic generator is valid for both Stokes and Navier-Stokes formulations. To switch the formulation set the "_formulation_" variable to either "_Stokes_" or "_NavierStokes_"
 *  Linearisation settings: "_FullNR_" considers the convective velocity as _"v-vmesh"_, implying that _v_ is taken into account in the differenctiation of the **LHS** and **RHS**. On the contrary "_Picard_" option (_a.k.a. QuasiNR_) defines the convective velocity as "a", so it is considered as a constant in the differenctiation of the **LHS** and **RHS**. Note that this option is only meaningful in the Navier-Stokes formulation as there is no convective term in the Stokes case.
 *  Divide by rho: If _divide_by_rho_ flag is set to "_True_" the mass conservation equation is divided by the density for the sake of having a better conditioned system of equations. It is therefore advised to keep it activated.
 *  Artificial compressiblity: If set to "_True_", the time derivative of the density is introduced in the mass conservation equation and rearanged to the pressure time derivative by using the simple state equation $$\partial p/\partial\rho = c^{2}$$, being "_c_" the fluid speed of sound. This adds some extra terms to the usual **Navier-Stokes** equations that are intended to act as a soft artificial compressibility. Such artificial compressibility is controlled by the value of "_c_", meaning that the artificial compressibility terms vanishes as "_c_" tends to infinite. So far this option is only available in the Navier-Stokes formulation. Although the generator is already prepared to include it in the Stokes element, slight modifications would be required in the header and source template files.
+* Darcy term: If set to "_True_" a darcy term is included in the momentum conservation equation. Such Darcy term is governed by a resistance constant (1/L^2) retrieved from the elemental database.
 
 ## INSTRUCTIONS
 Run:
 ~~~py
-python generate_navier_stokes_element.py
+python generate_weakly_compressible_navier_stokes_element.py
 ~~~
 Then, depending on the selected formulation, a file "_symbolic_stokes.cpp_" or "_weakly_compressible_navier_stokes.cpp_" is automatically generated. Such file muest be copied within the "_custom_elements_" folder of the **FluidDynamicsApplication**. The corresponding header file ("symbolic_stokes.h" or "weakly_compressible_navier_stokes.h") is already stored in such folder.

applications/FluidDynamicsApplication/symbolic_generation/navier_stokes/generate_navier_stokes_element.py → applications/FluidDynamicsApplication/python_scripts/symbolic_generation/weakly_compressible_navier_stokes/generate_weakly_compressible_navier_stokes_element.py

@@ -28,12 +28,14 @@
 formulation = "WeaklyCompressibleNavierStokes" # Element type. Options: "WeaklyCompressibleNavierStokes", "Stokes"
 
 if formulation == "WeaklyCompressibleNavierStokes":
+    darcy_term = True
     convective_term = True
     artificial_compressibility = True
     linearisation = "Picard" # Convective term linearisation type. Options: "Picard", "FullNR"
     output_filename = "weakly_compressible_navier_stokes.cpp"
     template_filename = "weakly_compressible_navier_stokes_cpp_template.cpp"
 elif formulation == "Stokes":
+    darcy_term = False
     convective_term = False
     artificial_compressibility = False
     output_filename = "symbolic_stokes.cpp"
@@ -52,14 +54,14 @@
 print(info_msg)
 
 #TODO: DO ALL ELEMENT TYPES FOR N-S TOO
-if formulation == "NavierStokes" or formulation == "WeaklyCompressibleNavierStokes":
-    if (dim_to_compute == "2D"):
+if formulation == "WeaklyCompressibleNavierStokes":
+    if dim_to_compute == "2D":
         dim_vector = [2]
         nnodes_vector = [3] # tria
-    elif (dim_to_compute == "3D"):
+    elif dim_to_compute == "3D":
         dim_vector = [3]
         nnodes_vector = [4] # tet
-    elif (dim_to_compute == "Both"):
+    elif dim_to_compute == "Both":
         dim_vector = [2, 3]
         nnodes_vector = [3, 4] # tria, tet
 elif formulation == "Stokes":
@@ -81,9 +83,9 @@
 
 for dim, nnodes in zip(dim_vector, nnodes_vector):
 
-    if(dim == 2):
+    if dim == 2:
         strain_size = 3
-    elif(dim == 3):
+    elif dim == 3:
         strain_size = 6
 
     impose_partion_of_unity = False
@@ -124,13 +126,17 @@
     stress = DefineVector('stress',strain_size)
 
     ## Other simbols definition
-    dt  = sympy.Symbol('dt', positive = True)         # Time increment
-    nu  = sympy.Symbol('nu', positive = True)         # Kinematic viscosity (mu/rho)
-    mu  = sympy.Symbol('mu', positive = True)         # Dynamic viscosity
+    dt = sympy.Symbol('dt', positive = True)         # Time increment
+    mu = sympy.Symbol('mu', positive = True)         # Dynamic viscosity
     h = sympy.Symbol('h', positive = True)
     dyn_tau = sympy.Symbol('dyn_tau', positive = True)
     stab_c1 = sympy.Symbol('stab_c1', positive = True)
     stab_c2 = sympy.Symbol('stab_c2', positive = True)
+    gauss_weight = sympy.Symbol('gauss_weight', positive = True)
+
+    ## Symbols for Darcy term
+    sigma = sympy.Symbol('sigma', positive = True) if darcy_term else 0.0 # Resistance (permeability inverse, 1/m^2)
+    stab_c3 = sympy.Symbol('stab_c3', positive = True) if darcy_term else 0.0 # Darcy term stabilization constant
 
     ## Backward differences coefficients
     bdf0 = sympy.Symbol('bdf0')
@@ -158,11 +164,11 @@
         for i in range(0, dim):
             stab_norm_a += vconv_gauss[i]**2
         stab_norm_a = sympy.sqrt(stab_norm_a)
-        tau1 = 1.0/((rho*dyn_tau)/dt + (stab_c2*rho*stab_norm_a)/h + (stab_c1*mu)/(h*h)) # Stabilization parameter 1
-        tau2 = mu + (stab_c2*rho*stab_norm_a*h)/stab_c1                                  # Stabilization parameter 2
+        tau1 = 1.0/(rho*dyn_tau/dt + stab_c2*rho*stab_norm_a/h + stab_c1*mu/h**2 + stab_c3*sigma/h**2) # Stabilization parameter 1
+        tau2 = mu + (stab_c2*rho*stab_norm_a*h + stab_c3*sigma)/stab_c1                                # Stabilization parameter 2
     else:
-        tau1 = 1.0/((rho*dyn_tau)/dt + (stab_c1*mu)/(h*h)) # Stabilization parameter 1
-        tau2 = (h*h) / (stab_c1 * tau1)                    # Stabilization parameter 2
+        tau1 = 1.0/(rho*dyn_tau/dt + stab_c1*mu/h**2 + stab_c3*sigma/h**2) # Stabilization parameter 1
+        tau2 = mu + stab_c3*sigma/stab_c1                                  # Stabilization parameter 2
 
     ## Compute the rest of magnitudes at the Gauss points
     accel_gauss = (bdf0*v + bdf1*vn + bdf2*vnn).transpose()*N
@@ -189,7 +195,7 @@
 
     grad_sym_v = grad_sym_voigtform(DN,v)       # Symmetric gradient of v in Voigt notation
     grad_w_voigt = grad_sym_voigtform(DN,w)     # Symmetric gradient of w in Voigt notation
-    # Recall that the grad(w):sigma contraction equals grad_sym(w)*sigma in Voigt notation since sigma is a symmetric tensor.
+    # Recall that the grad(w):stress contraction equals grad_sym(w)*stress in Voigt notation since the stress is a symmetric tensor.
 
     # Convective term definition
     if convective_term:
@@ -198,8 +204,8 @@
 
     ## Compute galerkin functional
     # Navier-Stokes functional
-    if (divide_by_rho):
-        rv_galerkin = rho*w_gauss.transpose()*f_gauss - rho*w_gauss.transpose()*accel_gauss - grad_w_voigt.transpose()*stress + div_w*p_gauss - q_gauss*div_v
+    if divide_by_rho:
+        rv_galerkin = rho*w_gauss.transpose()*f_gauss - rho*w_gauss.transpose()*accel_gauss - grad_w_voigt.transpose()*stress + div_w*p_gauss - sigma*w_gauss.transpose()*v_gauss - q_gauss*div_v
         if artificial_compressibility:
             rv_galerkin -= (1/(rho*c*c))*q_gauss*pder_gauss
             if convective_term:
@@ -218,12 +224,12 @@
     ##  Stabilization functional terms
     # Momentum conservation residual
     # Note that the viscous stress term is dropped since linear elements are used
-    vel_residual = rho*f_gauss - rho*accel_gauss - grad_p
+    vel_residual = rho*f_gauss - rho*accel_gauss - grad_p - sigma*v_gauss
     if convective_term:
         vel_residual -= rho*convective_term_gauss.transpose()
 
     # Mass conservation residual
-    if (divide_by_rho):
+    if divide_by_rho:
         mas_residual = -div_v
         if artificial_compressibility:
             mas_residual -= (1/(rho*c*c))*pder_gauss
@@ -240,17 +246,18 @@
     mas_subscale = tau2*mas_residual
 
     # Compute the ASGS stabilization terms using the momentum and mass conservation residuals above
-    if (divide_by_rho):
+    if divide_by_rho:
         rv_stab = grad_q.transpose()*vel_subscale
     else:
         rv_stab = rho*grad_q.transpose()*vel_subscale
     if convective_term:
         rv_stab += rho*vconv_gauss.transpose()*grad_w*vel_subscale
         rv_stab += rho*div_vconv*w_gauss.transpose()*vel_subscale
+    rv_stab -= sigma*w_gauss.transpose()*vel_subscale
     rv_stab += div_w*mas_subscale
 
     ## Add the stabilization terms to the original residual terms
-    if (ASGS_stabilization):
+    if ASGS_stabilization:
         rv = rv_galerkin + rv_stab
     else:
         rv = rv_galerkin
@@ -259,7 +266,7 @@
     dofs = sympy.zeros(nnodes*(dim+1), 1)
     testfunc = sympy.zeros(nnodes*(dim+1), 1)
 
-    for i in range(0,nnodes):
+    for i in range(nnodes):
 
         # Velocity DOFs and test functions
         for k in range(0,dim):
@@ -273,22 +280,22 @@
     ## Compute LHS and RHS
     # For the RHS computation one wants the residual of the previous iteration (residual based formulation). By this reason the stress is
     # included as a symbolic variable, which is assumed to be passed as an argument from the previous iteration database.
-    print("Computing " + str(dim) + "D RHS Gauss point contribution\n")
+    print(f"Computing {dim}D{nnodes}N RHS Gauss point contribution\n")
     rhs = Compute_RHS(rv.copy(), testfunc, do_simplifications)
-    rhs_out = OutputVector_CollectingFactors(rhs, "rhs", mode)
+    rhs_out = OutputVector_CollectingFactors(gauss_weight*rhs, "rRHS", mode, assignment_op='+=')
 
     # Compute LHS (RHS(residual) differenctiation w.r.t. the DOFs)
     # Note that the 'stress' (symbolic variable) is substituted by 'C*grad_sym_v' for the LHS differenctiation. Otherwise the velocity terms
     # within the velocity symmetryc gradient would not be considered in the differenctiation, meaning that the stress would be considered as
     # a velocity independent constant in the LHS.
-    print("Computing " + str(dim) + "D LHS Gauss point contribution\n")
+    print(f"Computing {dim}D{nnodes}N LHS Gauss point contribution\n")
     SubstituteMatrixValue(rhs, stress, C*grad_sym_v)
     lhs = Compute_LHS(rhs, testfunc, dofs, do_simplifications) # Compute the LHS (considering stress as C*(B*v) to derive w.r.t. v)
-    lhs_out = OutputMatrix_CollectingFactors(lhs, "lhs", mode)
+    lhs_out = OutputMatrix_CollectingFactors(gauss_weight*lhs, "rLHS", mode, assignment_op='+=')
 
     ## Replace the computed RHS and LHS in the template outstring
-    outstring = outstring.replace("//substitute_lhs_" + str(dim) + 'D' + str(nnodes) + 'N', lhs_out)
-    outstring = outstring.replace("//substitute_rhs_" + str(dim) + 'D' + str(nnodes) + 'N', rhs_out)
+    outstring = outstring.replace(f"//substitute_lhs_{dim}D{nnodes}N", lhs_out)
+    outstring = outstring.replace(f"//substitute_rhs_{dim}D{nnodes}N", rhs_out)
 
 ## Write the modified template
 print("Writing output file \'" + output_filename + "\'")