code-saturne
Equivalent of Courant/Fourier criterion for transported scalars.
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.
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
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
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.
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
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).
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