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

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

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

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.

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

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