Equivalent of Courant/Fourier criterion for transported scalars.

[mathrack]

Currently, the listing contains the Courant and Fourier numbers for the momentum equation, and the combined Courant/Fourier criterion. It could be interesting to also have this combined Courant/Fourier criterion for transported scalars. This might be especially relevant for buoyant scalars.

[mathrack]

Attached is a first draft of the functionality written on the current development version. It probably does not cover all the options but could be used as a starting point. Feel free to comment it. EDIT: copy-paste patch as attachement is lost

Subject: [PATCH] Print combined Courant/Fourier for transported scalars.

---
 src/base/dttvar.f90 | 101 ++++++++++++++++++++++++++++++++++++++++++++++++++--
 1 file changed, 98 insertions(+), 3 deletions(-)

diff --git a/src/base/dttvar.f90 b/src/base/dttvar.f90
index 7762983..2828346 100644
--- a/src/base/dttvar.f90
+++ b/src/base/dttvar.f90
@@ -102,6 +102,7 @@ double precision ckupdc(6,ncepdp), smacel(ncesmp,nvar)
 ! Local variables
 
 character(len=8) :: cnom
+character(len=80) :: fname
 
 integer          ifac, iel, icfmax(1), icfmin(1), idiff0, iconv0, isym, flid
 integer          modntl
@@ -112,6 +113,7 @@ integer          nswrgp, imligp
 integer          f_id
 integer          nbrval, nclptr
 integer          ntcam1
+integer          ii
 
 double precision epsrgp, climgp, extrap
 double precision cfmax,cfmin, w1min, w2min, w3min
@@ -120,6 +122,7 @@ double precision xyzmax(3), xyzmin(3), vmin(1), vmax(1)
 double precision dtsdtm
 double precision hint
 double precision mult
+double precision prt
 
 double precision, allocatable, dimension(:) :: viscf, viscb
 double precision, allocatable, dimension(:) :: dam
@@ -128,11 +131,11 @@ double precision, allocatable, dimension(:) :: cofbft, coefbt, coefbr
 double precision, dimension(:,:,:), pointer :: coefbv, cofbfv
 double precision, allocatable, dimension(:,:) :: grad
 double precision, allocatable, dimension(:) :: w1, w2, w3, dtsdt0
-double precision, dimension(:), pointer :: imasfl, bmasfl
-double precision, dimension(:), pointer :: brom, crom
+double precision, dimension(:), pointer :: imasfl, bmasfl, imasflt, bmasflt
+double precision, dimension(:), pointer :: brom, crom, cromt
 double precision, dimension(:), pointer :: viscl, visct, cpro_cour, cpro_four
 
-type(var_cal_opt) :: vcopt_u, vcopt_p
+type(var_cal_opt) :: vcopt_u, vcopt_p, vcopt_t
 
 !===============================================================================
 
@@ -939,6 +942,98 @@ if (idtvar.ge.0) then
 
   endif
 
+
+!===============================================================================
+! 4.6 PRINT COMBINED COURANT/FOURIER NUMBER FOR TRANSPORTED SCALARS
+!===============================================================================
+
+  do ii = 1, nscal
+
+    ! Get field parameters
+    call field_get_key_struct_var_cal_opt(ivarfl(isca(ii)), vcopt_t)
+
+    if ( vcopt_t%idiff.ge.1 .or. vcopt_t%iconv.ge.1 ) then
+
+      ! Get mass fluxes
+      call field_get_key_int(ivarfl(isca(ii)), kimasf, iflmas)
+      call field_get_key_int(ivarfl(isca(ii)), kbmasf, iflmab)
+      call field_get_val_s(iflmas, imasflt)
+      call field_get_val_s(iflmab, bmasflt)
+
+      ! Get density
+      call field_get_key_int(ivarfl(isca(ii)), kromsl, flid)
+      if (flid.gt.-1) then
+        call field_get_val_s(flid, cromt)
+      else
+        call field_get_val_s(icrom, cromt)
+      endif
+
+      ! Compute diffusivity
+      if (vcopt_t%idiff.ge.1) then
+
+        call field_get_key_int(ivarfl(isca(ii)), kivisl, flid)
+        if (flid.gt.-1) then
+          call field_get_val_s(flid, viscl)
+        else
+          call field_get_val_s(iviscl, viscl)
+        endif
+        call field_get_key_double(ivarfl(isca(ii)), ksigmas, prt)
+        do iel = 1, ncel
+          w1(iel) = viscl(iel) + vcopt_t%idifft*visct(iel)/prt
+        enddo
+        call viscfa(imvisf, w1, viscf, viscb)
+
+      else
+
+        do ifac = 1, nfac
+          viscf(ifac) = 0.d0
+        enddo
+        do ifac = 1, nfabor
+          viscb(ifac) = 0.d0
+        enddo
+
+      endif
+
+      ! Boundary conditions for matrdt
+      do ifac = 1, nfabor
+        if (bmasflt(ifac).lt.0.d0) then
+
+          iel = ifabor(ifac)
+          hint = vcopt_t%idiff*( viscl(iel)                                 &
+                               + vcopt_t%idifft*visct(iel)/prt)/distb(ifac)
+          coefbt(ifac) = 0.d0
+          cofbft(ifac) = hint
+
+        else
+
+          coefbt(ifac) = 1.d0
+          cofbft(ifac) = 0.d0
+
+        endif
+      enddo
+
+      ! MATRICE A PRIORI NON SYMETRIQUE
+      isym = 1
+      if (vcopt_t%iconv.gt.0) isym = 2
+
+      call matrdt &
+      !==========
+ (vcopt_t%iconv, vcopt_t%idiff, isym, coefbt, cofbft, imasflt, bmasflt, &
+  viscf, viscb, dam)
+
+      do iel = 1, ncel
+        w1(iel) = dam(iel)/(cromt(iel)*volume(iel))
+      enddo
+
+      call field_get_name(ivarfl(isca(ii)), fname)
+      call log_iteration_add_array(fname(1:15), 'criterion',    &
+                                   MESH_LOCATION_CELLS, .true.,       &
+                                   1, w2)
+
+    endif
+
+  enddo
+
 !===============================================================================
 ! 5.   ALGORITHME STATIONNAIRE
 !===============================================================================
-- 
2.7.4

[mathrack]

Please note that a loop is missing: at the end, after computing w1 and before calling log_iteration_add_arry, please add this

    do iel = 1, ncel
      w2(iel) = w1(iel)*dt(iel)
    enddo

[YvanFournier]

Hello,

As the repository is a mirror, we'll need to check how to apply it, but a pull request could be interesting here (at least we can learn to better use this type of workflow). We can also check which extension will work here adding a patch generated by "git format-patch".

Since you have several scalars, An extension to what you do here would be to use a logic similar to scalar variable diffusivity (where we add a dependent field to the considered scalars, using a dedicated field keyword indicating the relation between a given scalar and its associated dependent field) to save the values equivalent to the Fourier number for selected scalars. This would allow not only logging min/max/mean values, but post-processing the associated field, just like we can do with the CFL and Fourier numbers., and would make the feature more "consistent" and complete.

This can be done in a second phase.

Could you try re-posting the patch as an addition or generate a merge request ? I can work with your pasted code, but it is a good time to try to improve our workflow with git.

[mathrack]

Hello Yvan,

I have send a pull request as '.patch' files are not supported as attachments. It should correspond to the first phase. I have not inserted the code showing the location of the minimum / maximum, it already appears several times in the subroutine and might be moved to a dedicated subroutine?

P.S.: in order to generate a merge request, one should fork, patch the fork, then a pull request can be generated from the fork.
P.S.2: it is probably more flexible to use the forum of Code_Saturne to exchange '.patch' files and discuss them

Best regards,
Cédric

[YvanFournier]

Hello,

I have pulled to code from you pull request, and am reviewing it. As this adds additional operations for logging, we will probably add a keyword to activate this only upon request.

I will check with the rest of the team whether a separate keyword, a var_cal_opt member, or an additional flag in a field's postprocessing option is best (I tend to prefer the latter option).

[mathrack]

Hello Yvan,

Thank you for the follow-up. I am not fully sure about the 'boundary conditions for matrdt' part of my patch, please double check it. Otherwise it should be functional.

Best regards,
Cédric