open TELEMAC-MASCARET The mathematically superior suite of solvers

Hello!

I want to know if it's possible to use the following Navier-Stokes and K-E in equations Telemac 3D (images attached - I couldn't see an option to inset equations). The idea is to model porosity as it varies vertically, representing different types of vegetation through the water column. I know porosity isn't an option in Telemac3D and so would need to do some programming myself.

The problem would require the inclusion of porosity terms in the Navier-Stokes and K-E equations. I can alter the source terms in the N-S equations using porosity values (i.e. head loss)but I need to include a varying porosity in the equation itself so as to avoid assuming a constant porosity - as a constant porosity isn't theoretically correct. I also need to alter the source terms in the K-E model with porosity, and so need to find out where they are calculated.

Is this possible within Telemac? I've previously modelled effects of porous media using just head loss but it's not sufficient. This approach would greatly help.

Many thanks!

Hello,

I personally do not think that the k-epsilon model is well suited for our applications which have large element sizes, too large to capture the real horizontal gradients. As a consequence, going into more sophistication may be deceiving.

Another point is that the porosity is treated as a multiplication factor at element level. It is well suited for classical finite elements which do variational formulations that result in integrals on elements. For other methods such as the distributive advection schemes, which are more and more used, this treatment is impossible, so I would not encourage people to develop new solutions based on porosity.

This being said, the principle is simple, you need to prepare an array giving a value between 0 and 1 for every element, this would be PHI, and when calling subroutine VECTOR or MATRIX you prescribe MSK=TRUE and give your array as the mask. It requires that you understand where the k-epsilon terms are computed, and possibly some terms will have to be split, as not all are multiplied by the porosity in your equations.

With best regards,

Jean-Michel Hervouet

Hello!

Thankyou for your reply.

So if I were to use, say, the second order PSi scheme it would be impossible to build the model on porosity for the argument above?

Would your method be used in the MASKEL subroutine?

Kind regards

Apologies for the second reply, but is there any documentation about the treatment of porosity at the element level? There's one section that describes porosity in Hervouet (2007), but I'm looking for recommendation particularly on the difficulty reconciling the above porosity treatment with distributive advection schemes. I'm not familiar with such schemes and want to get to grips with the theory.

Many thanks and kind regards.

Hello,

The distributive advection schemes derive all their computations from fluxes leaving points computed in a finite element way, thus taking into account initially the porosity. Then they just deal with fluxes between points and porosity is not considered, but as the initial fluxes are correct it could be OK. However no real test was done coupling distributive schemes and porosity and no theory was ever written. Moreover it is difficult to find quantitative test cases with porosity.

As for the treatment of porosity at element level the basic principle for dealing with porosity is that we consider that the available volume is reduced by porosity, hence the finite elements integrals computed on a triangle are multiplied by a coefficient which is the water in the volume divided by the total volume including obstacles. It means that if the porosity is constant the behaviour is unchanged, since the whole equation is multiplied by the same coefficient (as you know friction on obstacles has to be added by extra drag forces).

When adapting equations to porosity the difficulty is to find where to put the porosity coefficient, sometimes it will be inside derivatives, sometimes outside, and other terms like a prescribed discharge are not impacted. It is so true that the very first publications on porosity in shallow water equations were obviously wrong...

With best regards,

Jean-Michel Hervouet

Good afternoon!

Thankyou for the above. That clarifies a lot things.

I want to try porosity so as to demonstrate how Telemac deals with it, given the caveats. The initial part of the approach is easy enough, however I'm not hugely experienced at fortran coding so apologies for the novice questions.

After defining the array I need to call the SUBROUTINE MATRIX. Reading the programming guide not all the terms need to be used, so it it true to say that I only need to define MSK = 1 and MASKEL = array with PHI, or are there required terms?

The main part I'm unsure about is how to implement the mask so it affects the N-S terms. It's not clear to me what extra code is necessary to take the mask and apply it to the terms, and unfortunately I've not found any codes with similar approaches on the forum. From what I can see, the K-E source terms I would need to adapt are in VISCK (viscosity), KEPCL3 (KBOR, EBOR, and AUBOR), and SOUKEP (source terms). What is the step I'm missing here?

Kind regards

Hello,

This would be rather MSK=.TRUE., but MASKEL is of type BIEF_OBJ, if you want to create it locally you would have to declare it as such, and then to allocate it by calling BIEF_ALLVEC, giving its name, etc., and element 40 which corresponds to P0 prisms, which means one value per element (this allocation must be done once* and you should add a SAVE command if you are in a subroutine so that your array is kept for the next call). This mask is to be called with MATRIX and possibly also VECTOR, depending on your equations.

* possible implementation for allocating once :
TYPE(BIEF_OBJ) :: MASKEL
DATA ALREADY/.FALSE./

IF(.NOT.ALREADY) THEN
CALL BIEF_ALLVEC(1,MASKEL,'MASKEL',40,1,1,MESH)
ALREADY=.TRUE.
ENDIF


With best regards,

Jean-Michel Hervouet

Hello,

Thank-you for your reply!

I had assumed that MASKEL could just simply be prescribed PHI (e.g. MASKEL =0.13), but I see that's not the case? In this case, I could write a code for a simple mass matrix such as:

TYPE(BIEF_OBJ) :: MASKEL
DATA ALREADY/.FALSE./

IF(.NOT.ALREADY) THEN
CALL BIEF_ALLVEC(1,MASKEL,'MASKEL',40,1,1,MESH)
ALREADY=.TRUE.
ENDIF

FORMUL=MATMAS
CALL MATRIX(AM1,'M=N ',FORMUL,IELMH,IELMH,
& SL1,TE3,S,S,S,S,S,MESH,MSK=.TRUE.,MASKEL)
SAVE

(from a Telemac 2D script - I guess I would have to alter it?)

If I wanted to have a mask that varies at different levels (e.g. unique values on 1-4th horizontal levels, then nothing above) then would I define MASKEL before the BIEF_ALLVEC is called? For example,

DO IPLAN = NPLAN-17,NPLAN-14
I3D = I+NPOIN2*(IPLAN-1)
IF (IPLAN==(NPLAN-17))THEN
MASKEL = PHI1
ELSE IF (IPLAN==(NPLAN-16))THEN
MASKEL = PHI2
ELSE IF (IPLAN==NPLAN-15) THEN
MASKEL = PHI3
ELSE IF (IPLAN==NPLAN-14) THEN
MASKEL = PHI4
ENDIF
END DO

and then use the above BIEF-OBj call code, or similar? I don't get any errors from my current code so Telemac is accepting it, so it's tricky to see what works but clearly I've been doing it incorrectly.

Hello,
See my comments and corrections in the code below,
With best regards,
Jean-Michel Hervouet



TYPE(BIEF_OBJ) :: MASKEL
DATA ALREADY/.FALSE./
SAVE ! to store MASKEL for further calls

IF(.NOT.ALREADY) THEN
CALL BIEF_ALLVEC(1,MASKEL,'MASKEL',40,1,1,MESH)
ALREADY=.TRUE.
ENDIF
!!!!!!!! YOU SHOULD DEFINE MASKEL HERE BEFORE THE CALL TO MATMAS !!!!!!!!!!!
OR EVEN WITHIN THE IF(.NOT.ALREADY)... ENDIF

initialisation of MASKEL outside NPLAN-17 NPLAN-14 ?
I suggest :
CALL OS('X=C ',X=MASKEL,C=1.D0)

then maskel is defined per element, so NELEM2 instead of NPOIN2

DO IPLAN = NPLAN-17,NPLAN-14
DO I=1,NELEM2
I3D = I+NELEM2*(IPLAN-1)
IF (IPLAN==NPLAN-17)THEN
MASKEL%R(I3D) = PHI1 ! here I suppose that PHI1 is a double precision number, not an array.
ELSE IF (IPLAN==NPLAN-16)THEN
MASKEL%R(I3D) = PHI2
ELSE IF (IPLAN==NPLAN-15) THEN
MASKEL%R(I3D) = PHI3
ELSE IF (IPLAN==NPLAN-14) THEN
MASKEL%R(I3D) = PHI4
ENDIF
ENDDO
ENDDO


SL1= ?????? coefficient to multiply the result, if not SL1=1.D0, but depending on what you do it could be 1.D0/DT as well

FORMUL=MATMAS
CALL MATRIX(AM1,'M=N ',FORMUL,41,41,SL1,S,S,S,S,S,S,MESH,MSK=.TRUE.,MASKEL)

Hello!

I've been running this version of the code but I keep getting a "RETURNING EXIT CODE: 2" after Telemac starts to run the simulaiton.

I'm pretty sure it's related to my declarations. The only one I can see being incorrect is MESH: in the programming guide MESH is stated as a BIEF_OBJECT. Is this correct for the code, or should ti rather be declared as something else? I've seen MESH declared under TELEMAC2D. but that's not appropriate here.

USE BIEF
USE DECLARATIONS_SPECIAL
USE DECLARATIONS_TELEMAC3D, ONLY: PRIVE2D, NPLAN
IMPLICIT NONE
INTEGER :: I , IPLAN, I3D, NELEM, NELEM2
DOUBLE PRECISION:: PHI1, PHI2, PHI3, PHI4, SL1
TYPE(BIEF_OBJ) :: MASKEL
TYPE(BIEF_MESH) :: MESH
LOGICAL :: ALREADY
DATA ALREADY/.FALSE./
IF(.NOT.ALREADY) THEN
CALL BIEF_ALLVEC(1,MASKEL,'MASKEL',40,1,1,MESH)
ALREADY=.TRUE.
ENDIF

CALL OS('X=C ',X=MASKEL,C=1.D0)
SL1=1.D0
PHI1= 0.13
PHI2= 0.47
PHI3= 0.76
PHI4= 0.91

DO IPLAN = NPLAN-17,NPLAN-14
DO I=1,NELEM2
I3D = I+NELEM2*(IPLAN-1)
IF (IPLAN==NPLAN-17)THEN
MASKEL%R(I3D) = PHI1
ELSE IF (IPLAN==NPLAN-16)THEN
MASKEL%R(I3D) = PHI2
ELSE IF (IPLAN==NPLAN-15) THEN
MASKEL%R(I3D) = PHI3
ELSE IF (IPLAN==NPLAN-14) THEN
MASKEL%R(I3D) = PHI4
ENDIF
ENDDO
ENDDO
RETURN
END