### MMS for poisson equation with periodic boundary conditions

80

views

0

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!

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

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)
```

Please login to add an answer/comment or follow this question.

_{D}was affecting convergence. Also thanks for the alternative method for entering expressions and calculating the L^{2}norm.