!------------------------------------------------------------------------------- ! Code_Saturne version 2.0.6 ! -------------------------- ! This file is part of the Code_Saturne Kernel, element of the ! Code_Saturne CFD tool. ! Copyright (C) 1998-2009 EDF S.A., France ! contact: saturne-support@edf.fr ! The Code_Saturne Kernel 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. ! The Code_Saturne Kernel 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 the Code_Saturne Kernel; if not, write to the ! Free Software Foundation, Inc., ! 51 Franklin St, Fifth Floor, ! Boston, MA 02110-1301 USA !------------------------------------------------------------------------------- subroutine usproj & !================ ( idbia0 , idbra0 , & ndim , ncelet , ncel , nfac , nfabor , nfml , nprfml , & nnod , lndfac , lndfbr , ncelbr , & nvar , nscal , nphas , & nbpmax , nvp , nvep , nivep , ntersl , nvlsta , nvisbr , & nideve , nrdeve , nituse , nrtuse , & ifacel , ifabor , ifmfbr , ifmcel , iprfml , maxelt , lstelt , & ipnfac , nodfac , ipnfbr , nodfbr , itepa , & idevel , ituser , ia , & xyzcen , surfac , surfbo , cdgfac , cdgfbo , xyznod , volume , & dt , rtpa , rtp , propce , propfa , propfb , & coefa , coefb , & ettp , ettpa , tepa , statis , stativ , tslagr , parbor , & rdevel , rtuser , ra ) !=============================================================================== ! Purpose: ! ------- ! User subroutine. ! Called at end of each time step, very general purpose ! (i.e. anything that does not have another dedicated user subroutine) ! Several examples are given here: ! - compute a thermal balance ! (if needed, see note below on adapting this to any scalar) ! - compute global efforts on a subset of faces ! - arbitrarily modify a calculation variable ! - extract a 1 d profile ! - print a moment ! - examples on using parallel utility functions ! These examples are valid when using periodicity (iperio .gt. 0) ! and in parallel (irangp .ge. 0). ! The thermal balance compution also illustates a few other features, ! including the required precautions in parallel or with periodicity): ! - gradient calculation ! - computation of a value depending on cells adjacent to a face ! (see synchronization of Dt and Cp) ! - computation of a global sum in parallel (parsom) ! Cells, boundary faces and interior faces identification ! ======================================================= ! Cells, boundary faces and interior faces may be identified using ! the subroutines 'getcel', 'getfbr' and 'getfac' (respectively). ! getfbr(string, nelts, eltlst): ! - string is a user-supplied character string containing selection criteria; ! - nelts is set by the subroutine. It is an integer value corresponding to ! the number of boundary faces verifying the selection criteria; ! - lstelt is set by the subroutine. It is an integer array of size nelts ! containing the list of boundary faces verifying the selection criteria. ! string may contain: ! - references to colors (ex.: 1, 8, 26, ...) ! - references to groups (ex.: inlet, group1, ...) ! - geometric criteria (ex. x < 0.1, y >= 0.25, ...) ! These criteria may be combined using logical operators ('and', 'or') and ! parentheses. ! Example: '1 and (group2 or group3) and y < 1' will select boundary faces ! of color 1, belonging to groups 'group2' or 'group3' and with face center ! coordinate y less than 1. ! Similarly, interior faces and cells can be identified using the 'getfac' ! and 'getcel' subroutines (respectively). Their syntax are identical to ! 'getfbr' syntax. ! For a more thorough description of the criteria syntax, it can be referred ! to the user guide. !------------------------------------------------------------------------------- ! Arguments !__________________.____._____.________________________________________________. ! name !type!mode ! role ! !__________________!____!_____!________________________________________________! ! idbia0 ! i ! <-- ! number of first free position in ia ! ! idbra0 ! i ! <-- ! number of first free position in ra ! ! ndim ! i ! <-- ! spatial dimension ! ! ncelet ! i ! <-- ! number of extended (real + ghost) cells ! ! ncel ! i ! <-- ! number of cells ! ! nfac ! i ! <-- ! number of interior faces ! ! nfabor ! i ! <-- ! number of boundary faces ! ! nfml ! i ! <-- ! number of families (group classes) ! ! nprfml ! i ! <-- ! number of properties per family (group class) ! ! nnod ! i ! <-- ! number of vertices ! ! lndfac ! i ! <-- ! size of nodfac indexed array ! ! lndfbr ! i ! <-- ! size of nodfbr indexed array ! ! ncelbr ! i ! <-- ! number of cells with faces on boundary ! ! nvar ! i ! <-- ! total number of variables ! ! nscal ! i ! <-- ! total number of scalars ! ! nphas ! i ! <-- ! number of phases ! ! nbpmax ! i ! <-- ! max. number of particles allowed ! ! nvp ! i ! <-- ! number of particle-defined variables ! ! nvep ! i ! <-- ! number of real particle properties ! ! nivep ! i ! <-- ! number of integer particle properties ! ! ntersl ! i ! <-- ! number of return coupling source terms ! ! nvlsta ! i ! <-- ! number of Lagrangian statistical variables ! ! nvisbr ! i ! <-- ! number of boundary statistics ! ! nideve, nrdeve ! i ! <-- ! sizes of idevel and rdevel arrays ! ! nituse, nrtuse ! i ! <-- ! sizes of ituser and rtuser arrays ! ! ifacel(2, nfac) ! ia ! <-- ! interior faces -> cells connectivity ! ! ifabor(nfabor) ! ia ! <-- ! boundary faces -> cells connectivity ! ! ifmfbr(nfabor) ! ia ! <-- ! boundary face family numbers ! ! ifmcel(ncelet) ! ia ! <-- ! cell family numbers ! ! iprfml ! ia ! <-- ! property numbers per family ! ! (nfml, nprfml) ! ! ! lstelt ! ! maxelt ! i ! <-- ! max number of cells and faces (int/boundary) ! ! lstelt(maxelt) ! ia ! --- ! work array ! ! ipnfac(nfac+1) ! ia ! <-- ! interior faces -> vertices index (optional) ! ! nodfac(lndfac) ! ia ! <-- ! interior faces -> vertices list (optional) ! ! ipnfbr(nfabor+1) ! ia ! <-- ! boundary faces -> vertices index (optional) ! ! nodfbr(lndfbr) ! ia ! <-- ! boundary faces -> vertices list (optional) ! ! itepa ! ia ! <-- ! integer particle attributes ! ! (nbpmax, nivep) ! ! ! (containing cell, ...) ! ! idevel(nideve) ! ia ! <-- ! integer work array for temporary development ! ! ituser(nituse) ! ia ! <-- ! user-reserved integer work array ! ! ia(*) ! ia ! --- ! main integer work array ! ! xyzcen ! ra ! <-- ! cell centers ! ! (ndim, ncelet) ! ! ! ! ! surfac ! ra ! <-- ! interior faces surface vectors ! ! (ndim, nfac) ! ! ! ! ! surfbo ! ra ! <-- ! boundary faces surface vectors ! ! (ndim, nfabor) ! ! ! ! ! cdgfac ! ra ! <-- ! interior faces centers of gravity ! ! (ndim, nfac) ! ! ! ! ! cdgfbo ! ra ! <-- ! boundary faces centers of gravity ! ! (ndim, nfabor) ! ! ! ! ! xyznod ! ra ! <-- ! vertex coordinates (optional) ! ! (ndim, nnod) ! ! ! ! ! volume(ncelet) ! ra ! <-- ! cell volumes ! ! dt(ncelet) ! ra ! <-- ! time step (per cell) ! ! rtp, rtpa ! ra ! <-- ! calculated variables at cell centers ! ! (ncelet, *) ! ! ! (at current and previous time steps) ! ! propce(ncelet, *)! ra ! <-- ! physical properties at cell centers ! ! propfa(nfac, *) ! ra ! <-- ! physical properties at interior face centers ! ! propfb(nfabor, *)! ra ! <-- ! physical properties at boundary face centers ! ! coefa, coefb ! ra ! <-- ! boundary conditions ! ! (nfabor, *) ! ! ! ! ! ettp, ettpa ! ra ! <-- ! particle-defined variables ! ! (nbpmax, nvp) ! ! ! (at current and previous time steps) ! ! tepa ! ra ! <-- ! real particle properties ! ! (nbpmax, nvep) ! ! ! (statistical weight, ... ! ! statis ! ra ! <-- ! statistic means ! ! (ncelet, nvlsta)! ! ! ! ! stativ(ncelet, ! ra ! <-- ! accumulator for variance of volume statisitics ! ! nvlsta -1)! ! ! ! ! tslagr ! ra ! <-- ! Lagrangian return coupling term ! ! (ncelet, ntersl)! ! ! on carrier phase ! ! parbor ! ra ! <-- ! particle interaction properties ! ! (nfabor, nvisbr)! ! ! on boundary faces ! ! rdevel(nrdeve) ! ra ! <-> ! real work array for temporary development ! ! rtuser(nrtuse) ! ra ! <-- ! user-reserved real work array ! ! ra(*) ! ra ! --- ! main real work array ! !__________________!____!_____!________________________________________________! ! Type: i (integer), r (real), s (string), a (array), l (logical), ! and composite types (ex: ra real array) ! mode: <-- input, --> output, <-> modifies data, --- work array !=============================================================================== implicit none !=============================================================================== ! Common blocks !=============================================================================== include "dimfbr.h" include "paramx.h" include "pointe.h" include "numvar.h" include "optcal.h" include "cstphy.h" include "cstnum.h" include "entsor.h" include "lagpar.h" include "lagran.h" include "parall.h" include "period.h" include "ppppar.h" include "ppthch.h" include "ppincl.h" !=============================================================================== ! Arguments integer idbia0 , idbra0 integer ndim , ncelet , ncel , nfac , nfabor integer nfml , nprfml integer nnod , lndfac , lndfbr , ncelbr integer nvar , nscal , nphas integer nbpmax , nvp , nvep , nivep integer ntersl , nvlsta , nvisbr integer nideve , nrdeve , nituse , nrtuse integer ifacel(2,nfac) , ifabor(nfabor) integer ifmfbr(nfabor) , ifmcel(ncelet) integer iprfml(nfml,nprfml) integer maxelt, lstelt(maxelt) integer ipnfac(nfac+1), nodfac(lndfac) integer ipnfbr(nfabor+1), nodfbr(lndfbr) integer itepa(nbpmax,nivep) integer idevel(nideve), ituser(nituse) integer ia(*) double precision xyzcen(ndim,ncelet) double precision surfac(ndim,nfac), surfbo(ndim,nfabor) double precision cdgfac(ndim,nfac), cdgfbo(ndim,nfabor) double precision xyznod(ndim,nnod), volume(ncelet) double precision dt(ncelet), rtp(ncelet,*), rtpa(ncelet,*) double precision propce(ncelet,*) double precision propfa(nfac,*), propfb(ndimfb,*) double precision coefa(ndimfb,*), coefb(ndimfb,*) double precision ettp(nbpmax,nvp) , ettpa(nbpmax,nvp) double precision tepa(nbpmax,nvep) double precision statis(ncelet,nvlsta), stativ(ncelet,nvlsta-1) double precision tslagr(ncelet,ntersl) double precision parbor(nfabor,nvisbr) double precision rdevel(nrdeve), rtuser(nrtuse), ra(*) ! Local variables integer idebia, idebra integer iel , ielg , ifac , ifacg , ivar integer iel1 , iel2 , ieltsm integer iortho , impout integer ifinia , ifinra integer igradx , igrady , igradz integer itravx , itravy , itravz , itreco integer inc , iccocg integer nswrgp , imligp , iphydp , iwarnp integer iutile , iphas , iclvar , iii integer ipcrom , ipcvst , iflmas , iflmab , ipccp, ipcvsl integer idimte , itenso , iscal integer ii , nbr , irangv , irang1 , npoint integer imom , ipcmom , idtcm integer itab(3), iun integer ncesmp , icesmp , ismacp , itpsmp integer ilelt , nlelt,isca_num double precision xrtpa , xrtp double precision xbilan , xbilvl , xbilpa , xbilpt double precision xbilsy , xbilen , xbilso , xbildv double precision xbilmi , xbilma double precision epsrgp , climgp , extrap double precision xfluxf , xgamma double precision diipbx, diipby, diipbz, surfbn, distbr double precision visct, flumab , xcp , xvsl, cp0iph, rrr, ctb1, ctb2 double precision xfor(3), xyz(3), xabs, xu, xv, xw, xk, xeps double precision flux_rj, flux_total, flux_entree1, flux_entree2 double precision xc,yc,zc,temp,dscal,vitesse_sortie,double_nlelt double precision vitesse_mag,flux1,flux2,flux3,vx,vy,vz,vitesse_entree,rho_s double precision rho_e !=============================================================================== !=============================================================================== ! 1. Initialization !=============================================================================== ! Memory management idebia = idbia0 idebra = idbra0 !=============================================================================== ! 2. Example: compute energy balance relative to temperature ! ------------------------------------------------------- ! We assume that we want to compute balances (convective and diffusive) ! at the boundaries of the calculation domain represented below ! (with boundaries marked by colors). ! The scalar considered if the temperature. We will also use the ! specific heat (to obtain balances in Joules) ! Domain and associated boundary colors ! ------------------------------------- ! 6 ! -------------------------- ! | | ! | | ! 7 | 1 | 5 ! | ^ | ! | | | ! -------------------------- ! 2 3 4 ! 2, 4, 7 : adiabatic walls ! 6 : wall with fixed temperature ! 3 : inlet ! 5 : outlet ! 1 : symmetry !------------------------------------------------------------------------------- ! To ensure calculations have physical meaning, it is best to use ! a spatially unifor time step (idtvar = 0 or 1). ! In addition, when restarting a calculation, the balance is ! incorrect if inpdt0 = 1 (visct not initialized and t(n-1) not known) !------------------------------------------------------------------------------- ! Temperature variable: iflux= flux +rtp(iel,iu(1))*propce(iel,ipproc(irom(1)))*rtp(iel,isca(1))*surfbo(1,ifac)var = isca(iscalt) (use rtp(iel, ivar)) !------------------------------------------------------------------------------- ! The balance at time step n is equal to: ! n iel=ncelet n-1 ! balance = sum { volume(iel)*cp*rom(iel)*(rtpa(iel,ivar)-rtp(iel,ivar)) } ! iel=1 ! ifac=nfabor ! + sum { ! ifac=1 ! surfbn(ifac)*dt(ifabor(ifac))*cp ! * [visls0(iscalt) + visct(ifabor(ifac))/sigmas(iscalt) ] ! / distbr(ifac) ! * [ coefa(ifac,iclvar) ! + (coefb(ifac,iclvar)-1.d0)*rtp(ifabor(ifac,ivar))] ! } ! ifac=nfabor ! + sum { ! ifac=1 ! dt(ifabor(ifac))*cp ! * rtp(ifabor(ifac,ivar))*(-flumab(ifac)) ! } ! The first term is negative if the amount of energy in the volume ! has decreased (it is 0 in a steady regime). ! The other terms (convection, diffusion) are positive if the amount ! of energy in the volume has increased due to boundary conditions. ! In a steady regime, a positive balance thus indicates an energy gain. !------------------------------------------------------------------------------- ! With 'rom' calculated using the density law from the usphyv subroutine, ! for example: ! n-1 ! rom(iel) = p0 / [rr * (rtpa(iel,ivar) + tkelv)] !------------------------------------------------------------------------------- ! Cp and lambda/Cp may be variable !------------------------------------------------------------------------------- ! Adaptation to an arbitrary scalar ! --------------------------------- ! The approach may be used for the balance of any other scalar (but the ! balances are not in Joules and the specific heat is not used) ! In this case: ! - replace iscalt(iphas) by the number iscal of the required scalar, ! iscal having an allowed range of 1 to nscal. ! - set ipccp to 0 independently of the value of icp(iphas) and assign ! 1 to cp0iph (instead of cp0(iphas)). !=============================================================================== ! The balance is not valid if inpdt0=1 flux_rj = 0. flux1=0. flux2=0. flux3=0. flux_total = 0. flux_entree1 = 0. flux_entree2 = 0. vitesse_sortie=0. vitesse_entree=0. double_nlelt=0. dscal=0. temp=0. vx=0.0 vy=0.0 rho_s=0. vz=0.0 call getfbr('sortie', nlelt, lstelt) !========== isca_num=1 do ilelt = 1, nlelt ifac = lstelt(ilelt) iel = ifabor(ifac) ! associated boundary cell vx=rtp(iel,iu(1))*(rtp(iel,iu(1))) vy=rtp(iel,iv(1))*(rtp(iel,iv(1))) vz=rtp(iel,iw(1))*(rtp(iel,iw(1))) !scalaire 2 dans le xml temp=rtp(iel,isca(isca_num)) ! flux_rj = flux_rj +rtp(iel,iu(1))*propce(iel,ipproc(irom(1)))*rtp(iel,isca(2))*surfbo(1,ifac) flux1=rtp(iel,iu(1))*propce(iel,ipproc(irom(1)))*surfbo(1,ifac) flux2=rtp(iel,iv(1))*propce(iel,ipproc(irom(1)))*surfbo(2,ifac) flux3=rtp(iel,iw(1))*propce(iel,ipproc(irom(1)))*surfbo(3,ifac) flux_rj = flux_rj + rtp(iel,isca(isca_num))*(flux1+flux2+flux3) dscal=max(temp,dscal) rho_s=rho_s+propce(iel,ipproc(irom(1))) vitesse_sortie=vitesse_sortie+sqrt(vx+vy+vz) flux_total = flux_total +flux1+flux2+flux3 if(ilelt.eq.nlelt) then double_nlelt =(dble(nlelt)) vitesse_sortie=vitesse_sortie/double_nlelt rho_s=rho_s/double_nlelt endif enddo flux1=0. flux2=0. flux3=0. vitesse_mag=0. vitesse_entree=0. vx=0.0 vy=0.0 vz=0.0 rho_e=0.0 call getfbr('entree', nlelt, lstelt) !========== do ilelt = 1, nlelt ifac = lstelt(ilelt) iel = ifabor(ifac) ! associated boundary cell vx=rtp(iel,iu(1))*(rtp(iel,iu(1))) vy=rtp(iel,iv(1))*(rtp(iel,iv(1))) vz=rtp(iel,iw(1))*(rtp(iel,iw(1))) vitesse_entree=vitesse_entree+sqrt(vx+vy+vz) flux3=rtp(iel,iw(1))*propce(iel,ipproc(irom(1)))*surfbo(3,ifac) flux2=rtp(iel,iv(1))*propce(iel,ipproc(irom(1)))*surfbo(2,ifac) flux1=rtp(iel,iu(1))*propce(iel,ipproc(irom(1)))*surfbo(1,ifac) rho_e=rho_e+propce(iel,ipproc(irom(1))) flux_entree1 = flux_entree1 +flux1+flux2+flux3 if(ilelt.eq.nlelt) then double_nlelt =(dble(nlelt)) vitesse_entree=vitesse_entree/double_nlelt rho_e=rho_e/double_nlelt endif enddo vx=0.0 vy=0.0 vz=0.0 nbr=1 iel = 1 vx=rtp(iel,iu(1))*(rtp(iel,iu(1))) vy=rtp(iel,iv(1))*(rtp(iel,iv(1))) vz=rtp(iel,iw(1))*(rtp(iel,iw(1))) xyz(1) = sqrt(vx+vy+vz) do iel = 1,ncel if(iel.eq.1) then write (nfecra, *) 'nb = ',ncel endif vitesse_mag= sqrt(vx+vy+vz) xyz(1) = max(xyz(1),vitesse_mag) enddo if (irangp.ge.0) then call parsom(flux_rj) call parsom(dscal) call parsom(flux_total) call parsom(flux_entree1) call parsom(vitesse_sortie) call parsom(vitesse_entree) call parsom(rho_s) call parsom(rho_e) call parrmn(nbr,xyz) endif write (nfecra, *) 'usproj - Flux_polluant = ',flux_rj write (nfecra, *) 'usproj - dscal = ',dscal write (nfecra, *) 'usproj - Flux_total = ',flux_total write (nfecra, *) 'usproj - Flux_entree1 = ', -flux_entree1 write (nfecra, *) 'usproj - vitesse_sortie = ', vitesse_sortie write (nfecra, *) 'usproj - vitesse_entree = ', vitesse_entree write (nfecra, *) 'usproj - masse vol sortie= ', rho_s write (nfecra, *) 'usproj - masse vol entree = ', rho_e write (nfecra, *) 'usproj - vitesse_max = ', -xyz(1) return end subroutine