open TELEMAC-MASCARET The mathematically superior suite of solvers

Good morning Telemac users,

I would like to run a parameter estimation of the bottom friction coefficient for my tidal model, but I want it to be done with only one time serie of free surface at one point of my mesh.

I tried to modify mesures.f but I am not really good with Fortran, so I wasn't able to do it.

Does anyone know what is the Fortran code I have to implement in order to achieve this ?

Thank you :)

Hello,

If you have only one point of measurement (possibly depending on iteration ITER) what you can do in subroutine mesures.f is :

First cancelling all values and weights:

CALL OS( 'X=0 ' , X=HD )
CALL OS( 'X=0 ' , X=UD )
CALL OS( 'X=0 ' , X=VD )
CALL OS( 'X=0 ' , X=ALPHA1 )
CALL OS( 'X=0 ' , X=ALPHA2 )
CALL OS( 'X=0 ' , X=ALPHA3 )

and keep the integral of test functions (we put weights proportionnal to the area around points):
CALL VECTOR(T1,'=','MASBAS ',
& HD%ELM,1.D0,T3,T3,T3,T3,T3,
& T3,MESH,MSK,MASKEL)

Then suppose you know that the depth of point 325 is 3.2 m:

HD%R(325)=3.2D0 (here it may also depend on the iteration ITER)
ALPHA1%R(325)=T1%R(325)

That should do it (not tested...)

Regards,

Jean-Michel Hervouet

Good morning
Could you explain more clealy which subroutine have to modified please
Thanks you

Hello,

See how example "estimation" works. The difference with your case is that in the example "estimation" the values are read on a file, while in your case they are all with weights equal to 0, except on one point (but they vary in time).

So you need to program the subroutine mesures.f, replacing the existing implementation by:

CALL OS( 'X=0 ' , X=HD )
CALL OS( 'X=0 ' , X=UD )
CALL OS( 'X=0 ' , X=VD )
CALL OS( 'X=0 ' , X=ALPHA1 )
CALL OS( 'X=0 ' , X=ALPHA2 )
CALL OS( 'X=0 ' , X=ALPHA3 )

Then suppose you know that the depth of point 325 is 3.2 m:

HD%R(325)=3.2D0 (but in your case a function of time which is TT in the subroutine)
ALPHA1%R(325)=1.D0 (I put here a weight of 1.)

Having one time series at one point, I would start by estimating a single value of friction, so only one zone to start with.

I hope this helps,

With best regards,

Jean-Michel Hervouet

Good morning
I have try with estimation cas of telemac2d.
I try modified mesures.f like follow

! ******************
SUBROUTINE MESURES
! ******************
!
&(ITER,TT)

USE BIEF
USE DECLARATIONS_TELEMAC2D
!
IMPLICIT NONE
INTEGER LNG,LU
COMMON/INFO/LNG,LU
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
INTEGER, INTENT(IN) :: ITER
DOUBLE PRECISION, INTENT(IN) :: TT
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
DOUBLE PRECISION TPS,C
LOGICAL OKH,OKU,OKV
INTEGER I
!
!

---

!
IF(T2D_FILES(T2DREF)%NAME(1:1).NE.' ') THEN
!
!

---

!
! WHEN MEASUREMENTS ARE IN A SELAFIN FILE
!
CALL FIND_IN_SEL(HD,TEXTE(4)(1:16),T2D_FILES(T2DREF)%LU,
& W,OKH,RECORD=ITER,TIME=TPS)
CALL FIND_IN_SEL(UD,TEXTE(1)(1:16),T2D_FILES(T2DREF)%LU,
& W,OKU,RECORD=ITER,TIME=TPS)
CALL FIND_IN_SEL(VD,TEXTE(2)(1:16),T2D_FILES(T2DREF)%LU,
& W,OKV,RECORD=ITER,TIME=TPS)

IF(.NOT.OKH.OR..NOT.OKU.OR..NOT.OKV) THEN
IF(LNG.EQ.1) THEN
WRITE(LU,*) 'MESURES : PROBLEME DE LECTURE DE HD, UD OU VD'
ENDIF
IF(LNG.EQ.2) THEN
WRITE(LU,*) 'MESURES : PROBLEM WHEN READIND HD, UD, OR VD'
ENDIF
CALL PLANTE(1)
STOP
ENDIF
IF(ABS(TT-TPS).GT.1.D-3) THEN
IF(LNG.EQ.1) THEN
WRITE(LU,*) 'MESURES : PROBLEME DE LECTURE DU TEMPS'
ENDIF
IF(LNG.EQ.2) THEN
WRITE(LU,*) 'MESURES : PROBLEM WHEN READIND TIME'
ENDIF
CALL PLANTE(1)
STOP
ENDIF
! UD AND VD MAY BE QUASI-BUBBLE
!(BUT ALPHA2 AND ALPHA3 WILL BE SET TO 0)
IF(UD%ELM.EQ.12) THEN
CALL CHGDIS(UD,11,12,MESH)
CALL CHGDIS(VD,11,12,MESH)
ENDIF
!
!

---

!
ELSE
!
!

---

!
! CASE TO BE IMPLEMENTED BY USER HERE (OTHER FILE FORMAT, ETC.)
!
IF(LNG.EQ.1) WRITE(LU,*) 'MESURES A PROGRAMMER DANS MESURES'
IF(LNG.EQ.2) WRITE(LU,*) 'MEASUREMENTS TO IMPLEMENT IN MESURES'

CALL OS( 'X=0 ' , X=HD )
CALL OS( 'X=0 ' , X=UD )
CALL OS( 'X=0 ' , X=VD )
CALL OS( 'X=0 ' , X=ALPHA1 )
CALL OS( 'X=0 ' , X=ALPHA2 )
CALL OS( 'X=0 ' , X=ALPHA3 )

IF(TT.eq.0)THEN
ALPHA1%R(264)=1.D0
HD%R(264)=0.5D0
ELSEIF(TT.eq.200)THEN
ALPHA1%R(264)=1.D0
HD%R(264)=0.6625382304191589D0
ELSEIF(TT.eq.400)THEN
ALPHA1%R(264)=1.D0
HD%R(264)=0.7598909735679626D0
ELSEIF(TT.eq.600)THEN
ALPHA1%R(264)=1.D0
HD%R(264)=0.7875543832778931D0
ELSEIF(TT.eq.800)THEN
ALPHA1%R(264)=1.D0
HD%R(264)=0.7979275584220886D0
ELSEIF(TT.eq.1000)THEN
ALPHA1%R(264)=1.D0
HD%R(264)=0.8018307685852051D0
ELSEIF(TT.eq.1200)THEN
ALPHA1%R(264)=1.D0
HD%R(264)=0.8032307028770447D0
ELSEIF(TT.eq.1400)THEN
ALPHA1%R(264)=1.D0
HD%R(264)=0.8032264113426208D0
ELSEIF(TT.eq.1600)THEN
ALPHA1%R(264)=1.D0
HD%R(264)=0.8037620186805725D0
ELSEIF(TT.eq.1800)THEN
ALPHA1%R(264)=1.D0
HD%R(264)=0.8043234944343567D0
ELSEIF(TT.eq.2000)THEN
ALPHA1%R(264)=1.D0
HD%R(264)=0.8034544587135315D0
ENDIF
print*,ALPHA1%R(264)
print*,HD%R(264)
pause
CALL PLANTE(1)


STOP
!
!

---

!
ENDIF
!
!

---

!
! WEIGHT FUNCTIONS FOR ALL THE TIMESTEPS
!
CALL VECTOR(T1,'=','MASBAS ',
& HD%ELM,1.D0,T3,T3,T3,T3,T3,
& T3,MESH,MSK,MASKEL)
CALL OS( 'X=Y ' , ALPHA1 , T1 , T1 , C )

! CASE OF QUASI-BUBBLE ELEMENT FOR UD
IF(HD%ELM.NE.UD%ELM) THEN
CALL VECTOR(T1,'=','MASBAS ',
& UD%ELM,1.D0,T3,T3,T3,T3,T3,
& T3,MESH,MSK,MASKEL)
ENDIF

CALL OS( 'X=Y ' , ALPHA2 , T1 , T1 , C )
CALL OS( 'X=Y ' , ALPHA3 , T1 , T1 , C )
!
! CANCELS WEIGHTS FOR QUASI-BUBBLE POINTS
!
IF(UD%ELM.EQ.12) THEN
DO I=NPOIN+1,NPOIN+NELEM
ALPHA2%R(I)=0.D0
ALPHA3%R(I)=0.D0
ENDDO
ENDIF


!
!

---

!
RETURN
END

It don't run
Could you explain.
Thanks you very much

Hello,

I think that your tests IF(TT.EQ.200) will never work, because 200 is an integer and TT is a double precision number, that may be in fact 199.9999999999,
so it would be better to write :

IF(ABS(TT-200.D0).LT.1.D-4) THEN


for example, allowing a truncation error.

With best regards,

Jean-Michel Hervouet

Good morning
Could you send us your example Loire Dampiere with estimation of friction please?
Thanks you

Excuse me,
I have tried many times our real case, and Telemac always give us the sensibility infini, grad(J)= xxx E+50.
Could you help us

Hello,

This Loire Dampierre case was a real study and has not been maintained, however it is exactly like the estimation case in the examples. Rather post your Fortran and steering file, so that we can have a look at it.

Regards,

JMH

Good morning
My steering file is quiet identic to the exemple 023_estimation :

/

---

/ TELEMAC2D Version v6p1 Mar 12, 2014
/ TELEMAC 2D : TEST DE CANAL$
/

---

/

---

/ ENTREES-SORTIES, FICHIERS
/

---

FICHIER DES RESULTATS = 'r2d_estimation_3000_1zone.slf'

FICHIER DES CONDITIONS AUX LIMITES ='cas.conlim'

FICHIER DES FRONTIERES LIQUIDES = 'cas.liq'

/FICHIER DE REFERENCE ='f2d_estimation.slf'

/CONDITIONS INITIALES : VOIR CONDIN
FICHIER FORTRAN ='t2d_estimation.f'

FICHIER DE GEOMETRIE ='luoi'

FICHIER DES PARAMETRES ='t2d_estimation.cas'

FICHIER DU CALCUL PRECEDENT = 'CI_3000_70.slf'

/ HP C3700 compilateur HP : 60 s 5.5 26/11/2004
/ HP C3700 compilateur HP : 57 s 5.7 18/04/2007
/ Dell 2.8 MHz compilateur pgi : 43 s 5.7 19/04/2007
/ HP C3700 compilateur HP : 59 s 5.8 13/11/2007
/ HP C3700 compilateur Nag : 147 s 5.8 20/12/2007
/ Dell 2.8 MHz compilateur pgi : 44 s 5.8 19/12/2007
/ HP C3700 compilateur HP : 57 s 5.9 15/10/2008
/ HP C3700 compilateur Nag : 148 s 5.9 16/10/2008
/ Dell 2.8 MHz compilateur pgi : 43 s 5.9 15/10/2008
/ HP C3700 compilateur HP : 55 s 6.0 24/11/2009
/ HP C3700 compilateur Nag : 152 s 6.0 26/11/2009
/ Dell 2.8 MHz compilateur Intel : 32 s 6.0 26/11/2009
/ EN MODE ESTIMATION

ESTIMATION DE PARAMETRE ='FROTTEMENT, PERMANENT'

FICHIER DE RESULTATS FORMATE ='rfo_estimation.txt'

FICHIER DE RESULTATS BINAIRE ='rbi_estimation.slf'


/

---

/ ENTREES-SORTIES, GENERALITES
/

---

SUITE DE CALCUL = OUI

VALIDATION =NON

TITRE ='TELEMAC 2D : TEST DE CANAL$'


/

---

/ ENTREES-SORTIES, GRAPHIQUES ET LISTING
/

---

PERIODE DE SORTIE LISTING =1000

VARIABLES POUR LES SORTIES GRAPHIQUES = 'U,V,H,B,S'

PERIODE POUR LES SORTIES GRAPHIQUES =1000

BILAN DE MASSE =OUI

INFORMATIONS SUR LE SOLVEUR =OUI


/

---

/ EQUATIONS
/

---

CONVECTION =OUI

COEFFICIENT DE FROTTEMENT =100.

COEFFICIENT DE DIFFUSION DES VITESSES =0.

/ REFERENCE : 35.
/LOI DE FROTTEMENT SUR LE FOND = 3 COEFFICIENT DE FROTTEMENT = 35.
LOI DE FROTTEMENT SUR LE FOND =3


/

---

/ EQUATIONS, CONDITIONS INITIALES
/

---

CONDITIONS INITIALES ='COTE CONSTANTE'

COTE INITIALE =3.66


/

---

/ EQUATIONS, CONDITIONS LIMITES
/

---

DEBITS IMPOSES =0.; 50.

COTES IMPOSEES =0.5; 0.


/

---

/ PARAMETRES NUMERIQUES
/

---

PRECISIONS POUR L'IDENTIFICATION =1.E-3;1.E-3;1.E-3;1.E-5

PRODUIT MATRICE-VECTEUR =1

DEFINITION DE ZONES =OUI

NOMBRE DE PAS DE TEMPS =36000

FORME DE LA CONVECTION =1;5

/ 1 : GRADIENT 2 : GRADIENT CONJUGUE 3 : INTERPOLATION DE LAGRANGE
METHODE D'IDENTIFICATION =2

DISCRETISATIONS EN ESPACE =11; 11

PRECONDITIONNEMENT =2

PAS DE TEMPS =0.5

OPTION DE SUPG =1;2

/

---

/---/
STOCKAGE DES MATRICES =1

FONCTION COUT =1

MAXIMUM D'ITERATIONS POUR L'IDENTIFICATION =20


/

---

/ PARAMETRES NUMERIQUES, SOLVEUR
/

---

SOLVEUR =7

PRECISION DU SOLVEUR =1.E-4

OPTION DE TRAITEMENT DES BANCS DECOUVRANTS = 2

CLIPPING DE H = OUI

VALEUR MINIMUM DE H = 0.1
/

---

/ PARAMETRES NUMERIQUES, VITESSE-CELERITE-HAUTEUR
/

---

ORDRE DU TIR INITIAL POUR H =1

IMPLICITATION POUR LA HAUTEUR =0.55

IMPLICITATION POUR LA VITESSE =0.55