Define the forward Euler

1

Hi everyone,
I have read the FENICS course at here and it really helps.

But I have a very easy question just about defining the forward Euler. I am solving the same heat equation at the course and get good result at the backward Euler and the mid-point Euler. I thinks things only goes differently at the stiffness matrix. So I only change something at the stiffness matrix.But the results goes wrong.

Just wonder is there anything wrong at defining the weak form(forward Euler method =2).

#include FEniCS headers
from __future__ import print_function
from fenics import *
from ufl import nabla_div, nabla_grad, div
import time

#Problem parameters variables
H = 1.00   #Plate Height
W = 1.00   #Plate Width
k_bar=385  #watt.m−1K−1
value = 3.8151e6 #N.m−2K−1
T=1000.0
num_steps=100  #time steps
dt = T / num_steps  #delta t
method=2    #method=1 backward method=2 forward method=3 mid-point

#Create mesh
#The points represent the diagonal ends of the box, followed by the number of elments along each dimension
mesh = RectangleMesh(Point(0, 0), Point(W, H), 20, 10) #2D mesh

#Define function space.
#This creates a H1 functiPon space from which we can construct the solution (trail) and test function spaces.
V = FunctionSpace(mesh, 'P', 1)  #With degree 1, we get the linear Lagrange element

# Define needed u (trial space), w (test space) function spaces
u = TrialFunction(V)
w = TestFunction(V)

# Define boundary condition
g_left= Constant(300.0)          #along x = 0 m (left edge)  K
g_right= Constant(310.0)         #along x = 1 m (right edge) K


# Define function to identify the left boundary at 300K
def on_left(x, on_boundary):
    return on_boundary and x[0] < DOLFIN_EPS
# Create the boundary condition
bc_left = DirichletBC(V, g_left, on_left) 

# Define function to identify the right boundary at 310K
def on_right(x, on_boundary):
    return on_boundary and x[0] > W - DOLFIN_EPS 
# Create the boundary condition
bc_right = DirichletBC(V, g_right, on_right) 


bc = [ bc_left, bc_right]

# Define expressions for heat conduction

def j(u):
    return -k_bar*nabla_grad(u) 

# Define initial value

u_0 = Expression('x[0] < 0.5  ? 300 : 300+20*(x[0]-0.5)', degree=1)
u_n=interpolate(u_0,V) 



vtkfile = File('heat/heat.pvd')

#Define the weak form of the probleminner
f= Constant(0.0)
t=0
if method==1:
      F =value*u*w*dx - dt*dot(j(u),nabla_grad(w))*dx - (value*u_n + dt*f)*w*dx    #define the weak form(backward) 
elif method==2:
      F =value*u*w*dx - dt*dot(j(u_n),nabla_grad(w))*dx - (value*u_n + dt*f)*w*dx    #define the weak form(forward)
elif method==3:
      F =value*u*w*dx - 0.5* dt*dot(j(u),nabla_grad(w))*dx - 0.5* dt*dot(j(u_n),nabla_grad(w))*dx- (value*u_n + dt*f)*w*dx    #define the weak form(mid-point)
a, L = lhs(F), rhs(F)

#solve the problem 
u = Function(V)

#vtkfile << (u_n, t)
for n in range(num_steps):

    # Update current time
    t += dt
  
    # Compute solution
    solve(a == L, u, bc)

    # Save to file and plot solution
    vtkfile << (u, t)

    # Update previous solution
    u_n.assign(u)
2

The forward Euler method is only conditionally stable. It requires a very small time step for parabolic PDEs like the heat equation. Look up the term “CFL condition”. For the heat equation, you’d need a time step that scales like some dimensionless constant times the square of the element size, divided by k_bar, which will be much smaller than what you’re using.

1 Like