open TELEMAC-MASCARET The mathematically superior suite of solvers

Hi,

I made a small sketch to discuss the way how boundary conditions are implemented in the Finite Volume Methods. On the left had side I defined a simple inflow boundary (red line).




The sketch shows 3 possible control volumes that belong to the boundary (BCV 1, BCV 2, BCV 3). I've plotted also two different ways how the segments for the inflow might be defined and the outward normal vector. First I thought that the inflow is separated into two boundary segments to handle cases such as BCV 2 and BCV3. But I couldn't find this in the code.

I hope that the control volumes at the boundary are drawn correctly?


If I look in the HLLC solver (I think it is not different in the other FVM schemes), the boundary flux is computed in a subroutine called CDL_HLLC. The length of a boundary segment (VNL) is defined via the unit outward components of a boundary node.



The fluxes are computed via the HLLC solver and the result is then multiplied by VNL



Looking at the two possibilities how the segment length might be defined, there can occur different problems. For case A, the segment length must not be identical with the segment length that I expect for the BCV 1, so the computed inflow will be wrong. For case B, the length will be correct.
However, both cases can not resolve BCV 3. Case B might be able to work for BCV 2.

I'm not sure if I made a mistake in my sketch?

With best regards,
Leo

If I print out VNL and compare it with the measured lengths in my mesh it looks like case B is implemented.

However, I could not figure out how the length (VBL) at the corner is computed.

But now I think I figured it out

Sorry for posting so often, I think the BCV's of my sketch were wrong. The CV goes normal through the boundary and I've computed a wrong segment length.

For that case, the flow between to BCV's across the boundary is zero for a straight boundary (flux normal into the domain).

Maybe I can improve the sketch and include the vectors at BCV 2 and BCV 3.

Best Regards,

Leo

Locking back at the previous example. I looked at the normalized normals at the boundary nodes near the corner (XNN,YNN in CDL_HLLC). The corner node has a component in x and y direction.




Also the velocity at the corner has a component in the y direction.

Hello Leo

You will find in the attached file the subroutine debimp.f that fixes the bug on the imposed discharge.
The problem was in the definition of the mask that tells to subroutine vector which segments have to be considered in the integration by part. Unlike FE, first and last segments are added to the mask.
The definition of the outward normal in the corner can be done with several manners and for Telemac FV, we (first developpers) have chosen to consider the mean of the normals of the two adjacent segments weighted by the length of these latters. For FE, other definition is taken into account.
A final remark, the accuracy of the imposed discharge is linked to the velocity profile you have chosen. The most accurate option (for complex sheped boundaries) is option 2 which implies a given distribution in the boundary conditions file. This conclusion is based on some numerical tests I have achieved on the "confluence" test case. More in deep investigations are necessary to have a final conclusion and, for sure, to improve the accuracy.
This file will be commited to the svn trunk as soon as possible.

With my best regards

Riadh

Hi Riadh,

thank you very much!

Leo

Hi Riadh,

this means that the actual implemented Finite Volume Schemes are limited to problems where no inflow boundary must be imposed via q.


Most (if not all) finite volume models for the SWEs use cell centered methods with triangles or rectangles. This means that the problem of defining q on the boundary is much simpler (e.g. via ghostcells). As far as I know from a friend, the Riemann invariants of the one-dimensional SWEs are usually used to compute the fluxes on the boundary.

This is also done in the subroutine FLUSEW when H or V is imposed. For an imposed q it might be necessary to iteratively solve for h and u at the boundary (see Song et al. "A robust well-balanced finite volume model for shallow water flows with wetting and drying over irregular terrain", Advances in Water Resources).

This leads to two questions:

Is it necessary to replace the calculations done in debimp.f with an iterative procedure to impose q correctly?

Is it mathematically proven that it is possible to handle corner values (all possible angles) as it is actually done? I'm not sure since the boundary itself is no longer 1D.


With my best regards

Leo

Hello Leo

Grateful for these nice discussions !

The ghost particle approach is not exclusively used with cell centered FV, we use exactly the same way here by considering the sum of the fluxes through all the edges linked to a boundary node (stored in vector CE (see routine cdl_hllc. or cdl.)) and then we use a fictitious (ghost) particle to which we give a state vector Ue=(h,hu,hv)_e that ensures the imposition of the desired boundary condition.
Furthermore, in Telemac also we use Riemann invariants to impose boundary conditions. Even in FE, the thompson BC is based on the characteristic theory which assumes mainly that the Riemann invariants is conserved on the incoming (or outgoing) characteristics.
For FV, only kinetic schemes handle exhaustively all the possible cases (see subroutine cdl.f). For the remaining schemes, only imposed discharge and elevation are implemented (see subroutine flusew.f and cdl_hllc.f for instance). and as I have notified it in the header of flusew.f several cases need to be implemented, and I'm working on it these days.
For 2D case, the use of Riemann invariants do not lead to an iterative process unlike the 1D case. This is due to the use of conservative variables (Q,S: discharge and wetted area) for 1D model.
A final remark, about corners, as far as I know, there is no mathematical study dedicated to this question. Actually, the way to define the normal is a full responsability of the code developer, and I think that there is no major problem while the correct discharge is entered (or exited) to (from) the domain. Only very local disturbance can be observed in the vicinity of the corner depending on the shape of the liquid boundary(straight, bended, with 90° angle or with higher or lower angle) and on the velocity profile chosen by the user.

I hope that I was clear enough, otherwise, do not hesitate to go on with this topic :)

With my best regards,

Riadh

Hi Riadh,

oh I haven't heard about the ghost particle approach. I was confused because when the mesh is discretized the "Telemac way" a node must have more than one ghost cell. But I think I can understand this now.

It is no problem when the velocity distribution is not 100% perfect as long as the mass that enters the domain is correct. Actually the inflow (mass) is still not 100% correct in my model. But the corner problem is not longer a real problem for me

I've learned again that it is hard to try to define Neumann like boundary conditions for the SWEs and that it is necessary to accept that it does not work as simple/accurate as for other equations (Richards-equation, heat-equations).


The Riemann invariants for the discharge are used by the 2D HLLC models. As far as I understood it can be handled like the rotation invariance in the HLLC scheme (we always solve a 1D problem). I haven't finished reading the papers but it is quite similar to that what in flusew.f happens when h or v is imposed. However, in flusew.f only h is set when the inflow q is given. At this point the other models use an iterative method to compute h and v,u from a given q via the Riemann invariants.

In Telemac the velocity is given (comes from debimp.f) and only h is computed. I do not understand the following part in flusew.f (R/R1 come form the invariants)


...


UI and VI are computed but not used anymore? Inserting R and RLAMBO to compute HI leads to HI = HJ. At this point the other models do the iteration and compute HI,UI, and VI from q. Maybe I understand this tomorrow. But now I'm too tired, sorry for spelling errors


Thank you,

Leo

I looked in the paper, and the condition for a given H at the boundary is implemented as written in the paper.

HI is the value at the boundary (taken from AMINF), R1 is computed with the variable of h,u,v inside the domain.



For an inflow boundary q is given at the boundary


and inserted into the following formula


which is then solved iteratively.


However, in Telemac the corresponding velocities of a given q are computed in debimb. (No iterative routine)

I think there is an error in the formula that is used to compute h from the given velocities. Looking in the paper, RLAMBO the velocity at the boundary should be used to compute RLAMBO. I will test this later.