From 5b993c8ad1042fe1678db7ed09f94e3a9e3b4ddf Mon Sep 17 00:00:00 2001 From: Alice MAISON Date: Tue, 8 Jun 2021 14:53:06 +0200 Subject: [PATCH] Upload src --- Large_canyon/cs_user_physical_properties.f90 | 280 +++++++++++++++++++ 1 file changed, 280 insertions(+) create mode 100644 Large_canyon/cs_user_physical_properties.f90 diff --git a/Large_canyon/cs_user_physical_properties.f90 b/Large_canyon/cs_user_physical_properties.f90 new file mode 100644 index 0000000..05770f0 --- /dev/null +++ b/Large_canyon/cs_user_physical_properties.f90 @@ -0,0 +1,280 @@ +!------------------------------------------------------------------------------- + +!VERS + +! 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_physical_properties.f90 +!> +!> \brief Definition of physical variable laws. +!> +!> See \subpage physical_properties for examples. +!> + +!=============================================================================== +!> \brief Definition of physical variable laws. +!> +!> \section Warning +!> +!> It is \b forbidden to modify turbulent viscosity \c visct here +!> (a specific subroutine is dedicated to that: \ref usvist) +!> +!> - icp = 0 must have been specified +!> in \ref usipsu if we wish to define a variable specific heat +!> cpro_cp (otherwise: memory overwrite). +!> +!> - the kivisl field integer key (diffusivity_id) +!> must have been specified +!> in \ref usipsu if we wish to define a variable viscosity +!> \c viscls. +!> +!> +!> \remarks +!> - This routine is called at the beginning of each time step +!> Thus, AT THE FIRST TIME STEP (non-restart case), the only +!> values initialized before this call are those defined +!> - in the GUI or \ref usipsu (cs_user_parameters.f90) +!> - density (initialized at \c ro0) +!> - viscosity (initialized at \c viscl0) +!> - in the GUI or \ref cs_user_initialization +!> - calculation variables (initialized at 0 by defaut +!> or to the value given in the GUI or in \ref cs_user_initialization) +!> +!> - We may define here variation laws for cell properties, for: +!> - density: rom kg/m3 +!> - density at boundary faces: romb kg/m3) +!> - molecular viscosity: cpro_viscl kg/(m s) +!> - specific heat: cpro_cp J/(kg degrees) +!> - diffusivities associated with scalars: cpro_vscalt kg/(m s) +!> +!> \b Warning: if the scalar is the temperature, cpro_vscalt corresponds +!> to its conductivity (Lambda) in W/(m K) +!> +!> +!> The types of boundary faces at the previous time step are available +!> (except at the first time step, where arrays \c itypfb and \c itrifb have +!> not been initialized yet) +!> +!> It is recommended to keep only the minimum necessary in this file +!> (i.e. remove all unused example code) +!> +!> +!> \section usphyv_cell_id Cells identification +!> +!> Cells may be identified using the \ref getcel subroutine. +!> The syntax of this subroutine is described in the +!> \ref cs_user_boundary_conditions subroutine, +!> but a more thorough description can be found in the user guide. +! +!------------------------------------------------------------------------------- + +!------------------------------------------------------------------------------- +! Arguments +!______________________________________________________________________________. +! mode name role ! +!______________________________________________________________________________! +!> \param[in] nvar total number of variables +!> \param[in] nscal total number of scalars +!> \param[in] mbrom indicator of filling of romb array +!> \param[in] dt time step (per cell) +!_______________________________________________________________________________ + +subroutine usphyv & + ( nvar , nscal , & + mbrom , & + dt ) + +!=============================================================================== + +!=============================================================================== +! Module files +!=============================================================================== + +use atincl +use paramx +use pointe +use numvar +use optcal +use cstphy +use cstnum +use entsor +use parall +use period +use ppppar +use ppthch +use ppincl +use field +use mesh +use lagran +use cs_c_bindings + +!=============================================================================== + +implicit none + +! Arguments + +integer nvar , nscal + +integer mbrom + +double precision dt(ncelet) + +! Local variables + +integer ivart, iel + +double precision xvart, rhum, rscp, pp, zent +double precision lrhum, lrscp +double precision qsl, deltaq +double precision qliq, qwt, tliq, dum + +double precision, dimension(:), pointer :: brom, crom +double precision, dimension(:), pointer :: cvar_vart, cvar_totwt +double precision, dimension(:), pointer :: cpro_tempc, cpro_liqwt +double precision, dimension(:), pointer :: cpro_beta + +logical activate + +!=============================================================================== +! 0. INITIALISATIONS A CONSERVER +!=============================================================================== + +activate = .false. + +! Initialize variables to avoid compiler warnings + +ivart = -1 + +if (idilat.eq.0) then + call field_get_val_s_by_name("thermal_expansion", cpro_beta) +endif + +! This routine computes the density and the thermodynamic temperature. +! The computations require the pressure profile which is here taken from +! the meteo file. If no meteo file is used, the user should +! give the laws for RHO and T in cs_user_physical_properties.f90 + +if (imeteo.ne.0) return + +!=============================================================================== + +! Positions des variables, coefficients +! ------------------------------------- + +! --- Numero de variable thermique +! (et de ses conditions limites) +! (Pour utiliser le scalaire utilisateur 2 a la place, ecrire +! IVART = ISCA(2) + +if (iscalt.gt.0) then + ivart = isca(iscalt) +else + write(nfecra,9010) iscalt + call csexit (1) +endif + +! --- Masse volumique + +call field_get_val_s(icrom, crom) +call field_get_val_s(ibrom, brom) + +call field_get_val_s(itempc,cpro_tempc) +call field_get_val_s(ivarfl(ivart), cvar_vart) + +! From potential temperature, compute: +! - Temperature in Celsius +! - Density +! ---------------------- + +! Computes the perfect gaz constants according to the physics + +rhum = rair +rscp = rair/cp0 + +lrhum = rair +lrscp = rair/cp0 + +do iel = 1, ncel + + xvart = cvar_vart(iel) ! The thermal scalar is potential temperature + + zent = xyzcen(3,iel) + + pp = p0 + + ! Temperature in Celsius in cell centers: + ! --------------------------------------- + ! law: T = theta * (p/ps) ** (Rair/Cp0) + + cpro_tempc(iel) = xvart*(pp/ps)**lrscp + cpro_tempc(iel) = cpro_tempc(iel) - tkelvi + + ! Density in cell centers: + ! ------------------------ + ! law: RHO = P / ( Rair * T(K) ) + + if (idilat.eq.0) then + crom(iel) = ro0 + ! "delta rho = - beta rho0 delta theta" gives + ! "beta = 1 / theta" + cpro_beta(iel) = 1.d0 / xvart + else + crom(iel) = pp/(lrhum*xvart)*(ps/pp)**lrscp + endif + +enddo + +!=============================================================================== +! FORMATS +!---- + +9010 format( & +'@ ',/,& +'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,& +'@ ',/,& +'@ @@ ATTENTION : ARRET LORS DU CALCUL DES GRANDEURS PHYSIQUES',/,& +'@ ========= ',/,& +'@ APPEL A csexit DANS LE SOUS PROGRAMME atphyv ',/,& +'@ ',/,& +'@ La variable dont dependent les proprietes physiques ne ',/,& +'@ semble pas etre une variable de calcul. ',/,& +'@ En effet, on cherche a utiliser la temperature alors que',/,& +'@ ISCALT = ',I10 ,/,& +'@ Le calcul ne sera pas execute. ',/,& +'@ ',/,& +'@ Verifier le codage de cs_user_physical_properties,', /,& +'@ (et le test lors de la definition de IVART).' ,/,& +'@ Verifier la definition des variables de calcul dans ',/,& +'@ usipsu. Si un scalaire doit jouer le role de la ',/,& +'@ temperature, verifier que ISCALT a ete renseigne. ',/,& +'@ ',/,& +'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,& +'@ ',/) + +!---- +! FIN +!---- + +return + +end subroutine usphyv -- GitLab