Heat equation with explicit euler

General
1

Hello everyone,
I wanted to use an explicit Euler as time discretization for the heat equation from the tutorial. When i run the code every second time step is zero, when i review the xdmf.file. Has anyone an idea how to fix this?

import numpy as np

from mpi4py import MPI
from petsc4py import PETSc

from dolfinx import fem, mesh, io, plot

# Define temporal parameters
t = 0 # Start time
T = 1.0 # Final time
num_steps = 100     
dt = T / num_steps # time step size

# Define mesh
nx, ny = 5, 5
domain = mesh.create_rectangle(MPI.COMM_WORLD, [np.array([-2, -2]), np.array([2, 2])], 
                               [nx, ny], mesh.CellType.triangle)
V = fem.FunctionSpace(domain, ("CG", 1))

# Create initial condition
def initial_condition(x, a=5):
    return np.exp(-a*(x[0]**2+x[1]**2))
u_n = fem.Function(V)
u_n.name = "u_n"
u_n.interpolate(initial_condition)

# Create boundary condition
fdim = domain.topology.dim - 1
boundary_facets = mesh.locate_entities_boundary(
    domain, fdim, lambda x: np.full(x.shape[1], True, dtype=bool))
bc = fem.dirichletbc(PETSc.ScalarType(0), fem.locate_dofs_topological(V, fdim, boundary_facets), V)

xdmf = io.XDMFFile(domain.comm, "explicit.xdmf", "w")
xdmf.write_mesh(domain)

# Define solution variable, and interpolate initial solution for visualization in Paraview
uh = fem.Function(V)
uh.name = "uh"
uh.interpolate(initial_condition)
xdmf.write_function(uh, t)

import ufl
u, v = ufl.TrialFunction(V), ufl.TestFunction(V)
f = fem.Constant(domain, PETSc.ScalarType(0))
a = u * v * ufl.dx 
L = (dt*ufl.dot(ufl.grad(u_n), ufl.grad(v))-u_n*v-dt*f*v)*ufl.dx

bilinear_form = fem.form(a)
linear_form = fem.form(L)

A = fem.petsc.assemble_matrix(bilinear_form, bcs=[bc])
A.assemble()
b = fem.petsc.create_vector(linear_form)

solver = PETSc.KSP().create(domain.comm)
solver.setOperators(A)
solver.setType(PETSc.KSP.Type.PREONLY)
solver.getPC().setType(PETSc.PC.Type.LU)

for i in range(num_steps):
    t += dt

    # Update the right hand side reusing the initial vector
    with b.localForm() as loc_b:
        loc_b.set(0)
    fem.petsc.assemble_vector(b, linear_form)
    
    # Apply Dirichlet boundary condition to the vector
    fem.petsc.apply_lifting(b, [bilinear_form], [[bc]])
    b.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
    fem.petsc.set_bc(b, [bc])

    # Solve linear problem
    solver.solve(b, uh.vector)

    # Update solution at previous time step (u_n)
    u_n.x.array[:] = uh.x.array

    # Write solution to file
    xdmf.write_function(uh, t)
xdmf.close()
2

There is a mistake here, this should be
L = (-dt*ufl.dot(ufl.grad(u_n), ufl.grad(v))+u_n*v-dt*f*v)*ufl.dx

You need to make the signs consistent with a=L, not F=0

2 Likes
3

Hello Fenics users,

I want to solve the problem using a nonlinear solver to solve F = 0.

But I do not know which couple lines of code I should use to save my solution in xdmf file and then visualise it with paraview.

Here is my code:

#https://fenicsproject.discourse.group/t/heat-equation-with-explicit-euler/10987

import numpy as np

from mpi4py import MPI
from petsc4py import PETSc

from dolfinx import fem, mesh, io, plot
from dolfinx.io import XDMFFile, VTXWriter

from dolfinx.fem.petsc import assemble_vector, assemble_matrix, create_vector, apply_lifting, set_bc
from dolfinx.mesh import CellType, create_rectangle, locate_entities_boundary
from dolfinx.fem import (Constant, dirichletbc, Function, FunctionSpace, assemble_scalar,
                         form, locate_dofs_geometrical, locate_dofs_topological, set_bc)
from dolfinx.nls.petsc import NewtonSolver

from ufl import dx, grad, inner, TrialFunction, derivative, dot
from dolfinx.fem.petsc import NonlinearProblem, LinearProblem

# Define temporal parameters
t = 0 # Start time
T = 1.0 # Final time
num_steps = 100     
dt = T / num_steps # time step size

# Define mesh
nx, ny = 16, 16
domain = mesh.create_rectangle(MPI.COMM_WORLD, [np.array([0, 0]), np.array([1, 1])], 
                               [nx, ny], mesh.CellType.triangle)
V = fem.FunctionSpace(domain, ("Lagrange", 1))

# Create initial condition
# Interpolate initial condition


#def initial_condition(x, a=0):
#    return np.exp(-a*(x[0]**2+x[1]**2))
u_n = fem.Function(V)
u_n.name = "u_n"
#u_n.interpolate(initial_condition)
u_n.interpolate(lambda x: np.full(x.shape[1], PETSc.ScalarType(0.), dtype=np.float64)) 

#We create two python functions for determining the facets to apply boundary conditions to



def right(x):
    return np.isclose(x[0], 1)  #right (x = 1)


def top(x):
    return np.isclose(x[1], 1)  #top (y = 1)

# 
###################################################
######## Dirichlet boundary conditions #############
###################################################
# 

boundary_facets_right = locate_entities_boundary(domain, domain.topology.dim - 1, right)
boundary_dofs_right = locate_dofs_topological(V, domain.topology.dim - 1, boundary_facets_right)
bc_m_right = dirichletbc(PETSc.ScalarType(1.), boundary_dofs_right, V)

boundary_facets_top = locate_entities_boundary(domain, domain.topology.dim - 1, top)
boundary_dofs_top = locate_dofs_topological(V, domain.topology.dim - 1, boundary_facets_top)
bc_m_top = dirichletbc(PETSc.ScalarType(1.), boundary_dofs_top, V)

bcs = [bc_m_right, bc_m_top]  #dirichlet boundary conditions

#xdmf = io.XDMFFile(domain.comm, "explicit_30th_june_2024.xdmf", "w")
#xdmf.write_mesh(domain)

file = XDMFFile(domain.comm, "implicit_30th_june_2024_Nonlinear_solver.xdmf", "w")
file.write_mesh(domain)

# Define solution variable, and interpolate initial solution for visualization in Paraview
uh = fem.Function(V)
uh.name = "uh"
#uh.interpolate(initial_condition)
uh.interpolate(lambda x: np.full(x.shape[1], PETSc.ScalarType(0.), dtype=np.float64)) 
file.write_function(uh, t)

import ufl
u, v = fem.Function(V), ufl.TestFunction(V)
f = fem.Constant(domain, PETSc.ScalarType(0))
a = (u * v + dt*ufl.dot(ufl.grad(u), ufl.grad(v))) * ufl.dx 
L = (u_n*v + dt*f*v)*ufl.dx
F = a - L

# Non linear problem definition
du  = TrialFunction(V)
#Compute the jacobian
Jac = derivative(F, u, du)
problem = NonlinearProblem(F, u, bcs = bcs, J = Jac)

  
  
solver = NewtonSolver(MPI.COMM_WORLD, problem)
solver.convergence_criterion = "incremental"      #residual    #incremental
# Absolute tolerance
solver.atol = 5e-10
solver.rtol = 1e-12 #0.0 #1e-14 #0.0 #1e-12
#solver.max_it = 70 
solver.report = True
#solver.max_it = 100        #To increase the number of iterations
solver.error_on_nonconvergence = False
solver.max_it = 3

info = solver.solve(u)
print(info)

Thank you.

4

Please look at the tutorial
https://jsdokken.com/dolfinx-tutorial/chapter2/diffusion_code.html#time-dependent-output

5

Thank you @dokken. This is my benchmark tutorial. It uses linear solver to solve a=L. But I am struggling to modify

in my own that’s why I use nonlinear solver to solve F=0.
But I do not know how to save and visualise my results with Paraview.

6

This is my attempt:

for i in range(num_steps):
    t += dt

    # Update the right hand side reusing the initial vector
    #with b.localForm() as loc_b:
        #loc_b.set(0)
    #fem.petsc.assemble_vector(b, linear_form)
    
    # Apply Dirichlet boundary condition to the vector
    #fem.petsc.apply_lifting(b, [bilinear_form], [[bc]])
    #b.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
    #fem.petsc.set_bc(b, [bc])

    # Solve nonlinear problem
    #solver.solve(b, uh.vector)
    solver.solve(u)

    # Update solution at previous time step (u_n)
    u_n.x.array[:] = uh.x.array
    u_n.x.scatter_forward()

But the solution is the same at every second time step.

7

for i in range(num_steps):
t += dt

# Update the right hand side reusing the initial vector
#with b.localForm() as loc_b:
    #loc_b.set(0)
#fem.petsc.assemble_vector(b, linear_form)

# Apply Dirichlet boundary condition to the vector
#fem.petsc.apply_lifting(b, [bilinear_form], [[bc]])
#b.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
#fem.petsc.set_bc(b, [bc])

# Solve nonlinear problem
#solver.solve(b, uh.vector)
solver.solve(u)

# Update solution at previous time step (u_n)
u_n.x.array[:] = uh.x.array
u_n.x.scatter_forward()
#u = uh

# Write solution to file
file.write_function(uh, t)

file.close()

8

So your issue is not storing the solution, But updating the solution over time?

As you haven’t provided a code that reproduces the behavior it is hard to give any specific pointers.

9

Yes @dokken , this is my issue.

Here is my code for approximating the 2D heat equation using FEM. The initial condition u(t=0, x, y) = 0 and boundaries conditions (top and right) are equal to 1.

Here is my code :

#https://fenicsproject.discourse.group/t/heat-equation-with-explicit-euler/10987

import numpy as np

from mpi4py import MPI
from petsc4py import PETSc

from dolfinx import fem, mesh, io, plot
from dolfinx.io import XDMFFile, VTXWriter

from dolfinx.fem.petsc import assemble_vector, assemble_matrix, create_vector, apply_lifting, set_bc
from dolfinx.mesh import CellType, create_rectangle, locate_entities_boundary
from dolfinx.fem import (Constant, dirichletbc, Function, FunctionSpace, assemble_scalar,
                         form, locate_dofs_geometrical, locate_dofs_topological, set_bc)
from dolfinx.nls.petsc import NewtonSolver

from ufl import dx, grad, inner, TrialFunction, derivative, dot
from dolfinx.fem.petsc import NonlinearProblem, LinearProblem

# Define temporal parameters
t = 0 # Start time
T = 1.0 # Final time
num_steps = 100     
dt = T / num_steps # time step size

# Define mesh
nx, ny = 16, 16
domain = mesh.create_rectangle(MPI.COMM_WORLD, [np.array([0, 0]), np.array([1, 1])], 
                               [nx, ny], mesh.CellType.triangle)
V = fem.FunctionSpace(domain, ("Lagrange", 1))

# Create initial condition
# Interpolate initial condition


#def initial_condition(x, a=0):
#    return np.exp(-a*(x[0]**2+x[1]**2))
u_n = fem.Function(V)
u_n.name = "u_n"
#u_n.interpolate(initial_condition)
u_n.interpolate(lambda x: np.full(x.shape[1], PETSc.ScalarType(0.), dtype=np.float64)) 

#We create two python functions for determining the facets to apply boundary conditions to



def right(x):
    return np.isclose(x[0], 1)  #right (x = 1)


def top(x):
    return np.isclose(x[1], 1)  #top (y = 1)

# 
###################################################
######## Dirichlet boundary conditions #############
###################################################
# 

boundary_facets_right = locate_entities_boundary(domain, domain.topology.dim - 1, right)
boundary_dofs_right = locate_dofs_topological(V, domain.topology.dim - 1, boundary_facets_right)
bc_m_right = dirichletbc(PETSc.ScalarType(1.), boundary_dofs_right, V)

boundary_facets_top = locate_entities_boundary(domain, domain.topology.dim - 1, top)
boundary_dofs_top = locate_dofs_topological(V, domain.topology.dim - 1, boundary_facets_top)
bc_m_top = dirichletbc(PETSc.ScalarType(1.), boundary_dofs_top, V)

bcs = [bc_m_right, bc_m_top]  #dirichlet boundary conditions

#xdmf = io.XDMFFile(domain.comm, "explicit_30th_june_2024.xdmf", "w")
#xdmf.write_mesh(domain)

file = XDMFFile(domain.comm, "implicit_30th_june_2024_Nonlinear_solver.xdmf", "w")
file.write_mesh(domain)

# Define solution variable, and interpolate initial solution for visualization in Paraview
uh = fem.Function(V)
uh.name = "uh"
#uh.interpolate(initial_condition)
uh.interpolate(lambda x: np.full(x.shape[1], PETSc.ScalarType(0.), dtype=np.float64)) 
file.write_function(uh, t)

import ufl
u, v = fem.Function(V), ufl.TestFunction(V)
f = fem.Constant(domain, PETSc.ScalarType(0))
a = (u * v + dt*ufl.dot(ufl.grad(u), ufl.grad(v))) * ufl.dx 
L = (u_n*v + dt*f*v)*ufl.dx
F = a - L

# Non linear problem definition
du  = TrialFunction(V)
#Compute the jacobian
Jac = derivative(F, u, du)
problem = NonlinearProblem(F, u, bcs = bcs, J = Jac)


#Solve variational problem
# Create nonlinear problem and Newton solver

  
solver = NewtonSolver(MPI.COMM_WORLD, problem)
solver.convergence_criterion = "incremental"      #residual    #incremental
# Absolute tolerance
solver.atol = 5e-10
solver.rtol = 1e-12 #0.0 #1e-14 #0.0 #1e-12

solver.report = True

solver.error_on_nonconvergence = False
solver.max_it = 3

info = solver.solve(u)
print(info)
10

But in this code you have not attempted to write the solution to file. Could you please provide a reproducible example were you store multiple time steps.
In the code above you only call:

once, for t=0

11

Yes, here is:

for i in range(num_steps):
    t += dt

    # Update the right hand side reusing the initial vector
    #with b.localForm() as loc_b:
        #loc_b.set(0)
    #fem.petsc.assemble_vector(b, linear_form)
    
    # Apply Dirichlet boundary condition to the vector
    #fem.petsc.apply_lifting(b, [bilinear_form], [[bc]])
    #b.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
    #fem.petsc.set_bc(b, [bc])

    # Solve nonlinear problem
    #solver.solve(b, uh.vector)
    solver.solve(u)

    # Update solution at previous time step (u_n)
    u_n.x.array[:] = uh.x.array
    u_n.x.scatter_forward()
    #u = uh

    # Write solution to file
    file.write_function(uh, t)
file.close()