My Project
programmer's documentation
Calculation of a local Nusselt number

This function is called at the end of each time step, and has a very general purpose (i.e. anything that does not have another dedicated user subroutine)

cs_f_user_extra_operations subroutine

Local variables

! Local variables
integer iscal, ivar, iortho, iprev
integer inc, iccocg, ilelt, irangv
integer ifac, iel, npoint, iel1, irang1, ii, iun, impout, nlelt, neltg
double precision diipbx, diipby, diipbz
double precision xtbulk, xubulk, xyz(3), xtb, xub, tfac, lambda, xab
character*19 namfil
integer, allocatable, dimension(:) :: lstelt
double precision, dimension(:), pointer :: coefap, coefbp, cofafp
double precision, allocatable, dimension(:,:) :: grad
double precision, allocatable, dimension(:) :: treco,treglo,treloc
double precision, allocatable, dimension(:) :: xnusselt
double precision, allocatable, dimension(:) :: xabs, xabsg
double precision, dimension(:,:), pointer :: vel
double precision, dimension(:), pointer :: cvar, viscl
double precision :: height, prandtl, qwall
parameter(height = 1.0d0)
parameter(prandtl = 1.0d0)
parameter(qwall = 1.0d0)

Calculation of the Nusselt number

if (ntcabs.eq.ntmabs) then
allocate(lstelt(max(ncel,nfac,nfabor)))
call field_get_val_v(ivarfl(iu), vel)
iscal = iscalt ! temperature scalar number
ivar = isca(iscal) ! temperature variable number
call field_get_val_s(iviscl, viscl)
call field_get_coefa_s(ivarfl(ivar), coefap)
call field_get_coefb_s(ivarfl(ivar), coefbp)
call field_get_val_s(ivarfl(ivar), cvar)

Compute value reconstructed at I' for boundary faces

allocate(treco(nfabor))
iortho = 0

General cas (for non-orthogonal meshes)

if (iortho.eq.0) then

Allocate a work array for the gradient calculation

allocate(grad(3,ncelet))

Compute gradient

iprev = 0
inc = 1
iccocg = 1
call field_gradient_scalar(ivarfl(ivar), iprev, imrgra, inc, &
iccocg, &
grad)

Compute reconstructed value in boundary cells

do ifac = 1, nfabor
iel = ifabor(ifac)
diipbx = diipb(1,ifac)
diipby = diipb(2,ifac)
diipbz = diipb(3,ifac)
treco(ifac) = coefap(ifac) + coefbp(ifac)*(cvar(iel) &
+ diipbx*grad(1,iel) &
+ diipby*grad(2,iel) &
+ diipbz*grad(3,iel))
enddo

Free memory

deallocate(grad)

Case of orthogonal meshes

else

Compute reconstructed value

Note
Here, we assign the non-reconstructed value.
do ifac = 1, nfabor
iel = ifabor(ifac)
treco(ifac) = coefap(ifac) + coefbp(ifac)*cvar(iel)
enddo
endif
impout = impusr(1)
if (irangp.le.0) then
open(file="Nusselt.dat",unit=impout)
endif
call getfbr('normal[0,-1,0,0.1] and y < 0.01',nlelt,lstelt)
neltg = nlelt
if (irangp.ge.0) then
call parcpt(neltg)
endif
allocate(xabs(nlelt))
allocate(treloc(nlelt))
allocate(xnusselt(neltg))
if (irangp.ge.0) then
allocate(xabsg(neltg))
allocate(treglo(neltg))
endif
do ilelt = 1, nlelt
ifac = lstelt(ilelt)
xabs(ilelt) = cdgfbo(1,ifac)
treloc(ilelt) = treco(ifac)
enddo
if (irangp.ge.0) then
call paragv(nlelt, neltg, xabs, xabsg)
call paragv(nlelt, neltg, treloc, treglo)
endif
do ilelt = 1,neltg

Calculation of the bulk temperature

if (irangp.ge.0) then
xab = xabsg(ilelt)
else
xab = xabs(ilelt)
endif
xtbulk = 0.0d0
xubulk = 0.0d0
npoint = 200
iel1 = -999
do ii = 1, npoint
xyz(1) = xab
xyz(2) = float(ii-1)/float(npoint-1)
xyz(3) = 0.d0
call findpt(ncelet, ncel, xyzcen, xyz(1), xyz(2), xyz(3), iel, irangv)
!==========
if ((iel.ne.iel1).or.(irangv.ne.irang1)) then
iel1 = iel
irang1 = irangv
if (irangp.eq.irangv) then
xtb = volume(iel)*cvar(iel)*vel(1,iel)
xub = volume(iel)*vel(1,iel)
xtbulk = xtbulk + xtb
xubulk = xubulk + xub
lambda = cp0 * viscl(iel) / prandtl
endif
endif
enddo
if (irangp.ge.0) then
iun = 1
call parall_bcast_r(irangv, lambda)
call parsom(xtbulk)
call parsom(xubulk)
endif
xtbulk = xtbulk/xubulk
if (irangp.ge.0) then
tfac = treglo(ilelt)
else
tfac = treloc(ilelt)
endif
xnusselt(ilelt) = qwall * 2.d0 * height / lambda / (tfac - xtbulk)
enddo
if (irangp.eq.-1) then
do ii = 1, neltg
write(impout,'(2E17.9)') xabs(ii)*10, xnusselt(ii)/(0.023d0*30000.d0**(0.8d0)*0.71d0**(0.4d0))
enddo
endif
if (irangp.eq.0) then
call sortc2(xabsg, xnusselt, neltg)
do ii = 1, neltg
write(impout,'(2E17.9)') xabsg(ii)*10, xnusselt(ii)/(0.023d0*30000.d0**(0.8d0)*0.71d0**(0.4d0))
enddo
endif
close(impout)
deallocate(treco)
deallocate(xabs)
deallocate(xnusselt)
deallocate(lstelt)
deallocate(treloc)
if (irangp.ge.0) then
deallocate(xabsg)
deallocate(treglo)
endif
endif
return

sortc2 subroutine

Local variables

implicit none
integer n
double precision tab1(n), tab2(n)
integer ns, ii, jj, kk
double precision tabmin, t2

Body

ns = 1
do ii = 1, n-1
tabmin = 1.d20
kk = 0
do jj = ns,n
if (tabmin.gt.tab1(jj)) then
tabmin = tab1(jj)
kk = jj
endif
enddo
t2 = tab2(kk)
tab1(kk) = tab1(ns)
tab2(kk) = tab2(ns)
tab2(ns) = t2
tab1(ns) = tabmin
ns = ns + 1
enddo
return
end subroutine sortc2
getfbr
subroutine getfbr(fstr, facnb, faces)
Build the list of boundary faces matching a criteria string.
Definition: cs_selector_f2c.f90:111
lambda
Definition: cs_field_pointer.h:112
findpt
subroutine findpt(ncelet, ncel, xyzcen, xx, yy, zz, node, ndrang)
Definition: findpt.f90:57
cs_f_user_extra_operations
subroutine cs_f_user_extra_operations(nvar, nscal, dt)
Definition: cs_user_extra_operations.f90:53