1!-------------------------------------------------------------------------------
2
3! This file is part of Code_Saturne, a general-purpose CFD tool.
4!
5! Copyright (C) 1998-2021 EDF S.A.
6!
7! This program is free software; you can redistribute it and/or modify it under
8! the terms of the GNU General Public License as published by the Free Software
9! Foundation; either version 2 of the License, or (at your option) any later
10! version.
11!
12! This program is distributed in the hope that it will be useful, but WITHOUT
13! ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
14! FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more
15! details.
16!
17! You should have received a copy of the GNU General Public License along with
18! this program; if not, write to the Free Software Foundation, Inc., 51 Franklin
19! Street, Fifth Floor, Boston, MA 02110-1301, USA.
20
21!-------------------------------------------------------------------------------
22
23!> \file dvvpst.f90
24!> \brief Standard output of variables on post-processing meshes
25!> (called after \ref cs_user_extra_operations).
26!>
27!------------------------------------------------------------------------------
28
29!------------------------------------------------------------------------------
30! Arguments
31!------------------------------------------------------------------------------
32!   mode          name          role
33!------------------------------------------------------------------------------
34!> \param[in]     nummai        post-processing mesh number
35!> \param[in]     numtyp        post-processing type number
36!>                               - -1: volume
37!>                               - -2: edge
38!>                               - default: nummai
39!> \param[in]     nvar          total number of variables
40!> \param[in]     ncelps        post-processing mesh cells number
41!> \param[in]     nfbrps        number of boundary faces
42!> \param[in]     lstcel        post-processing mesh cell numbers
43!> \param[in]     lstfbr        post-processing mesh boundary faces numbers
44!> \param[in,out] tracel        post processing cell real values
45!> \param[in,out] trafbr        post processing boundary faces real values
46!______________________________________________________________________________
47
48subroutine dvvpst &
49 ( nummai , numtyp ,                                              &
50   nvar   ,                                                       &
51   ncelps , nfbrps ,                                              &
52   lstcel , lstfbr ,                                              &
53   tracel , trafbr )
54
55!===============================================================================
56! Module files
57!===============================================================================
58
59use paramx
60use pointe
61use entsor
62use cstnum
63use cstphy
64use optcal
65use numvar
66use parall
67use period
68use lagran
69use ppppar
70use ppthch
71use ppincl
72use cplsat
73use mesh
74use field
75use field_operator
76use post
77use cs_f_interfaces
78use rotation
79use turbomachinery
80
81!===============================================================================
82
83implicit none
84
85! Arguments
86
87integer          nummai , numtyp
88integer          nvar
89integer          ncelps , nfbrps
90
91integer          lstcel(ncelps), lstfbr(nfbrps)
92
93double precision tracel(ncelps*3)
94double precision trafbr(nfbrps*3)
95
96! Local variables
97
98character(len=80) :: name80
99
100logical          ientla, ivarpr
101integer          inc   , iccocg
102integer          ifac  , iloc  , ivar
103integer          ipp   , idimt , kk   , ll, iel
104integer          fldid, fldprv, keycpl, iflcpl
105integer          ifcsii, iflpst, itplus, iprev, f_id
106
107double precision rbid(1)
108double precision vr(3)
109double precision visls_0
110
111double precision, allocatable, dimension(:,:) :: grad
112double precision, dimension(:), pointer :: tplusp
113double precision, dimension(:), pointer :: valsp, coefap, coefbp
114double precision, dimension(:,:), pointer :: valvp, cofavp, cofbvp
115double precision, dimension(:,:,:), pointer :: cofbtp
116double precision, dimension(:), pointer :: crom
117double precision, dimension(:,:), pointer :: vel
118double precision, dimension(:), pointer :: cpotr, cpoti, cvisii
119double precision, dimension(:), pointer :: cvar_pr
120
121!===============================================================================
122
123! Initialize variables to avoid compiler warnings
124
125ipp = 0
126
127!===============================================================================
128! Fluid domain
129!===============================================================================
130
131if (numtyp .eq. -1) then
132
133  ! Map field arrays
134  call field_get_val_s(ivarfl(ipr), cvar_pr)
135  call field_get_val_v(ivarfl(iu), vel)
136
137  !  Automatic additional variables
138  !  ------------------------------
139
140  ! Relative pressure and velocity in case of turbomachinery
141
142  if (iturbo.ne.0) then
143
144    call field_get_val_s(icrom, crom)
145
146    idimt = 1
147    ientla = .true.
148    ivarpr = .false.
149
150    do iloc = 1, ncelps
151
152      iel = lstcel(iloc)
153      if (irotce(iel).gt.0) then
154        call rotation_velocity(irotce(iel), xyzcen(:,iel), vr)
155      else
156        vr(1) = 0
157        vr(2) = 0
158        vr(3) = 0
159      endif
160
161      tracel(iloc) =   cvar_pr(iel) &
162                     - crom(iel)*0.5d0*(vr(1)**2 + vr(2)**2 + vr(3)**2)
163
164    enddo
165
166    call post_write_var(nummai, 'Rel Pressure', idimt, ientla, ivarpr,  &
167                        ntcabs, ttcabs, tracel, rbid, rbid)
168
169    idimt = 3
170    ientla = .true.
171    ivarpr = .false.
172
173    do iloc = 1, ncelps
174
175      iel = lstcel(iloc)
176      if (irotce(iel).gt.0) then
177        call rotation_velocity(irotce(iel), xyzcen(:,iel), vr)
178      else
179        vr(1) = 0
180        vr(2) = 0
181        vr(3) = 0
182      endif
183
184      tracel(1 + (iloc-1)*idimt) = vel(1,iel) - vr(1)
185      tracel(2 + (iloc-1)*idimt) = vel(2,iel) - vr(2)
186      tracel(3 + (iloc-1)*idimt) = vel(3,iel) - vr(3)
187
188    enddo
189
190    call post_write_var(nummai, 'Rel Velocity', idimt, ientla, ivarpr,  &
191                        ntcabs, ttcabs, tracel, rbid, rbid)
192
193  endif
194
195  ! Absolute pressure and velocity in case of relative coordinate system
196
197  if (icorio.eq.1) then
198
199    call field_get_val_s(icrom, crom)
200
201    idimt = 1
202    ientla = .true.
203    ivarpr = .false.
204
205    do iloc = 1, ncelps
206
207      iel = lstcel(iloc)
208      call rotation_velocity(0, xyzcen(:,iel), vr)
209
210      tracel(iloc) = cvar_pr(iel) + &
211             0.5d0*crom(iel)*(vr(1)**2 + vr(2)**2 + vr(3)**2)
212
213    enddo
214
215    call post_write_var(nummai, 'Abs Pressure', idimt, ientla, ivarpr,  &
216                        ntcabs, ttcabs, tracel, rbid, rbid)
217
218    idimt = 3
219    ientla = .true.
220    ivarpr = .false.
221
222    do iloc = 1, ncelps
223
224      iel = lstcel(iloc)
225      call rotation_velocity(0, xyzcen(:,iel), vr)
226
227      tracel(1 + (iloc-1)*idimt) = vel(1,iel) + vr(1)
228      tracel(2 + (iloc-1)*idimt) = vel(2,iel) + vr(2)
229      tracel(3 + (iloc-1)*idimt) = vel(3,iel) + vr(3)
230
231    enddo
232
233    call post_write_var(nummai, 'Abs Velocity', idimt, ientla, ivarpr,  &
234                        ntcabs, ttcabs, tracel, rbid, rbid)
235
236  endif
237
238!===============================================================================
239! Boundary
240!===============================================================================
241
242else if (numtyp .eq. -2) then
243
244  !  Projection of variables at boundary with no reconstruction
245  !  ----------------------------------------------------------
246
247  call field_get_key_id('coupled', keycpl)
248
249  fldprv = -1
250
251  do ivar = 1, nvar  ! Loop on main cell-based variables
252
253    fldid = ivarfl(ivar)
254
255    if (fldid .eq. fldprv) cycle ! already output for multiple components
256
257    fldprv = fldid
258
259    call field_get_key_int(fldid, keyvis, iflpst)
260
261    if (iand(iflpst, POST_BOUNDARY_NR) .eq. 0) cycle ! nothing to do
262
263    call field_get_dim (fldid, idimt)
264    call field_get_name(fldid, name80(4:80))
265    name80(1:3) = 'bc_'
266
267    !  Compute non-reconstructed values at boundary faces
268
269    if (idimt.ne.1) then
270      call field_get_key_int(fldid, keycpl, iflcpl)
271    else
272      iflcpl = 0
273    endif
274
275    if (idimt.eq.1) then  ! Scalar
276
277      call field_get_val_s(fldid, valsp)
278      call field_get_coefa_s(fldid, coefap)
279      call field_get_coefb_s(fldid, coefbp)
280
281      do iloc = 1, nfbrps
282
283        ifac = lstfbr(iloc)
284        iel = ifabor(ifac)
285
286        trafbr(iloc) =   coefap(ifac) + coefbp(ifac)*valsp(iel)
287
288      enddo
289
290    else if (iflcpl.eq.0) then  ! Uncoupled vector or tensor
291
292      call field_get_val_v(fldid, valvp)
293      call field_get_coefa_v(fldid, cofavp)
294      call field_get_coefb_uv(fldid, cofbvp)
295
296      do kk = 0, idimt-1
297
298        do iloc = 1, nfbrps
299
300          ifac = lstfbr(iloc)
301          iel = ifabor(ifac)
302
303          trafbr(kk + (iloc-1)*idimt + 1)                      &
304               =   cofavp(kk+1,ifac)                           &
305                 + cofbvp(kk+1,ifac)*valvp(kk+1,iel)
306
307        enddo
308
309      enddo
310
311    else ! Coupled vector or tensor
312
313      call field_get_val_v(fldid, valvp)
314      call field_get_coefa_v(fldid, cofavp)
315      call field_get_coefb_v(fldid, cofbtp)
316
317      do kk = 0, idimt-1
318
319        do iloc = 1, nfbrps
320
321          ifac = lstfbr(iloc)
322          iel = ifabor(ifac)
323
324          trafbr(kk + (iloc-1)*idimt + 1) = cofavp(kk+1,ifac)
325
326          do ll = 1, idimt
327            trafbr(kk + (iloc-1)*idimt + 1)                    &
328               =   trafbr(kk + (iloc-1)*idimt + 1)             &
329                 + cofbtp(kk+1,ll,ifac)*valvp(ll,iel)
330          enddo
331
332        enddo
333
334      enddo
335
336    endif ! test on field dimension and interleaving
337
338    ientla = .true.  ! interleaved result values
339    ivarpr = .false. ! defined on work array
340
341    call post_write_var(nummai, trim(name80), idimt, ientla, ivarpr,  &
342                        ntcabs, ttcabs, rbid, rbid, trafbr)
343
344  enddo ! End of loop on variables
345
346  ! Handle stresses at boundary
347  ! ---------------------------
348
349  if (iand(ipstdv(ipstfo), 1) .ne. 0) then
350
351    ! Compute variable values on boundary faces
352
353    call post_stress(nfbrps, lstfbr, trafbr)
354
355    idimt = 3        ! variable dimension
356    ientla = .true.  ! interleaved values
357    ivarpr = .false. ! defined on work array
358
359    call post_write_var(nummai, 'Stress', idimt, ientla, ivarpr,  &
360                        ntcabs, ttcabs, rbid, rbid, trafbr)
361
362  endif
363
364  if (iand(ipstdv(ipstfo), 2) .ne. 0) then
365
366    ! Compute variable values on boundary faces
367
368    call post_stress_tangential(nfbrps, lstfbr, trafbr)
369
370    idimt = 3        ! variable dimension
371    ientla = .true.  ! interleaved values
372    ivarpr = .false. ! defined on work array
373
374    call post_write_var(nummai, 'Shear Stress', idimt, ientla, ivarpr,  &
375                        ntcabs, ttcabs, rbid, rbid, trafbr)
376
377  endif
378
379  if (iand(ipstdv(ipstfo), 4) .ne. 0) then
380
381    ! Calcul des valeurs de la variable sur les faces de bord
382
383    call post_stress_normal(nfbrps, lstfbr, trafbr)
384
385    idimt = 1        ! variable dimension
386    ientla = .true.  ! interleaved values
387    ivarpr = .false. ! defined on work array
388
389    call post_write_var(nummai, 'Normal Stress', idimt, ientla, ivarpr,  &
390                        ntcabs, ttcabs, rbid, rbid, trafbr)
391
392  endif
393
394  ! T+ near the boundary
395  ! --------------------
396
397  if (ipstdv(ipsttp).ne.0) then
398
399    call field_get_id_try('tplus', itplus)
400
401    if (itplus.ge.0) then
402
403      call field_get_val_s(itplus, tplusp)
404
405      idimt = 1        ! variable dimension
406      ientla = .true.  ! interleaved values
407      ivarpr = .true.  ! defined on parent array
408
409      if (itherm .eq. 1) then
410        name80 = 'Tplus'
411      else if (itherm .eq. 2) then
412        name80 = 'Hplus'
413      else if (itherm .eq. 3) then
414        name80 = 'Eplus'
415      else
416        return
417      endif
418
419      call post_write_var(nummai, name80, idimt, ientla, ivarpr,  &
420                          ntcabs, ttcabs, rbid, rbid, tplusp)
421
422    endif ! end of test on presence ot T+
423
424  endif ! end of test on output of y+
425
426  ! Thermal flux at boundary
427  ! ------------------------
428  !  If working with enthalpy, compute an enthalpy flux
429
430  if (ipstdv(ipstft).ne.0) then
431
432    if (iscalt.gt.0) then
433
434      call post_boundary_thermal_flux(nfbrps, lstfbr, trafbr)
435
436      idimt = 1        ! variable dimension
437      ientla = .true.  ! interleaved values
438      ivarpr = .false. ! defined on work array
439
440      call post_write_var(nummai, 'Input thermal flux', idimt, ientla, ivarpr,  &
441                          ntcabs, ttcabs, rbid, rbid, trafbr)
442
443    endif
444
445  endif
446
447  ! Nusselt at the boundary
448  ! -----------------------
449
450  if (ipstdv(ipstnu).ne.0) then
451
452    idimt = 1        ! variable dimension
453    ientla = .true.  ! interleaved values
454    ivarpr = .false. ! defined on work array
455
456    ! Compute variable on boundary faces
457
458    call post_boundary_nusselt(nfbrps, lstfbr, trafbr)
459
460    call post_write_var(nummai, 'Dimensionless heat flux', idimt, ientla, ivarpr,  &
461                        ntcabs, ttcabs, rbid, rbid, trafbr)
462
463  endif ! end of test on output of Nusselt
464
465endif ! end of test on postprocessing mesh number
466
467!===============================================================================
468! Electric module variables
469!===============================================================================
470
471if (numtyp.eq.-1) then
472
473  if (     ippmod(ieljou).ge.1                                      &
474      .or. ippmod(ielarc).ge.1) then
475
476    allocate(grad(3,ncelet))
477
478    ! For Joule Heating by direct conduction:
479    !   gradient of the imaginary component of the potential
480
481    if (ippmod(ieljou).eq.2 .or. ippmod(ieljou).eq.4) then
482
483      call field_get_id('elec_pot_i', f_id)
484
485      inc = 1
486      iprev = 0
487      iccocg = 1
488
489      call field_gradient_scalar(f_id, iprev, 0, inc, iccocg, grad)
490
491      idimt  = 3
492      ientla = .true.
493      ivarpr = .true.
494
495      call post_write_var(nummai, 'Pot_Gradient_Im', idimt, ientla, ivarpr,  &
496                          ntcabs, ttcabs, grad, rbid, rbid)
497
498    endif
499
500    ! For Joule heating by direct conduction:
501    !   imaginary component of the current density
502
503    if (ippmod(ieljou).eq.2 .or. ippmod(ieljou).eq.4) then
504
505      call field_get_id('elec_pot_i', f_id)
506
507      ! As in elflux
508
509      inc = 1
510      iprev = 0
511      iccocg = 1
512
513      call field_gradient_scalar(f_id, iprev, 0, inc, iccocg, grad)
514
515      call field_get_key_int (f_id, kivisl, ifcsii)
516      if (ifcsii .ge. 0) then
517        call field_get_val_s(ifcsii, cvisii)
518        do iloc = 1, ncelps
519          iel = lstcel(iloc)
520          tracel(1 + (iloc-1)*idimt) = -cvisii(iel)*grad(1,iel)
521          tracel(2 + (iloc-1)*idimt) = -cvisii(iel)*grad(2,iel)
522          tracel(3 + (iloc-1)*idimt) = -cvisii(iel)*grad(3,iel)
523        enddo
524      else
525        call field_get_key_double(f_id, kvisl0, visls_0)
526        do iloc = 1, ncelps
527          iel = lstcel(iloc)
528          tracel(1 + (iloc-1)*idimt) = -visls_0*grad(1,iel)
529          tracel(2 + (iloc-1)*idimt) = -visls_0*grad(2,iel)
530          tracel(3 + (iloc-1)*idimt) = -visls_0*grad(3,iel)
531        enddo
532      endif
533
534      idimt  = 3
535      ientla = .true.
536      ivarpr = .false.
537
538      call post_write_var(nummai, 'Current_Im', idimt, ientla, ivarpr,       &
539                          ntcabs, ttcabs, tracel, rbid, rbid)
540
541    endif
542
543    ! Calculation of Module and Argument of the complex potential if IELJOU = 4
544
545    if (ippmod(ieljou).eq.4) then
546
547      ivar = 0
548
549      call field_get_val_s_by_name('elec_pot_r', cpotr)
550      call field_get_val_s_by_name('elec_pot_i', cpoti)
551
552      do iloc = 1, ncelps
553        iel = lstcel(iloc)
554        tracel(iloc) = sqrt(cpotr(iel)*cpotr(iel) + cpoti(iel)*cpoti(iel))
555      enddo
556
557      idimt  = 1
558      ientla = .true.
559      ivarpr = .false.
560
561      call post_write_var(nummai, 'Pot_Module', idimt, ientla, ivarpr,  &
562                          ntcabs, ttcabs, tracel, rbid, rbid)
563
564      do iloc = 1, ncelps
565
566        iel = lstcel(iloc)
567
568        if (cpotr(iel) .ne. 0.d0) then
569          if (cpotr(iel) .ge. 0.d0) then
570            tracel(iloc) = atan(cpoti(iel)/cpotr(iel))
571          else
572            if (cpoti(iel) .gt. 0.d0) then
573              tracel(iloc) = 4.d0*atan(1.d0)                      &
574                             + atan(cpoti(iel) / cpotr(iel))
575            else
576              tracel(iloc) = -4.d0*atan(1.d0)                     &
577                             + atan(cpoti(iel) / cpotr(iel))
578            endif
579          endif
580        else
581          tracel(iloc) = 2.d0*atan(1.d0)
582        endif
583
584        if (tracel(iloc) .lt. 0.d0) then
585          tracel(iloc) = tracel(iloc) + 8.d0**atan(1.d0)
586        endif
587
588      enddo
589
590      idimt  = 1
591      ientla = .true.
592      ivarpr = .false.
593
594      call post_write_var(nummai, 'Pot_Arg', idimt, ientla, ivarpr,  &
595                          ntcabs, ttcabs, tracel, rbid, rbid)
596
597    endif
598
599    ! Free memory
600    deallocate(grad)
601
602  endif
603
604endif
605
606!----
607! End
608!----
609
610return
611end subroutine
612