unstable for 1D advection equation by CG

Question

Hi
I am trying to do the time depend problem by explicit RK4 method and starting by 1D advection equation in periodic UnitInterval. but it is very unstable even if I try very tiny time step size to full fill CFL number limitation. actually I tried to do it manually in matlab, which is very stable while I choose the same mesh size and time step size.

here is code:

from dolfin import *
mesh = UnitInterval(128)
V = FunctionSpace(mesh,'CG',1)
c0=Expression('exp(-pow((x[0]-0.5),2)/0.025)')

# Sub domain for Periodic boundary condition
class PeriodicBoundary(SubDomain):

    def inside(self, x, on_boundary):
        return bool(x[0] < DOLFIN_EPS and x[0] > -DOLFIN_EPS and on_boundary)

    def map(self, x, y):
        y[0] = x[0] - 1.0

# Create periodic boundary condition
pbc = PeriodicBoundary()
bc = PeriodicBC(V, pbc)

c=TrialFunction(V) 
v=TestFunction(V)

a_K=-nabla_grad(v)*c*dx
# a_K=-nabla_grad(c)*v*dx  
a_M=c*v*dx

M=assemble(a_M)
bc.apply(M)
K=assemble(a_K)

c=Function(V)

c_1=interpolate(c0,V)
c.assign(c_1)

k=[]
for i in range(4):
    k.append(Function(V))

T = 4
t = 0
dt =0.001

while t<=T:
    b1=K*c.vector()
    bc.apply(b1)
    solve(M,k[0].vector(),b1)

    b2=K*(c.vector()+0.5*dt*k[0].vector())
    bc.apply(b2)
    solve(M,k[1].vector(),b2)

    b3=K*(c.vector()+0.5*dt*k[1].vector())
    bc.apply(b2)
    solve(M,k[2].vector(),b3)

    b4=K*(c.vector()+dt*k[2].vector())
    bc.apply(b2)
    solve(M,k[3].vector(),b4)

    c.vector()[:]= c.vector()+dt*(k[0].vector()+ 2 * k[1].vector() + 2 * k[2].vector() +k[3].vector())/6.0

    t += dt
    plot(c,rescale=False)    

Comments

Simplify your code to increase the chances of getting and answer.

---

In addition to Gath's comment, edit your question using the code environment, please. Your post is currently unreadable because some syntax symbol (for example `*') has been parsed by this website engine.

---

Hi
sorry for that. I am a fresh man here not used to that.
actually I tried both weak form and strong form. both are not stable.

---

Are you aware of a mathematical reason why this should work without added diffusion/stabilization?

Did you test that bc works? You could also try to use newer approach for periodic BCs implemented by directly constraining function space in 1.2.0 and development version.

---

I am not very sure. but one thing strange is that I compared the mass and derivative matrix for both Fenics and my maltab code. they are the same, and eigen value of operator inv(M)*K are the same as well . here is the code I used for compute the eigen value.

A=M.array()
A=np.linalg.inv(A)
A=np.dot(A,K.array())
ee,ev=np.linalg.eig(B)
re=np.real(ee)
im=np.imag(ee)
import matplotlib.pyplot as plt
fig, ax = plt.subplots()
ax.plot(a, b, 'o')

except the results.. that is very strange. I don't know if it could be the algebra solver?

---

Well, I can reproduce it with FEniCS 1.1.0 with PETSc backend (see figures) and with FEniCS 1.0.0 with uBLAS backend.

---

I compared the mass and derivative matrix for both Fenics and my maltab code.

Did you compared them after an application of BCs?

---

yes. I checked both. and for eigen value I use the mass matrix with BC while derivative one not. everything is the same in matlab and Fenics.

---

There's just some Chinese placeholder image in place of the figures in your last post.

Can you convert the code for the new implementation of periodic BC and run in 1.2.0 or development version? The implementation of BC is maybe crucial for stability.

Answer 1

Problem is probably in an implementation of periodic BC in DOLFIN 1.0.0 and 1.1.0 which seems being very susceptible to instability. Your code updated to new implementation which directly restricts function space works fine in 1.2.0 or development version (see figures).

from dolfin import *
mesh = UnitInterval(128)
c0=Expression('exp(-pow((x[0]-0.5),2)/0.025)')

# Sub domain for Periodic boundary condition
class PeriodicBoundary(SubDomain):

    def inside(self, x, on_boundary):
        return bool(x[0] < DOLFIN_EPS and x[0] > -DOLFIN_EPS and on_boundary)

    def map(self, x, y):
        y[0] = x[0] - 1.0

# Create function space constrained by periodic BC
pbc = PeriodicBoundary()
V = FunctionSpace(mesh, 'CG', 1, constrained_domain=pbc)

c=TrialFunction(V)
v=TestFunction(V)

a_K = -v.dx(0)*c*dx
a_M = c*v*dx

M=assemble(a_M)
K=assemble(a_K)

c=Function(V)

c_1=interpolate(c0,V)
c.assign(c_1)

k=[]
for i in range(4):
    k.append(Function(V))

T = 4
t = 0
dt =0.001

while t<=T:
    b1=K*c.vector()
    solve(M,k[0].vector(),b1)

    b2=K*(c.vector()+0.5*dt*k[0].vector())
    solve(M,k[1].vector(),b2)

    b3=K*(c.vector()+0.5*dt*k[1].vector())
    solve(M,k[2].vector(),b3)

    b4=K*(c.vector()+dt*k[2].vector())
    solve(M,k[3].vector(),b4)

    c.vector()[:]= c.vector()+dt*(k[0].vector()+ 2 * k[1].vector() + 2 * k[2].vector() +k[3].vector())/6.0

    t += dt
    plot(c,rescale=False)

Comments

Thanks. sounds great!
could it be the problem of mass matrix when apply this BC? because I tried spetrum element in myself. and it works fine

---

Let's have DOF1 and DOF2 subject to periodic BC. Former approach I guess worked like: set row of a matrix corresponding to DOF1 to (0, 0, ..., -1, 0, ..., 1, 0, ...) with -1 at column DOF2 and 1 at diagonal entry; set RHS at DOF1 entry to 0. Current approach: eliminate one DOF2 directly from function space. Maybe there's some good numerical reason why latter works better. You should ask the authors of the implementation.

---

yes, the former one sounds working as what you mentioned. it seems one can't check how the later affect directly. but the thing is that I apply the exact the same one as the first method in my matlab code. as eigen value of is the same for both FEnics and Matlab. That is to say the instability of algebra equations should be the same. then how to explain so?

/Wei

---

eigen value of is the same

Which one?

---

OK, I can finnaly handle the figures. I plot eigen value of both matlab and Fenics codes while mass matrix is applied the BC as you mentioned before

---

it seems one can't check how the later affect directly

This is hidden somewhere in dolfin/fem/DofMapBuilder.cpp.

---

Fenics 1.1.0 and my matlab code give same eigen value of inv(M)*K. while M is already enforce the periodic bc by add last colume to the first. and set the last one as (1,0,0,0...,-1)

---

I guess that particular implementation of LU solver may affect whether numerical error in respective matrix subspace is created to cause a growth of instability. What is condition number of $M^{-1}K$?

---

very large,,,

cond(A)
ans = 6.1966e+16