/usr/share/code_saturne/user_examples/cs_user_source_terms-scalar_in_a_channel.f90 is in code-saturne-data 4.3.3+repack-1build1.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 | !===============================================================================
! 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 4.3.3
! --------------------------
! This file is part of Code_Saturne, a general-purpose CFD tool.
!
! Copyright (C) 1998-2016 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.
!-------------------------------------------------------------------------------
!===============================================================================
! Purpose:
! -------
!> \file cs_user_source_terms-scalar_in_a_channel.f90
!>
!> \brief User source terms for a scalar in a channel example.
!>
!> 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)
!>
!> \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).
!>
!> 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] ncssmp 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
!===============================================================================
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(ncepdp,6), smacel(ncesmp,nvar)
double precision crvexp(3,ncelet), crvimp(3,3,ncelet)
! Local variables
character*80 chaine
integer iel
double precision ckp, qdm
double precision, dimension(:), pointer :: cpro_rom
!===============================================================================
!===============================================================================
! 1. Initialization
!===============================================================================
if (iwarni(ivar).ge.1) then
call field_get_label(ivarfl(ivar), chaine)
write(nfecra,1000) chaine(1:8)
endif
call field_get_val_s(icrom, cpro_rom)
!===============================================================================
! 2. Example of arbitrary source term for component u:
! S = A * u + B
! appearing in the equation under the form:
! rho*du/dt = S (+ standard Navier-Stokes terms)
!In the following example:
! A = -rho*CKP
! B = XMMT
!
!with:
! CKP = 1.d0 [1/s ] (return term on velocity)
! MMT = 100.d0 [kg/m2/s2] (momentum production by volume and time unit)
!
!which yields:
! crvimp(1, 1, iel) = volume(iel)* A = - volume(iel)*(rho*CKP )
! crvexp(1, iel) = volume(iel)* B = volume(iel)*(XMMT)
! ----------------------------------------------
ckp = 0.d0
qdm = 1.d0
do iel = 1, ncel
crvimp(1, 1, iel) = - volume(iel)*cpro_rom(iel)*ckp
enddo
do iel = 1, ncel
crvexp(1, iel) = volume(iel)*qdm
enddo
!--------
! Formats
!--------
1000 format(' User source termes for variable ',A8,/)
!----
! End
!----
return
end subroutine ustsnv
!===============================================================================
!===============================================================================
! Purpose:
! -------
! User subroutine.
!> \file cs_user_source_terms-scalar_in_a_channel.f90
!> \brief Additional right-hand side source terms for scalar equations (user
!> scalars and specific physics scalars).
!>
!>
!> 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 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] ncssmp 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 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 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(ncepdp,6), smacel(ncesmp,nvar)
double precision crvexp(ncelet), crvimp(ncelet)
!< [loc_var_scal]
! Local variables
integer iel
double precision, dimension(:), pointer :: imasfl
double precision, dimension(:,:), pointer :: vel
double precision ubulk, surf, xx
integer iflmas, ifac, nn
double precision ubulkm
data ubulkm /0.d0/
save ubulkm
!< [loc_var_scal]
!===============================================================================
if (iscal.ne.iscalt) return
!===============================================================================
! 1. Initialization
!===============================================================================
!< [init_scal]
ubulk= 0.d0
surf = 0.d0
!< [init_scal]
! Map field arrays
!< [map_field_scal]
call field_get_val_v(ivarfl(iu), vel)
call field_get_key_int(ivarfl(iu), kimasf, iflmas)
call field_get_val_s(iflmas, imasfl)
!< [map_field_scal
nn= 0
do ifac = 1, nfac
xx=cdgfac(1,ifac)
if(abs(xx).lt.1.d-5) then
nn = nn + 1
ubulk = ubulk + imasfl(ifac)
surf = surf + surfac(1,ifac)
endif
enddo
!< [map_field_scal]
! Parallelism
!< [parallelism_scal]
if (irangp.ge.0) then
call parsom(ubulk)
call parsom(surf)
endif
ubulk = abs(ubulk)/abs(surf)
ubulk = max(ubulk,1d-12)
do iel = 1, ncel
crvimp(iel) = 0.d0
crvexp(iel) = volume(iel)*vel(1,iel)/ubulk
enddo
!< [parallelism_scal]
!--------
! Formats
!--------
!< [format_scal]
1000 format(' User source terms for variable ',A8,/)
!< [format_scal]
!----
! End
!----
!< [end_scal]
return
end subroutine ustssc
!< [end_scal]
|