/usr/share/code_saturne/user_examples/cs_user_source_terms-richards_decay.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 | !-------------------------------------------------------------------------------
! 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-richards_decay.f90
!>
!> \brief User source terms 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.
!
!-------------------------------------------------------------------------------
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
use darcy_module
!===============================================================================
implicit none
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)
character*80 chaine
integer iiscvr, iel
integer ilelt, nlelt
double precision lambda
integer, allocatable, dimension(:) :: lstelt
double precision, dimension(:), pointer :: cpro_rom
integer delay_id
double precision, dimension(:), pointer :: delay, saturation
character*80 fname
!===============================================================================
! In groundwater module flow, the radioactive decay of solute is treated as a
! source term in the transport equation
allocate(lstelt(ncel))
call field_get_label(ivarfl(isca(iscal)), chaine)
iiscvr = iscavr(iscal)
if (iwarni(isca(iscal)).ge.1) then
write(nfecra,1000) chaine(1:8)
endif
!< [richards_decay]
! Set the first order decay coefficient
lambda = 1.d-2
! Set radioactive decay for the first solute
if (isca(iscal).eq.1) then
do iel = 1, ncel
crvimp(iel) = -volume(iel)*lambda
enddo
endif
!< [richards_decay]
1000 format(' User source terms for variable ',A8,/)
deallocate(lstelt)
return
end subroutine ustssc
|