Open mathrack opened 5 years ago
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
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.