Extract a 1D profile
This is an example of cs_user_extra_operations which performs 1D profile.
Local variables to be added
integer iel
integer iel1
integer impout
integer ii , irangv , irang1 , npoint
integer iun
double precision xyz(3), xabs, xu, xv, xw, xk, xeps
double precision, dimension(:,:), pointer :: cvar_vel
double precision, dimension(:), pointer :: cvar_k, cvar_ep, cvar_omg
double precision, dimension(:), pointer :: cvar_r11, cvar_r22, cvar_r33
Body
We seek here to extract the profile of U, V, W, k and epsilon on an arbitrary 1D curve based on a curvilear abscissa. The profile is described in the 'profile.dat' file (do not forget to define it as user data in the run script).
- the curve used here is the segment: [(0;0;0),(0;0.1;0)], but the generalization to an arbitrary curve is simple.
- the routine handles parallelism an periodicity, as well as the different turbulence models.
- the 1D curve is discretized into 'npoint' points. For each of these points, we search for the closest cell center and we output the variable values at this cell center. For better consistency, the coordinate which is output is that of the cell center (instead of the initial point).
- we avoid using the same cell multiple times (in case several points an the curve are associated with the same cell).
Here is the corresponding code:
if (ntcabs.eq.ntmabs) then
call field_get_val_v(ivarfl(iu), cvar_vel)
if (itytur.eq.2 .or. iturb.eq.50 &
.or. iturb.eq.60) then
call field_get_val_s(ivarfl(ik), cvar_k)
elseif (itytur.eq.3) then
call field_get_val_s(ivarfl(ir11), cvar_r11)
call field_get_val_s(ivarfl(ir22), cvar_r22)
call field_get_val_s(ivarfl(ir33), cvar_r33)
endif
if ( itytur.eq.2 .or. itytur.eq.3 &
.or. iturb.eq.50) then
call field_get_val_s(ivarfl(iep), cvar_ep)
elseif (iturb.eq.60) then
call field_get_val_s(ivarfl(ik), cvar_k)
call field_get_val_s(ivarfl(iomg), cvar_omg)
endif
impout = impusr(1)
if (irangp.le.0) then
open(impout,file='profile.dat')
write(impout,*) &
'# z(m) U(m/s) V(m/s) W(m/s) k(m2/s2) eps(m2/s3)'
endif
npoint = 200
iel1 = -999
irang1 = -999
do ii = 1, npoint
xyz(1) = 0.d0
xyz(2) = float(ii-1)/float(npoint-1)*0.1d0
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
xabs = xyzcen(2,iel)
xu = cvar_vel(1,iel)
xv = cvar_vel(2,iel)
xw = cvar_vel(3,iel)
xk = 0.d0
xeps = 0.d0
if ( itytur.eq.2 .or. iturb.eq.50 &
.or. iturb.eq.60) then
xk = cvar_k(iel)
elseif (itytur.eq.3) then
xk = ( cvar_r11(iel) + cvar_r22(iel) &
+ cvar_r33(iel)) / 2.d0
endif
if ( itytur.eq.2 .or. itytur.eq.3 &
.or. iturb.eq.50) then
xeps = cvar_ep(iel)
elseif (iturb.eq.60) then
xeps = cmu*cvar_k(iel)*cvar_omg(iel)
endif
else
xabs = 0.d0
xu = 0.d0
xv = 0.d0
xw = 0.d0
xk = 0.d0
xeps = 0.d0
endif
if (irangp.ge.0) then
call parall_bcast_r(irangv, xabs)
call parall_bcast_r(irangv, xu)
call parall_bcast_r(irangv, xv)
call parall_bcast_r(irangv, xw)
call parall_bcast_r(irangv, xk)
call parall_bcast_r(irangv, xeps)
endif
if (irangp.le.0) write(impout,99) xabs, xu, xv, xw, xk, xeps
99 format(6g17.9)
endif
enddo
if (irangp.le.0) close(impout)
endif