Periodic boundary condition not giving periodic solution

General
1

Hello All,
I am trying to solve the following non linear equation with Periodic boundary condition.

u" = u + 25/(1+exp(u)) on x ∈ [0, 2*pi].

I want the solution to be periodic with periodicity pi. Here is the code for the same

from fenics import *
import matplotlib.pyplot as plt
import numpy as np

L = 2*np.pi
c0 = 1
Pe = 25
D = 1
gamma = 1

l = L/(2*np.pi)
n = (l**2)*gamma
U = 25/l
stress_const = U*np.sqrt(gamma*n)

mesh_size = 1000
mesh = IntervalMesh(mesh_size, 0, L)
plot(mesh)

class PeriodicBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return bool(x[0] > -DOLFIN_EPS and x[0] <= L+DOLFIN_EPS)
    def map(self, x, y):
        y[0] = x[0] - L/2

pbc = PeriodicBoundary()


V = FunctionSpace(mesh, "CG", 1, constrained_domain=pbc)

u = project(Constant(2), V)
v = TestFunction(V)

F = dot(grad(u), grad(v))*dx - (1/l**2)*u*v*dx - (1/(D*n))*stress_const*(1 - ((exp(u)/(1+exp(u)))))*v*dx 
J = derivative(F,u)

problem = NonlinearVariationalProblem(F, u, J=J)

solver = NonlinearVariationalSolver(problem)
stype = 'newton'
solver.parameters['nonlinear_solver']=stype
sprms = solver.parameters[stype+'_solver']

sprms['maximum_iterations'] = 10000


solver.solve()

Please help me figure out the issue

2

You should just solve your problem in the domain [0,\pi] if your problem is periodic, i.e.

from fenics import *
import numpy as np

L = 2*np.pi
c0 = 1
Pe = 25
D = 1
gamma = 1
l = L/(2*np.pi)
n = (l**2)*gamma
U = 25/l
stress_const = U*np.sqrt(gamma*n)

mesh_size = 1000
mesh = IntervalMesh(mesh_size, 0, np.pi)

class PeriodicBoundary(SubDomain):
    def inside(self, x, on_boundary):
        """Left boundary"""
        return near(x[0], 0)
    def map(self, x, y):
        """ map from right to left boundary"""
        y[0] = x[0] - np.pi

pbc = PeriodicBoundary()


V = FunctionSpace(mesh, "CG", 1, constrained_domain=pbc)

u = project(Constant(2), V)
v = TestFunction(V)

F = dot(grad(u), grad(v))*dx - (1/l**2)*u*v*dx - (1/(D*n))*stress_const*(1 - ((exp(u)/(1+exp(u)))))*v*dx 
J = derivative(F,u)

problem = NonlinearVariationalProblem(F, u, J=J)

solver = NonlinearVariationalSolver(problem)
stype = 'newton'
solver.parameters['nonlinear_solver']=stype
sprms = solver.parameters[stype+'_solver']

sprms['maximum_iterations'] = 10000


solver.solve()
File("u.pvd") << u

You should just solve on the half interval, as dolfin relies on the mesh spacing macthing your mapping perfectly, which is really hard when you use numbers such as pi to describe its length.

1 Like
3

I don’t think the domain necessarily needs to be [0,\pi] explicitly. I think the problem is more that

y[0] = x[0] - L/2

was being used instead of

y[0] = x[0] - L

when the domain was defined as [0,L]. Perhaps the domain should be [0,L/2].
(of course, L/2=\pi so it’s an equivalent resolution)

1 Like
4

Thank you for your valuable input @dokken and @dparsons. I just have one question.

In the code you have mapped left boundary to the right boundary. Shouldn’t we map the a period of the function to the domain?

I mean something like this

class PeriodicBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return bool(x[0] > -DOLFIN_EPS and x[0] <= L+DOLFIN_EPS)
    def map(self, x, y):
        y[0] = x[0] - L/2

Where L is the domain.

5

I dont understand why you would do that. You are saying that your solution is periodic, that means that any point x\in[0,\pi] maps to y=x+pi\in[\pi,2\pi]. This means that you only need to solve an equation on [0,\pi] but that u(0)=u(\pi)

1 Like
6

Thanks a lot for the clarification. I misunderstood the implementation of PBC in fenics which led to this confusion.