/*============================================================================ * This function is called each time step to define physical properties *============================================================================*/ /* VERS */ /* This file is part of Code_Saturne, a general-purpose CFD tool. Copyright (C) 1998-2021 EDF S.A. This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version. This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. */ /*----------------------------------------------------------------------------*/ #include "cs_defs.h" /*---------------------------------------------------------------------------- * Standard C library headers *----------------------------------------------------------------------------*/ #include #include /*---------------------------------------------------------------------------- * PLE library headers *----------------------------------------------------------------------------*/ #include /*---------------------------------------------------------------------------- * Local headers *----------------------------------------------------------------------------*/ #include "cs_headers.h" /*----------------------------------------------------------------------------*/ BEGIN_C_DECLS /*----------------------------------------------------------------------------*/ /*! * \file cs_user_source_terms-momentum.c * * \brief Base examples for additional right-hand side source terms for * momentum equations. * * See the reference \ref cs_user_source_terms.c for documentation. */ /*----------------------------------------------------------------------------*/ /*============================================================================ * User function definitions *============================================================================*/ /*----------------------------------------------------------------------------*/ /*! * \brief Function called at each time step to define source terms. * * \param[in, out] domain pointer to a cs_domain_t structure * \param[in] f_id field id of the variable * \param[out] st_exp explicit source term * \param[out] st_imp implicit part of the source term */ /*----------------------------------------------------------------------------*/ void cs_user_source_terms(cs_domain_t *domain, int f_id, cs_real_t *st_exp, cs_real_t *st_imp) { CS_NO_WARN_IF_UNUSED(domain); /*! [st_meta] */ /* field structure */ const cs_field_t *f = cs_field_by_id(f_id); /* mesh quantities */ /* total number of cells*/ const cs_lnum_t n_cells = cs_glob_mesh->n_cells; /* cell volume*/ const cs_real_t *cell_f_vol = cs_glob_mesh_quantities->cell_vol; /* cell centre coefficient */ const cs_real_3_t *cell_cen = (const cs_real_3_t *)cs_glob_mesh_quantities->cell_cen; /* reference density */ const cs_real_t ro0 = cs_glob_fluid_properties->ro0; /**********************************************************************/ /* density grad*/ /* inicialize a mesh variable m */ const cs_mesh_t *m = cs_glob_mesh; /* inicializing gradp */ cs_real_3_t *gradp; /* Allocate memeory for gradp*/ BFT_MALLOC(gradp,m->n_cells_with_ghosts, cs_real_3_t); /* alloting hyd_p_flag */ int hyd_p_flag = cs_glob_velocity_pressure_param->iphydr; cs_real_3_t *f_ext = (hyd_p_flag ==1) ? (cs_real_3_t *)cs_field_by_name_try("volume_forces"):NULL; /* incialize and allocate values*/ bool use_previous_t = false; /* true if values at previous time step ae used, false if otherwise */ /* inc indicates incriments. 0 for increment and 1 means do not increment */ int inc = 1; /*1 or 0 ; for recompute COCG or not */ int recompute_cocg = 1; /* calculate gradient of PRESSURE and assign it to gradp*/ cs_field_gradient_potential(CS_F_(p), use_previous_t, inc, recompute_cocg, hyd_p_flag, f_ext, gradp); /**********************************************************************/ /* ! [rotational terms in heat equation] ** Warning : * It is assumed here that the thermal scalar is an enthalpy. * If the scalar is a temperature, PWatt does not need to be divided * by Cp because Cp is put outside the diffusion term and multiplied * in the temperature equation as follows: * rho*cell_f_vol*dh/dt + .... =+S for h * rho*Cp*cell_f_vol*dT/dt + .... =S for t * */ if (f == CS_F_(h)) /* enthalpy */ { /* Velocity */ cs_real_3_t *cvar_vel = (cs_real_3_t *)CS_F_(vel)->val; /* logging */ const cs_equation_param_t *eqp = cs_field_get_equation_param_const(f); for (cs_lnum_t i = 0; i < n_cells; i++) { double VelL = cvar_vel[i][0]*gradp[i][0]+cvar_vel[i][2]*gradp[i][2]+cvar_vel[i][1]*gradp[i][1]; st_exp[i] = cell_f_vol[i]*VelL; } } else if (f == CS_F_(t)) /* Temperature */ { /* Velocity */ cs_real_3_t *cvar_vel = (cs_real_3_t *)CS_F_(vel)->val; /* logging */ const cs_equation_param_t *eqp = cs_field_get_equation_param_const(f); for (cs_lnum_t i = 0; i < n_cells; i++) { double VelL = cvar_vel[i][0]*gradp[i][0]+cvar_vel[i][2]*gradp[i][2]+cvar_vel[i][1]*gradp[i][1]; st_exp[i] = cell_f_vol[i]*VelL; } } /* ! [rotational terms in heat equation] */ /*----------------------------------------------------------------------------*/ BFT_FREE(gradp); } /*----------------------------------------------------------------------------*/ END_C_DECLS