From 5fb281773ceedfb5984173b802169363b7072973 Mon Sep 17 00:00:00 2001 From: Weston Ortiz Date: Thu, 25 Jan 2024 11:55:20 -0700 Subject: [PATCH] tmp --- CMakeLists.txt | 12 ++ include/mm_as_structs.h | 19 +++ include/mm_fill_stabilization.h | 3 + include/mm_mp_const.h | 1 + include/mm_names.h | 54 ++++++- include/rf_bc_const.h | 3 + include/rf_fem_const.h | 6 +- src/bc_colloc.c | 20 +++ src/bc_contact.c | 2 + src/load_field_variables.c | 168 +++++++++++++++++++++- src/main.c | 6 +- src/mm_as_alloc.c | 6 + src/mm_fill.c | 65 +++++++-- src/mm_fill_momentum.c | 241 +++++++++++++++++++++++++++++--- src/mm_fill_ptrs.c | 21 +++ src/mm_fill_stabilization.c | 107 +++++++++++++- src/mm_flux.c | 18 +++ src/mm_input.c | 10 ++ src/mm_input_bc.c | 3 + src/mm_input_mp.c | 6 +- src/mm_input_util.c | 4 + src/mm_prob_def.c | 2 +- src/mm_unknown_map.c | 72 ++++++++++ src/mm_viscosity.c | 19 ++- 24 files changed, 826 insertions(+), 42 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 750a0ed46..51bf8054b 100644 --- a/CMakeLists.txt +++ b/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 diff --git a/include/mm_as_structs.h b/include/mm_as_structs.h index 2f7ec6f72..c6959e0b6 100644 --- a/include/mm_as_structs.h +++ b/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. */ diff --git a/include/mm_fill_stabilization.h b/include/mm_fill_stabilization.h index 8bfb6b873..497f047cd 100644 --- a/include/mm_fill_stabilization.h +++ b/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 diff --git a/include/mm_mp_const.h b/include/mm_mp_const.h index 6a6bf2ec6..6a41b282d 100644 --- a/include/mm_mp_const.h +++ b/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 diff --git a/include/mm_names.h b/include/mm_names.h index 0eb76df28..8c6695d9f 100644 --- a/include/mm_names.h +++ b/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); diff --git a/include/rf_bc_const.h b/include/rf_bc_const.h index f67a2ec15..37e108826 100644 --- a/include/rf_bc_const.h +++ b/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 diff --git a/include/rf_fem_const.h b/include/rf_fem_const.h index 46cdcd5fe..8f3eebbe9 100644 --- a/include/rf_fem_const.h +++ b/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 diff --git a/src/bc_colloc.c b/src/bc_colloc.c index 26bc6d24c..8670846ba 100644 --- a/src/bc_colloc.c +++ b/src/bc_colloc.c @@ -26,6 +26,7 @@ #include #include +#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; diff --git a/src/bc_contact.c b/src/bc_contact.c index 3251c03f2..56d2059d7 100644 --- a/src/bc_contact.c +++ b/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 diff --git a/src/load_field_variables.c b/src/load_field_variables.c index a1f503738..07c891e1d 100644 --- a/src/load_field_variables.c +++ b/src/load_field_variables.c @@ -820,6 +820,22 @@ int load_fv(void) stateVector[EDDY_NU] = fv->eddy_nu; } + if (pdgv[TURB_K]) { + v = TURB_K; + scalar_fv_fill(esp->turb_k, esp_dot->turb_k, esp_old->turb_k, bf[v]->phi, + ei[upd->matrix_index[v]]->dof[v], &(fv->turb_k), &(fv_dot->turb_k), + &(fv_old->turb_k)); + stateVector[TURB_K] = fv->turb_k; + } + + if (pdgv[TURB_OMEGA]) { + v = TURB_OMEGA; + scalar_fv_fill(esp->turb_omega, esp_dot->turb_omega, esp_old->turb_omega, bf[v]->phi, + ei[upd->matrix_index[v]]->dof[v], &(fv->turb_omega), &(fv_dot->turb_omega), + &(fv_old->turb_omega)); + stateVector[TURB_OMEGA] = fv->turb_omega; + } + if (pdgv[LIGHT_INTP]) { v = LIGHT_INTP; scalar_fv_fill(esp->poynt[0], esp_dot->poynt[0], esp_old->poynt[0], bf[v]->phi, @@ -1753,10 +1769,10 @@ int load_fv(void) if (upd->turbulent_info->use_internal_wall_distance) { fv->wall_distance = 0.; if (pdgv[pd->ShapeVar]) { - dofs = ei[pg->imtrx]->dof[pd->ShapeVar]; + dofs = ei[upd->matrix_index[pd->ShapeVar]]->dof[pd->ShapeVar]; for (i = 0; i < dofs; i++) { fv->wall_distance += - upd->turbulent_info->wall_distances[ei[pg->imtrx]->gnn_list[pd->ShapeVar][i]] * + upd->turbulent_info->wall_distances[ei[upd->matrix_index[pd->ShapeVar]]->gnn_list[pd->ShapeVar][i]] * bf[pd->ShapeVar]->phi[i]; } } @@ -3191,6 +3207,40 @@ int load_fv_grads(void) } } + if (pd->gv[TURB_K]) { + v = TURB_K; + bfn = bf[v]; + dofs = ei[upd->matrix_index[v]]->dof[v]; + for (p = 0; p < VIM; p++) { + fv->grad_turb_k[p] = 0.0; + for (i = 0; i < dofs; i++) { + fv->grad_turb_k[p] += *esp->turb_k[i] * bfn->grad_phi[i][p]; + } + } + } else if (zero_unused_grads && upd->vp[pg->imtrx][TURB_K] == -1) { + for (p = 0; p < VIM; p++) { + fv->grad_turb_k[p] = 0.0; + } + } + + if (pd->gv[TURB_OMEGA]) { + v = TURB_OMEGA; + bfn = bf[v]; + dofs = ei[upd->matrix_index[v]]->dof[v]; + for (p = 0; p < VIM; p++) { + fv->grad_turb_omega[p] = 0.0; + fv_old->grad_turb_omega[p] = 0.0; + for (i = 0; i < dofs; i++) { + fv->grad_turb_omega[p] += *esp->turb_omega[i] * bfn->grad_phi[i][p]; + fv_old->grad_turb_omega[p] += *esp_old->turb_omega[i] * bfn->grad_phi[i][p]; + } + } + } else if (zero_unused_grads && upd->vp[pg->imtrx][TURB_OMEGA] == -1) { + for (p = 0; p < VIM; p++) { + fv->grad_turb_omega[p] = 0.0; + } + } + /* * grad(APR) */ @@ -4892,6 +4942,120 @@ int load_fv_mesh_derivs(int okToZero) memset(&(fv->d_grad_eddy_nu_dmesh[0][0][0]), 0, siz); } + if (pd->gv[TURB_K]) { + v = TURB_K; + bfv = bf[v]; + vdofs = ei[upd->matrix_index[v]]->dof[v]; +#ifdef DO_NOT_UNROLL + siz = sizeof(double) * DIM * DIM * MDE; + memset(&(fv->d_grad_turb_k_dmesh[0][0][0]), 0, siz); + for (i = 0; i < vdofs; i++) { + T_i = *esp->turb_k[i]; + for (p = 0; p < dimNonSym; p++) { + for (b = 0; b < dim; b++) { + for (j = 0; j < mdofs; j++) { + fv->d_grad_turb_k_dmesh[p][b][j] += T_i * bfv->d_grad_phi_dmesh[i][p][b][j]; + } + } + } + } +#else + for (j = 0; j < mdofs; j++) { + T_i = *esp->turb_k[0]; + + fv->d_grad_turb_k_dmesh[0][0][j] = T_i * bfv->d_grad_phi_dmesh[0][0][0][j]; + fv->d_grad_turb_k_dmesh[1][1][j] = T_i * bfv->d_grad_phi_dmesh[0][1][1][j]; + fv->d_grad_turb_k_dmesh[1][0][j] = T_i * bfv->d_grad_phi_dmesh[0][1][0][j]; + fv->d_grad_turb_k_dmesh[0][1][j] = T_i * bfv->d_grad_phi_dmesh[0][0][1][j]; + + if (dimNonSym == 3) { + fv->d_grad_turb_k_dmesh[2][2][j] = T_i * bfv->d_grad_phi_dmesh[0][2][2][j]; + fv->d_grad_turb_k_dmesh[2][0][j] = T_i * bfv->d_grad_phi_dmesh[0][2][0][j]; + fv->d_grad_turb_k_dmesh[2][1][j] = T_i * bfv->d_grad_phi_dmesh[0][2][1][j]; + fv->d_grad_turb_k_dmesh[0][2][j] = T_i * bfv->d_grad_phi_dmesh[0][0][2][j]; + fv->d_grad_turb_k_dmesh[1][2][j] = T_i * bfv->d_grad_phi_dmesh[0][1][2][j]; + } + + for (i = 1; i < vdofs; i++) { + T_i = *esp->turb_k[i]; + + fv->d_grad_turb_k_dmesh[0][0][j] += T_i * bfv->d_grad_phi_dmesh[i][0][0][j]; + fv->d_grad_turb_k_dmesh[1][1][j] += T_i * bfv->d_grad_phi_dmesh[i][1][1][j]; + fv->d_grad_turb_k_dmesh[1][0][j] += T_i * bfv->d_grad_phi_dmesh[i][1][0][j]; + fv->d_grad_turb_k_dmesh[0][1][j] += T_i * bfv->d_grad_phi_dmesh[i][0][1][j]; + + if (dimNonSym == 3) { + fv->d_grad_turb_k_dmesh[2][2][j] += T_i * bfv->d_grad_phi_dmesh[i][2][2][j]; + fv->d_grad_turb_k_dmesh[2][0][j] += T_i * bfv->d_grad_phi_dmesh[i][2][0][j]; + fv->d_grad_turb_k_dmesh[2][1][j] += T_i * bfv->d_grad_phi_dmesh[i][2][1][j]; + fv->d_grad_turb_k_dmesh[0][2][j] += T_i * bfv->d_grad_phi_dmesh[i][0][2][j]; + fv->d_grad_turb_k_dmesh[1][2][j] += T_i * bfv->d_grad_phi_dmesh[i][1][2][j]; + } + } + } +#endif + } else if (upd->vp[pg->imtrx][TURB_K] != -1) { + siz = sizeof(double) * DIM * DIM * MDE; + memset(&(fv->d_grad_turb_k_dmesh[0][0][0]), 0, siz); + } + + if (pd->gv[TURB_OMEGA]) { + v = TURB_OMEGA; + bfv = bf[v]; + vdofs = ei[upd->matrix_index[v]]->dof[v]; +#ifdef DO_NOT_UNROLL + siz = sizeof(double) * DIM * DIM * MDE; + memset(&(fv->d_grad_turb_omega_dmesh[0][0][0]), 0, siz); + for (i = 0; i < vdofs; i++) { + T_i = *esp->turb_omega[i]; + for (p = 0; p < dimNonSym; p++) { + for (b = 0; b < dim; b++) { + for (j = 0; j < mdofs; j++) { + fv->d_grad_turb_omega_dmesh[p][b][j] += T_i * bfv->d_grad_phi_dmesh[i][p][b][j]; + } + } + } + } +#else + for (j = 0; j < mdofs; j++) { + T_i = *esp->turb_omega[0]; + + fv->d_grad_turb_omega_dmesh[0][0][j] = T_i * bfv->d_grad_phi_dmesh[0][0][0][j]; + fv->d_grad_turb_omega_dmesh[1][1][j] = T_i * bfv->d_grad_phi_dmesh[0][1][1][j]; + fv->d_grad_turb_omega_dmesh[1][0][j] = T_i * bfv->d_grad_phi_dmesh[0][1][0][j]; + fv->d_grad_turb_omega_dmesh[0][1][j] = T_i * bfv->d_grad_phi_dmesh[0][0][1][j]; + + if (dimNonSym == 3) { + fv->d_grad_turb_omega_dmesh[2][2][j] = T_i * bfv->d_grad_phi_dmesh[0][2][2][j]; + fv->d_grad_turb_omega_dmesh[2][0][j] = T_i * bfv->d_grad_phi_dmesh[0][2][0][j]; + fv->d_grad_turb_omega_dmesh[2][1][j] = T_i * bfv->d_grad_phi_dmesh[0][2][1][j]; + fv->d_grad_turb_omega_dmesh[0][2][j] = T_i * bfv->d_grad_phi_dmesh[0][0][2][j]; + fv->d_grad_turb_omega_dmesh[1][2][j] = T_i * bfv->d_grad_phi_dmesh[0][1][2][j]; + } + + for (i = 1; i < vdofs; i++) { + T_i = *esp->turb_omega[i]; + + fv->d_grad_turb_omega_dmesh[0][0][j] += T_i * bfv->d_grad_phi_dmesh[i][0][0][j]; + fv->d_grad_turb_omega_dmesh[1][1][j] += T_i * bfv->d_grad_phi_dmesh[i][1][1][j]; + fv->d_grad_turb_omega_dmesh[1][0][j] += T_i * bfv->d_grad_phi_dmesh[i][1][0][j]; + fv->d_grad_turb_omega_dmesh[0][1][j] += T_i * bfv->d_grad_phi_dmesh[i][0][1][j]; + + if (dimNonSym == 3) { + fv->d_grad_turb_omega_dmesh[2][2][j] += T_i * bfv->d_grad_phi_dmesh[i][2][2][j]; + fv->d_grad_turb_omega_dmesh[2][0][j] += T_i * bfv->d_grad_phi_dmesh[i][2][0][j]; + fv->d_grad_turb_omega_dmesh[2][1][j] += T_i * bfv->d_grad_phi_dmesh[i][2][1][j]; + fv->d_grad_turb_omega_dmesh[0][2][j] += T_i * bfv->d_grad_phi_dmesh[i][0][2][j]; + fv->d_grad_turb_omega_dmesh[1][2][j] += T_i * bfv->d_grad_phi_dmesh[i][1][2][j]; + } + } + } +#endif + } else if (upd->vp[pg->imtrx][TURB_OMEGA] != -1) { + siz = sizeof(double) * DIM * DIM * MDE; + memset(&(fv->d_grad_turb_omega_dmesh[0][0][0]), 0, siz); + } + if (pd->gv[LIGHT_INTP]) { v = LIGHT_INTP; bfv = bf[v]; diff --git a/src/main.c b/src/main.c index 5077678a8..b8f03ef2c 100644 --- a/src/main.c +++ b/src/main.c @@ -315,9 +315,6 @@ int main(int argc, char **argv) /********************** BEGIN EXECUTION ***************************************/ -#ifdef FP_EXCEPT - feenableexcept((FE_OVERFLOW | FE_DIVBYZERO | FE_INVALID)); -#endif /* assume number of commands is less than or equal to the number of * arguments in the command line minus 1 (1st is program name) */ @@ -802,6 +799,9 @@ int main(int argc, char **argv) wr_dpi(DPI_ptr, ExoFileOut); } } +#ifdef FP_EXCEPT + feenableexcept((FE_OVERFLOW | FE_DIVBYZERO | FE_INVALID)); +#endif /***********************************************************************/ /***********************************************************************/ diff --git a/src/mm_as_alloc.c b/src/mm_as_alloc.c index 5d544606e..bcfc87751 100644 --- a/src/mm_as_alloc.c +++ b/src/mm_as_alloc.c @@ -1232,6 +1232,12 @@ int assembly_alloc(Exo_DB *exo) if (Num_Var_In_Type[imtrx][EDDY_NU]) { esp->eddy_nu = (dbl **)alloc_ptr_1(MDE); } + if (Num_Var_In_Type[imtrx][TURB_K]) { + esp->turb_k = (dbl **)alloc_ptr_1(MDE); + } + if (Num_Var_In_Type[imtrx][TURB_OMEGA]) { + esp->turb_omega = (dbl **)alloc_ptr_1(MDE); + } } /* End of loop over matrices */ diff --git a/src/mm_fill.c b/src/mm_fill.c index fe5c863f9..8857463bd 100644 --- a/src/mm_fill.c +++ b/src/mm_fill.c @@ -25,6 +25,7 @@ /* GOMA include files */ +#include "ad_turbulence.h" #include "ac_stability.h" #include "ac_stability_util.h" #include "bc/rotate.h" @@ -1129,6 +1130,10 @@ Revised: Summer 1998, SY Tam (UNM) err = load_fv_grads(); GOMA_EH(err, "load_fv_grads"); +#ifdef GOMA_ENABLE_SACADO + fill_ad_field_variables(); +#endif + if (pd->gv[R_MESH1]) { err = load_fv_mesh_derivs(1); GOMA_EH(err, "load_fv_mesh_derivs"); @@ -1393,6 +1398,9 @@ Revised: Summer 1998, SY Tam (UNM) err = load_fv_grads(); GOMA_EH(err, "load_fv_grads"); +#ifdef GOMA_ENABLE_SACADO + fill_ad_field_variables(); +#endif if (pd->gv[R_MESH1]) { err = load_fv_mesh_derivs(1); GOMA_EH(err, "load_fv_mesh_derivs"); @@ -2408,7 +2416,11 @@ Revised: Summer 1998, SY Tam (UNM) } if (pde[R_EDDY_NU]) { +#ifdef GOMA_ENABLE_SACADO + err = ad_assemble_spalart_allmaras(time_value, theta, delta_t, &pg_data); +#else err = assemble_spalart_allmaras(time_value, theta, delta_t, &pg_data); +#endif GOMA_EH(err, "assemble_spalart_allmaras"); #ifdef CHECK_FINITE err = CHECKFINITE("assemble_spalart_allmaras"); @@ -2417,6 +2429,34 @@ Revised: Summer 1998, SY Tam (UNM) #endif } + if (pde[R_TURB_K]) { +#ifdef GOMA_ENABLE_SACADO + err = ad_assemble_turb_k(time_value, theta, delta_t, &pg_data); +#else + GOMA_EH(-1, "TURB_K requires Sacado for assembly"); +#endif + GOMA_EH(err, "assemble_turb_k"); +#ifdef CHECK_FINITE + err = CHECKFINITE("assemble_spalart_allmaras"); + if (err) + return -1; +#endif + } + + if (pde[R_TURB_OMEGA]) { +#ifdef GOMA_ENABLE_SACADO + err = ad_assemble_turb_omega(time_value, theta, delta_t, &pg_data); +#else + GOMA_EH(-1, "TURB_OMEGA requires Sacado for assembly"); +#endif + GOMA_EH(err, "assemble_turb_omega"); +#ifdef CHECK_FINITE + err = CHECKFINITE("assemble_spalart_allmaras"); + if (err) + return -1; +#endif + } + if (pde[R_MOMENT0] || pde[R_MOMENT1] || pde[R_MOMENT2] || pde[R_MOMENT3]) { err = assemble_moments(time_value, theta, delta_t, &pg_data); GOMA_EH(err, "assemble_moments"); @@ -4087,6 +4127,10 @@ int matrix_fill_stress(struct GomaLinearSolverData *ams, GOMA_EH(err, "load_fv_mesh_derivs"); } +#ifdef GOMA_ENABLE_SACADO + fill_ad_field_variables(); +#endif + computeCommonMaterialProps_gp(time_value); /* @@ -4989,8 +5033,8 @@ static void load_lec(Exo_DB *exo, /* ptr to EXODUS II finite element mesh db */ int Print_Zeroes = TRUE; lec_it++; - sprintf(lec_name, "lec_dump_%d_%d.txt", DPI_ptr->elem_index_global[ei[pg->imtrx]->ielem], ProcID); - sprintf(ler_name, "ler_dump_%d_%d.txt", DPI_ptr->elem_index_global[ei[pg->imtrx]->ielem], ProcID); + sprintf(lec_name, "lec_dump_%d.txt", ProcID); + sprintf(ler_name, "ler_dump_%d.txt", ProcID); llll = fopen(lec_name, "a"); rrrr = fopen(ler_name, "a"); fprintf(rrrr, "------------------------------------------------------\n"); @@ -4998,6 +5042,11 @@ static void load_lec(Exo_DB *exo, /* ptr to EXODUS II finite element mesh db */ fprintf(rrrr, "lec_it = %d\n", lec_it); fprintf(rrrr, "global element = %d\n", DPI_ptr->elem_index_global[ei[pg->imtrx]->ielem]); fprintf(rrrr, "\nGlobal_NN Proc_NN Equation idof Proc_SolnNum ResidValue\n"); + fprintf(llll, "------------------------------------------------------\n"); + fprintf(llll, "local element = %d, Proc = %d\n", ei[pg->imtrx]->ielem, ProcID); + fprintf(llll, "lec_it = %d\n", lec_it); + fprintf(llll, "global element = %d\n", DPI_ptr->elem_index_global[ei[pg->imtrx]->ielem]); + fprintf(llll, "\nGlobal_NN Proc_NN Equation idof Proc_SolnNum ResidValue\n"); #endif if (strcmp(Matrix_Format, "epetra") == 0) { @@ -5047,7 +5096,7 @@ static void load_lec(Exo_DB *exo, /* ptr to EXODUS II finite element mesh db */ { if (fabs(lec->R[LEC_R_INDEX(MAX_PROB_VAR + ke, i)]) > DBL_SMALL || Print_Zeroes) { - fprintf(rrrr, "%7d %7d MF%-3d %9d - %12d - %10.3g\n", + fprintf(rrrr, "%7d %7d MF%-3d %9d - %12d - %.10f\n", DPI_ptr->node_index_global[gnn], gnn, ke, i, ie, lec->R[LEC_R_INDEX(MAX_PROB_VAR + ke, i)]); } @@ -5096,7 +5145,7 @@ static void load_lec(Exo_DB *exo, /* ptr to EXODUS II finite element mesh db */ { if (fabs(lec->J[LEC_J_INDEX(pe, pv, i, j)]) > DBL_SMALL || Print_Zeroes) { - fprintf(llll, "%9d %9d %9d %9d - %12d - %10.3g\n", pe, pv, i, j, + fprintf(llll, "%9d %9d %9d %9d - %12d - %.10f\n", pe, pv, i, j, ja, lec->J[LEC_J_INDEX(pe, pv, i, j)]); } } @@ -5125,7 +5174,7 @@ static void load_lec(Exo_DB *exo, /* ptr to EXODUS II finite element mesh db */ { if (fabs(lec->J[LEC_J_INDEX(pe, pv, i, j)]) > DBL_SMALL || Print_Zeroes) { - fprintf(llll, "%9d %9d %9d %9d - %12d - %10.3g\n", pe, pv, i, j, + fprintf(llll, "%9d %9d %9d %9d - %12d - %.10f\n", pe, pv, i, j, ja, lec->J[LEC_J_INDEX(pe, pv, i, j)]); } } @@ -5156,7 +5205,7 @@ static void load_lec(Exo_DB *exo, /* ptr to EXODUS II finite element mesh db */ { if (fabs(lec->R[LEC_R_INDEX(pe, i)]) > DBL_SMALL || Print_Zeroes) { - fprintf(rrrr, "%9d %9d - %12d - %10.3g\n", pe, i, ie, + fprintf(rrrr, "%9d %9d - %12d - %.10f\n", pe, i, ie, lec->R[LEC_R_INDEX(pe, i)]); } } @@ -5205,7 +5254,7 @@ static void load_lec(Exo_DB *exo, /* ptr to EXODUS II finite element mesh db */ { if (fabs(lec->J[LEC_J_INDEX(pe, pv, i, j)]) > DBL_SMALL || Print_Zeroes) { - fprintf(llll, "%9d %9d %9d %9d - %12d - %10.3g\n", pe, pv, i, j, + fprintf(llll, "%9d %9d %9d %9d - %12d - %.10f\n", pe, pv, i, j, ja, lec->J[LEC_J_INDEX(pe, pv, i, j)]); } } @@ -5243,7 +5292,7 @@ static void load_lec(Exo_DB *exo, /* ptr to EXODUS II finite element mesh db */ { if (fabs(lec->J[LEC_J_INDEX(pe, pv, i, j)]) > DBL_SMALL || Print_Zeroes) { - fprintf(llll, "%9d %9d %9d %9d - %12d - %10.3g\n", pe, pv, i, j, ja, + fprintf(llll, "%9d %9d %9d %9d - %12d - %.10f\n", pe, pv, i, j, ja, lec->J[LEC_J_INDEX(pe, pv, i, j)]); } } diff --git a/src/mm_fill_momentum.c b/src/mm_fill_momentum.c index 54695be55..f17888230 100644 --- a/src/mm_fill_momentum.c +++ b/src/mm_fill_momentum.c @@ -10,7 +10,7 @@ /* GOMA include files */ #define GOMA_MM_FILL_MOMENTUM_C #include "mm_fill_momentum.h" - +#include "ad_turbulence.h" #include "ac_particles.h" #include "az_aztec.h" #include "bc_colloc.h" @@ -291,7 +291,7 @@ int assemble_momentum(dbl time, /* current time */ rho = density(d_rho, time); if (supg != 0.) { - tau_momentum_shakib(&supg_terms, dim, dt, FALSE); + ad_tau_momentum_shakib(&supg_terms, dim, dt, FALSE); } /* end Petrov-Galerkin addition */ @@ -1189,6 +1189,183 @@ int assemble_momentum(dbl time, /* current time */ } } + if (pdv[TURB_K]) { + var = TURB_K; + pvar = upd->vp[pg->imtrx][var]; + for (j = 0; j < ei[pg->imtrx]->dof[var]; j++) { + phi_j = bf[var]->phi[j]; + + mass = 0.; + advection = 0.; + source = 0.0; + if (supg != 0.) { + dbl d_wt_func = 0; + for (p = 0; p < dim; p++) { + d_wt_func += supg * supg_terms.d_tau_dturb_k[j] * v[p] * bfm->grad_phi[i][p]; + } + if (transient_run) { + if (mass_on) { + mass = v_dot[a] * rho; + mass *= -d_wt_func * d_area; + mass *= mass_etm; + } + + if (porous_brinkman_on) { + mass /= por; + } + + if (particle_momentum_on) { + mass *= ompvf; + } + } + + if (advection_on) { +#ifdef DO_NO_UNROLL + for (p = 0; p < WIM; p++) { + advection += (v[p] - x_dot[p]) * grad_v[p][a]; + } +#else + advection += (v[0] - x_dot[0]) * grad_v[0][a]; + advection += (v[1] - x_dot[1]) * grad_v[1][a]; + if (WIM == 3) + advection += (v[2] - x_dot[2]) * grad_v[2][a]; +#endif + + if (upd->PSPG_advection_correction) { + advection -= pspg[0] * grad_v[0][a]; + advection -= pspg[1] * grad_v[1][a]; + if (WIM == 3) + advection -= pspg[2] * grad_v[2][a]; + } + + advection *= rho; + advection *= -d_wt_func * d_area; + advection *= advection_etm; + + if (porous_brinkman_on) { + por2 = por * por; + advection /= por2; + } + + if (particle_momentum_on) { + advection *= ompvf; + } + } + + /* + * Source term... + */ + + if (source_on) { + source += f[a]; + source *= d_wt_func * d_area; + source *= source_etm; + } + } + + diffusion = 0.; + if (diffusion_on) { + for (p = 0; p < VIM; p++) { + for (q = 0; q < VIM; q++) { + diffusion += grad_phi_i_e_a[p][q] * d_Pi->turb_k[q][p][j]; + } + } + diffusion *= -d_area; + diffusion *= pd->etm[pg->imtrx][eqn][(LOG2_DIFFUSION)]; + } + + lec->J[LEC_J_INDEX(peqn, pvar, ii, j)] += mass + advection + diffusion + source; + } + } + if (pdv[TURB_OMEGA]) { + var = TURB_OMEGA; + pvar = upd->vp[pg->imtrx][var]; + for (j = 0; j < ei[pg->imtrx]->dof[var]; j++) { + phi_j = bf[var]->phi[j]; + + mass = 0.; + advection = 0.; + source = 0.0; + if (supg != 0.) { + dbl d_wt_func = 0; + for (p = 0; p < dim; p++) { + d_wt_func += supg * supg_terms.d_tau_dturb_omega[j] * v[p] * bfm->grad_phi[i][p]; + } + if (transient_run) { + if (mass_on) { + mass = v_dot[a] * rho; + mass *= -d_wt_func * d_area; + mass *= mass_etm; + } + + if (porous_brinkman_on) { + mass /= por; + } + + if (particle_momentum_on) { + mass *= ompvf; + } + } + + if (advection_on) { +#ifdef DO_NO_UNROLL + for (p = 0; p < WIM; p++) { + advection += (v[p] - x_dot[p]) * grad_v[p][a]; + } +#else + advection += (v[0] - x_dot[0]) * grad_v[0][a]; + advection += (v[1] - x_dot[1]) * grad_v[1][a]; + if (WIM == 3) + advection += (v[2] - x_dot[2]) * grad_v[2][a]; +#endif + + if (upd->PSPG_advection_correction) { + advection -= pspg[0] * grad_v[0][a]; + advection -= pspg[1] * grad_v[1][a]; + if (WIM == 3) + advection -= pspg[2] * grad_v[2][a]; + } + + advection *= rho; + advection *= -d_wt_func * d_area; + advection *= advection_etm; + + if (porous_brinkman_on) { + por2 = por * por; + advection /= por2; + } + + if (particle_momentum_on) { + advection *= ompvf; + } + } + + /* + * Source term... + */ + + if (source_on) { + source += f[a]; + source *= d_wt_func * d_area; + source *= source_etm; + } + } + + diffusion = 0.; + if (diffusion_on) { + for (p = 0; p < VIM; p++) { + for (q = 0; q < VIM; q++) { + diffusion += grad_phi_i_e_a[p][q] * d_Pi->turb_omega[q][p][j]; + } + } + diffusion *= -d_area; + diffusion *= pd->etm[pg->imtrx][eqn][(LOG2_DIFFUSION)]; + } + + lec->J[LEC_J_INDEX(peqn, pvar, ii, j)] += mass + advection + diffusion + source; + } + } + /* * J_m_F */ @@ -2619,23 +2796,6 @@ void ve_polymer_stress(double gamma[DIM][DIM], } } -static bool is_evss_f_model(int model) { - switch (model) { - case EVSS_F: - case EVSS_GRADV: - case LOG_CONF: - case LOG_CONF_GRADV: - case LOG_CONF_TRANSIENT: - case LOG_CONF_TRANSIENT_GRADV: - case CONF: - case SQRT_CONF: - return true; - - default: - return false; - } -} - /* * Calculate the total stress tensor for a fluid at a single gauss point * This includes the diagonal pressure contribution @@ -2911,6 +3071,28 @@ void fluid_stress(double Pi[DIM][DIM], STRESS_DEPENDENCE_STRUCT *d_Pi) { } } + var = TURB_K; + if (d_Pi != NULL && pd->v[pg->imtrx][var]) { + for (p = 0; p < VIM; p++) { + for (q = 0; q < VIM; q++) { + for (j = 0; j < ei[pg->imtrx]->dof[var]; j++) { + d_mu->turb_k[j] = mu_num * d_mu->turb_k[j]; + } + } + } + } + + var = TURB_OMEGA; + if (d_Pi != NULL && pd->v[pg->imtrx][var]) { + for (p = 0; p < VIM; p++) { + for (q = 0; q < VIM; q++) { + for (j = 0; j < ei[pg->imtrx]->dof[var]; j++) { + d_mu->turb_omega[j] = mu_num * d_mu->turb_omega[j]; + } + } + } + } + var = PRESSURE; if (d_Pi != NULL && pd->v[pg->imtrx][var]) { for (j = 0; j < ei[pg->imtrx]->dof[var]; j++) { @@ -3465,6 +3647,27 @@ void fluid_stress(double Pi[DIM][DIM], STRESS_DEPENDENCE_STRUCT *d_Pi) { } } + var = TURB_K; + if (d_Pi != NULL && pd->v[pg->imtrx][var]) { + for (p = 0; p < VIM; p++) { + for (q = 0; q < VIM; q++) { + for (j = 0; j < ei[pg->imtrx]->dof[var]; j++) { + d_Pi->turb_k[p][q][j] = d_mu->turb_k[j] * gamma[p][q]; + } + } + } + } + var = TURB_OMEGA; + if (d_Pi != NULL && pd->v[pg->imtrx][var]) { + for (p = 0; p < VIM; p++) { + for (q = 0; q < VIM; q++) { + for (j = 0; j < ei[pg->imtrx]->dof[var]; j++) { + d_Pi->turb_omega[p][q][j] = d_mu->turb_omega[j] * gamma[p][q]; + } + } + } + } + var = PRESSURE; if (d_Pi != NULL && pd->v[pg->imtrx][var]) { for (p = 0; p < VIM; p++) { diff --git a/src/mm_fill_ptrs.c b/src/mm_fill_ptrs.c index 4c9c75c29..d88c1ef9c 100644 --- a/src/mm_fill_ptrs.c +++ b/src/mm_fill_ptrs.c @@ -1865,6 +1865,16 @@ int load_elem_dofptr(const int ielem, load_varType_Interpolation_ptrs(eqn, esp->eddy_nu, esp_old->eddy_nu, esp_dot->eddy_nu); } + eqn = R_TURB_K; + if (upd->ep[pg->imtrx][eqn] >= 0) { + load_varType_Interpolation_ptrs(eqn, esp->turb_k, esp_old->turb_k, esp_dot->turb_k); + } + + eqn = R_TURB_OMEGA; + if (upd->ep[pg->imtrx][eqn] >= 0) { + load_varType_Interpolation_ptrs(eqn, esp->turb_omega, esp_old->turb_omega, esp_dot->turb_omega); + } + eqn = R_STRESS11; if (upd->ep[pg->imtrx][eqn] >= 0) { /* This should loop through all the stress variables @@ -2756,6 +2766,17 @@ int load_elem_dofptr_all(const int ielem, const Exo_DB *exo) { load_varType_Interpolation_ptrs_mat(imtrx, eqn, esp->eddy_nu, esp_old->eddy_nu, esp_dot->eddy_nu); } + eqn = R_TURB_K; + if (upd->ep[imtrx][eqn] >= 0) { + load_varType_Interpolation_ptrs_mat(imtrx, eqn, esp->turb_k, esp_old->turb_k, + esp_dot->turb_k); + } + + eqn = R_TURB_OMEGA; + if (upd->ep[imtrx][eqn] >= 0) { + load_varType_Interpolation_ptrs_mat(imtrx, eqn, esp->turb_omega, esp_old->turb_omega, + esp_dot->turb_omega); + } eqn = R_STRESS11; if (upd->ep[imtrx][eqn] >= 0) { diff --git a/src/mm_fill_stabilization.c b/src/mm_fill_stabilization.c index 207bcfc3e..2b4a78c37 100644 --- a/src/mm_fill_stabilization.c +++ b/src/mm_fill_stabilization.c @@ -22,6 +22,23 @@ #include "std.h" #include "user_mp.h" +bool is_evss_f_model(int model) { + switch (model) { + case EVSS_F: + case EVSS_GRADV: + case LOG_CONF: + case LOG_CONF_GRADV: + case LOG_CONF_TRANSIENT: + case LOG_CONF_TRANSIENT_GRADV: + case CONF: + case SQRT_CONF: + return true; + + default: + return false; + } +} + static dbl yzbeta1(dbl scale, int dim, dbl Y, @@ -148,10 +165,16 @@ void tau_momentum_shakib(momentum_tau_terms *tau_terms, int dim, dbl dt, int psp DENSITY_DEPENDENCE_STRUCT *d_rho = &d_rho_struct; VISCOSITY_DEPENDENCE_STRUCT d_mu_struct; VISCOSITY_DEPENDENCE_STRUCT *d_mu = &d_mu_struct; + VISCOSITY_DEPENDENCE_STRUCT d_mup_struct; + VISCOSITY_DEPENDENCE_STRUCT *d_mup = &d_mup_struct; if (pspg_scale) { dbl rho = density(d_rho, dt); - inv_rho = 1.0 / rho; + if (rho > 0) { + inv_rho = 1.0 / rho; + } else { + inv_rho = 1.0; + } } int interp_eqn = VELOCITY1; @@ -171,7 +194,77 @@ void tau_momentum_shakib(momentum_tau_terms *tau_terms, int dim, dbl dt, int psp } mu = viscosity(gn, gamma, d_mu); + dbl mup = 0.0; + if (pd->gv[POLYMER_STRESS11] && is_evss_f_model(vn->evssModel)) { + for (int mode = 0; mode < vn->modes; mode++) { + /* get polymer viscosity */ + mup = viscosity(ve[mode]->gn, gamma, d_mup); + mu += mup; + + int var = VELOCITY1; + int a, j; + if (pd->v[pg->imtrx][var]) { + for (a = 0; a < WIM; a++) { + for (j = 0; j < ei[pg->imtrx]->dof[var]; j++) { + d_mu->v[a][j] += d_mup->v[a][j]; + } + } + } + + var = MESH_DISPLACEMENT1; + if (pd->v[pg->imtrx][var]) { + for (a = 0; a < dim; a++) { + for (j = 0; j < ei[pg->imtrx]->dof[var]; j++) { + d_mu->X[a][j] += d_mup->X[a][j]; + } + } + } + + var = TEMPERATURE; + if (pd->v[pg->imtrx][var]) { + for (j = 0; j < ei[pg->imtrx]->dof[var]; j++) { + d_mu->T[j] += d_mup->T[j]; + } + } + var = BOND_EVOLUTION; + if (pd->v[pg->imtrx][var]) { + for (j = 0; j < ei[pg->imtrx]->dof[var]; j++) { + d_mu->nn[j] += d_mup->nn[j]; + } + } + + var = RESTIME; + if (pd->v[pg->imtrx][var]) { + for (j = 0; j < ei[pg->imtrx]->dof[var]; j++) { + d_mu->degrade[j] += d_mup->degrade[j]; + } + } + + var = FILL; + if (pd->v[pg->imtrx][var]) { + for (j = 0; j < ei[pg->imtrx]->dof[var]; j++) { + d_mu->F[j] += d_mup->F[j]; + } + } + + var = PRESSURE; + if (pd->v[pg->imtrx][var]) { + for (j = 0; j < ei[pg->imtrx]->dof[var]; j++) { + d_mu->P[j] += d_mup->P[j]; + } + } + + var = MASS_FRACTION; + if (pd->v[pg->imtrx][var]) { + for (int w = 0; w < pd->Num_Species_Eqn; w++) { + for (j = 0; j < ei[pg->imtrx]->dof[var]; j++) { + d_mu->C[w][j] += d_mup->C[w][j]; + } + } + } + } // for mode + } dbl coeff = (12.0 * mu * mu); dbl diff_g_g = 0; for (int i = 0; i < dim; i++) { @@ -260,6 +353,18 @@ void tau_momentum_shakib(momentum_tau_terms *tau_terms, int dim, dbl dt, int psp inv_rho * -0.5 * (d_mu->eddy_nu[k] * d_diff_g_g_dmu) * supg_tau_cubed; } } + if (pd->e[pg->imtrx][TURB_K]) { + for (int k = 0; k < ei[pg->imtrx]->dof[TURB_K]; k++) { + tau_terms->d_tau_dturb_k[k] = + inv_rho * -0.5 * (d_mu->turb_k[k] * d_diff_g_g_dmu) * supg_tau_cubed; + } + } + if (pd->e[pg->imtrx][TURB_OMEGA]) { + for (int k = 0; k < ei[pg->imtrx]->dof[TURB_OMEGA]; k++) { + tau_terms->d_tau_dturb_omega[k] = + inv_rho * -0.5 * (d_mu->turb_omega[k] * d_diff_g_g_dmu) * supg_tau_cubed; + } + } if (pd->e[pg->imtrx][MASS_FRACTION]) { for (int w = 0; w < pd->Num_Species_Eqn; w++) { for (int k = 0; k < ei[pg->imtrx]->dof[MASS_FRACTION]; k++) { diff --git a/src/mm_flux.c b/src/mm_flux.c index 76e7790b5..45752c44c 100644 --- a/src/mm_flux.c +++ b/src/mm_flux.c @@ -7389,6 +7389,24 @@ static int load_fv_sens(void) { } } + v = TURB_K; + fv_sens->turb_k = 0.; + if (pd->v[pg->imtrx][v]) { + dofs = ei[pg->imtrx]->dof[v]; + for (i = 0; i < dofs; i++) { + fv_sens->turb_k += *esp_old->turb_k[i] * bf[v]->phi[i]; + } + } + + v = TURB_OMEGA; + fv_sens->turb_omega = 0.; + if (pd->v[pg->imtrx][v]) { + dofs = ei[pg->imtrx]->dof[v]; + for (i = 0; i < dofs; i++) { + fv_sens->turb_omega += *esp_old->turb_omega[i] * bf[v]->phi[i]; + } + } + /* * Acoustic Pressure */ diff --git a/src/mm_input.c b/src/mm_input.c index 1eb328977..6cd1a6888 100644 --- a/src/mm_input.c +++ b/src/mm_input.c @@ -8454,6 +8454,10 @@ void rd_eq_specs(FILE *ifp, char *input, const int mn) { ce = set_eqn(R_CUR_STRAIN, mtrx_index0, pd_ptr); } else if (!strcasecmp(ts, "eddy_visc")) { ce = set_eqn(R_EDDY_NU, mtrx_index0, pd_ptr); + } else if (!strcasecmp(ts, "turb_k")) { + ce = set_eqn(R_TURB_K, mtrx_index0, pd_ptr); + } else if (!strcasecmp(ts, "turb_omega")) { + ce = set_eqn(R_TURB_OMEGA, mtrx_index0, pd_ptr); } else if (!strcasecmp(ts, "shell_diff_flux")) { ce = set_eqn(R_SHELL_DIFF_FLUX, mtrx_index0, pd_ptr); pd_ptr->Do_Surf_Geometry = 1; @@ -8955,6 +8959,10 @@ void rd_eq_specs(FILE *ifp, char *input, const int mn) { cv = set_var(GRAD_S_V_DOT_N3, mtrx_index0, pd_ptr); } else if (!strcasecmp(ts, "EDDY_NU")) { cv = set_var(EDDY_NU, mtrx_index0, pd_ptr); + } else if (!strcasecmp(ts, "TURB_K")) { + cv = set_var(TURB_K, mtrx_index0, pd_ptr); + } else if (!strcasecmp(ts, "TURB_OMEGA")) { + cv = set_var(TURB_OMEGA, mtrx_index0, pd_ptr); } else if (!strcasecmp(ts, "APR")) { cv = set_var(ACOUS_PREAL, mtrx_index0, pd_ptr); } else if (!strcasecmp(ts, "API")) { @@ -9765,6 +9773,8 @@ void rd_eq_specs(FILE *ifp, char *input, const int mn) { case R_EM_H2_IMAG: case R_EM_H3_IMAG: case R_EDDY_NU: + case R_TURB_K: + case R_TURB_OMEGA: if (fscanf(ifp, "%lf %lf %lf %lf %lf", &(pd_ptr->etm[mtrx_index0][ce][(LOG2_MASS)]), &(pd_ptr->etm[mtrx_index0][ce][(LOG2_ADVECTION)]), diff --git a/src/mm_input_bc.c b/src/mm_input_bc.c index fe9c28e3d..9a8b061dc 100644 --- a/src/mm_input_bc.c +++ b/src/mm_input_bc.c @@ -341,6 +341,7 @@ void rd_bc_specs(FILE *ifp, char *input) { case LS_SOLID_FLUID_BC: overlap_bc = TRUE; /* fall through */ + case SA_WALL_FUNC_BC: case PSPG_BC: case KIN_ELECTRODEPOSITION_BC: /* RSL 5/27/02 */ case VNORM_ELECTRODEPOSITION_BC: /* RSL 5/30/02 */ @@ -719,6 +720,8 @@ void rd_bc_specs(FILE *ifp, char *input) { case EM_CONT_IMAG_BC: case SHELL_TFMP_SAT_BC: case EDDY_NU_BC: + case TURB_K_BC: + case TURB_OMEGA_BC: if (fscanf(ifp, "%lf", &BC_Types[ibc].BC_Data_Float[0]) != 1) { sprintf(err_msg, "%s: Expected 1 flt for %s.", yo, BC_Types[ibc].desc->name1); diff --git a/src/mm_input_mp.c b/src/mm_input_mp.c index 3865d578f..0dee7ca0e 100644 --- a/src/mm_input_mp.c +++ b/src/mm_input_mp.c @@ -1470,6 +1470,8 @@ void rd_mp_specs(FILE *imp, char input[], int mn, char *echo_file) ConstitutiveEquation = TURBULENT_SA; } else if (!strcmp(model_name, "TURBULENT_SA_DYNAMIC")) { ConstitutiveEquation = TURBULENT_SA_DYNAMIC; + } else if (!strcmp(model_name, "TURBULENT_K_OMEGA")) { + ConstitutiveEquation = TURBULENT_K_OMEGA; } else { GOMA_EH(GOMA_ERROR, "Unrecognizable Constitutive Equation"); } @@ -1483,7 +1485,9 @@ void rd_mp_specs(FILE *imp, char input[], int mn, char *echo_file) /* read in constants for constitutive equation if they are input */ - if ((ConstitutiveEquation == NEWTONIAN) || (ConstitutiveEquation == TURBULENT_SA)) { + if ((ConstitutiveEquation == NEWTONIAN) || (ConstitutiveEquation == TURBULENT_SA) || + (ConstitutiveEquation == TURBULENT_SA_DYNAMIC) || + (ConstitutiveEquation == TURBULENT_K_OMEGA)) { model_read = look_for_mat_proptable( imp, "Viscosity", &(mp_glob[mn]->ViscosityModel), &(mp_glob[mn]->viscosity), &(mp_glob[mn]->u_viscosity), &(mp_glob[mn]->len_u_viscosity), diff --git a/src/mm_input_util.c b/src/mm_input_util.c index ec0d0a565..1c07d05db 100644 --- a/src/mm_input_util.c +++ b/src/mm_input_util.c @@ -731,6 +731,10 @@ int variable_string_to_int(const char *input, const char *err_string) var = WSTAR; else if (!strcmp(input, "EDDY_NU")) var = EDDY_NU; + else if (!strcmp(input, "TURB_K")) + var = TURB_K; + else if (!strcmp(input, "TURB_OMEGA")) + var = TURB_OMEGA; /* * Kluge to break up large if block. Problems with HP compiler! diff --git a/src/mm_prob_def.c b/src/mm_prob_def.c index 9e85128ad..abbcfbb9e 100644 --- a/src/mm_prob_def.c +++ b/src/mm_prob_def.c @@ -282,7 +282,7 @@ int setup_pd(void) { (ce == R_EM_H2_REAL) || (ce == R_EM_H3_REAL) || (ce == R_EM_H1_IMAG) || (ce == R_EM_H2_IMAG) || (ce == R_EM_H3_IMAG) || (ce == R_SHELL_FILMP) || (ce == R_SHELL_FILMH) || (ce == R_SHELL_PARTC) || (ce == R_SHELL_ENERGY) || - (ce == R_EDDY_NU) || (ce == R_MOMENT0) || (ce == R_MOMENT1) || + (ce == R_EDDY_NU) || (ce == R_TURB_K) || (ce == R_TURB_OMEGA) || (ce == R_MOMENT0) || (ce == R_MOMENT1) || (ce == R_MOMENT2) || (ce == R_MOMENT3)) { if (pd_glob[mn]->etm[imtrx][ce][(LOG2_MASS)] != 0.) { pd_glob[mn]->e[imtrx][ce] |= T_MASS; diff --git a/src/mm_unknown_map.c b/src/mm_unknown_map.c index 73b88df47..4c8c3c2bd 100644 --- a/src/mm_unknown_map.c +++ b/src/mm_unknown_map.c @@ -1585,6 +1585,12 @@ static void set_interaction_masks(Exo_DB *exo) eqn_var_mask[imtrx][e][v] = 1; v = EDDY_NU; + if (Num_Var_In_Type[imtrx][v]) + eqn_var_mask[imtrx][e][v] = 1; + v = TURB_K; + if (Num_Var_In_Type[imtrx][v]) + eqn_var_mask[imtrx][e][v] = 1; + v = TURB_OMEGA; if (Num_Var_In_Type[imtrx][v]) eqn_var_mask[imtrx][e][v] = 1; @@ -2493,6 +2499,12 @@ static void set_interaction_masks(Exo_DB *exo) if (Num_Var_In_Type[imtrx][v]) eqn_var_mask[imtrx][e][v] = 1; v = EDDY_NU; + if (Num_Var_In_Type[imtrx][v]) + eqn_var_mask[imtrx][e][v] = 1; + v = TURB_K; + if (Num_Var_In_Type[imtrx][v]) + eqn_var_mask[imtrx][e][v] = 1; + v = TURB_OMEGA; if (Num_Var_In_Type[imtrx][v]) eqn_var_mask[imtrx][e][v] = 1; break; @@ -3172,6 +3184,66 @@ static void set_interaction_masks(Exo_DB *exo) if (Num_Var_In_Type[imtrx][v]) eqn_var_mask[imtrx][e][v] = 1; + v = MESH_DISPLACEMENT1; + if (Num_Var_In_Type[imtrx][v]) + eqn_var_mask[imtrx][e][v] = 1; + v = MESH_DISPLACEMENT2; + if (Num_Var_In_Type[imtrx][v]) + eqn_var_mask[imtrx][e][v] = 1; + v = MESH_DISPLACEMENT3; + if (Num_Var_In_Type[imtrx][v]) + eqn_var_mask[imtrx][e][v] = 1; + + v = VELOCITY1; + if (Num_Var_In_Type[imtrx][v]) + eqn_var_mask[imtrx][e][v] = 1; + v = VELOCITY2; + if (Num_Var_In_Type[imtrx][v]) + eqn_var_mask[imtrx][e][v] = 1; + v = VELOCITY3; + if (Num_Var_In_Type[imtrx][v]) + eqn_var_mask[imtrx][e][v] = 1; + + break; + case R_TURB_K: + v = TURB_K; + if (Num_Var_In_Type[imtrx][v]) + eqn_var_mask[imtrx][e][v] = 1; + v = TURB_OMEGA; + if (Num_Var_In_Type[imtrx][v]) + eqn_var_mask[imtrx][e][v] = 1; + + + v = MESH_DISPLACEMENT1; + if (Num_Var_In_Type[imtrx][v]) + eqn_var_mask[imtrx][e][v] = 1; + v = MESH_DISPLACEMENT2; + if (Num_Var_In_Type[imtrx][v]) + eqn_var_mask[imtrx][e][v] = 1; + v = MESH_DISPLACEMENT3; + if (Num_Var_In_Type[imtrx][v]) + eqn_var_mask[imtrx][e][v] = 1; + + v = VELOCITY1; + if (Num_Var_In_Type[imtrx][v]) + eqn_var_mask[imtrx][e][v] = 1; + v = VELOCITY2; + if (Num_Var_In_Type[imtrx][v]) + eqn_var_mask[imtrx][e][v] = 1; + v = VELOCITY3; + if (Num_Var_In_Type[imtrx][v]) + eqn_var_mask[imtrx][e][v] = 1; + + break; + case R_TURB_OMEGA: + v = TURB_OMEGA; + if (Num_Var_In_Type[imtrx][v]) + eqn_var_mask[imtrx][e][v] = 1; + v = TURB_K; + if (Num_Var_In_Type[imtrx][v]) + eqn_var_mask[imtrx][e][v] = 1; + + v = MESH_DISPLACEMENT1; if (Num_Var_In_Type[imtrx][v]) eqn_var_mask[imtrx][e][v] = 1; diff --git a/src/mm_viscosity.c b/src/mm_viscosity.c index 11ff10568..e0c605590 100644 --- a/src/mm_viscosity.c +++ b/src/mm_viscosity.c @@ -23,7 +23,7 @@ #include /* GOMA include files */ - +#include "ad_turbulence.h" #include "density.h" #include "el_elm.h" #include "mm_as.h" @@ -130,7 +130,8 @@ double viscosity(struct Generalized_Newtonian *gn_local, if ((gn_local->ConstitutiveEquation == NEWTONIAN) || (gn_local->ConstitutiveEquation == TURBULENT_SA) || - (gn_local->ConstitutiveEquation == TURBULENT_SA_DYNAMIC)) { + (gn_local->ConstitutiveEquation == TURBULENT_SA_DYNAMIC) || + (gn_local->ConstitutiveEquation == TURBULENT_K_OMEGA)) { if (mp->ViscosityModel == USER) { err = usr_viscosity(mp->u_viscosity); mu = mp->viscosity; @@ -308,9 +309,12 @@ double viscosity(struct Generalized_Newtonian *gn_local, /* Calculate contribution from turbulent viscosity */ if ((gn_local->ConstitutiveEquation == TURBULENT_SA) || (gn_local->ConstitutiveEquation == TURBULENT_SA_DYNAMIC)) { + #ifdef GOMA_ENABLE_SACADO + mu = ad_sa_viscosity(gn_local, d_mu); + #else dbl scale = 1.0; DENSITY_DEPENDENCE_STRUCT d_rho; - if (TURBULENT_SA_DYNAMIC) { + if (gn_local->ConstitutiveEquation == TURBULENT_SA_DYNAMIC) { scale = density(&d_rho, tran->time_value); } int negative_mu_e = FALSE; @@ -340,6 +344,15 @@ double viscosity(struct Generalized_Newtonian *gn_local, } } } + #endif + } + + if (gn_local->ConstitutiveEquation == TURBULENT_K_OMEGA) { +#ifdef GOMA_ENABLE_SACADO + mu = ad_turb_k_omega_sst_viscosity(d_mu); +#else + GOMA_EH(GOMA_ERROR, "TURBULENT_K_OMEGA requires Sacado"); +#endif } } /* end Newtonian section */