0
$\begingroup$

The original post was on stackoverflow : I transfert it here.

I have to solve numerically the advection equation with periodic boundaries conditions : u(t,0) = u(t,L) with L the length of system to solve.

I start also with u(0,x) = uexacte(0,x) = sin(2*pi*x/L)

Here the main part of the code with loop time (we use here FTCS scheme) :

V=1
L=1
# analytical solution --------------------------
def uexacte(t,x):
    return sin(2*pi*(x-V*t)/L)

# 1. Centre FTCS (Forward Time Centered Space)
cfl = 0.25
nx = 10
tend = 1

#
dx = L/(nx-1.)
dt = cfl*dx/V
nt = int(tend/dt)+1
print "CFL=%5.2f tend=%4.1f --> %i iterations en temps"%(cfl,tend,nt)

# Arrays
x = linspace(0,L,nx)
# Bounadry condition
u0 = uexacte(0,x)

# Starting solution
t=0.0 ; u=copy(u0)

# Time loop
for i in range(1,nt):
    # FTCS
    #u[1:nx-1] = u[1:nx-1] - cfl/2*(u[2:nx] - u[0:nx-2])
    # Using roll
    u = u + - cfl/2*(roll(u,-1)- roll(u,1))
    # Update time
    t = t+dt

I don't understand the solution given by teacher who uses the python function roll in this way :

# Using roll
u = u - cfl/2*(roll(u,-1)- roll(u,1))

One says that with the using of roll, we are sure to respect the periodic boundary conditions but I don't understand why ?

Indeed, my first approach was to do :

u[0] = u[nx-1]
u[1:nx-1] = u[1:nx-1] - cfl/2*(u[2:nx] - u[0:nx-2])

but this doesn't work and I don't know how to implement this periodic conditions in this way (without using roll function).

If someone could explain this matter and the trick with roll function, this would be nice to tell it.

UPDATE 1 :

I tried with classical approach (simple recurrence formula ) like this :

# Time loop
for i in range(1,nt):
    # FTCS
    u[1:nx-1] = u[1:nx-1] - cfl/2*(u[2:nx] - u[0:nx-2])

    # Try to impose periodic boundary conditions but without success
    u[0] = u[0] - cfl/2*(u[0] - u[nx-1])

    # Update time
    t = t+dt

Indeed, the result is bad (values for each side are not the same). I could impose the theorical values at each step but in practise, we don't know always the analytical solution.

What is the trick to impose this periodic boundary condition on the numerical solution at each step of time ?

The problem with periodic boundary conditions is formulated as :

enter image description here

UPDATE 2 :

The solution given by LonelyProf below works fine. Thanks

$\endgroup$

closed as off-topic by David Ketcheson, Anton Menshov, nicoguaro Dec 11 '18 at 16:10

This question appears to be off-topic. The users who voted to close gave this specific reason:

If this question can be reworded to fit the rules in the help center, please edit the question.

  • $\begingroup$ I have a general query, unrelated to my answer (which addressed only the original question, concerning periodic boundaries). You seem to be using the FTCS scheme for the advection equation. Isn't this unstable? See en.m.wikipedia.org/wiki/FTCS_scheme $\endgroup$ – LonelyProf Dec 3 '18 at 21:38
  • $\begingroup$ yes, you're right but I think if I choose a little CFL, I have, at least at the beginning, the numerical solution closed to analytical solution $\endgroup$ – youpilat13 Dec 3 '18 at 21:39
  • $\begingroup$ With the updates, this question has definitely moved into the realm of "debug my code" and doesn't seem to be applicable to others any longer. $\endgroup$ – David Ketcheson Dec 11 '18 at 11:43
2
$\begingroup$

Similar issues arise in other languages, so this is not Python-specific.

The NumPy function roll applies a cyclic shift to the array given as its first argument.

  • roll(u,-1) returns an array with elements u[1],u[2],u[3]...,u[nx-1],u[0]
  • roll(u,1) returns an array with elements u[nx-1],u[0],u[1]...,u[nx-3],u[nx-2].

This is exactly what you want on the right hand side of your assignment, accessing the elements on either side of every element of u, using periodic boundaries.

You could do this by array slices, but then you would need to cope with the end effects by hand. So you can write a single statement to update all of u[1:nx-2] in terms of u[0:nx-3] and u[2:nx-1], and write separate statements to update u[0] in terms of u[nx-1] and u[1], and to update u[nx-1] in terms of u[nx-2] and u[0].

Currently, you are making mistakes in the way you try to do this by hand (using an illegal index nx, for instance, or simply writing the wrong indices on the right).

Another issue is that using a function such as roll enables you to do the update effectively simultaneously for all the elements of u. In some applications, perhaps yours, it may be an error (and is certainly different) to update any of the elements, and then use the updated value on the right-hand side to update other elements. You can avoid this by defining temporary variables to hold the updated values, if you are sticking with the "by hand" approach. Or define an array du to hold the change in u, calculated from the current values, and actually perform the update u=u+du in a separate statement, all at once. Using the roll function (or similar) allows the code to be more compact and simpler.


[EDIT following OP comments]

Here is an example of what I mean by using temporary variables.

utemp1=u[0] - cfl/2*(u[nx-1] - u[1]) utemp2=u[nx-1] - cfl/2*(u[nx-2] - u[0]) u[1:nx-2] = u[1:nx-2] - cfl/2*(u[0:nx-3] - u[2:nx-1]) u[0] = utemp1 u[nx-1] = utemp2

$\endgroup$
  • $\begingroup$ thanks for your answer, could you take a look please at my UPDATE 2, I don't know where is my error $\endgroup$ – youpilat13 Dec 3 '18 at 20:16
  • $\begingroup$ You seem to have extended the question somewhat, @youpilat13 . I am not sure that this site is intended to help debug your code. Can I first ask whether you understood the last paragraph of my answer? You are updating values u[1] and u[nx-2], before using them in the update for u[0] and u[nx-1]. What happens if you use the roll function? $\endgroup$ – LonelyProf Dec 3 '18 at 21:31
  • $\begingroup$ -@LonelyProf I am near to the solution but what do you mean by So you can write a single statement to update all of u[1:nx-2] in terms of u[0:nx-3] and u[2:nx-1], and write separate statements to update u[0] in terms of u[nx-1] and u[1], and to update u[nx-1] in terms of u[nx-2] and u[0]. ?, I mean what's the traduction in terms of python code ? $\endgroup$ – youpilat13 Dec 3 '18 at 21:35
  • $\begingroup$ Basically, I mean what you have written. Except that you should write the update equations for u[0] and u[nx-1] first, storing the results in temporary variables; then update the array slice; and only afterwards, actually update u[0] and u[nx-1], using the temporary variables just calculated. $\endgroup$ – LonelyProf Dec 3 '18 at 21:43
  • $\begingroup$ I have edited my answer in an attempt to clarify it. I hope that it is helpful. Since this site is definitely not intended to help with coding problems in specific languages, I think that's as far as I can go, in that direction. $\endgroup$ – LonelyProf Dec 3 '18 at 22:03

Not the answer you're looking for? Browse other questions tagged or ask your own question.