/*============================================================================ * Additional user-defined source terms for variable equations. *============================================================================*/ /* VERS */ /* This file is part of code_saturne, a general-purpose CFD tool. Copyright (C) 1998-2022 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.c * * \brief Additional source terms for variable equations. * * See \ref user_source_terms for examples. */ /*----------------------------------------------------------------------------*/ /*============================================================================ * User function definitions *============================================================================*/ /*----------------------------------------------------------------------------*/ /*! * \brief Additional user-defined source terms for variable equations * (momentum, scalars, turbulence...). * * This function is called at each time step, for each relevant field. * It is therefore necessary to * test the value of the field id or name to separate * the treatments of the different variables. * * The additional source term is decomposed into an explicit part (st_exp) and * an implicit part (st_imp) that must be provided here. * The resulting equation solved by the code for a scalar f is: * * \f[ \rho*volume*\frac{df}{dt} + .... = st\_imp*f + st\_exp \f] * * Note that st_exp and st_imp are defined after the Finite Volume integration * over the cells, so they include the "volume" term. More precisely: * - st_exp is expressed in kg.[var]/s, where [var] is the unit of the * variable. * Its dimension is the one of the variable (3 for vectors) * - st_imp is expressed in kg/s. * Its dimension is 1 for scalars, 3x3 for vectors. * * The st_exp and st_imp arrays are already initialized to 0 (or a value * defined through the GUI or defined by a model) before entering * the function. It is generally not useful to reset them here. * * For stability reasons, code_saturne will not add -st_imp directly to the * diagonal of the matrix, but Max(-st_imp,0). This way, the st_imp term is * treated implicitely only if it strengthens the diagonal of the matrix. * However, when using the second-order in time scheme, this limitation cannot * be done anymore and -st_imp is added directly. The user should therefore * check for the negativity of st_imp. * * When using the second-order in time scheme, one should supply: * - st_exp at time n * - st_imp at time n+1/2 * * \warning * \parblock * * If the variable is a temperature, the resulting equation solved is: * * rho*Cp*volume*dT/dt + .... = st_imp*T + st_exp * * \endparblock * * Note that st_exp and st_imp are defined after the Finite Volume integration * over the cells, so they include the "volume" term. More precisely: * - st_exp is expressed in W * - st_imp is expressed in W/K * * \par Steep source terms * \parblock * * In case of a complex, non-linear source term, say F(f), for variable f, the * easiest method is to implement the source term explicitly. * * df/dt = .... + F(f(n)) * where f(n) is the value of f at time tn, the beginning of the time step. * * This yields: * st_exp = volume*F(f(n)) * st_imp = 0 * * However, if the source term is potentially steep, this fully explicit * method will probably generate instabilities. It is therefore wiser to * partially implicit the term by writing: * * df/dt = .... + dF/df*f(n+1) - dF/df*f(n) + F(f(n)) * * This yields: * st_exp = volume*( F(f(n)) - dF/df*f(n) ) * st_imp = volume*dF/df * * \endparblock * * \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) { const cs_field_t *fl = cs_field_by_id(f_id); cs_real_3_t *_st_exp = (cs_real_3_t *)st_exp; const cs_lnum_t n_cells = domain->mesh->n_cells; const cs_lnum_t n_i_faces = domain->mesh->n_i_faces; const cs_real_t *cell_vol = domain->mesh_quantities->cell_vol; const cs_real_t *restrict i_face_surf = domain->mesh_quantities->i_face_surf; const cs_real_3_t *i_fac_cog = (const cs_real_3_t *) domain->mesh_quantities->i_face_cog; const cs_real_3_t *i_face_normal = (const cs_real_3_t *)domain->mesh_quantities->i_face_normal; const int nt_cur = cs_glob_time_step->nt_cur; const int nt_max = cs_glob_time_step->nt_max; const int nt_prev = cs_glob_time_step->nt_prev; const int kimasf = cs_field_key_id("INLET"); const cs_real_t *i_massflux = cs_field_by_id(cs_field_get_key_int(fl, kimasf))->val; static cs_real_t mflowa = -1; static cs_real_t dp = -1; cs_real_t mflow = 0; cs_real_t tot_sur = 0; /* Computation of mass flow rate and surface */ for (cs_lnum_t face_id = 0; face_id < n_i_faces; face_id++) { const cs_real_t x_abs = i_fac_cog[face_id][0]; if (x_abs < 0.0) { const cs_real_t x_s = i_face_surf[face_id]; const cs_real_t x_n = i_face_normal[face_id][0]; const cs_real_t x_mf = i_massflux[face_id]; tot_sur += x_s; mflow += x_mf*x_n/x_s; } } /* Bulk on the volume */ cs_real_t xtbul = 0; cs_real_t xubul = 0; cs_real_t xvtot = 0; for (cs_lnum_t cell_id = 0; cell_id < n_cells; cell_id ++) { xtbul += CS_F_(vel)->val[3*cell_id + 0] * cell_vol[cell_id] * CS_F_(t)->val[cell_id]; xubul += CS_F_(vel)->val[3*cell_id + 0] * cell_vol[cell_id]; xvtot += cell_vol[cell_id]; } cs_real_t xdh = cs_notebook_parameter_value_by_name("D"); cs_real_t xre = cs_notebook_parameter_value_by_name("Re"); cs_real_t xmu = cs_notebook_parameter_value_by_name("_mu"); cs_real_t xrho = cs_notebook_parameter_value_by_name("_rho"); cs_real_t ubulk_ref = xre * xmu / (xrho * xdh ); cs_real_t mflow_ref = ubulk_ref*tot_sur; cs_real_t bbuf[5] = {tot_sur, mflow, xtbul, xubul, xvtot}; cs_parall_sum(5, CS_REAL_TYPE, bbuf); tot_sur = bbuf[0], mflow = bbuf[1]; xtbul = bbuf[2], xubul = bbuf[3], xvtot = bbuf[4]; /* For the velocity * =============== */ if (f_id == CS_F_(vel)->id) { bft_printf("surtot %e \n", tot_sur); /* Get the mass flow rate and the pressure gradient * from file in case of restart */ if (nt_cur == 1) { mflow = 0.99*mflow_ref; mflowa = 1.01*mflow_ref; } else if (nt_cur == nt_prev + 1) { if (cs_glob_rank_id <= 0) { FILE *file = fopen("restart/MassFlow.dat", "r"); fscanf(file,"%le %le \n", &mflowa, &dp); fclose(file); bbuf[0] = mflowa; bbuf[1] = dp; } cs_parall_bcast(0, 2, CS_REAL_TYPE, bbuf); mflowa = bbuf[0]; dp = bbuf[1]; } /* Compute the pressure gradient to impose */ const cs_real_t dtref = cs_glob_time_step->dt_ref; dp -= 0.8*(mflow_ref - 2.0*mflow + mflowa )/(tot_sur * dtref); if (nt_cur == nt_max) { if (cs_glob_rank_id <= 0) { FILE *file = fopen("checkpoint/MassFlow.dat", "w"); fprintf(file,"%e %e \n", mflow, dp); fclose(file); } } cs_real_t tot_vol = 0; for (cs_lnum_t i = 0; i < n_cells; i++) { const cs_real_t c_v = cell_vol[i]; tot_vol += c_v; _st_exp[i][0] = -c_v*dp; } cs_parall_sum(1, CS_REAL_TYPE, &tot_vol); domain->mesh_quantities->tot_vol = tot_vol; const cs_real_t ratio = 0.8; const cs_real_t err = (ratio*(mflow_ref-mflow) + (1.0-ratio)*(mflow_ref-mflowa))/tot_vol; bft_printf("target: %e actual: %e old: %e \n error %e\n", mflow_ref, mflow, mflowa, err); mflowa = mflow; } } /*----------------------------------------------------------------------------*/ END_C_DECLS