Accumulation at inlet with no sediment influx

System Modules GAIA / SISYPHE
1

Hi [name removed],

I have run into a perplexing issue that runs counter to what I would’ve expected. In essence, I’m trying to determine the armoured GSD from a known, pre-armoured GSD and so I’m running the model without any sediment input (fluxes = 0.,0.). I am using GAIA, [name removed] and [name removed], and bedload only for this simulation. My GSD runs from over 56 mm to under 1 mm (gravel stream).

My problem is that the bed at the inlet boundary is beginning to rise in ripples perpendicular to the flow direction. This accumulation is on the order of 7 cm in 35 minutes. This doesn’t sound like much, but it does cause a numerical instability problem where at the scale of my model (flow depth ≤ 30 cm), it causes supercritical inflow and since I don’t have an H at the inlet, it halt-and-catch-fires. Obviously, I could bandage the model by setting a known H, but this only moves the problem back a step.

Fundamentally, what I don’t understand is why would the bed rise, and why would sediment build up, at the inlet when there is no sediment influx? I would expect the opposite - scouring at the inlet - rather than a buildup. Maybe I’m missing something obvious, but this makes no physical sense to me. I have checked the cumulative bed evolution on the boundary itself and it is > 0 cm, so it’s not digging downward on the boundary to accumulate sediment at the next nodes.

Could it be worth it to try with a different set of equations? [name removed] seems most suitable to my application, so I’d rather not try others for now, but I can try MP&M as well.

I have attached a screenshot of the phenomenon as well as my steering files (apologies for the mess; I have no idea how best to structure a .cas). If needed, I can provide .cli and .slf, but I’d prefer not to if possible.

Thanks for any insight on why GAIA would accumulate sediment!

[name removed]

Ripples.jpg
Ripples.jpg1787×1200 45.1 KB

hdyn-coup-01h00.cas (4.23 KB)
gaia-wc.cas (2.51 KB)

2

Hi [name removed],

Apologies for the double post. I tried with an influx of 8 kg/min (0.1375 kg/s) and this time, the same conditions (accumulation of 7 cm, supercritical crash) were observed after 55 minutes of simulation time instead of 37 minutes. This is most perplexing…

3

[name removed],

how did you set your open inlet boundary (one line from the cli-file would be enough)?

What happens to the open boundary nodes. Is there an accumulation? From the file I would think there is no evolution directly at the boundary. What about the shear stress at the boundary? Is it high enough that the sediments should move?

Two ideas to smooth instable results especially at the boundaries:

You can try FINITE VOLUMES = YES and play around with upwind criteria (UPWINDING FOR BEDLOAD default is 0.5. If you use up to 1.0 it is more diffusive but also more stable.

Furthermore sometimes it can be helpful to start with a rigid bed at the inflow boundary. Typically, this also smoothes results because it decouples the morphodynamics from the hydrodynamics for the uncertain boundary part.

[name removed]

4

Hi [name removed],

Thanks for your reply with many avenues to pursue!

I will try the FINITE VOLUMES / UPWINDING FOR BEDLOAD keywords to begin with for sure. This will take a while to run, assuming it does run to begin with :laughing:

For the rigid bed, would I have to look into creating my own noerod.f or is there a way to define these areas in the .slf (eg. Blue Kenue)? If the former, which examples would be best to look at, if you know of any in particular?

For shear stress on the upstream boundary, it ranges from 13 to 16 N/m[sup]2[/sup] which is plenty to move this granulometry. Cumulative bed evolution on the upstream boundary ranges from 0 to 18 mm upwards; as mentioned, nodes 1–2 rows from downstream reach 70 mm.

Here are two lines from my .cli file, one for a downstream BC and one for a closed node:

5 4 4  0.000 0.000 0.000 0.000  4  0.000 0.000 0.000           1           1   # downstream (689 - 1)

2 2 2  0.000 0.000 0.000 0.000  2  0.000 0.000 0.000          11           2   #

I modified it as per this article.

Thank you again for your helpful advice!

[name removed]

5

Hi [name removed],

for the upstream boundary in gaia and a prescribed solid discharge you need 1.0 in the 4th column. I didn’t see that in the website you refer to. But maybe I just did not found it. (see validation example flume_bc-t2d)

4 5 5 1.0 0.0 0.0 0.0 4

In gaia there is no noerod.f anymore. You can define the rigid bed by the layer thickness(es). If all layer-thicknesses are zero then you have a rigid bed. You can do this in user_bed_init.f (see the hippodrom-t2d validation test case).

Or you create a previous sedimentological computation file and changes the active layer thicknesses via blue-kenue.

[name removed]

6

Hi [name removed],

Thanks for these ideas, I have adjusted my .cli accordingly. Someday I will understand the CONLIM structure :laughing:

A few updates.

First, I tried the FINITE VOLUMES / UPWINDING FOR BEDLOAD with values of 0.50, 0.75, and 1.00 for the latter keyword. The 0.50 completed a 1h00 simulation, and the 0.75 and 1.00 both failed after just under 0h30 of simulated time. Moreover, all simulations still have rippling at the inlet, and the rippling is somehow more pronounced than it was! That is to say, it seems to propagate further downstream…

Now, I am trying to make the first metre or so non-erodible to see if this will work. Unfortunately, this is turning out to be more complicated than I thought (and it looked complicated already!). I grabbed a list of elements from the mesh that I wanted to make non-erodible, put them in a text file, and tried reading the file as part of user_bed_init.f (attached).

I have debugged a few points already, including the need to manually copy the noerod.csv into the working directory (scratch directory on an HPC cluster) and apparently needing to transform points from global to local (I was getting segfaults if I just fed it the following).

! DO loop here where I feed it a node number IPOIN

            ESTRATUM(1,IPOIN) = 10000.D0

! loop...

It should be noted that I am doing this computation in parallel, as it is a rather heavy simulation and I am not sure how to run only the user_bed_init step in series.

My latest problem is at the compilation stage. Here is the relevant code:

OPEN(UNIT=10, FILE="noerod.csv", [name removed]", ACTION="read")

      DO

            READ(10,*,END=99) IPOIN

            ESTRATUM(1,GLOBAL_TO_LOCAL_POINT(IPOIN,MESH)) = 0.D0

      END DO

99    CONTINUE

! close the file, return, etc.

This apparently leads to memory errors and I’m not sure why; at this stage, the errors become:

================================================================================

 ITERATION        0    TIME:   0.0000 S

 PRERES: MAXIMUM COURANT NUMBER:    0.1021956



                            BALANCE OF OMEGA (UNIT: NA * M3)

     INITIAL QUANTITY OF TRACER   :    0.000000

 TELEMAC2D COUPLED WITH:

 GAIA



 TELEMAC-2D: INTERNAL COUPLING WITH GAIA

 INBIEF (BIEF): NOT A VECTOR MACHINE (ACCORDING TO YOUR DATA)



 LIQUID BOUNDARIES:           2



free(): invalid pointer                    GAIA INITIAL MASS OF SEDIMENTS:



 TOTAL MASS OF CLASS  1 =    2629.891    ( KG )

 TOTAL MASS OF CLASS  2 =    840.7029    ( KG )

 TOTAL MASS OF CLASS  3 =    2910.125    ( KG )

 TOTAL MASS OF CLASS  4 =    2220.318    ( KG )

 TOTAL MASS OF CLASS  5 =    1767.632    ( KG )

 TOTAL MASS OF CLASS  6 =    711.3640    ( KG )

 TOTAL MASS OF CLASS  7 =    603.5815    ( KG )

 TOTAL MASS OF CLASS  8 =    1530.510    ( KG )

 TOTAL MASS OF CLASS  9 =    1853.858    ( KG )

 TOTAL MASS OF CLASS 10 =    6488.501    ( KG )

free(): invalid pointer

free(): invalid pointer

free(): invalid pointer

free(): invalid pointer

free(): invalid pointer

free(): invalid pointer

free(): invalid pointer

free(): invalid pointer (and so on and so forth)

I realize this is now a Fortran question rather than a GAIA question, but how can I get this code working, especially with parallelization? Or can just this subroutine be run in serial? Would it be better to somehow import all the points in noerod.csv into user_bed_init.f? I have attached the full [group removed] file.

The previous computation file sounds like an interesting option, but in my mind it would still be simpler to have that as part of the geometry variables in the .slf object :slight_smile: Maybe I will try this, run just one time step and then use the output as previous computation.

Again, thank you so much for your help!

[name removed]

user_bed_init.f (2 KB)

7

Hi

You should read the telemac user manual on input/output files and use available keyword to manage correctly the copy/partition of user files and then avoid management of opening/close (see formated file 1 …)

8

Hi [name removed],

Thanks for these hints. Based on the user and reference manual as well as the developer guide, I tried changing my code in user_bed_init.f to the following, as well as using the FORMATTED DATA FILE 1 keyword in the .cas file:

DO

         READ(T2D_FILES(T2DFO1)%LU,END=99) IPOIN

         ESTRATUM(1,GLOBAL_TO_LOCAL_POINT(IPOIN,MESH)) = 0.D0

      END DO

99    CONTINUE

However, at runtime, it gave an error for every thread:

forrtl: severe (256): unformatted I/O to unit open for formatted transfers, unit 13

It sounds like the FORMATTED keyword is, for want of a better term, creating unrealistic expectations for what the I/O is supposed to see. As a reminder, my input file is simply a newline-separated list of integers referring to node IDs (Unix format on Unix machine), nothing else.

Since this failed, I tried the BINARY DATA FILE 1 keyword instead with BINARY DATA FILE 1 FORMAT set to ‘BIN’ since this is listed as an option in the T2D v8p4 reference manual. However, this gave a PLANTE(2) saying there are only 3 options (SERAFIN, SERAFIND, MED). It would appear there is a small discrepancy between the reference manual and the state of the code… so I have gone back to trying with FORMATTED instead.

Anyway, maybe I shouldn’t be trying any of this, as I am very much inexperienced with Fortran and the Telemac architecture hanging over what seems to me to simple operations. However, it seems to be the only “true” way so far to create a spatially heterogeneous active layer thickness (which occurs quite frequently in nature!), and I have really struggled to find any examples of a successful GAIA case where a file is read in parallel without errors.

Further advice is therefore most welcome, although I will certainly continue hacking away in the meantime to try and get this working. When I do, I will document it as breadcrumbs for the next poor soul trying this :stuck_out_tongue_winking_eye:

[name removed]

9

Hi

You could try a dos2unix for the formatted file …

Otherwise, you could have a look at the telemac2d pluie testcase particularly those with cn_geo.

THe cn value is directly implemented in the geometry file so the code is adapted to cope with such data.

Hope this helps

10

Hi [name removed],

in case you skip the fortran definition you can also try using a PREVIOUS SEDIMENTOLOGICAL COMPUTATION FILE. Just run your simulation only 1 time step and create a gaia-result file which should contain all layer thicknesses. Then you can set the layer thicknesses near the inlet to zero via blueKenue (I hope this is possible). Afterwards you start from that prev. sed. computation file.

Another idea we do sometimes is to set a special roughness coefficient (e.g. 0.0301 instead of 0.0300 otherwise which would not disturb your hydrodynamics, CHESTR is the roughness coefficient in telemac, you can include it in your subroutine by the “USE” command) near the inlet. Then the coding in user_bed_init is more easy:

USE DECLARATIONS_TELEMAC2D, ONLY : CHESTR

…

DO I=1,NPOIN

IF ([name removed]) THEN

ESTRATUM(1,I) = 0.D0

ESTRATUM(2,I) = 0.D0

ENDIF

END DO

[name removed]

1 Like
11

Hi [name removed],

Thanks for the ideas I tried the dos2unix utility and it appears to already have been in Unix format, so that’s unlikely to be the issue.

I then looked at the pluie case, as well as #39507 about setting a spatially variable bed (even though it may not work). My code now looks like this:

USE DECLARATIONS_TELEMAC2D, ONLY: MESH,PRIVE

! ...

      DO IPOIN=1,NPOIN

        ESTRATUM(ISTRAT,IPOIN) = PRIVE%ADR(1)%P%R(IPOIN)

      END DO

In the steering file, with [group removed] the variable I added to the geometry selafin:

NUMBER OF PRIVATE ARRAYS = 1

NAMES OF PRIVATE VARIABLES = 'ALTHICKNESS'

(Question, is the NAMES OF PRIVATE VARIABLES keyword necessary? It isn’t present in the pluie_cn case.)

And of course, it’s a new error this time, a good old segfault for each thread:

forrtl: severe (174): SIGSEGV, segmentation fault occurred

I feel like it must be so simple to take the data from a selafin variable and apply it to the layer thickness, and I must be missing something, but I cannot see for the life of me what it is.

My last resort will be to try setting every point as discussed further in the #39507 thread linked above, even though this seems like a dreadful approach.

I forgot to press Submit last night, so some updates:

I removed the NAMES OF PRIVATE VARIABLES keyword and the simulation ran, although I suspect it didn’t actually change anything as the evolution seems virtually identical to previous tests. Moreover, when I look at LAYER 1 THICKNESS in the results file, it appears that the layer thickness is uniform at 0.02 m.

Here are the relevant lines from my steering file; is there anything that could be causing the AL thickness to change back from my prescribed values in the user_bed_init subroutine?

ACTIVE LAYER THICKNESS = 0.020

CONSTANT ACTIVE LAYER THICKNESS = NO

NUMBER OF LAYERS FOR INITIAL STRATIFICATION = 3

LAYERS INITIAL THICKNESS = 0.020;0.020;0.020

As well, I’ve taken note of [name removed]'s reply. These both seem like helpful ideas, I will try them! Thank you so much for your help!

12

As a follow-up, I can confirm that [name removed]'s strategy of using the bottom roughness as a map for the active layer thickness is successful! I went into Blue Kenue (other preprocessors are available) and created a buffer zone at the inlet with [name removed]'s roughness = 0.02501, just off from 0.025.

For anyone who is also down this rabbit hole, I have attached my Fortran file (the _man is for [name removed]). The loop of interest is here:

DO ISTRAT=1,NUMSTRAT

        DO IPOIN=1,NPOIN

          IF ((ABS(CHESTR%R(IPOIN)-0.02501).LT.0.000001).OR.

     &    (ABS(CHESTR%R(IPOIN)-0.01001).LT.0.000001)) THEN

            ESTRATUM(ISTRAT,IPOIN) = 0.D0

          END IF

        END DO

      END DO

I have it set up to loop through all the layers (NUMSTRAT) and all the points (NPOIN). Maybe it would be best to only define the thickness of the first (top?) layer? If so, you can replace

DO ISTRAT=1,NUMSTRAT

with

ISTRAT=1

and remove the corresponding END DO statement at the end.

Note as well that this is set up to look for two bottom thickness values: 0.02501 (off from 0.025, main roughness) and 0.01001 (off from roughness of a smooth, non-erodible feature in my waterway, which I needed to delimit anyway!). Because of [name removed]'s clever code involving ABS(), if my understanding is correct, this will only activate the ALT modification for these two values of CHESTR and no others. If you have more roughnesses to add, follow the same pattern (don’t forget the & in column 6!):

IF ((ABS(CHESTR%R(IPOIN)-some_value).LT.0.000001).OR.

     &    (ABS(CHESTR%R(IPOIN)-new_value).LT.0.000001)).OR.

     &    (ABS(CHESTR%R(IPOIN)-another_value).LT.0.000001)

Remember to add this line to the Telemac steering file:

FORTRAN FILE='user_bed_init_man.f'

All this to say, the original problem of this topic has been solved* – no more ripples at the entrance – and I’ve figured out how to define a spatially variable bed, which I needed to learn anyway. It takes some extra preparation in the preprocessing stage, but it’s not much more complicated than it would be to create a Selafin variable with ALT values.

Thank you again, [name removed] and [name removed], for your willingness to help a poor student reach numerical stability!

  • That said, there is still a bit of strangeness at the limit of my non-erodible zone (about 1 m long). This is somewhat to be expected and I’ll just have to exclude this area from my analysis, which could be problematic because of how close this area is to my test subject… but that’s a problem for future me :notes:

user_bed_init_man.f (1.99 KB)

13

Renault

Hello, how did you finally solve this problem? Is it convenient to share the solution on the forum? I also met your problem in the process of simulation, but I didn’t get a solution. I asked you for advice and guidance.

[name removed]

14

Hello

Just to contribute to the topic, following [name removed]'s suggestion of using friction coefficient, this is how it worked in case of [name removed], version 9.0.

I have a “geo.slf” file with friction factor 0.09001 m in the riverbed and 0.09 m in the floodplain. I want different bottom layers in river/floodplain. The following code is inserted in “user_bed_init.f”, for 2 classes of sediment and 2 estrata:

USE DECLARATIONS_TELEMAC3D, ONLY : RUGOF



        DO IPOIN=1,NPOIN

           IF (RUGOF%R(IPOIN).GT.0.090005) THEN

                ESTRATUM(1,IPOIN) = 5.D0

                RATIO_INIT(1,1,IPOIN) = 0.D0

                RATIO_INIT(2,1,IPOIN) = 1.D0

                ESTRATUM(2,IPOIN) = 5.D0

                RATIO_INIT(1,2,IPOIN) = 1.D0

                RATIO_INIT(2,2,IPOIN) = 0.D0

           ELSE

                ESTRATUM(1,IPOIN) = 0.5D0

                RATIO_INIT(1,1,IPOIN) = 1.D0

                RATIO_INIT(2,1,IPOIN) = 0.D0   

                ESTRATUM(2,IPOIN) = 9.5D0

                RATIO_INIT(1,2,IPOIN) = 0.D0

                RATIO_INIT(2,2,IPOIN) = 1.D0

           ENDIF

          ENDDO