Culverts

This issue was created automatically from an original CUE issue. Further discussion may take place here.

---

The joined example illustrate a problem in buse.f causing NaN
In this case the crash occur at line 261 of buse.f but it could also occur on line 300 for the opposite flow.
There is no verification that S1-(RD2+H2) is positive...
By the way, in the formula, is this really S1-(RD2+H2) that we should have in this case?

bug_buse.zip

assigned to @chi-tuan.pham

added module:telemac-2d type:bug labels

changed milestone to %v8p2

changed milestone to %v8p3

Hello Christophe,

Happy New Year and best wishes for 2021, it is my first cue message this year! :),

Sorry for the delay, after the release of tag v8p2r0 at the beginning of December, lots of work to do before my holidays mid December and this bug seemed to need relevant time to be solved as I am not so familiar with this feature (Agnes used to deal with it a few years ago).

Anyway, you are right, S1-(RD2+H2) can be negative in some cases. Looking at the TELEMAC-3D theory guide, in that configuration, it is not exactly the same formula implemented in the code. It seems that the modification was done when deleting siphon formulation and changing it in a culvert feature in 2017.

What about adding a new test if S1-(RD2+H2) is negative and then give a flowrate formula with S1-S2 rather than S1-(RD2+H2), the remaining source code still unchanged?

IF(S1.GE.(RD1+H1).AND.S2.LT.(RD2+H2)) THEN
                  Q = SECBUS(N)*SQRT(2.D0*GRAV*(S1-(RD2+H2))
     &              /(1.D0+CE1+L))
                ELSE
                  Q = SECBUS(N)*SQRT(2.D0*GRAV*(S1-S2)/(L+CS2+CE1))
                ENDIF

changed to

IF(S1.GE.(RD1+H1).AND.S2.LT.(RD2+H2)) THEN
                  IF(S1.GE.(RD2+H2)) THEN
                    Q = SECBUS(N)*SQRT(2.D0*GRAV*(S1-(RD2+H2))
     &                /(1.D0+CE1+L))
                  ELSE
                    Q = SECBUS(N)*SQRT(2.D0*GRAV*(S1-S2)
     &                /(1.D0+CE1+L))
                  ENDIF
                ELSE
                  Q = SECBUS(N)*SQRT(2.D0*GRAV*(S1-S2)/(L+CS2+CE1))
                ENDIF

or

IF(S1.GE.(RD1+H1).AND.S2.LT.(RD2+H2).AND.
     &             S1.GE.(RD2+H2)) THEN
                  Q = SECBUS(N)*SQRT(2.D0*GRAV*(S1-(RD2+H2))
     &              /(1.D0+CE1+L))
                ELSE
                  Q = SECBUS(N)*SQRT(2.D0*GRAV*(S1-S2)/(L+CS2+CE1))
                ENDIF

My knowledge of culverts is quite limited, I do not know which one would be the best (my first idea would be to choose the 2nd one). Perhaps you may have a better opinion?

Chi-Tuan

Hello,

We proprose with @lucie.mesquita and Thomas Viard to fix this problem with a verification of the formula.

We will work on this branch https://gitlab.pam-retd.fr/otm/telemac-mascaret/-/commits/culvert_fix

Matthieu

closed with commit 6ccfbd5b

Only the fix for possible negative values under square roots has been pushed in SHA 6ccfbd5b

Other fixes from branch culvert_fix are postponed at the moment.

mentioned in issue #761 (closed)

mentioned in commit 6ccfbd5b

mentioned in commit 5f166a24

mentioned in commit 6ea83f83

mentioned in commit c4230db2

mentioned in commit 96e49325

mentioned in commit a44db17f

mentioned in commit 7e81da70

mentioned in merge request !84 (merged)

mentioned in commit 19e2afa4

mentioned in commit 15cc708b