Periodic boundary conditions for compressible Euler equations
 

Hi,

I am currently building up a 2D hydrodynamics code. Specifically, I am numerically solving the 2D Euler equations in order to simulate Kelvin-Helmholtz instability. The boundary conditions are periodic, however I am running into issues when enforcing such conditions.

Given the governing equations






I define Q to be the left-hand side of the system of PDEs. As such in discrete form its a (nx,ny,4) array. Given we have periodic boundary conditions I set

Q(1,1:ny,1:4) = Q(nx,1:ny,1:4) (Q(1) = Q(n) in x)
Q(1:nx,1,1:4) = Q(1:nx,ny,1:4) (Q(1) = Q(n) in y)

With this implementation I am not getting the desired effects of the Kelvin-Helmholtz. I conclude the error is at the boundary due to the plot. I start my arrays at 1 so my boundaries are at 1 and nx.

I am assuming a uniform grid and am using the Rusanov (Local Lax Friedrichs) flux discretization, so there shouldnt be any special treatment of ghost points needed. I have done periodic boundary conditions for non-systems of PDEs and the implementation is straight forward however, I am suspecting that I way of doing PBCs is not correct given I have a system.

NOTE 1: I can post my Fortran code if it will help.
NOTE 2: I have tested this code on various 2D Riemann problems and I get correct results.

"I define Q to be the left-hand side of the system of PDEs." :confused:

Maybe the RHS? Have you tried to set periodicity on the fluxes?

Periodic BC
 

Hi Denaro,

Sorry for not being more precise. I meant to say that Q is the left-hand side (time-derivative terms) of the hyperbolic system. In terms of applying the boundary conditions to the fluxes, I have tried that with little luck.

The flux is a simple Rusanov flux.

As I wrote, you have to set the same fluxes

If you are resetting the boundary values after the flux evaluation, then as FMDenaro points out, your flux values won't match the Q values. So you can either reset the flux values or call your BC routine prior to the flux evaluation.

Periodic BC
 

Hi Guys,

Thanks for the comments. I will try it out when I get to a computer.

Results of PBC
 

Hi,

So if I'm understanding what you said, I tried called my BC subroutine before the flux evaluations:

Is this what you were suggesting? With that said in my main time loop, do I need to call the BCs a second time? For my 2D Riemann problem which I had non-periodic BC, I did that.

Immagine a flux vector fluxes_species(1:N) and the periodicity being between 1 and N and where fluxes_species(i) is at the left side of the components species(i). So the update of the species(1) depends on the difference
fluxes_species(2)-fluxes_species(1)

where fluxes_species(1) lies on the left boundary and must be computed.
When you arrive to compute the update of species(N) you should compute the difference
fluxes_species(N+1)-fluxes_species(N)

that setting the periodicity becomes
fluxes_species(1)-fluxes_species(N)

You want the Q values to match, prior to the evaluation of the fluxes F and G.

I do not understand the vector you pass to the subroutine BC. You send RHS which is the sum of the fluxes for each element, it is equivalent to say dQ1/dt = dQn/dt, the update of conservative variables in the first and last cell are equal, which is not true in general.

As FMDenaro suggested, you have to enforce periodicity on the fluxes of the boundary of the domain, not on the sum of fluxes of bounding elements (RHS). What goes out on one side enters the other side.

If you use ghost cells you can very easily enforce periodicity on the conservative variables before evaluating the fluxes the normal way. Q_in(0,:,:)=Q_in(nx,:,:) (ghost cell left = last layer of real cells), the same for cell layers i=nx+1, j=0, j=ny+1.

Have a look to the book of Blazek (COMPUTATIONAL FLUID DYNAMICS:
PRINCIPLES AND APPLICATIONS), section 8.7 periodic boundary conditions

Periodic BC
 

Hi all,

Thank you for the helpful replies. I have opted to try and implement the periodic boundary conditions via ghost points. With that said I have some question:

1. When using ghost points in this context my arrays are (nx,ny) for Q, (nx-1,ny,4) for F_half, and (nx,ny-1,4) for G_half. Does that mean My arrays for Q have to be (-1:nx+1,-1:ny+1,4)? the -1 in the array is to account for the -1 ghost point.

In terms of calling my periodic BC routine prior to the total flux evaluation (RHS), I have tried doing the following

Fh(1,:,: ) = Fh(nx-1,:,: )
Fh(:,1,: ) = Fh(:,ny,: )

Gh(1,:,: ) = Gh(nx,:,: )
Gh(:,1,: ) = Gh(:,ny-1,: )

I choose nx-1 and ny-1 for Fh and Gh respectively as they are sizes (nx-1,ny,4) and (nx,ny-1,4). Yet, I still do not get the correct results.

If you are using colocated grid, the things are the same in x and y direction. If the periodicity in one general direction is between node 1 and node N then you have for a general function the links

f(0)=f(N-1)
f(-1)=f(N-2)
....

and

f(N+1)=f(2)
f(N+2)=f(3)
....

Periodic BC
 

Hi Dr. Denaro,

That clears up a lot of questions. I'm working on a uniform grid currently. So given I have a first order accurate discretization (in space) I have points

which means in terms of periodicity I have

which means my arrays will need to be structured

Periodic BC
 

Hi everyone,

After some odd days I am still sheepishly coming back with no luck. I feel like I have tried almost everything. Attached is my main time loop and my RHS with my most recent attempt. The two code snippets are where I am enforcing boundary conditions. For purposes of clarity I am using the Shu TVD Runge-Kutta 3.

Main time loop for advancing the solutions

Sorry to say that I think you have not programmed correctly the BC.s.
I give you a very simple 1D example. Consider the domain going from node 1 (x=0) to node N (x=L). Therefore, the periodicity is such that any function f has f(1)=f(N). As a consequence, in terms of ghost nodes you have f(0)=f(N-1) and f(N+1)=f(2).

Let us work with a first order FD discretization of the linear wave equation (u>0):

fold will be prescribed for i=1,N

fnew(1)=fold(1) -u*dt*(fold(1)-fold(N-1))/dx
do i=2,N-1
fnew(i)=fold(i) -u*dt*(fold(i)-fold(i-1))/dx
end do
fnew(N)=fnew(1)

then swap the vector fnew to fold and continue the cycle in t.

A finite volume discretization is slightly different. You have to work on the flux function to be computed on each face and inserted in the equation

fav_new(i)=fav_old(i) -dt*(flux(i+1/2)-flux(i-1/2))/dx

so that you prescribe the periodicity on the vector flux

Periodic BC
 

Hi FM Denaro,

Thank you for your patience on this thread. I will digest this and try to implement this in my solver.

EDIT: So what I'm doing wrong is that I need to be evaluating the flux "manually" i.e. what you did for fnew outside the interior loop?

If you do not want to treat explicitly the equation wher the BC enters into, then you have to use ghost points. Therefore, you start at t0 prescribing fold for i=0:N+1 then compute fnew(1) using the general scheme in the interior. Do not forget to explicitly set the BC for fnew(0), fnew(N+1).

Ghost points
 

Hi again,

I have implemented your suggestion on ghost points, but again I am getting incorrect results. To be sure I have posted my code.

The RHS function with the flux evaluated at the boundary

In order to compensate for the ghost points, r, u, v, E, p, H, c, Q_in, RHS, Fh, Gh, Q1, Q2, Qn, and Q are all (0:nx+1,0:ny+1,4). Given nx+1 and 0 are technically not in my domain, I am not taking those indices into account when looping. I am quite stumped on what is possibly the issue at this point as I have tried everything suggested on the threat (unless I am missing something?)

I have uploaded a screenshot at T = 1.0. As you can see on the right, there is build-up when it should just be periodic. The density is on a 128 x 128 resolution. Its not an optimal resolution but I'm doing it for speed/debugging purposes.

EDIT: So I removed the boundary conditions and I get the same result. Furthermore, the arrays are unchanged when I remove the BC.

I think you are mixing up things. As explained in the earlier answers, you have two ways of enforcing periodicity, either on the fluxes OR on the state variables through ghost points.

If I understood correctly your indices, your physical domain is decomposed into finite volumes from 1 to nx and from 1 to ny. Your ghost cells are then located at 0, nx+1, 0, ny+1. And your faces go from 1 to nx+1, 1 to ny in the x direction and from 1 to nx, 1 to ny+1 in the y direction.

Beforce calling the routine RHS, set the ghost cell boundary condition on the vector Q. Not Q1 ! That is to say on Qold and not Qnew.

Q(0,:,: )=Q(nx,:,: )
Q(nx+1,:,: )=Q(1,:,: )
Q(:,0,: )=Q(:,ny,: )
Q(:,ny+1,: )=Q(:,1,: )

That's it.

Then you construct your fluxes from face 1 to nx+1, 1 to ny in the x direction and from 1 to nx, 1 to ny+1 in the y direction. That way, in the x direction the flux computed at the face i=1 uses Q(0,:,: ) and Q(1,:,: ) and the flux computed at the last face nx+1 uses Q(nx,:,: ) and Q(nx+1,:,: ). And you do not need to set periodicity on Fh and Gh.

Have a look to the book I mentionned, it is pretty well explained.

Pbc
 

Hi Francois,

You are indeed correct. 0, nx+1, and ny+1 are the ghost points. Upon reading the literature reference you gave me and the book by Leveque I implemented what I believe to be the correct way to set the boundary conditions.

The problem I'm seeing when I print out the arrays and through the plots is that the boundary conditions are not being satisfied, i.e. when I remove the boundary conditions of Q I get the exact same result. I believe something else maybe wrong in my code.

I've restructured some of my solver similar to what Leveque does for boundary conditions. Thank you for the patience and time.