tmp

CMakeLists.txt

@@ -108,6 +108,16 @@ if(${stratimikos_package_index} GREATER_EQUAL 0)
   endif()
 endif()
 
+list(FIND Trilinos_PACKAGE_LIST Sacado sacado_package_index)
+if(${sacado_package_index} GREATER_EQUAL 0)
+  option(ENABLE_SACADO "ENABLE_SACADO" ON)
+  if(ENABLE_SACADO)
+    message(STATUS "TRILINOS: Sacado found, enabling in Goma")
+    list(APPEND GOMA_COMPILE_DEFINITIONS GOMA_ENABLE_SACADO)
+  endif()
+endif()
+
+
 # Exodus
 find_package(SEACASExodus REQUIRED HINTS ${Trilinos_DIR}/../SEACASExodus)
 message(STATUS "SEACASExodus_DIR = ${SEACASExodus_DIR}")
@@ -416,6 +426,7 @@ set(GOMA_INCLUDES
     include/ac_stability.h
     include/ac_stability_util.h
     include/ac_update_parameter.h
+    include/ad_turbulence.h
     include/bc_colloc.h
     include/bc_contact.h
     include/bc_curve.h
@@ -571,6 +582,7 @@ set(GOMA_SOURCES
     src/ac_stability.c
     src/ac_stability_util.c
     src/ac_update_parameter.c
+    src/ad_turbulence.cpp
     src/bc_colloc.c
     src/bc_contact.c
     src/bc_curve.c

include/mm_as_structs.h

@@ -570,6 +570,8 @@ struct Element_Variable_Pointers {
   dbl *sh_sat_3[MDE];   /* Porous shell saturation layer 3 */
 
   dbl *eddy_nu[MDE]; /* Eddy viscosity for turbulent flow */
+  dbl *turb_k[MDE]; /* Eddy viscosity for turbulent flow */
+  dbl *turb_omega[MDE]; /* Eddy viscosity for turbulent flow */
 };
 
 /*___________________________________________________________________________*/
@@ -697,6 +699,8 @@ struct Element_Stiffness_Pointers {
   dbl **sh_sat_3; /* Porous shell saturation layer 3 */
 
   dbl **eddy_nu; /* Eddy viscosity for turbulent flow */
+  dbl **turb_k; /* Eddy viscosity for turbulent flow */
+  dbl **turb_omega; /* Eddy viscosity for turbulent flow */
 
   /*
    * These are for debugging purposes...
@@ -1730,6 +1734,8 @@ struct Field_Variables {
   dbl sh_sat_3;   /* Porous shell saturation layer 3 */
 
   dbl eddy_nu;       /* Eddy viscosity for turbulent flow */
+  dbl turb_k;       /* Eddy viscosity for turbulent flow */
+  dbl turb_omega;       /* Eddy viscosity for turbulent flow */
   dbl wall_distance; /* Distance to nearest wall */
 
   /*
@@ -1780,6 +1786,8 @@ struct Field_Variables {
   dbl grad_sh_sat_3[DIM];  /* Gradient of porous shell saturation layer 3 */
 
   dbl grad_eddy_nu[DIM];       /* Gradient of Eddy viscosity */
+  dbl grad_turb_k[DIM];       /* Gradient of Eddy viscosity */
+  dbl grad_turb_omega[DIM];       /* Gradient of Eddy viscosity */
   dbl grad_wall_distance[DIM]; /* Distance to nearest wall */
 
   /*
@@ -1965,6 +1973,8 @@ struct Field_Variables {
   dbl d_max_strain_dmesh[DIM][MDE];
   dbl d_cur_strain_dmesh[DIM][MDE];
   dbl d_grad_eddy_nu_dmesh[DIM][DIM][MDE];
+  dbl d_grad_turb_k_dmesh[DIM][DIM][MDE];
+  dbl d_grad_turb_omega_dmesh[DIM][DIM][MDE];
   dbl d_grad_restime_dmesh[DIM][DIM][MDE];
   /*
    * Values at surfaces for integrated boundary conditions
@@ -2118,6 +2128,8 @@ struct Diet_Field_Variables {
   dbl sh_sat_3;   /* Porous shell saturation layer 3 */
 
   dbl eddy_nu; /* Eddy viscosity for turbulent flow */
+  dbl turb_k; /* Eddy viscosity for turbulent flow */
+  dbl turb_omega; /* Eddy viscosity for turbulent flow */
 
   dbl grad_em_er[DIM][DIM]; /* EM wave Fields */
   dbl grad_em_ei[DIM][DIM]; /* EM wave Fields */
@@ -2146,6 +2158,7 @@ struct Diet_Field_Variables {
   dbl grad_tfmp_sat[DIM];  /* Gradient of the thin-film multi-phase lubrication saturation */
 
   dbl grad_n[DIM][DIM]; /* Normal to level set function OR shell normal */
+  dbl grad_turb_omega[DIM];
   dbl div_n;            /* Divergence of LS normal field */
 
   /* Material tensors used at old time values */
@@ -3028,6 +3041,8 @@ struct stress_dependence {
   double pf[DIM][DIM][MAX_PHASE_FUNC][MDE];
   double degrade[DIM][DIM][MDE];
   double eddy_nu[DIM][DIM][MDE];
+  double turb_k[DIM][DIM][MDE];
+  double turb_omega[DIM][DIM][MDE];
 };
 typedef struct stress_dependence STRESS_DEPENDENCE_STRUCT;
 
@@ -3077,6 +3092,8 @@ struct viscosity_dependence {
   double pf[MAX_PHASE_FUNC][MDE]; /* phase function */
   double degrade[MDE];            /* amount of degradation */
   double eddy_nu[MDE];            /* Turbulent viscosity */
+  double turb_k[MDE];            /* Turbulent k */
+  double turb_omega[MDE];            /* Turbulent omega */
 };
 typedef struct viscosity_dependence VISCOSITY_DEPENDENCE_STRUCT;
 
@@ -3138,6 +3155,8 @@ struct pspg_dependence {
   double v[DIM][DIM][MDE]; /* velocity dependence. */
   double T[DIM][MDE];      /* temperature dependence. */
   double eddy_nu[DIM][MDE];
+  double turb_k[DIM][MDE];
+  double turb_omega[DIM][MDE];
   double P[DIM][MDE];           /* pressure dependence. */
   double C[DIM][MAX_CONC][MDE]; /* conc dependence. */
   double X[DIM][DIM][MDE];      /* mesh dependence. */

include/mm_fill_stabilization.h

@@ -21,6 +21,8 @@ typedef struct {
   dbl d_tau_dF[MDE];           /* FILL dependence. */
   dbl d_tau_dnn[MDE];          /* bond concentration dependence */
   dbl d_tau_dEDDY_NU[MDE];     /* Turbulent viscosity */
+  dbl d_tau_dturb_k[MDE];     /* Turbulent viscosity */
+  dbl d_tau_dturb_omega[MDE];     /* Turbulent viscosity */
 } momentum_tau_terms;
 
 void supg_tau_shakib(SUPG_terms *supg_terms, int dim, dbl dt, dbl diffusivity, int interp_eqn);
@@ -83,4 +85,5 @@ int calc_cont_gls /* mm_fill_terms.c                           */
      dbl, /* current time                              */
      const PG_DATA *);
 
+bool is_evss_f_model(int model);
 #endif // GOMA_MM_FILL_STABILIZATION_H

include/mm_mp_const.h

@@ -335,6 +335,7 @@ extern int Num_Var_Init_Mat[MAX_NUMBER_MATLS]; /* number of variables to overwri
 /* Turbulent viscosity models for Reynolds Averaged NS */
 #define TURBULENT_SA         52 /* Spallart Allmaras */
 #define TURBULENT_SA_DYNAMIC 53 /* Spallart Allmaras */
+#define TURBULENT_K_OMEGA    54 
 
 /*
  *  Heat source modeling

include/mm_names.h

@@ -477,6 +477,16 @@ struct BC_descriptions BC_Desc[] = {
      {1, 1, 1, 1, 1, 1, 1, 1, 1, 1},
      SINGLE_PHASE,
      DVI_SINGLE_PHASE_DB},
+    {"SA_WALL_FUNC",
+     "SA_WALL_FUNC_BC",
+     COLLOCATE_SURF,
+     SA_WALL_FUNC_BC,
+     R_EDDY_NU,
+     SCALAR,
+     NO_ROT,
+     {1, 1, 1, 1, 1, 1, 1, 1, 1, 1},
+     SINGLE_PHASE,
+     DVI_SINGLE_PHASE_DB},
     {"P",
      "P_BC",
      DIRICHLET,
@@ -7228,6 +7238,36 @@ struct BC_descriptions BC_Desc[] = {
       0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1},
      SINGLE_PHASE,
      DVI_SINGLE_PHASE_DB},
+    {"TURB_K",
+     "TURB_K_BC",
+     DIRICHLET,
+     TURB_K_BC,
+     R_TURB_K,
+     SCALAR,
+     NO_ROT,
+     {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+      0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+      0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+      0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+      0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+      0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1},
+     SINGLE_PHASE,
+     DVI_SINGLE_PHASE_DB},
+    {"TURB_OMEGA",
+     "TURB_OMEGA_BC",
+     DIRICHLET,
+     TURB_OMEGA_BC,
+     R_TURB_OMEGA,
+     SCALAR,
+     NO_ROT,
+     {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+      0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+      0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+      0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+      0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+      0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1},
+     SINGLE_PHASE,
+     DVI_SINGLE_PHASE_DB},
 
 };
 
@@ -7501,6 +7541,8 @@ struct Equation_Names EQ_Name[] = {
     {"R_VSTAR", "VSTAR", R_VSTAR},
     {"R_WSTAR", "WSTAR", R_WSTAR},
     {"R_EDDY_NU", "EDDY_NU", R_EDDY_NU},
+    {"R_TURB_K", "TURB_K", R_TURB_K},
+    {"R_TURB_OMEGA", "TURB_OMEGA", R_TURB_OMEGA},
 
     /*
      *  Note -> these entries must remain until we get rid
@@ -7801,9 +7843,11 @@ struct Equation_Names Var_Name[] = {
     {"VSTAR", "USY", VSTAR},
     {"WSTAR", "USZ", WSTAR},
     {"EDDY_NU", "EDDY_NU", EDDY_NU}, // 214
+    {"TURB_K", "TURB_K", TURB_K}, // 215
+    {"TURB_OMEGA", "TURB_OMEGA", TURB_OMEGA}, // 215
 
     {"MESH_POSITION1", "X", MESH_POSITION1},
-    {"MESH_POSITION2", "Y", MESH_POSITION2}, /* 216 */
+    {"MESH_POSITION2", "Y", MESH_POSITION2}, /* 218 */
     {"MESH_POSITION3", "Z", MESH_POSITION3},
 
     {"VEL_NORM", "VN", VEL_NORM},
@@ -7817,14 +7861,14 @@ struct Equation_Names Var_Name[] = {
 
     {"D_X1_DT", "XDOT", D_X1_DT},
     {"D_X2_DT", "YDOT", D_X2_DT},
-    {"D_X3_DT", "ZDOT", D_X3_DT}, /* 227 */
+    {"D_X3_DT", "ZDOT", D_X3_DT}, /* 229 */
     {"D_S_DT", "SDOT", D_S_DT},
 
     {"D_P_DT", "PDOT", D_P_DT},
 
     {"SOLID_POSITION1", "X_RS", SOLID_POSITION1},
     {"SOLID_POSITION2", "Y_RS", SOLID_POSITION2},
-    {"SOLID_POSITION3", "Z_RS", SOLID_POSITION3} /* 232 */
+    {"SOLID_POSITION3", "Z_RS", SOLID_POSITION3} /* 234 */
 };
 
 int Num_Var_Names = sizeof(Var_Name) / sizeof(struct Equation_Names);
@@ -8072,6 +8116,8 @@ struct Equation_Names Exo_Var_Names[] = {
     {"V Int.", "USY", VSTAR},
     {"W Int.", "USZ", WSTAR},
     {"Eddy Turbulence Viscosity.", "EDDY_NU", EDDY_NU},
+    {"Turbulent K", "TURB_K", TURB_K},
+    {"Turbulent OMEGA", "TURB_OMEGA", TURB_OMEGA},
 };
 
 int Num_Exo_Var_Names = sizeof(Exo_Var_Names) / sizeof(struct Equation_Names);
@@ -8379,6 +8425,8 @@ struct Equation_Names Var_Units[] = {
     {"VSTAR", "[1]", VSTAR},
     {"WSTAR", "[1]", WSTAR},
     {"EDDY_NU", "[1]", EDDY_NU},
+    {"TURB_K", "[1]", TURB_K},
+    {"TURB_OMEGA", "[1]", TURB_OMEGA},
 };
 
 int Num_Var_Units = sizeof(Var_Units) / sizeof(struct Equation_Names);

include/rf_bc_const.h

@@ -800,6 +800,7 @@
 #define DOUBLE_RAD_BC              961123500
 #define FEATURE_ROLLON_BC          961223400
 #define ROLL_FLUID_BC              961124500
+#define SA_WALL_FUNC_BC              961124501
 #define TENSION_SHEET_BC           96210200
 #define MOVING_PLANE_BC            96110000
 #define MOVING_PLANE_ETCH_BC       96115000
@@ -903,6 +904,8 @@
 #define VELO_SLIP_ROT_FLUID_BC 964910000
 
 #define EDDY_NU_BC 966666666
+#define TURB_K_BC 966666667
+#define TURB_OMEGA_BC 966666668
 
 /* Structural Shells */
 #define SH_K_BC             970000000

include/rf_fem_const.h

@@ -519,6 +519,8 @@
 #define VSTAR              212
 #define WSTAR              213
 #define EDDY_NU            214
+#define TURB_K            215
+#define TURB_OMEGA            216
 /*
  * define a variable to hold an external field which will be
  * held fixed in the problem but parametered by the basis functions
@@ -925,7 +927,9 @@
 #define R_VSTAR        212
 #define R_WSTAR        213
 #define R_EDDY_NU      214
-#define V_LAST         215
+#define R_TURB_K            215
+#define R_TURB_OMEGA            216
+#define V_LAST         217
 
 /* MMH
  * This is used for those parts of the code that want to ensure

src/bc_colloc.c

@@ -26,6 +26,7 @@
 #include <stdio.h>
 #include <string.h>
 
+#include "ad_turbulence.h"
 #include "ac_stability.h"
 #include "ac_stability_util.h"
 #include "bc/rotate_coordinates.h"
@@ -182,6 +183,8 @@ int apply_point_colloc_bc(double resid_vector[], /* Residual vector for the curr
       err = load_bf_mesh_derivs();
       GOMA_EH(err, "load_bf_mesh_derivs");
 
+      fill_ad_field_variables();
+
       /* calculate the shape functions and their gradients */
 
       /* calculate the determinant of the surface jacobian  and the normal to
@@ -363,6 +366,13 @@ int apply_point_colloc_bc(double resid_vector[], /* Residual vector for the curr
                               BC_Types[bc_input_id].BC_Data_Float[2], mp_2);
             } break;
 
+            case SA_WALL_FUNC_BC:
+              memset(kfunc, 0, DIM * sizeof(double));
+              ad_sa_wall_func(kfunc, d_kfunc);
+              func = kfunc[0];
+              d_func[EDDY_NU] = 1.0;
+              break;
+
             case PLANEX_BC:
             case PLANEY_BC:
             case PLANEZ_BC:
@@ -2587,6 +2597,16 @@ int load_variable(double *x_var,       /* variable value */
     var = EDDY_NU;
     *d_x_var = 1.;
     break;
+  case TURB_K:
+    *x_var = fv->turb_k;
+    var = TURB_K;
+    *d_x_var = 1.;
+    break;
+  case TURB_OMEGA:
+    *x_var = fv->turb_omega;
+    var = TURB_OMEGA;
+    *d_x_var = 1.;
+    break;
   case LIGHT_INTP:
     *x_var = fv->poynt[0];
     var = LIGHT_INTP;

src/bc_contact.c

@@ -30,6 +30,7 @@
 
 /* GOMA include files */
 
+#include "ad_turbulence.h"
 #include "ac_stability.h"
 #include "ac_stability_util.h"
 #include "bc_colloc.h"
@@ -2484,6 +2485,7 @@ void setup_shop_at_point(int ielem, double *xi, const Exo_DB *exo)
   err = load_fv_grads();
   GOMA_EH(err, "load_fv_grads");
 
+  fill_ad_field_variables();
   /*
    *  Just as in the main element assembly, we ensure that the current element
    *  actually has mesh equations associated with it before calculation