MMS for poisson equation with periodic boundary conditions


80
views
0
6 weeks ago by
I am trying to solve the poisson equation with periodic boundary conditions on the top and bottom boundaries of a unit square. When I try to use the method of manufactured solutions (MMS) my solution, u, varies significantly from the prescribed solution, u_D. My solution does appear to be periodic.

from dolfin import *
%matplotlib inline

# Define mesh
mesh = UnitSquareMesh(32, 32)

# Define classes for Dirichlet section (left and right sides) of the boundary
class DirichletBoundary(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], 0.0) or near(x[0], 1.0) and on_boundary


#Create and map bottom boundary to top boundary for periodic conditions
class PeriodicBoundary(SubDomain):
def inside(self, x, on_boundary):
return near(x[1], 0.0) and on_boundary # Creates bottom boundary

# Maps bottom to top boundary
def map(self, x, y):
y[0] = x[0]
y[1] = x[1] - 1

# Define input data
u_D = Expression("sin(x[0]*x[1]*DOLFIN_PI)", degree=2)
f = Expression("x[1]*x[1]*DOLFIN_PI*DOLFIN_PI*sin(x[0]*x[1]*DOLFIN_PI) \
+ x[0]*x[0]*DOLFIN_PI*DOLFIN_PI*sin(x[0]*x[1]*DOLFIN_PI)", degree=2)

# Define function space and basis functions
V = FunctionSpace(mesh, "CG", 2, constrained_domain=PeriodicBoundary())
u = TrialFunction(V)
v = TestFunction(V)

# Define boundary conditions at Dirichlet boundaries
bcs = DirichletBC(V, u_D, DirichletBoundary())

# Define variational form
F = (inner(grad(u), grad(v))*dx - f*v*dx)

# Separate left and right hand sides of equation
a, L = lhs(F), rhs(F)

# Solve problem
u = Function(V)
solve(a == L, u, bcs)

plot(u)
plot(mesh)

# Compute error in L2 norm
error_L2 = errornorm(u_D, u, 'L2')

# Compute maximum error at vertices
vertex_values_u_D = u_D.compute_vertex_values(mesh)
vertex_values_u = u.compute_vertex_values(mesh)
import numpy as np
error_max = np.max(np.abs(vertex_values_u_D - vertex_values_u))

# Print errors
print('error_L2 =', error_L2)
print('error_max =', error_max)


I am new to FEniCS, python, and solving PDEs so any insights as to why I have such high error would be very much appreciated.

Thanks!
Community: FEniCS Project

1 Answer


3
6 weeks ago by
The difficulty is that the choice of \(u_D = \operatorname{sin}\left(\pi x_0x_1\right)\) does not have the desired periodicity.  (Consider fixing the value of \(x_0\) at, say \(x_0=0.5\), and looking at the behavior of the resulting function of \(x_1\).)  If you choose a \(u_D\) with periodicity in \(x_1\) for every fixed \(x_0\), you will see convergence, e.g.,

from dolfin import *

# Define mesh
NEL = 32
mesh = UnitSquareMesh(NEL,NEL)

# Define classes for Dirichlet section (left and right sides) of the boundary
class DirichletBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], 0.0) or near(x[0], 1.0) and on_boundary


#Create and map bottom boundary to top boundary for periodic conditions
class PeriodicBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[1], 0.0) and on_boundary # Creates bottom boundary

    # Maps bottom to top boundary
    def map(self, x, y):
        y[0] = x[0]
        y[1] = x[1] - 1
    
# Define input data
#u_D = Expression("sin(x[0]*x[1]*DOLFIN_PI)", degree=2)
#f = Expression("x[1]*x[1]*DOLFIN_PI*DOLFIN_PI*sin(x[0]*x[1]*DOLFIN_PI) + x[0]*x[0]*DOLFIN_PI*DOLFIN_PI*sin(x[0]*x[1]*DOLFIN_PI)", degree=2)
x = SpatialCoordinate(mesh)
u_D = cos(x[0]*2.0*pi)*sin(x[1]*2.0*pi)
f = -div(grad(u_D))

# Define function space and basis functions
V = FunctionSpace(mesh, "CG", 2, constrained_domain=PeriodicBoundary())
u = TrialFunction(V)
v = TestFunction(V)

# Define boundary conditions at Dirichlet boundaries
bcs = DirichletBC(V, u_D, DirichletBoundary())

# Define variational form
F = (inner(grad(u), grad(v))*dx - f*v*dx)

# Separate left and right hand sides of equation
a, L = lhs(F), rhs(F)

# Solve problem
u = Function(V)
solve(a == L, u, bcs)

import matplotlib.pyplot as plt
plot(u)
plot(mesh)
plt.show()

# Compute error in L2 norm
#error_L2 = errornorm(u_D, u, 'L2')
import math
error_L2 = math.sqrt(assemble(((u-u_D)**2)*dx))

# Compute maximum error at vertices
#vertex_values_u_D = u_D.compute_vertex_values(mesh)
#vertex_values_u = u.compute_vertex_values(mesh)
#import numpy as np
#error_max = np.max(np.abs(vertex_values_u_D - vertex_values_u))

# Print errors
print('error_L2 =', error_L2)
#print('error_max =', error_max)
​
Thanks for the reply! I understand now why my uD was affecting convergence. Also thanks for the alternative method for entering expressions and calculating the L2 norm.
written 5 weeks ago by Michael Mazza  
Please login to add an answer/comment or follow this question.

Similar posts:
Search »