diff --git a/NEWS.txt b/NEWS.txt index 190511a2c6f4015faaacac863f4d2857254938e7..0f5cc967810e7b1336ca993ad7161413a5c649ad 100644 --- a/NEWS.txt +++ b/NEWS.txt @@ -1,6 +1,15 @@ Latest changes ============== +TELEMAC-2D: better time step managment in the finite volume kernel. +- allowed constant time step +- allowed variable time step with CFL > 0.9 (before: was forced to 0.9); +- warning if CFL > 1 with advised time step; +- correction of the second order time step limitation; +- Add new keyword ADAPTATION OF TIME STEP TO GRAPHIC PRINTOUTS, + If NO results at tn+1 are written if tn<trec<tn+1; + If YES the time step dtn is adapted so that tn+dtn=trec (old behavior); + TELEMAC-2D/TELEMAC-3D/TOMAWAC/GAIA: add Nearest Neighbor search method for Time Series files diff --git a/documentation/telemac2d/user/latex/chap7.tex b/documentation/telemac2d/user/latex/chap7.tex index cf78a00da9ba435d1cbbe5cb16fe8d4f75085b5b..01b8b0efd0a5eaa5f58addae0ef8bc900d5dee19 100644 --- a/documentation/telemac2d/user/latex/chap7.tex +++ b/documentation/telemac2d/user/latex/chap7.tex @@ -114,6 +114,8 @@ linearization is to be performed, by using the keyword \section{Numerical schemes} \subsection{Finite elements} +\label{FE_schemes} + Finite elements resolution is based on the primitive equations. It is possible to replace the original equations by a generalized wave equation obtained by eliminating the velocity from the continuity equation using a value @@ -463,6 +465,7 @@ The bigger the number is, the more conservative the scheme is, but the higher th computational costs are. \subsection{Finite volumes} +\label{FV_schemes} When using the finite volume scheme, the primitive equations written in conservative form are solved. The finite volume schemes use control volumes @@ -505,12 +508,9 @@ The algorithm of finite volume schemes is explicit which means that the Courant number must be inferior or equal to 1 to ensure stability. It is however recommended to set the keyword \telkey{DESIRED COURANT NUMBER} to 0.9. This should be combined with the keyword \telkey{VARIABLE TIME-STEP} -set to YES (default: NO). -\telemac{2D} then adjusts the calculation time step so as to satisfy +set to YES (default: NO), in which case +\telemac{2D} adjusts the calculation time step so as to satisfy this Courant number criterion. -The graphic printout is handled with the time step given in the steering file. -However, it should be noted that this leads to irregular sampling -from the control listing. The keyword \telkey{OPTION OF THE HYDROSTATIC RECONSTRUCTION} enables to choose the option for hydrostatic reconstruction for the @@ -892,29 +892,43 @@ It is not activated with the wave equation (\telkey{TREATMENT OF THE LINEAR SYSTEM} = 2). -\section{Courant number management} +\section{Time step management and Courant number} -During a model simulation, the Courant number value (number of grid cells +During a simulation, the Courant number value (number of grid cells crossed by a water particle during a time step) considerably influences -the quality of the results. -Irrespective of numerical schemes with a stability condition on the Courant -number, experience shows that result quality decreases if the Courant number is -above 7 or 8. -Yet it is not so easy to estimate the value of the Courant number -- especially in sea models with a large tidal range. -To help, \telemac{2D} allows the user to check the Courant number during -computation: -the software automatically executes intermediate time steps -so that the Courant number keeps below a given value. - -This function is activated using the keyword \telkey{VARIABLE TIME-STEP} -(default value = NO) and the maximum Courant number value can be specified using -the keyword \telkey{DESIRED COURANT NUMBER} (default value = 1). - -It should be stressed that when a variable time step is used, +the quality of the results. As a consequence, the determination of the time step +must be done with great caution. + +\subsection{Finite elements} + +With \telkey{EQUATIONS='SAINT-VENANT FE'}, time step is constant and determined by +the user via the \telkey{TIME STEP} keyword. Some FE numerical schemes (distributive schemes) +perform sub-iterations in some cases to check numerical properties, +such as the positivity of the water depth. See the section \ref{FE_schemes} for more info. +However there are no guarantee that the CFL condition of the chosen numerical scheme is fullfilled. + +\subsection{Finite volumes} + +With \telkey{EQUATIONS='SAINT-VENANT FV'}, time step can either be defined as a constant (as with the +FE schemes, althought not recommended) or variable by using \telkey{VARIABLE TIME-STEP}=YES +(default value = NO). +In this case the time step is computed with the CFL condition and +the Courant number value can be specified using +the keyword \telkey{DESIRED COURANT NUMBER} (default value = 1.). + +\ + +It should be stressed that when \telkey{VARIABLE TIME-STEP}=YES, sampling from the results file and control listing is no longer regular in time, as it depends directly on the time step value. - +Note that specifying a \telkey{TIME STEP} is still usefull as +the graphic printout is handled with the time step given in the steering file. +The keyword \telkey{ADAPTATION OF TIME STEP TO GRAPHIC PRINTOUTS} (default value = YES) +is used to specify the way the result are recorded. +If set to YES the time step is adjusted to graphic printouts, in +which case $dt^n$ is modified so that $t^{n}+dt^n = t^{n+1} = t^{rec}$. This option is not +recommended as the results may varry depending on the period of graphic printouts. +If set to NO for a given record time $t^{rec}$, results at $t^{n+1}$ are written if $t^{n} < t^{rec} \leq t^{n+1}$. \section{Tidal flats} \textbf{This section mainly concerns the Finite Elements method.} diff --git a/examples/telemac2d/pluie/f2d_rain_CN_fv.slf b/examples/telemac2d/pluie/f2d_rain_CN_fv.slf index 7c252bacf7a48f22d2eecf8b82bb2d7a9c8466ad..834a57f16b8a09dae1dad6a06b29001a55efcaa0 100644 --- a/examples/telemac2d/pluie/f2d_rain_CN_fv.slf +++ b/examples/telemac2d/pluie/f2d_rain_CN_fv.slf @@ -1,3 +1,3 @@ version https://git-lfs.github.com/spec/v1 -oid sha256:eed4042fe92ddddd1f0aba8845f73c2fcb99b9ac064aa0594787b627eb4f095b +oid sha256:a3026280421b027ba3ff36b293c797d390d52908af580c77c76f69e06da00cef size 212368 diff --git a/examples/telemac2d/pluie/f2d_rain_uniform_fc_fv_h.slf b/examples/telemac2d/pluie/f2d_rain_uniform_fc_fv_h.slf index 14a5eca14e0cc7d44a93eb48ff3a4dc666449807..d028cac42c49b0fad8c3ec6aba3c9ae5a24639ca 100644 --- a/examples/telemac2d/pluie/f2d_rain_uniform_fc_fv_h.slf +++ b/examples/telemac2d/pluie/f2d_rain_uniform_fc_fv_h.slf @@ -1,3 +1,3 @@ version https://git-lfs.github.com/spec/v1 -oid sha256:b3b96b13754fa05283070e6692b73dee88bd1939331f0989dec67af344dad0f2 +oid sha256:bf6d8090775f4920d83569cb5b3035cdd91ac2df733931ea99d373328d65c8f0 size 212520 diff --git a/examples/telemac2d/pluie/f2d_rain_uniform_ks_fv_ga.slf b/examples/telemac2d/pluie/f2d_rain_uniform_ks_fv_ga.slf index 1b43873c0816b911a8cb2ce825be103bbb28bc21..6adb1904a999316037bc3b75002460e421f78559 100644 --- a/examples/telemac2d/pluie/f2d_rain_uniform_ks_fv_ga.slf +++ b/examples/telemac2d/pluie/f2d_rain_uniform_ks_fv_ga.slf @@ -1,3 +1,3 @@ version https://git-lfs.github.com/spec/v1 -oid sha256:66be10e48c2bf74aa22db0af858d7f35a2abd9a6950bb81d9312e4029eea2275 +oid sha256:f10d05406a27199c4ce902c1417928e962c4c4d662e35169c60959fdcb3f3f27 size 212520 diff --git a/examples/telemac2d/pluie/vnv_pluie_fv.py b/examples/telemac2d/pluie/vnv_pluie_fv.py index 43e107dacf16de575987e90958cc6b4678793d4b..157e26e032ce019bbb5e803ed179fd2630044f12 100644 --- a/examples/telemac2d/pluie/vnv_pluie_fv.py +++ b/examples/telemac2d/pluie/vnv_pluie_fv.py @@ -76,34 +76,20 @@ class VnvStudy(AbstractVnvStudy): Post-treatment processes """ # Comparison with the last time frame of the reference file. - self.check_epsilons('vnv_1:T2DRES', - 'f2d_rain_uniform_fv.slf', - eps=[1.E-15]) - + self.check_epsilons('vnv_1:T2DRES', 'f2d_rain_uniform_fv.slf', eps=[1.E-8]) # Comparison seq/par. - self.check_epsilons('vnv_1:T2DRES', - 'vnv_2:T2DRES', - eps=[1.E-15]) + self.check_epsilons('vnv_1:T2DRES', 'vnv_2:T2DRES', eps=[1.E-8]) # Comparison with the last time frame of the reference file. - self.check_epsilons('vnv_3:T2DRES', - 'f2d_rain_CN_fv.slf', - eps=[1.E-15]) - + self.check_epsilons('vnv_3:T2DRES', 'f2d_rain_CN_fv.slf', eps=[1.E-8]) # Comparison seq/par. - self.check_epsilons('vnv_3:T2DRES', - 'vnv_4:T2DRES', - eps=[1.E-15]) + self.check_epsilons('vnv_3:T2DRES', 'vnv_4:T2DRES', eps=[1.E-8]) # Comparison with the last time frame of the reference file. - self.check_epsilons('vnv_5:T2DRES', - 'f2d_rain_CN_hyetograph_fv.slf', - eps=[1.E-11]) - + self.check_epsilons('vnv_5:T2DRES', 'f2d_rain_CN_hyetograph_fv.slf', eps=[1.E-4]) # Comparison seq/par. - self.check_epsilons('vnv_5:T2DRES', - 'vnv_6:T2DRES', - eps=[1.E-11]) + self.check_epsilons('vnv_5:T2DRES', 'vnv_6:T2DRES', eps=[1.E-4]) + def _post(self): """ diff --git a/examples/telemac2d/pluie/vnv_pluie_fv_ga.py b/examples/telemac2d/pluie/vnv_pluie_fv_ga.py index a30ec468c4613a081eac76a8528d98a517dcde4e..af23685adcbd8f0eb8da7d4466fefd63d68b42ce 100644 --- a/examples/telemac2d/pluie/vnv_pluie_fv_ga.py +++ b/examples/telemac2d/pluie/vnv_pluie_fv_ga.py @@ -76,34 +76,19 @@ class VnvStudy(AbstractVnvStudy): Post-treatment processes """ # Comparison with the last time frame of the reference file. - self.check_epsilons('vnv_1:T2DRES', - 'f2d_rain_uniform_fv_ga.slf', - eps=[1.E-15]) - + self.check_epsilons('vnv_1:T2DRES', 'f2d_rain_uniform_fv_ga.slf', eps=[1.E-8]) # Comparison seq/par. - self.check_epsilons('vnv_1:T2DRES', - 'vnv_2:T2DRES', - eps=[1.E-15]) + self.check_epsilons('vnv_1:T2DRES', 'vnv_2:T2DRES', eps=[1.E-8]) # Comparison with the last time frame of the reference file. - self.check_epsilons('vnv_3:T2DRES', - 'f2d_rain_uniform_ks_fv_ga.slf', - eps=[1.E-12]) - + self.check_epsilons('vnv_3:T2DRES', 'f2d_rain_uniform_ks_fv_ga.slf', eps=[1.E-8]) # Comparison seq/par. - self.check_epsilons('vnv_3:T2DRES', - 'vnv_4:T2DRES', - eps=[1.E-12]) + self.check_epsilons('vnv_3:T2DRES', 'vnv_4:T2DRES', eps=[1.E-8]) # Comparison with the last time frame of the reference file. - self.check_epsilons('vnv_5:T2DRES', - 'f2d_hyetograph_ks_geo_fv_ga.slf', - eps=[1.E-13]) - + self.check_epsilons('vnv_5:T2DRES', 'f2d_hyetograph_ks_geo_fv_ga.slf', eps=[1.E-4]) # Comparison seq/par. - self.check_epsilons('vnv_5:T2DRES', - 'vnv_6:T2DRES', - eps=[1.E-12]) + self.check_epsilons('vnv_5:T2DRES', 'vnv_6:T2DRES', eps=[1.E-4]) def _post(self): diff --git a/examples/telemac2d/pluie/vnv_pluie_fv_h.py b/examples/telemac2d/pluie/vnv_pluie_fv_h.py index 27be4f424f9b340f97c1a6d6e5cab15c5c37b238..cbdda7768a736cdd38d57ace0512f74d17a61bf2 100644 --- a/examples/telemac2d/pluie/vnv_pluie_fv_h.py +++ b/examples/telemac2d/pluie/vnv_pluie_fv_h.py @@ -76,34 +76,19 @@ class VnvStudy(AbstractVnvStudy): Post-treatment processes """ # Comparison with the last time frame of the reference file. - self.check_epsilons('vnv_1:T2DRES', - 'f2d_rain_uniform_fv_h.slf', - eps=[1.E-15]) - + self.check_epsilons('vnv_1:T2DRES', 'f2d_rain_uniform_fv_h.slf', eps=[1.E-8]) # Comparison seq/par. - self.check_epsilons('vnv_1:T2DRES', - 'vnv_2:T2DRES', - eps=[1.E-15]) + self.check_epsilons('vnv_1:T2DRES', 'vnv_2:T2DRES', eps=[1.E-8]) # Comparison with the last time frame of the reference file. - self.check_epsilons('vnv_3:T2DRES', - 'f2d_rain_uniform_fc_fv_h.slf', - eps=[1.E-12]) - + self.check_epsilons('vnv_3:T2DRES', 'f2d_rain_uniform_fc_fv_h.slf', eps=[1.E-8]) # Comparison seq/par. - self.check_epsilons('vnv_3:T2DRES', - 'vnv_4:T2DRES', - eps=[1.E-12]) + self.check_epsilons('vnv_3:T2DRES', 'vnv_4:T2DRES', eps=[1.E-8]) # Comparison with the last time frame of the reference file. - self.check_epsilons('vnv_5:T2DRES', - 'f2d_hyetograph_fc_geo_fv_h.slf', - eps=[1.E-13]) - + self.check_epsilons('vnv_5:T2DRES', 'f2d_hyetograph_fc_geo_fv_h.slf', eps=[1.E-4]) # Comparison seq/par. - self.check_epsilons('vnv_5:T2DRES', - 'vnv_6:T2DRES', - eps=[1.E-12]) + self.check_epsilons('vnv_5:T2DRES', 'vnv_6:T2DRES', eps=[1.E-4]) def _post(self): diff --git a/sources/telemac2d/caldt.f b/sources/telemac2d/caldt.f index 5b6b90d5610e02756c4f50b538fd39dfd30e3cd2..f29df313c3b6ba26116543264c26058502e6979d 100644 --- a/sources/telemac2d/caldt.f +++ b/sources/telemac2d/caldt.f @@ -57,11 +57,12 @@ !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! USE BIEF_DEF, ONLY:NCSIZE - USE DECLARATIONS_TELEMAC2D,ONLY:DTINI,LEOPRD,PTINIG,GRAV,TORDER, + USE DECLARATIONS_TELEMAC2D,ONLY:LEOPRD,PTINIG,GRAV,TORDER, & LEO_TRAC,NPOIN,ICIN,SORDER,LT,HN, & U,V,CFLWTD,AT,TMAX,MESH,DTVARI, & ENTET,DTCAS,LIMPRO,IKLE,MVISUV, - & VISC,NTRAC,DIFNU,MVIST,DTMIN_FV + & VISC,NTRAC,DIFNU,MVIST,DTMIN_FV, + & EPS_DIV0,ADAPT_DT_TO_LEOPRD USE INTERFACE_TELEMAC2D, EX_CALDT => CALDT USE DECLARATIONS_SPECIAL USE INTERFACE_PARALLEL, ONLY : P_MIN @@ -76,72 +77,54 @@ ! INTEGER IS,ITRAC LOGICAL DEJA,THEEND - DOUBLE PRECISION RA3,EPSL,SIGMAX,UA2,UA3,UNORM,DTT - DOUBLE PRECISION RESTE,GPRDTIME,DT2 + DOUBLE PRECISION RA3,SIGMAX,UA2,UA3,UNORM,DTT,DTV + DOUBLE PRECISION RESTE0,RESTE1,GPRDTIME,DT2 ! INTRINSIC MIN,CEILING ! !----------------------------------------------------------------------- ! +! INIT OF GRAPHICAL PRINTOUTS PARAMETERS THEEND=.FALSE. IF(LEO.OR.LT.EQ.1) LEO_TRAC = .FALSE. LEO =.FALSE. -! +++++++++++++++++++++++++ -! VARIABLE TIME STEP - IF(DTVARI)THEN -! +++++++++++++++++++++++++ +! +! GRAPHICAL PRINTOUTS TIME DELTA + GPRDTIME = LEOPRD*DTCAS +! +! ****************************************************************** +! COMPUTE TIME STEP +! ****************************************************************** DEJA=.FALSE. DT = DTCAS - EPSL = 0.01D0 -! -! KINETIC SCHEME -! ============== - IF(ICIN.EQ.1)THEN - - RA3 = SQRT(1.5D0*GRAV) - DO IS=1,NPOIN - IF(HN%R(IS).LT.0.D0.AND.ENTET.AND..NOT.DEJA)THEN - WRITE(LU,*) 'CALDT WARNING : NEGATIVE WATER DEPTH' - WRITE(LU,*) ' SEE NODE:',IS,' FOR EXAMPLE' - DEJA = .TRUE. - ELSE - SIGMAX = ABS(HN%R(IS)) - UA2 = U%R(IS) - UA3 = V%R(IS) - UNORM=SQRT(UA2*UA2 + UA3*UA3) - SIGMAX= MAX(EPSL, RA3*SQRT(SIGMAX) +UNORM ) - DT = MIN(DT, CFLWTD*MESH%DTHAUT%R(IS)/SIGMAX) - ENDIF - ENDDO ! -! SCHEMES OF ROE, ZOKAGOA, TCHAMEN, HLLC AND WAF -! ============================================== - ELSEIF(ICIN.EQ.0.OR.ICIN.EQ.2.OR.ICIN.EQ.3.OR. - & ICIN.EQ.4.OR.ICIN.EQ.5) THEN -! - DO IS=1,NPOIN - IF(HN%R(IS).LT.0.D0.AND.ENTET.AND..NOT.DEJA)THEN - WRITE(LU,*) 'CALDT WARNING : NEGATIVE WATER DEPTH' - WRITE(LU,*) ' SEE NODE:',IS,' FOR EXAMPLE' - DEJA = .TRUE. +! COMPUTE DT FOR ADVECTION (CFL CONDITION) + DO IS=1,NPOIN + IF(HN%R(IS).LT.0.D0.AND.ENTET.AND..NOT.DEJA)THEN + WRITE(LU,*) 'CALDT WARNING : NEGATIVE WATER DEPTH' + WRITE(LU,*) ' SEE NODE:',IS,' FOR EXAMPLE' + DEJA = .TRUE. + ELSE + UA2 = U%R(IS) + UA3 = V%R(IS) + UNORM=SQRT(UA2*UA2 + UA3*UA3) +! KINETIC SCHEME + IF(ICIN.EQ.1)THEN + RA3 = SQRT(1.5D0*GRAV) + SIGMAX= MAX(EPS_DIV0, RA3*SQRT(ABS(HN%R(IS))) + UNORM) +! OTHER SCHEMES ELSE - UA2 = U%R(IS) - UA3 = V%R(IS) - UNORM=SQRT(UA2*UA2 + UA3*UA3) - SIGMAX= MAX(EPSL,SQRT(GRAV*ABS(HN%R(IS)))+UNORM) -! DTHAUT=|Ci|/Sum(Lij) - DT = MIN(DT, CFLWTD*MESH%DTHAUT%R(IS)/SIGMAX) - ENDIF - ENDDO -! - ELSE - WRITE(LU,4020) ICIN -4020 FORMAT(1X,'CALDT: ERROR IN THE CHOICE OF OPTFV: ',1I6) - CALL PLANTE(1) - STOP - ENDIF + SIGMAX= MAX(EPS_DIV0, SQRT(GRAV*ABS(HN%R(IS))) + UNORM) + ENDIF + DT = MIN(DT, CFLWTD*MESH%DTHAUT%R(IS)/SIGMAX) +! DTHAUT=|Ci|/Sum(Lij) + ENDIF + ENDDO ! +! UPDATE DT FROM BOUNDARY CONDITIONS CALL CDL_DT(DT,LIMPRO%I) + +! UPDATE DT FROM DIFFUSION (FOURIER CONDITION) CALL DIFF_DT(DT,IKLE%I,MESH%NUBO%I,MESH%ELTSEG%I,MESH%IFABOR%I, & MESH%VNOIN%R,MESH%COORDR%R,MVISUV,VISC%R) IF(NTRAC.GT.0) THEN @@ -152,34 +135,37 @@ & DIFNU(ITRAC)) ENDDO ENDIF -! + ! FOR NEWMARK SCHEME -! ================== IF (TORDER.EQ.2) THEN - IF (LT.EQ.1) THEN - DTN = DT - ELSE + IF (LT.GT.1) THEN DT2 = 0.5D0*MIN(DT, DTN) - DTN = DT DT = DT2 ENDIF ENDIF ! -! +++++++++++++++++++++++++ +! FOR PARALLELISM + IF(NCSIZE.GT.1) THEN + DT = P_MIN(DT) ENDIF -! +++++++++++++++++++++++++ ! -! FOR PARALLELISM - IF(NCSIZE.GT.1) DT=P_MIN(DT) +! ****************************************************************** +! CONSTANT SET TIME STEP +! ****************************************************************** + IF(DTVARI.EQV..FALSE.) THEN + DTV = DT ! STORE DT FOR RECOMMENDATION IN WARNING + DT = DTCAS + ENDIF ! +! ****************************************************************** +! END OF COMPUTATION AND OTHER SANITY CHECKS +! ****************************************************************** IF(DTVARI) THEN -! - IF(TMAX.LT.DT)DT=TMAX !REALLY CRAZY CASES - DTT=TMAX-AT - IF(CFLWTD.GE.1.D0) DT=0.9D0*DT/CFLWTD + IF(TMAX.LT.DT) DT=TMAX !REALLY CRAZY CASES + DTT = TMAX-AT IF(DTT.LE.DT.AND.DTT.GT.0.D0) THEN !LAST TIME STEP - DT=DTT -! END OF COMUTATION WILL BE DETECTED IN RESOLU FOR ORDER 2 + DT = DTT +! END OF COMUTATION WILL BE DETECTED IN SECOND_ORDER FOR ORDER 2 IF(SORDER.EQ.1)THEEND=.TRUE. ENDIF IF(AT.GT.TMAX) THEN @@ -188,73 +174,100 @@ CALL PLANTE(1) STOP ENDIF -! - ELSE -! -! DT NOT VARIABLE - WRITE(LU,*) 'ERROR: FIXED TIME-STEP AND CFL NOT GIVEN!...! ' - WRITE(LU,*) 'TIME-STEP MAY NOT SATISFY CFL CONDITION ' - WRITE(LU,*) 'THIS OPTION IS NOT ACCEPTED WITH FINITE VOLUMES ' + ENDIF +! GPRDTIME TOO SMALL + IF(GPRDTIME.LT.1.D-12)THEN + WRITE(LU,*) 'CALDT: PROBLEM WITH PARAMETERS: DTCAS,LEOPRD', + & DTCAS,LEOPRD CALL PLANTE(1) STOP -! + ENDIF +! CASE WHERE TIME STEP IS BIGGER THEN GRAPHIC OUTPUT (GPRDTIME) + IF(GPRDTIME.LT.DT)THEN + DT = DTCAS + IF(ENTET)THEN + WRITE(LU,*) 'WARNING: GRAPHICAL OUTPUT NOT OPTIMIZED: ' + WRITE(LU,*) ' - INITIAL TIME STEP TOO SMALL ' + WRITE(LU,*) ' - AND/OR PERIOD OF GRAPHIC OUTPUT TOO SMALL' + WRITE(LU,*) 'THIS COULD REDUCE COMPUTATION TIME-STEP' + WRITE(LU,*) 'AND INCREASE CPU TIME' + ENDIF ENDIF ! -!************************************************************************* -! GRAPHIC OUTPUTS LIKE ASKED BY USER: DT ADAPTATION -!************************************************************************* - GPRDTIME=LEOPRD*DTINI - IF(GPRDTIME.LT.1.D-12)THEN - WRITE(LU,*) 'CALDT: PROBLEM WITH PARAMETERS: DTINI,LEOPRD', - & DTINI,LEOPRD - CALL PLANTE(1) - STOP -! CASE WHERE TIME STEP IS BIGGER THEN GRAPHIC OUTPUT (GPRDTIME) - ELSEIF(GPRDTIME.LT.DT)THEN - DT=DTINI - IF(ENTET)THEN - WRITE(LU,*) 'WARNING: GRAPHICAL OUTPUT NOT OPTIMIZED: ' - WRITE(LU,*) ' - INITIAL TIME STEP TOO SMALL ' - WRITE(LU,*) ' - AND/OR PERIOD OF GRAPHIC OUTPUT TOO SMALL' - WRITE(LU,*) 'THIS COULD REDUCE COMPUTATION TIME-STEP' - WRITE(LU,*) 'AND INCREASE CPU TIME' +! ****************************************************************** +! DETECTION OF GRAPHICAL PRINTOUTS (LEO = TRUE) +! ****************************************************************** +! +! ONLY FOR ORDER 1, ORDER 2 IS IN SECOND_ORDER + IF (SORDER.EQ.1) THEN +! +! VARIABLE DT + IF(DTVARI) THEN +! +! ADAPT DT TO TAKE INTO ACCOUNT GRAPHIC OUTPUT + IF (ADAPT_DT_TO_LEOPRD) THEN + IS = CEILING(AT/GPRDTIME) + RESTE1 = IS*GPRDTIME-AT + IF(THEEND.OR.(RESTE1.LE.DT + & .AND.RESTE1.GT.DTMIN_FV + & .AND.LT.GT.PTINIG)) THEN +! HERE THERE IS GRAPHICAL OUTPUTY + LEO = .TRUE. + DT = MIN(RESTE1,DT) + ENDIF +! +! DETECT CLOSEST GRAPHIC PRINTOUT TO CURRENT TIME (AT+DT) +! (DETECTION OF LEO AT TIME AT, BUT WRITES RESULTS AT TIME AT+DT) + ELSE + IS = CEILING((AT+DT)/GPRDTIME) + RESTE0 = AT+DT - (IS-1)*GPRDTIME + IF (LT.EQ.1) THEN + IF(THEEND.OR.(RESTE0.LT.DT .AND.LT.GT.PTINIG)) THEN + LEO = .TRUE. + ENDIF + ELSE + IF(THEEND.OR.(RESTE0.LE.DT.AND.LT.GT.PTINIG)) THEN + LEO = .TRUE. + ENDIF + ENDIF ENDIF - ENDIF -! ADAPT DT TO TAKE INTO ACCOUNT GRAPHIC OUTPUT -! ONLY FOR ORDER 1, ORDER 2 WILL BE DONE IN SECOND_ORDER - IF(SORDER.EQ.1)THEN - IS=CEILING(AT/GPRDTIME) - RESTE=IS*GPRDTIME-AT -! - IF(THEEND.OR. - & (RESTE.LE.DT.AND.RESTE.GT.DTMIN_FV.AND.LT.GT.PTINIG))THEN -! HERE THERE IS GRAPHICAL OUTPUT +! +! CONSTANT DT + ELSE + IF (THEEND.OR.(MOD(LT,LEOPRD)==0)) THEN LEO = .TRUE. -! IF RESTE < DTMIN_VF : NO MODIFICATION OF DT -! TODO: DO NOT MODIFY DT BY DEFAULT: ADD INTERPOLATION INSTEAD -! TODO: ADD NEW KEYWORD FOR ADAPTATION OF DT (TO KEEP THIS FUNCTIONALITY) -! NOTE: SAME CODE IN SECOND ORDER - DT=MIN(RESTE,DT) ENDIF ENDIF ! + ENDIF ! SORDER1 ! -!************************************************************************* -! ERROR TREATMENT AND LISTING OUTPUTS -!************************************************************************* -! -! CASE DT <=0 +! ****************************************************************** +! ERRORS AND WARNINGS +! ****************************************************************** ! +! CASE DT <= 0 IF(DT.LE.0.D0) THEN - WRITE(LU,*) 'NEGATIVE (OR NIL) TIME-STEP: ',DT + WRITE(LU,*) 'NEGATIVE (OR NIL) TIME-STEP: ', DT CALL PLANTE(1) STOP ENDIF ! +! CASE CFL >= 1 IF(CFLWTD.GE.1.D0) THEN IF(ENTET) THEN - WRITE(LU,*) 'WARNING: CFL NOT GIVEN OR >1 !...! ' - WRITE(LU,*) 'TIME-STEP (WITH CFL = 0.9): ',DT + WRITE(LU,*) 'WARNING: CFL NOT GIVEN OR >=1 !...! ' + WRITE(LU,*) 'RECOMMENDED TIME-STEP (WITH CFL = 0.9): ', + & 0.9D0*DT/CFLWTD + ENDIF + ENDIF +! +! CASE CONSTANT TIME STEP + IF((DTVARI.EQV..FALSE.).AND.(SORDER.EQ.1)) THEN + IF(ENTET) THEN + WRITE(LU,*) 'WARNING: FIXED TIME-STEP (NOT RECOMMENDED)' + WRITE(LU,*) 'TIME-STEP MAY NOT SATISFY CFL CONDITION, ' + WRITE(LU,*) 'RECOMMENDED TIME-STEP (WITH CFL = 0.9): ', + & 0.9D0*DTV/CFLWTD ENDIF ENDIF ! diff --git a/sources/telemac2d/declarations_telemac2d.f b/sources/telemac2d/declarations_telemac2d.f index 6965e104a4621fd9d48551ac96749f7299f1f99e..21df9776b052cb2d304eab7ffa40da7e9262d13d 100644 --- a/sources/telemac2d/declarations_telemac2d.f +++ b/sources/telemac2d/declarations_telemac2d.f @@ -1860,12 +1860,20 @@ ! LOGICAL VIT_IN_T2DBB2 ! +! ADAPT FINITE VOLUME TIME STEP TO GRAPHICAL PRINTOUTS +! + LOGICAL ADAPT_DT_TO_LEOPRD + !----------------------------------------------------------------------- ! ! 6) REALS ! !----------------------------------------------------------------------- ! +! EPSILON TO AVOID DIVISION BY ZERO +! + DOUBLE PRECISION, PARAMETER :: EPS_DIV0 = 1.E-16 +! ! TIME STEP ! DOUBLE PRECISION, TARGET :: DT, DTCAS, DTN @@ -1874,10 +1882,6 @@ ! DOUBLE PRECISION, TARGET :: TMAX ! -! INITIAL TIME STEP -! - DOUBLE PRECISION DTINI -! ! MINIMUM TIME STEP IN FV TIME STEP ADAPTATION FOR GRAPHIC PRINTOUTS ! DOUBLE PRECISION, PARAMETER :: DTMIN_FV = 1.D-12 diff --git a/sources/telemac2d/init_fv.f b/sources/telemac2d/init_fv.f index 11b4de73ebee8ec8c137f2ae94e81211c1374765..e8ad26c87a797b7791891945c1e33140bd7e025a 100644 --- a/sources/telemac2d/init_fv.f +++ b/sources/telemac2d/init_fv.f @@ -36,7 +36,7 @@ USE BIEF USE INTERFACE_TELEMAC2D, EX_INIT_FV => INIT_FV USE DECLARATIONS_TELEMAC2D, ONLY: DIFVIT,LEOPRD,NPOIN,ITURB, - & NELMAX,IKLE,NELEM,DTINI,U,V,HN,X,Y, + & NELMAX,IKLE,NELEM,DTCAS,U,V,HN,X,Y, & AT,LT,SORDER,MXPTVS,NSEG,VERTIC ! !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+ @@ -62,7 +62,23 @@ ! !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+ ! - GPRDTIME=LEOPRD*DTINI +! +! ================================================================== +! CHECK SOLVER +! ================================================================== +! + IF((ICIN.GT.5).OR.(ICIN.LT.0)) THEN + WRITE(LU,4020) ICIN +4020 FORMAT(1X,'CALDT: ERROR IN THE CHOICE OF OPTFV: ',1I6) + CALL PLANTE(1) + STOP + ENDIF +! +! ================================================================== +! TIME DELTA OF GRAPHIC PRINTOUTS +! ================================================================== +! + GPRDTIME=LEOPRD*DTCAS ! ! ================================================================== ! INITIALIZE QU AND QV AND FLUX_OLD diff --git a/sources/telemac2d/lecdon_telemac2d.f b/sources/telemac2d/lecdon_telemac2d.f index 5c56a35971897083b3a159d5ad150e81bf9df8af..8f33a1dfb17856734090db0134c2958c2b05466a 100644 --- a/sources/telemac2d/lecdon_telemac2d.f +++ b/sources/telemac2d/lecdon_telemac2d.f @@ -1818,6 +1818,9 @@ MAX_PREV = MOTLOG(ADRESS(3,66)) ! VELOCITIES IN BINARY DATABASE 2 FOR TIDE VIT_IN_T2DBB2 = MOTLOG(ADRESS(3,67)) +! ADAPT FINITE VOLUME TIME STEP TO GRAPHICAL PRINTOUTS + ADAPT_DT_TO_LEOPRD = MOTLOG(ADRESS(3,21)) + ! IF(.NOT.DEFZON) NZONE=0 ! diff --git a/sources/telemac2d/majzz.f b/sources/telemac2d/majzz.f index bc0154ed1ec53aa78b12f317a244b292fd8324ec..230618a7ef6170688ce41e03c08731883cd07359 100644 --- a/sources/telemac2d/majzz.f +++ b/sources/telemac2d/majzz.f @@ -55,7 +55,7 @@ USE DECLARATIONS_TELEMAC2D, ONLY: NPOIN,LT,NPTFR,NTRAC,ICIN, & TORDER,MESH,RAIN,DT,GAMMA,V2DPAR,HN, & GRAV,PLUIE,H,FLUXT_OLD,SMTR,FLUXT, - & SMH,EPS_FV + & SMH,EPS_FV,DTN USE DECLARATIONS_TELEMAC, ONLY: KNEU USE DECLARATIONS_SPECIAL IMPLICIT NONE @@ -229,6 +229,10 @@ ! ENDIF ! +! UPDATE OF DTN +! ************* + DTN = DT +! !----------------------------------------------------------------------- ! RETURN diff --git a/sources/telemac2d/preres_telemac2d.f b/sources/telemac2d/preres_telemac2d.f index 11d600bea31cde4b09b6e9db749c0cb4ab3baa0c..13c7e1dff001f067285b53c52b3a8918d5507349 100644 --- a/sources/telemac2d/preres_telemac2d.f +++ b/sources/telemac2d/preres_telemac2d.f @@ -124,7 +124,7 @@ ! FVM IF(LT.GT.PTINIL)THEN ! LISTING OUTPUT - LPRDTIME=LISPRD*DTINI + LPRDTIME=LISPRD*DTCAS LTT=CEILING(AT/LPRDTIME) RESTE=(LTT*LPRDTIME-AT)/LPRDTIME IF(RESTE.LT.EPSS.OR.ABS(RESTE-1.D0).LT.EPSS.OR. diff --git a/sources/telemac2d/second_order.f b/sources/telemac2d/second_order.f index ba3b0cbb3d7549170780f89ac97f61721d69f951..ff327207ada0abc24fba0c6802d0fb3bb7958a20 100644 --- a/sources/telemac2d/second_order.f +++ b/sources/telemac2d/second_order.f @@ -50,7 +50,9 @@ & GRADIJ_FC,GRADJI_FC,DEJA_FC,IKLE, & NPTFR,ILIMHZ,ILIMUV,X,Y,V2DPAR,ZF, & GRAV,CFLWTD,MESH,EPS_FV,AT,TMAX, - & PTINIG,DTMIN_FV + & PTINIG,DTMIN_FV,EPS_DIV0,ICIN, + & ADAPT_DT_TO_LEOPRD,DTN,LEOPRD, + & DTVARI,DTCAS,ENTET USE DECLARATIONS_TELEMAC, ONLY: KDIR,KNEU USE DECLARATIONS_SPECIAL USE INTERFACE_PARALLEL, ONLY : P_MIN @@ -86,7 +88,7 @@ DOUBLE PRECISION RA3,RNN,SIGMAX,UAS11,UAS12,UAS21,UAS22 DOUBLE PRECISION UAS31,UAS32,UAS210,UAS220,UNORM,VNL,VNX DOUBLE PRECISION VNY,XNN,YNN,ZF1,ZF2,AMDS,AUX,AUX1,U,V - DOUBLE PRECISION CORR_TMP,DSH,RESTE + DOUBLE PRECISION CORR_TMP,DSH,RESTE0,RESTE1,DTV LOGICAL, ALLOCATABLE :: YESNO(:) DOUBLE PRECISION, ALLOCATABLE :: TMP1(:),TMP2(:) ! @@ -561,19 +563,27 @@ ! LIMITATION OF THE TIME STEP ! *************************** ! -! - SIGMAX=MAX( 1.D-2, RA3* SQRT(UAS11) + ABS(UAS21) ) - DTL = CFLWTD*AIRST(1,NSG)/(RNN*SIGMAX) - DTLL_FC(NUBO1) = MIN (DTL,DTLL_FC(NUBO1)) - DT = MIN(DT, DTL) -! - SIGMAX=MAX( 1.D-2, RA3* SQRT(UAS12) + ABS(UAS22) ) - DTL = CFLWTD*AIRST(2,NSG)/(RNN*SIGMAX) - DTLL_FC(NUBO2) = MIN (DTL,DTLL_FC(NUBO2)) - DT = MIN(DT, DTL) -! PARALLEL: TAKE MIN DT OF ALL SUBDOMAINS -! IF(NCSIZE.GT.1)DT = P_MIN(DT) !WILL BE PLACED AT THE END (SEE BELOW) -! +! KINETIC SCHEME + IF(ICIN.EQ.1)THEN + SIGMAX = MAX(EPS_DIV0, RA3*SQRT(UAS11) + ABS(UAS21)) +! OTHER SCHEMES + ELSE + SIGMAX = MAX(EPS_DIV0, SQRT(GRAV*UAS11) + ABS(UAS21)) + ENDIF + DTL = CFLWTD*AIRST(1,NSG)/(RNN*SIGMAX) + DTLL_FC(NUBO1) = MIN(DTL, DTLL_FC(NUBO1)) + DT = MIN(DT, DTL) +! +! KINETIC SCHEME + IF(ICIN.EQ.1)THEN + SIGMAX = MAX(EPS_DIV0, RA3*SQRT(UAS12) + ABS(UAS22)) +! OTHER SCHEMES + ELSE + SIGMAX = MAX(EPS_DIV0, SQRT(GRAV*UAS12) + ABS(UAS22)) + ENDIF + DTL = CFLWTD*AIRST(2,NSG)/(RNN*SIGMAX) + DTLL_FC(NUBO2) = MIN(DTL, DTLL_FC(NUBO2)) + DT = MIN(DT, DTL) ! 1234 CONTINUE ! @@ -584,7 +594,13 @@ ENDIF ENDDO ENDDO - +! + DEALLOCATE(YESNO) +! +! ****************************************************************** +! ADJUST TIME STEP +! ****************************************************************** +! IF(NPTFR.GT.0)THEN ! USEFUL FOR PARALLEL CASE DO K=1,NPTFR IS = MESH%NBOR%I(K) @@ -598,37 +614,111 @@ U = 0.D0 V = 0.D0 ENDIF - SIGMAX= SQRT(UA(1,IS)) UNORM=SQRT(U*U + V*V) - SIGMAX=MAX( 1.D-2, RA3*SIGMAX +UNORM ) +! KINETIC SCHEME + IF(ICIN.EQ.1)THEN + SIGMAX = MAX(EPS_DIV0, RA3*SQRT(UA(1,IS)) + UNORM) +! OTHER SCHEMES + ELSE + SIGMAX = MAX(EPS_DIV0, SQRT(GRAV*UA(1,IS)) + UNORM) + ENDIF DTL = CFLWTD*V2DPAR%R(IS)/(VNL*SIGMAX) AUX = DTL/DTLL_FC(IS) - AUX1 =AUX/(1.D0+AUX) - DT =MIN(DT, AUX1*DTLL_FC(IS)) + AUX1 = AUX/(1.D0+AUX) + DT = MIN(DT, AUX1*DTLL_FC(IS)) ENDDO ENDIF - -! FOR PARALLELISME +! +! FOR PARALLELISM IF(NCSIZE.GT.1) DT=P_MIN(DT) ! - DEALLOCATE(YESNO) +! ****************************************************************** +! CONSTANT SET TIME STEP +! ****************************************************************** + IF(DTVARI.EQV..FALSE.) THEN + DTV = DT ! STORE DT FOR RECOMMENDATION IN WARNING + DT = DTCAS + ENDIF ! +! ****************************************************************** +! END OF COMPUTATION +! ****************************************************************** THEEND =.FALSE. LEO =.FALSE. - IF((TMAX-AT).LE.DT)THEEND=.TRUE. !LAST TIME STEP -! ADAPT DT TO TAKE INTO ACCOUNT GRAPHIC OUTPUT - IS=CEILING(AT/GPRDTIME) - RESTE=IS*GPRDTIME-AT -! - IF(THEEND.OR. - & (RESTE.LE.DT.AND.RESTE.GT.DTMIN_FV.AND.LT.GT.PTINIG))THEN -! HERE THERE IS GRAPHICAL OUTPUT - LEO = .TRUE. - DT=MIN(RESTE,DT) + IF((TMAX-AT).LE.DT) THEEND=.TRUE. !LAST TIME STEP + DT = MIN(DT,TMAX-AT) +! +! ****************************************************************** +! DETECTION OF GRAPHICAL PRINTOUTS (LEO = TRUE) +! ****************************************************************** +! +! VARIABLE DT + IF(DTVARI) THEN +! +! ADAPT DT TO TAKE INTO ACCOUNT GRAPHIC OUTPUT + IF (ADAPT_DT_TO_LEOPRD) THEN + IS = CEILING(AT/GPRDTIME) + RESTE1 = IS*GPRDTIME-AT + IF(THEEND.OR.(RESTE1.LE.DT + & .AND.RESTE1.GT.DTMIN_FV + & .AND.LT.GT.PTINIG))THEN +! HERE THERE IS GRAPHICAL OUTPUT + LEO = .TRUE. + DT = MIN(RESTE1,DT) + ENDIF +! +! DETECT CLOSEST GRAPHIC PRINTOUT TO CURRENT TIME (AT+DT) + ELSE + IS = CEILING((AT+DT)/GPRDTIME) + RESTE0 = AT+DT - (IS-1)*GPRDTIME + IF (LT.EQ.1) THEN + IF(THEEND.OR.(RESTE0.LT.DT .AND.LT.GT.PTINIG)) THEN + LEO = .TRUE. + ENDIF + ELSE + IF(THEEND.OR.(RESTE0.LE.DT.AND.LT.GT.PTINIG)) THEN + LEO = .TRUE. + ENDIF + ENDIF + ENDIF +! +! CONSTANT DT + ELSE + IF (THEEND.OR.(MOD(LT,LEOPRD)==0)) THEN + LEO = .TRUE. + ENDIF ENDIF ! - DT = MIN(DT,TMAX-AT) +! ****************************************************************** +! ERRORS AND WARNINGS +! ****************************************************************** ! +! CASE DT <= 0 + IF(DT.LE.0.D0) THEN + WRITE(LU,*) 'NEGATIVE (OR NIL) TIME-STEP: ', DT + CALL PLANTE(1) + STOP + ENDIF +! +! CASE CFL >= 1 + IF(CFLWTD.GE.1.D0) THEN + IF(ENTET) THEN + WRITE(LU,*) 'WARNING: CFL NOT GIVEN OR >=1 !...! ' + WRITE(LU,*) 'RECOMMENDED TIME-STEP (WITH CFL = 0.9): ', + & 0.9D0*DT/CFLWTD + ENDIF + ENDIF +! +! CONSTANT SET TIME STEP + IF(DTVARI.EQV..FALSE.) THEN + IF(ENTET) THEN + WRITE(LU,*) 'WARNING: FIXED TIME-STEP (NOT RECOMMENDED)' + WRITE(LU,*) 'TIME-STEP MAY NOT SATISFY CFL CONDITION, ' + WRITE(LU,*) 'RECOMMENDED TIME-STEP (WITH CFL = 0.9): ', + & 0.9D0*DTV/CFLWTD + ENDIF + ENDIF +! !----------------------------------------------------------------------- ! RETURN diff --git a/sources/telemac2d/telemac2d.dico b/sources/telemac2d/telemac2d.dico index 2a8a4961b55f58d01e1fce7d2d6e980a917f6147..e35186114ff0a5e68db75f9844ccc7c9120efad8 100644 --- a/sources/telemac2d/telemac2d.dico +++ b/sources/telemac2d/telemac2d.dico @@ -7000,6 +7000,26 @@ AIDE1 = \item 2: Strong imposition. \end{itemize}' / +NOM = 'ADAPTATION DU PAS DE TEMPS AU SORTIES GRAPHIQUES' +NOM1 = 'ADAPTATION OF TIME STEP TO GRAPHIC PRINTOUTS' +TYPE = LOGICAL +INDEX = 21 +TAILLE = 1 +DEFAUT = OUI +DEFAUT1 = YES +MNEMO = 'ADAPT_DT_TO_LEOPRD' +RUBRIQUE = +'PARAMETRES NUMERIQUES';'AVANCEE';'' +RUBRIQUE1 = +'NUMERICAL PARAMETERS';'ADVANCED';'' +NIVEAU = 1 +AIDE = +'Adaptation du pas de temps a la periode de sortie graphique, +(seulement avec les volumes finis).' +AIDE1 = +'Adaptation of time step to graphic printouts, +(only with finite volumes).' +/ NOM = 'ETUDE DE CONVERGENCE' NOM1 = 'CONVERGENCE STUDY' TYPE = LOGICAL diff --git a/sources/telemac2d/telemac2d_init.F b/sources/telemac2d/telemac2d_init.F index 893010455e7975e7c9d1aff81979b69df99d4d99..f4c49cb669e826c49ef39501366b477887ca4658 100644 --- a/sources/telemac2d/telemac2d_init.F +++ b/sources/telemac2d/telemac2d_init.F @@ -696,7 +696,6 @@ ! IF(EQUA(1:15).EQ.'SAINT-VENANT VF') THEN ! - DTINI=DT CALL OS('X=YZ ', X=QU, Y=U, Z=H) CALL OS('X=YZ ', X=QV, Y=V, Z=H) ! @@ -713,9 +712,8 @@ ENDIF ENDIF ! - TMAX =DUREE+AT0 - DTINI=DT - EPS_FV=1.D-6 + TMAX = DUREE+AT0 + EPS_FV = 1.D-6 ! ! COMPUTE COORDINATES OF POINTS TO USE FOR FIELD ! RECONSTRUCTION IN FINITE VOLUME DIFFUSION SOLVER