Issue with using Quadratic(p=2) elements for a 1D advection diffusion problem

mesh
1

Hi,
I am trying to solve a model advection-diffusion problem in Chapter 2 of the book " Finite Element Methods for Flow Problems" by Donea and Huerta.
the system is

\begin {equation} - u_x - \nu u_{xx} = 1\, \forall\, x\in [0,1] \end {equation}

with the boundary conditons
u(0) = 0
u(1) = 0
I am using \nu=0.02 and 5 elements for my solution
I have been able to match my solutions when solving the problem with linear elements(p=1). I am unable to get any sensible result when I try to use quadratic (p=2) elements

I basically change the element order by making the function space order from

   V = FunctionSpace ( my_mesh, 'CG', 1)

to

  V = FunctionSpace ( my_mesh, 'CG', 2)

to get the quadatic mesh

However, when I solve the system and plot the results, I am getting a piecewise linear plot
image
the actual figure from the source is ( see Pe=5) which is a squiggly solution.

image
image688×613 33.6 KB

I also tried plotting the nodes, elements, and coordinates for the mesh, but I was surprised to find that they were still 2 noded elements used for the calculation even when the degree was 2 in the FunctionSpace. Since I am using 5 elements, the quadratic 1D mesh should have 11 elements, but I have only 6. (See output after the code)

Could someone guide me on how to get the correct solution?

My code

from __future__ import print_function
from fenics import *
from ufl import nabla_div
import pylab as plt
import numpy as np

## convection_diffusion simulates a 1D convection diffusion problem in Fig 2.1
#    - a*ux - (nu ux),x = f in Omega = the unit interval.
#      u(0) = 0
#      u(1) = 0
#      f=1
my_nu=0.02
n_ele=5
#Define the mesh
my_mesh = UnitIntervalMesh ( n_ele )

#  Set the function space.
V = FunctionSpace ( my_mesh, 'CG', 2)
#  Set the trial and test functions.
#
u = TrialFunction ( V )
v = TestFunction ( V )

#  Set the nu value
mu = Constant ( my_nu )
    
s = Constant ( 1.0 )
#
#  Define the bilinear form and right hand side
#
aa= Constant ( 1.0 )
a = ( aa* v* u.dx(0) + mu * u.dx(0) * v.dx(0) ) * dx
L = s * v * dx


#  Define the boundary condition.

def boundary ( x ):
   value = x[0] < DOLFIN_EPS or 1.0 - DOLFIN_EPS < x[0]
   return value

bc = DirichletBC ( V, 0, boundary )
#  Solve- use default guess
uh = Function ( V )
#  Solve the system.
solve ( a == L, uh, bc )

#  Plot the solution.
h=1/n_ele
fig = plt.figure (dpi=110 )
ax = plt.subplot ( 111 )
plot (uh, label = 'Computed Pe=%.2f'% (1*h/2/my_nu  ))
ax.legend ( )
ax.grid ( True )
ax.set_ylim([0,1.5])
plt.title ( 'convection_diffusion solutions, # ele %d' % ( n_ele ) )
plt.show()

#  Write NODES, ELEMENTS and VALUES 
#
node_num = my_mesh.num_vertices ( )
node_xy = my_mesh.coordinates ( )

print ( "Nodes")
for row in node_xy:
      print(row)
    


print ( "Elements")

#  Print ELEMENTS array

element_num = my_mesh.num_cells ( )
element_node = my_mesh.cells ( )

for row in element_node:
           print(row)

print ( " Nodal Values")
#  Create a VALUES array, and write it to a file.
#
u_nodal_values = uh.vector ( )
u_array = u_nodal_values.get_local ( )
value_num = my_mesh.num_vertices ( )
node_value = np.zeros ( node_num )

for i in range ( 0, value_num ):
    print(u_array[i] )

Output from Code

Nodes
[ 0.]
[ 0.2]
[ 0.4]
[ 0.6]
[ 0.8]
[ 1.]
Elements
[0 1]
[1 2]
[2 3]
[3 4]
[4 5]
Nodal Values
0.0
0.49944096704
1.12402169232
0.510899864052
0.769494036811
0.374829298033

1 Like
2

Standard plots in FEniCS evaluate the function at the vertices of the mesh. It’s a simple visualisation only since it’s cheap. For more complicated visualisations you should investigate more sophsticated plotting libraries/packages like pyvista or paraview. In this 1D case it’s pretty easy to plot the high order bases though. See the last few lines of the attached code.

from __future__ import print_function
from fenics import *
from ufl import nabla_div
import pylab as plt
import numpy as np

## convection_diffusion simulates a 1D convection diffusion problem in Fig 2.1
#    - a*ux - (nu ux),x = f in Omega = the unit interval.
#      u(0) = 0
#      u(1) = 0
#      f=1
my_nu=0.02
n_ele=5
#Define the mesh
my_mesh = UnitIntervalMesh ( n_ele )

#  Set the function space.
V = FunctionSpace ( my_mesh, 'CG', 2)
#  Set the trial and test functions.
#
u = TrialFunction ( V )
v = TestFunction ( V )

#  Set the nu value
mu = Constant ( my_nu )
    
s = Constant ( 1.0 )
#
#  Define the bilinear form and right hand side
#
aa= Constant ( 1.0 )
a = ( aa* v* u.dx(0) + mu * u.dx(0) * v.dx(0) ) * dx
L = s * v * dx


#  Define the boundary condition.

def boundary ( x ):
   value = x[0] < DOLFIN_EPS or 1.0 - DOLFIN_EPS < x[0]
   return value

bc = DirichletBC ( V, 0, boundary )
#  Solve- use default guess
uh = Function ( V )
#  Solve the system.
solve ( a == L, uh, bc )

x_plt = np.linspace(0, 1, 100)
uh_plt = np.fromiter(map(uh, x_plt), dtype=np.double)

plt.plot(x_plt, uh_plt)
plt.show()

image

3 Likes
3

@nate , thank you! So essesntially FENICS is solving the system with the correct element solution, but the plot tools use the end nodes of the elements for a basic visualization if i am correct?

I had also tried exporting the results as a Paraview pvd file to check the results. The code I used was

vtkfile = File('Pe.pvd")
vtkfile << uh

However, this sitll gives the old incorrect output. Is there a way to get the results correctly?

image
image810×825 9.7 KB

4

You should use XDMFFile and write checkpoint, as for instance showed here: When using CG degree 1, the result corresponds to degree 2 - #10 by dokken

1 Like