!=============================================================================== ! User source terms definition. ! ! 1) Momentum equation (coupled solver) ! 2) Species transport ! 3) Turbulence (k-epsilon, k-omega, Rij-epsilon, v2-f, Spalart-Allmaras) !=============================================================================== !------------------------------------------------------------------------------- ! Code_Saturne version 6.0 ! ------------------------ ! This file is part of Code_Saturne, a general-purpose CFD tool. ! ! Copyright (C) 1998-2019 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. !------------------------------------------------------------------------------- !=============================================================================== !> \file cs_user_source_terms.f90 !> !> \brief User subroutines for additional right-hand side source terms !> !> See \subpage cs_user_source_terms and !> \subpage cs_user_source_terms-scalar_in_a_channel for examples. !> !=============================================================================== !> \brief Additional right-hand side source terms for velocity components equation !> (Navier-Stokes) !> !> \deprecated Use \ref cs_user_source_terms instead. !> !> \section ustsnv_use Usage !> !> The additional source term is decomposed into an explicit part (\c crvexp) and !> an implicit part (\c crvimp) that must be provided here. !> The resulting equation solved by the code for a velocity is: !> \f[ !> \rho \norm{\vol{\celli}} \DP{\vect{u}} + .... !> = \tens{crvimp} \cdot \vect{u} + \vect{crvexp} !> \f] !> !> Note that \c crvexp and \c crvimp are defined after the Finite Volume integration !> over the cells, so they include the "volume" term. More precisely: !> - crvexp is expressed in kg.m/s2 !> - crvimp is expressed in kg/s !> !> The \c crvexp and \c crvimp arrays are already initialized to 0 !> before entering the !> the routine. It is not needed to do it in the routine (waste of CPU time). !> !> \remark The additional force on \f$ x_i \f$ direction is given by !> \c crvexp(i, iel) + vel(j, iel)* crvimp(j, i, iel). !> !> For stability reasons, Code_Saturne will not add -crvimp directly to the !> diagonal of the matrix, but Max(-crvimp,0). This way, the crvimp 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 -crvimp is added directly. The user should therefore test !> the negativity of crvimp by himself. !> !> When using the second-order in time scheme, one should supply: !> - crvexp at time n !> - crvimp at time n+1/2 !> !> The selection of cells where to apply the source terms is based on a !> \ref getcel command. For more info on the syntax of the \ref getcel command, !> refer to the user manual or to the comments on the similar command !> \ref getfbr in the routine \ref cs_user_boundary_conditions. ! !------------------------------------------------------------------------------- !------------------------------------------------------------------------------- ! Arguments !______________________________________________________________________________. ! mode name role ! !______________________________________________________________________________! !> \param[in] nvar total number of variables !> \param[in] nscal total number of scalars !> \param[in] ncepdp number of cells with head loss terms !> \param[in] ncesmp number of cells with mass source terms !> \param[in] ivar index number of the current variable !> \param[in] icepdc index number of cells with head loss terms !> \param[in] icetsm index number of cells with mass source terms !> \param[in] itypsm type of mass source term for each variable !> (see \ref cs_user_mass_source_terms) !> \param[in] dt time step (per cell) !> \param[in] ckupdc head loss coefficient !> \param[in] smacel value associated to each variable in the mass !> source terms or mass rate (see !> \ref cs_user_mass_source_terms) !> \param[out] crvexp explicit part of the source term !> \param[out] crvimp implicit part of the source term !_______________________________________________________________________________ subroutine ustsnv & ( nvar , nscal , ncepdp , ncesmp , & ivar , & icepdc , icetsm , itypsm , & dt , & ckupdc , smacel , & crvexp , crvimp ) !=============================================================================== !=============================================================================== ! Module files !=============================================================================== use paramx use numvar use entsor use optcal use cstphy use parall use period use mesh use field use field_operator use user_module ! Friess 20/4/21 use cs_c_bindings ! KD 030621 use pointe, only: itypfb !=============================================================================== implicit none ! Arguments integer nvar , nscal integer ncepdp , ncesmp integer ivar integer icepdc(ncepdp) integer icetsm(ncesmp), itypsm(ncesmp,nvar) double precision dt(ncelet) double precision ckupdc(6,ncepdp), smacel(ncesmp,nvar) double precision crvexp(3,ncelet), crvimp(3,3,ncelet) ! Local variables integer, allocatable, dimension(:) :: lstelt integer, allocatable, dimension(:) :: fgpa ! Friess 280521 integer, allocatable, dimension(:,:) :: finsta ! Friess 280521 integer, allocatable, dimension(:,:) :: fgam ! KD 030621 ! Friess 20/4/21 : useful pointers double precision, dimension(:,:), pointer :: vel double precision, dimension(:,:), pointer :: uifujf_moy double precision, dimension(:,:), pointer :: uif_moy double precision, dimension(:), pointer :: cvar_k double precision, dimension(:), pointer :: cvar_omg double precision, dimension(:), pointer :: xk_moy,robs,kres double precision, dimension(:), pointer :: omg_moy, fdes double precision, dimension(:), pointer :: psiomg, rho,rk integer :: i,j,iel, ilelt, nlelt double precision :: ru_L, xupiupi,xk,xe, ru_taur, fdesm1, xkrt, xemoy ! end Friess 20/4/21 : useful pointers !variables gam kd030621 !double precision, dimension(:), pointer :: coefark, coefagradrk, coefaduidf integer inc, iccocg, iprev, ifac double precision, allocatable, dimension(:) :: ktot, duidf, d2uidf2, laplacianf, coefa, coefb double precision xup, xxu1, xxu2, climgp, extrap, epsrgp double precision, allocatable, dimension(:,:) :: gradrk, gradduidf, Xgrad2rk, Ygrad2rk, Zgrad2rk double precision, allocatable, dimension(:) :: Agradrk, Bgradrk, Agradduidf, Bgradduidf, Agrad2rk, Bgrad2rk integer nswrgp, imligp, iwarnp, f_id type(var_cal_opt) :: vcopt_k !end variables gam kd030621 !=============================================================================== !=============================================================================== ! 1. Initialization !=============================================================================== ! Allocate a temporary array for cells selection allocate(lstelt(ncel)) ! Friess 280521 allocate(finsta(3,ncel)) allocate(fgpa(ncel)) allocate(fgam(3,ncel)) !KD 030621 allocate(gradrk(3,ncelet)) !KD allocate(gradduidf(3,ncelet)) !KD allocate(Xgrad2rk(3,ncelet)) !KD allocate(Ygrad2rk(3,ncelet)) allocate(Zgrad2rk(3,ncelet)) allocate(ktot(ncelet)) allocate(duidf(ncelet)) allocate(d2uidf2(ncelet)) allocate(laplacianf(ncel)) allocate(Agradrk(nfabor)) allocate(Bgradrk(nfabor)) allocate(Agradduidf(nfabor)) allocate(Bgradduidf(nfabor)) allocate(Agrad2rk(nfabor)) allocate(Bgrad2rk(nfabor)) ! end Friess 280521 ! Friess 20/4/21 ! terme source grey area mitigation ! appel des champs necessaires call field_get_val_prev_v(ivarfl(iu), vel) call field_get_val_prev_s(ivarfl(ik), cvar_k) call field_get_val_prev_s(ivarfl(iomg), cvar_omg) call field_get_val_v_by_name('uifujf_moy',uifujf_moy) call field_get_val_v_by_name('uif_moy',uif_moy) call field_get_val_s_by_name('k_m',xk_moy) call field_get_val_s_by_name('omega_m',omg_moy) call field_get_val_s_by_name('hybrid_blend',psiomg) call field_get_val_s_by_name('fdes',fdes) call field_get_val_s_by_name('robs',robs) call field_get_val_s_by_name('kres',kres) call field_get_val_s_by_name('hybrid_energy_ratio',rk) call field_get_val_s(icrom, rho) ! On met d'abord a jour les moyennes temporelles estimees ru_dtfilt=1.d1*ru_tphys do iel=1,ncel do j=1,3 uif_moy(j,iel) = uif_moy(j,iel) + dt(iel)/ru_dtfilt*(vel(j,iel)-uif_moy(j,iel)) enddo xk_moy(iel) = xk_moy(iel) + dt(iel)/ru_dtfilt*(cvar_k(iel)-xk_moy(iel)) omg_moy(iel) = omg_moy(iel) + dt(iel)/ru_dtfilt*(cvar_omg(iel)-omg_moy(iel)) do j=1,6 select case (j) case (1) uifujf_moy(j,iel) = uifujf_moy(j,iel) + dt(iel)/ru_dtfilt* & ( vel(1,iel)*vel(1,iel) - uifujf_moy(j,iel)) case (2) uifujf_moy(j,iel) = uifujf_moy(j,iel) + dt(iel)/ru_dtfilt* & ( vel(2,iel)*vel(2,iel) - uifujf_moy(j,iel)) case (3) uifujf_moy(j,iel) = uifujf_moy(j,iel) + dt(iel)/ru_dtfilt* & ( vel(3,iel)*vel(3,iel) - uifujf_moy(j,iel)) case (4) uifujf_moy(j,iel) = uifujf_moy(j,iel) + dt(iel)/ru_dtfilt* & ( vel(1,iel)*vel(2,iel) - uifujf_moy(j,iel)) case (5) uifujf_moy(j,iel) = uifujf_moy(j,iel) + dt(iel)/ru_dtfilt* & ( vel(2,iel)*vel(3,iel) - uifujf_moy(j,iel)) case (6) uifujf_moy(j,iel) = uifujf_moy(j,iel) + dt(iel)/ru_dtfilt* & ( vel(1,iel)*vel(3,iel) - uifujf_moy(j,iel)) end select enddo enddo do iel=1,ncel xupiupi = 0.d0 do j=1,3 xupiupi = xupiupi + uifujf_moy(j,iel) - uif_moy(j,iel) * uif_moy(j,iel) enddo kres(iel)=max(0.d0,0.5d0*xupiupi) robs(iel)= xk_moy(iel) / (xk_moy(iel)+kres(iel)) ktot(iel)=xk_moy(iel)+kres(iel) enddo ! Friess TS gradient de pression artificiel canal perio 120521 do iel=1,ncel fgpa(iel)=ru_utau**2/ru_Ly ! Forcage de Lundgren pour generer de l'instationnarite en debut de simu (-> 5 t_pass) ! Il peut etre necessaire de jouer sur ru_taur (petit -> booste mais risque de faire exploser le calcul) do j=1,3 ru_taur=2.d-1 !1000.d0*dt(iel) if ((rk(iel).le.0.99d0).and.(ttcabs.le.0.1d0*ru_tphys)) then finsta(j,iel)=(robs(iel)-rk(iel))/(max(1.d-1,1.d0-robs(iel)))/(2.d0*ru_taur)* & (vel(j,iel)-uif_moy(j,iel)) else finsta(j,iel)=0.d0 endif enddo enddo ! Terme force GAM 1D selon Y (pour canal pour le moment) KD020621 !calcul termes necessaires pour Cei do iel=1,ncelet call random_number(xxu1) call random_number(xxu2) xup=sqrt(-2.d0*log(xxu1))*cos(2.d0*pii*xxu2) ! Pour loi normale dans duidf et d2uidf2 duidf(iel)=-sqrt(ktot(iel)/(6.d0*(1.d0-rk(iel))))*xup d2uidf2(iel)=-sqrt(ktot(iel)/2.4d1*(1-rk(iel)**3))*xup enddo !arguments et calcul pour gradients iccocg = 1 inc = 1 imrgra = 1 nswrgp = vcopt_k%nswrgr epsrgp = vcopt_k%epsrgr imligp = vcopt_k%imligr iwarnp = vcopt_k%iwarni climgp = vcopt_k%climgr extrap = vcopt_k%extrag f_id = -1 do ifac=1,nfabor ! if (itypfb(ifac).eq.5) then ! iel=ifabor(ifac) Agradrk(ifac)=0.d0 Bgradrk(ifac)=1.d0 Agradduidf(ifac)=0.d0 Bgradduidf(ifac)=1.d0 ! else ! iel=ifabor(ifac) ! Agradrk(ifac)=rk(iel) ! Bgradrk(ifac)=0.d0 ! Agradduidf(ifac)=0.d0 ! Bgradduidf(ifac)=1.d0 ! endif enddo write(*,*) "flag2" ! call gradient_s & !( f_id , imrgra , inc , iccocg , nswrgp , imligp , & ! iwarnp , epsrgp , climgp , extrap , & ! rk , Agradrk , Bgradrk , & ! gradrk ) call field_gradient_scalar(19, iprev, 0, inc, iccocg, gradrk) write(*,*) "flag3" call gradient_s & ( f_id , imrgra , inc , iccocg , nswrgp , imligp , & iwarnp , epsrgp , climgp , extrap , & duidf , Agradduidf , Bgradduidf , & gradduidf ) write(*,*) "flag4" !coefA et coefB pour le grad(grad(rk(1,i))) do ifac=1,nfabor if (itypfb(ifac).eq.5) then iel=ifabor(ifac) Agrad2rk(ifac)=gradrk(1,iel) Bgrad2rk(ifac)=0.d0 else iel=ifabor(ifac) Agrad2rk(ifac)=gradrk(1,iel) Bgrad2rk(ifac)=0.d0 endif enddo write(*,*) "flag5" call gradient_s & ( f_id , imrgra , inc , iccocg , nswrgp , imligp , & iwarnp , epsrgp , climgp , extrap , & gradrk(1,:) , Agrad2rk , Bgrad2rk , & Xgrad2rk ) write(*,*) "flag6" !coefa et coefb pour le grad(grad(rk(2,i))) do ifac=1,nfabor if (itypfb(ifac).eq.5) then iel=ifabor(ifac) Agrad2rk(ifac)=gradrk(2,iel) Bgrad2rk(ifac)=0.d0 else iel=ifabor(ifac) Agrad2rk(ifac)=gradrk(2,iel) Bgrad2rk(ifac)=0.d0 endif enddo call gradient_s & ( f_id , imrgra , inc , iccocg , nswrgp , imligp , & iwarnp , epsrgp , climgp , extrap , & gradrk(2,:) , Agrad2rk , Bgrad2rk , & Ygrad2rk ) write(*,*) "flag7" !coefa et coefb pour le grad(grad(rk(3,i))) do ifac=1,nfabor if (itypfb(ifac).eq.5) then iel=ifabor(ifac) Agrad2rk(ifac)=gradrk(3,iel) Bgrad2rk(ifac)=0.d0 else iel=ifabor(ifac) Agrad2rk(ifac)=gradrk(3,iel) Bgrad2rk(ifac)=0.d0 endif enddo call gradient_s & ( f_id , imrgra , inc , iccocg , nswrgp , imligp , & iwarnp , epsrgp , climgp , extrap , & gradrk(3,:) , Agrad2rk , Bgrad2rk , & Zgrad2rk ) write(*,*) "flag8" !Calcul laplacien + Cei do iel=1,ncel laplacianf(iel)=Xgrad2rk(1,iel)+Ygrad2rk(2,iel)+Zgrad2rk(3,iel) fgam(2,iel)=duidf(iel)*(-viscl0* laplacianf(iel) )-viscl0* gradrk(2,iel) * gradrk(2,iel) * d2uidf2(iel) - & 2.d0*viscl0* gradrk(2,iel) * gradduidf(2,iel) !viscl0 a changer en nu0 par la suite enddo !end force GAM KD020621 call field_get_val_s(icrom, rho) do iel=1,ncel crvexp(1,iel)=rho(iel)*cell_f_vol(iel)*(fgpa(iel)+finsta(1,iel)) crvexp(2,iel)=rho(iel)*cell_f_vol(iel)*(finsta(2,iel)+fgam(2,iel)) !fgam que selon Y pour le canal crvexp(3,iel)=rho(iel)*cell_f_vol(iel)*finsta(3,iel) enddo ! end Friess TS gradient de pression canal perio 120521 !-------- ! Formats !-------- 1000 format(' User source terms for variable ', a8,/) !---- ! End !---- ! Deallocate the temporary array deallocate(lstelt) ! Friess 280521 deallocate(finsta) deallocate(fgpa) deallocate(gradrk) deallocate(gradduidf) deallocate(Xgrad2rk) deallocate(Ygrad2rk) deallocate(Zgrad2rk) deallocate(duidf) deallocate(d2uidf2) deallocate(laplacianf) ! end Friess 280521 return end subroutine ustsnv !=============================================================================== !=============================================================================== !> User subroutine. !> \brief Additional right-hand side source terms for scalar equations (user !> scalars and specific physics scalars). !> !> \deprecated Use \ref cs_user_source_terms instead. !> !> Usage !> ----- !> The routine is called for each scalar, user or specific physisc. It is !> therefore necessary to test the value of the scalar number iscal to separate !> the treatments of the different scalars (if (iscal.eq.p) then ....). !> !> The additional source term is decomposed into an explicit part (crvexp) and !> an implicit part (crvimp) that must be provided here. !> The resulting equation solved by the code for a scalar f is: !> !> \f[ \rho*volume*\frac{df}{dt} + .... = crvimp*f + crvexp \f] !> !> !> Note that crvexp and crvimp are defined after the Finite Volume integration !> over the cells, so they include the "volume" term. More precisely: !> - crvexp is expressed in kg.[scal]/s, where [scal] is the unit of the scalar !> - crvimp is expressed in kg/s !> !> !> The crvexp and crvimp arrays are already initialized to 0 before entering the !> the routine. It is not needed to do it in the routine (waste of CPU time). !> !> For stability reasons, Code_Saturne will not add -crvimp directly to the !> diagonal of the matrix, but Max(-crvimp,0). This way, the crvimp 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 -crvimp is added directly. The user should therefore test !> the negativity of crvimp by himself. !> !> When using the second-order in time scheme, one should supply: !> - crvexp at time n !> - crvimp at time n+1/2 !> !> !> The selection of cells where to apply the source terms is based on a getcel !> command. For more info on the syntax of the getcel command, refer to the !> user manual or to the comments on the similar command \ref getfbr in the routine !> \ref cs_user_boundary_conditions. !> WARNING: If scalar is the temperature, the resulting equation !> solved by the code is: !> !> rho*Cp*volume*dT/dt + .... = crvimp*T + crvexp !> !> !> Note that crvexp and crvimp are defined after the Finite Volume integration !> over the cells, so they include the "volume" term. More precisely: !> - crvexp is expressed in W !> - crvimp is expressed in W/K !> !> !> STEP SOURCE TERMS !>=================== !> In case of a complex, non-linear source term, say F(f), for scalar f, the !> easiest method is to implement the source term explicitely. !> !> df/dt = .... + F(f(n)) !> where f(n) is the value of f at time tn, the beginning of the time step. !> !> This yields : !> crvexp = volume*F(f(n)) !> crvimp = 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: !> crvexp = volume*( F(f(n)) - dF/df*f(n) ) !> crvimp = volume*dF/df !------------------------------------------------------------------------------- !------------------------------------------------------------------------------- ! Arguments !______________________________________________________________________________. ! mode name role !______________________________________________________________________________! !> \param[in] nvar total number of variables !> \param[in] nscal total number of scalars !> \param[in] ncepdp number of cells with head loss terms !> \param[in] ncesmp number of cells with mass source terms !> \param[in] iscal index number of the current scalar !> \param[in] icepdc index number of cells with head loss terms !> \param[in] icetsm index number of cells with mass source terms !> \param[in] itypsm type of mass source term for each variable !> \param[in] (see \ref cs_user_mass_source_terms) !> \param[in] dt time step (per cell) !> \param[in] ckupdc head loss coefficient !> \param[in] smacel value associated to each variable in the mass !> \param[in] source terms or mass rate (see \ref !> cs_user_mass_source_terms) !> \param[out] crvexp explicit part of the source term !> \param[out] crvimp implicit part of the source term !______________________________________________________________________________! subroutine ustssc & ( nvar , nscal , ncepdp , ncesmp , & iscal , & icepdc , icetsm , itypsm , & dt , & ckupdc , smacel , & crvexp , crvimp ) !=============================================================================== ! Module files !=============================================================================== use paramx use numvar use entsor use optcal use cstphy use parall use period use mesh use field !=============================================================================== implicit none ! Arguments integer nvar , nscal integer ncepdp , ncesmp integer iscal integer icepdc(ncepdp) integer icetsm(ncesmp), itypsm(ncesmp,nvar) double precision dt(ncelet) double precision ckupdc(6,ncepdp), smacel(ncesmp,nvar) double precision crvexp(ncelet), crvimp(ncelet) ! Local variables integer, allocatable, dimension(:) :: lstelt !=============================================================================== !=============================================================================== ! 1. Initialization !=============================================================================== ! Allocate a temporary array for cells selection allocate(lstelt(ncel)) !-------- ! Formats !-------- 1000 format(' User source terms for variable ',A8,/) !---- ! End !---- ! Deallocate the temporary array deallocate(lstelt) return end subroutine ustssc !=============================================================================== !=============================================================================== !> \brief Additional right-hand side source terms for vectorial equations !> (user vectors and specific physics vectors). !> !> \deprecated Use \ref cs_user_source_terms instead. !> !> Usage !> ----- !> The routine is called for each vector, user or specific physisc. It is !> therefore necessary to test the value of the vector number iscal to separate !> the treatments of the different vectors (if (iscal.eq.p) then ....). !> !> The additional source term is decomposed into an explicit part (crvexp) and !> an implicit part (crvimp) that must be provided here. !> The resulting equation solved by the code for a vector f is: !> !> \f[ \rho*volume*\frac{d\vect{f}}{dt} + .... = \tens{crvimp}*\vect{f} + !> \vect{crvexp} \f] !> !> !> Note that crvexp and crvimp are defined after the Finite Volume integration !> over the cells, so they include the "volume" term. More precisely: !> - crvexp is expressed in kg.[scal]/s, where [scal] is vector unit !> - crvimp is expressed in kg/s !> !> !> The crvexp and crvimp arrays are already initialized to 0 before entering the !> the routine. It is not needed to do it in the routine (waste of CPU time). !> !> For stability reasons, Code_Saturne will not add -crvimp directly to the !> diagonal of the matrix, but Max(-crvimp,0). This way, the crvimp 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 -crvimp is added directly. The user should therefore test !> the negativity of crvimp by himself. !> !> When using the second-order in time scheme, one should supply: !> - crvexp at time n !> - crvimp at time n+1/2 !> !> !> The selection of cells where to apply the source terms is based on a getcel !> command. For more info on the syntax of the getcel command, refer to the !> user manual or to the comments on the similar command \ref getfbr in the routine !> \ref cs_user_boundary_conditions. !> WARNING: If scalar is the temperature, the resulting equation !> solved by the code is: !> !> rho*Cp*volume*dT/dt + .... = crvimp*T + crvexp !> !> !> Note that crvexp and crvimp are defined after the Finite Volume integration !> over the cells, so they include the "volume" term. More precisely: !> - crvexp is expressed in W !> - crvimp is expressed in W/K !> !> !> STEEP SOURCE TERMS !>=================== !> In case of a complex, non-linear source term, say F(f), for scalar f, the !> easiest method is to implement the source term explicitely. !> !> df/dt = .... + F(f(n)) !> where f(n) is the value of f at time tn, the beginning of the time step. !> !> This yields : !> crvexp = volume*F(f(n)) !> crvimp = 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: !> crvexp = volume*( F(f(n)) - dF/df*f(n) ) !> crvimp = volume*dF/df !------------------------------------------------------------------------------- !------------------------------------------------------------------------------- ! Arguments !______________________________________________________________________________. ! mode name role !______________________________________________________________________________! !> \param[in] nvar total number of variables !> \param[in] nscal total number of scalars !> \param[in] ncepdp number of cells with head loss terms !> \param[in] ncesmp number of cells with mass source terms !> \param[in] iscal index number of the current scalar !> \param[in] icepdc index number of cells with head loss terms !> \param[in] icetsm index number of cells with mass source terms !> \param[in] itypsm type of mass source term for each variable !> \param[in] (see \ref cs_user_mass_source_terms) !> \param[in] dt time step (per cell) !> \param[in] ckupdc head loss coefficient !> \param[in] smacel value associated to each variable in the mass !> \param[in] source terms or mass rate !> (see \ref cs_user_mass_source_terms) !> \param[out] crvexp explicit part of the source term !> \param[out] crvimp implicit part of the source term !______________________________________________________________________________! subroutine ustsvv & ( nvar , nscal , ncepdp , ncesmp , & iscal , & icepdc , icetsm , itypsm , & dt , & ckupdc , smacel , & crvexp , crvimp ) !=============================================================================== ! Module files !=============================================================================== use paramx use numvar use entsor use optcal use cstphy use parall use period use mesh use field !=============================================================================== implicit none ! Arguments integer nvar , nscal integer ncepdp , ncesmp integer iscal integer icepdc(ncepdp) integer icetsm(ncesmp), itypsm(ncesmp,nvar) double precision dt(ncelet) double precision ckupdc(6,ncepdp), smacel(ncesmp,nvar) double precision crvexp(3,ncelet), crvimp(3,3,ncelet) ! Local variables integer, allocatable, dimension(:) :: lstelt !=============================================================================== !=============================================================================== ! 1. Initialization !=============================================================================== ! Allocate a temporary array for cells selection allocate(lstelt(ncel)) !-------- ! Formats !-------- 1000 format(' User source terms for variable ',A8,/) !---- ! End !---- ! Deallocate the temporary array deallocate(lstelt) return end subroutine ustsvv !=============================================================================== !=============================================================================== !> \brief Additional right-hand side source terms for turbulence models !> !> \deprecated Use \ref cs_user_source_terms instead. !> !> \section cs_user_turbulence_source_terms_use Usage !> !> The additional source term is decomposed into an explicit part (crvexp) and !> an implicit part (crvimp) that must be provided here. !> The resulting equations solved by the code are: !> \f[ !> \rho \norm{\vol{\celli}} \DP{\varia} + .... !> = \tens{crvimp} \varia + \vect{crvexp} !> \f] !> where \f$ \varia \f$ is the turbulence field of index \c f_id !> !> Note that crvexp, crvimp are defined after the Finite Volume !> integration over the cells, so they include the "volume" term. More precisely: !> - crvexp is expressed in kg.m2/s2 !> - crvimp is expressed in kg/s !> !> The crvexp, crvimp arrays are already initialized to 0 before !> entering the routine. It is not needed to do it in the routine (waste of CPU time). !> !> For stability reasons, Code_Saturne will not add -crvimp directly to the !> diagonal of the matrix, but Max(-crvimp,0). This way, the crvimp 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 -crvimp is added directly. The user should therefore test !> the negativity of crvimp by himself. !> !> When using the second-order in time scheme, one should supply: !> - crvexp at time n !> - crvimp at time n+1/2 !> !> The selection of cells where to apply the source terms is based on a getcel !> command. For more info on the syntax of the \ref getcel command, refer to the !> user manual or to the comments on the similar command \ref getfbr in the routine !> \ref cs_user_boundary_conditions. ! !------------------------------------------------------------------------------- !------------------------------------------------------------------------------- ! Arguments !______________________________________________________________________________. ! mode name role ! !______________________________________________________________________________! !> \param[in] nvar total number of variables !> \param[in] nscal total number of scalars !> \param[in] ncepdp number of cells with head loss terms !> \param[in] ncesmp number of cells with mass source terms !> \param[in] f_id field index of the current turbulent variable !> \param[in] icepdc index number of cells with head loss terms !> \param[in] icetsm index number of cells with mass source terms !> \param[in] itypsm type of mass source term for each variable !> (see \ref cs_user_mass_source_terms) !> \param[in] ckupdc head loss coefficient !> \param[in] smacel value associated to each variable in the mass !> source terms or mass rate (see !> \ref cs_user_mass_source_terms) !> \param[out] crvexp explicit part of the source term !> \param[out] crvimp implicit part of the source term !_______________________________________________________________________________ subroutine cs_user_turbulence_source_terms & ( nvar , nscal , ncepdp , ncesmp , & f_id , & icepdc , icetsm , itypsm , & ckupdc , smacel , & crvexp , crvimp ) !=============================================================================== !=============================================================================== ! Module files !=============================================================================== use paramx use numvar use entsor use optcal use cstphy use parall use period use mesh use field use cs_f_interfaces use cs_c_bindings !=============================================================================== implicit none ! Arguments integer nvar , nscal integer ncepdp , ncesmp integer f_id integer icepdc(ncepdp) integer icetsm(ncesmp), itypsm(ncesmp,nvar) double precision dt(ncelet) double precision ckupdc(6,ncepdp), smacel(ncesmp,nvar) double precision crvexp(ncelet), crvimp(ncelet) ! Local variables integer, allocatable, dimension(:) :: lstelt !=============================================================================== !=============================================================================== ! 1. Initialization !=============================================================================== ! Allocate a temporary array for cells selection allocate(lstelt(ncel)) !-------- ! Formats !-------- 1000 format(' User source terms for turbulence model',/) !---- ! End !---- ! Deallocate the temporary array deallocate(lstelt) return end subroutine cs_user_turbulence_source_terms !=============================================================================== !=============================================================================== !> \brief Additional right-hand side source terms for turbulence models and !>irijco =1 !> !> \deprecated Use \ref cs_user_source_terms instead. !> !> \section cs_user_rij_source_terms_use Usage !> !> The additional source term is decomposed into an explicit part (crvexp) and !> an implicit part (crvimp) that must be provided here. !> The resulting equations solved by the code are: !> \f[ !> \rho \norm{\vol{\celli}} \DP{\varia} + .... !> = \tens{crvimp} \varia + \vect{crvexp} !> \f] !> where \f$ \varia \f$ is the turbulence field of index \c f_id !> !> Note that crvexp, crvimp are defined after the Finite Volume !> integration over the cells, so they include the "volume" term. More precisely: !> - crvexp is expressed in kg.m2/s2 !> - crvimp is expressed in kg/s !> !> The crvexp, crvimp arrays are already initialized to 0 before !> entering the routine. It is not needed to do it in the routine (waste of CPU time). !> !> For stability reasons, Code_Saturne will not add -crvimp directly to the !> diagonal of the matrix, but Max(-crvimp,0). This way, the crvimp 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 -crvimp is added directly. The user should therefore test !> the negativity of crvimp by himself. !> !> When using the second-order in time scheme, one should supply: !> - crvexp at time n !> - crvimp at time n+1/2 !> !> The selection of cells where to apply the source terms is based on a getcel !> command. For more info on the syntax of the \ref getcel command, refer to the !> user manual or to the comments on the similar command \ref getfbr in the routine !> \ref cs_user_boundary_conditions. ! !------------------------------------------------------------------------------- !------------------------------------------------------------------------------- ! Arguments !______________________________________________________________________________. ! mode name role ! !______________________________________________________________________________! !> \param[in] nvar total number of variables !> \param[in] nscal total number of scalars !> \param[in] ncepdp number of cells with head loss terms !> \param[in] ncesmp number of cells with mass source terms !> \param[in] f_id field index of the current turbulent variable !> \param[in] icepdc index number of cells with head loss terms !> \param[in] icetsm index number of cells with mass source terms !> \param[in] itypsm type of mass source term for each variable !> (see \ref cs_user_mass_source_terms) !> \param[in] ckupdc head loss coefficient !> \param[in] smacel value associated to each variable in the mass !> source terms or mass rate (see !> \ref cs_user_mass_source_terms) !> \param[out] crvexp explicit part of the source term !> \param[out] crvimp implicit part of the source term !_______________________________________________________________________________ subroutine cs_user_turbulence_source_terms2 & ( nvar , nscal , ncepdp , ncesmp , & f_id , & icepdc , icetsm , itypsm , & ckupdc , smacel , & crvexp , crvimp ) !=============================================================================== !=============================================================================== ! Module files !=============================================================================== use paramx use numvar use entsor use optcal use cstphy use parall use period use mesh use field use cs_f_interfaces use cs_c_bindings !=============================================================================== implicit none ! Arguments integer nvar , nscal integer ncepdp , ncesmp integer f_id integer icepdc(ncepdp) integer icetsm(ncesmp), itypsm(ncesmp,nvar) double precision dt(ncelet) double precision ckupdc(6,ncepdp), smacel(ncesmp,nvar) double precision crvexp(6,ncelet), crvimp(6,6,ncelet) ! Local variables integer, allocatable, dimension(:) :: lstelt !=============================================================================== !=============================================================================== ! 1. Initialization !=============================================================================== ! Allocate a temporary array for cells selection allocate(lstelt(ncel)) !-------- ! Formats !-------- 1000 format(' User source terms for turbulence model',/) !---- ! End !---- ! Deallocate the temporary array deallocate(lstelt) return end subroutine cs_user_turbulence_source_terms2