Periodic Boundary Conditions for Poisson Problem

Question

Hello everybody!

I'm solving the cahn-hilliard-equation with periodic boundary conditions by a mixed finite element method. Here I have no problems with the periodic boundary conditions, everything works fine! Here is the code, just for comparison...

from dolfin import *
import numpy as np
import random


# Dolfin Parameters
parameters["form_compiler"]["optimize"]     = True
parameters["form_compiler"]["cpp_optimize"] = True
parameters["form_compiler"]["representation"] = "quadrature"


# Periodic Boundary Class
class PeriodicBC(SubDomain):
    def __init__(self, tolerance=DOLFIN_EPS, length = 1., length_scaling = 1.):
        SubDomain.__init__(self)
        self.tol = tolerance
        self.length = length
        self.length_scaling = length_scaling
    def inside(self, x, on_boundary):
        return (bool(near(x[1], 0.) and on_boundary) or bool(near(x[0], 0.) and on_boundary))
    def map(self, x, y):
        if x[1] > (self.length/self.length_scaling - self.tol):
            y[0] = x[0]
            y[1] = x[1] - self.length/self.length_scaling
        else:
            y[0] = x[0] - self.length/self.length_scaling
            y[1] = x[1]


# Initial Conditions
class InitialConditions(Expression):
    def __init__(self, MeanConcentration, PerturbationConcentration):
        random.seed(5)
        self.Concentration = MeanConcentration
        self.Perturbation = PerturbationConcentration
    def eval(self, values, x):
        values[0] = self.Concentration + self.Perturbation*(0.5 - random.random())
        values[1] = 0.
    def value_shape(self):
        return (2,)


# Parameters
Length = 8.0e-05
LengthScaling = 1.0e-04

TimeStepSize = 5.0e-02
TimeScaling = LengthScaling**2/1.0e-12
ScaledTimeStepSize = TimeStepSize/TimeScaling
ScaledFreeEnergy = 4.0
ScaledLineTension = 1.0e-12/LengthScaling**2


# Create Mesh
mesh = RectangleMesh(0.0, 0.0, Length/LengthScaling, Length/LengthScaling, \
             64, 64, "crossed")

# Build Function Space and Functions
PBC = PeriodicBC(length = Length, length_scaling = LengthScaling)
V  = FunctionSpace(mesh, "CG", 1, constrained_domain=PBC)             
MS = MixedFunctionSpace([V,V])

# Define Trial and Test Functions
du      = TrialFunction(MS)
wc, ws = TestFunctions(MS)

# Define Functions
u   = Function(MS)  # Current Solution
u0  = Function(MS)  # Old Solution

# Split mixed Functions
dc, ds = split(du)
c,  s  = split(u)
c0, s0 = split(u0)


# Define the Chemical Potential and Mobility
mu = 2.0*c*(1.0-c)*(1.0-c) - 2.0*c*c*(1.0-c)
M = 1.

# Weak Statement
F_1 = (c - c0)*wc*dx
F_2 = inner(M*grad(s), grad(wc))*dx
F_3 = s*ws*dx
F_4 = mu*ws*dx
F_5 = inner(grad(c), grad(ws))*dx

# Function F for Newton Method
F = F_1 + ScaledTimeStepSize*F_2 + (F_3 - ScaledFreeEnergy*F_4 - ScaledLineTension*F_5)       

# Derivative of F for Newton Method
J = derivative(F, u, du)

# Create Nonlinear Problem and Newton Solver
bcs =[]
Problem = NonlinearVariationalProblem(F, u, bcs, J)
Solver  = NonlinearVariationalSolver(Problem)


# Initial Conditions
TimeStep = 0
Time = 0.0
u.interpolate(InitialConditions(0.3, 0.01))
# Loop
while (TimeStep < 5):
    # Iterate
    TimeStep += 1
    Time += TimeStepSize
    # Store old Concentration
    u0.vector()[:] = u.vector()
    # Solve Nonlinear Equation
    Solver.solve()
    # Save Solution
    c, s = u.split(deepcopy = True) 
    plot(c)

Here is my problem: I want to solve the poisson problem on the same mesh and with the same boundary conditions, but in the following code the periodic boundary conditions do not work and I don't know why, so I need some help.

from dolfin import *
import numpy as np
import random

# Dolfin Parameters
parameters["form_compiler"]["optimize"]     = True
parameters["form_compiler"]["cpp_optimize"] = True
parameters["form_compiler"]["representation"] = "quadrature"


# Periodic Boundary Class
class PeriodicBC(SubDomain):
    def __init__(self, tolerance=DOLFIN_EPS, length = 1., length_scaling = 1.):
        SubDomain.__init__(self)
        self.tol = tolerance
        self.length = length
        self.length_scaling = length_scaling
    def inside(self, x, on_boundary):
        return (bool(near(x[1], 0.) and on_boundary) or bool(near(x[0], 0.) and on_boundary))
    def map(self, x, y):
        if x[1] > (self.length/self.length_scaling - self.tol):
            y[0] = x[0]
            y[1] = x[1] - self.length/self.length_scaling
        else:
            y[0] = x[0] - self.length/self.length_scaling
            y[1] = x[1]


# Initial Condition
class InitialConditions(Expression):
    def __init__(self, center=0.1, radius=0.08):
        self.center = center
        self.radius = radius
        return None
    def eval(self, values, x):
        if ((x[0]-self.center)**2)/((1.*self.radius)**2)+((x[1]-self.center)**2)/((self.radius)**2) < 1.0:
            values[0] = 1.
        else:
            values[0] = 0.                       

# Translation 
def translation(h):
        int_h = h *dx
        average_h = assemble(int_h)
        h_array = h.vector().array() - average_h
        return h_array

# Parameters
Length = 8.0e-05
LengthScaling = 1.0e-04


# Create Mesh
mesh = RectangleMesh(0.0, 0.0, Length/LengthScaling, Length/LengthScaling, 64, 64, "crossed")

# Build Function Space and Functions
PBC = PeriodicBC(length = Length, length_scaling = LengthScaling)
W = FunctionSpace(mesh, "CG", 1,  constrained_domain=PBC)
phi = Function(W)
c = Function(W)
v_phi = TrialFunction(W)
w_phi = TestFunction(W)
# Auxiliary Function 
h1 = Function(W)    

# Define Linear Variational Problem
a = (1./LengthScaling**2)*inner(grad(v_phi), grad(w_phi))*dx
L = h1*w_phi*dx

# Create Linear Problem and Solver
bcs = []
Problem = LinearVariationalProblem(a, L, phi, bcs) 
Solver = LinearVariationalSolver(Problem)

# Prepare RHS
c.interpolate(InitialConditions(0.1, 0.08))
h1.vector()[:] = translation(c)

# Solve Linear Equation, Plot Solution
Solver.solve()
plot(phi)
interactive()

I would be very grateful if anyone could help!

Comments

lisa, can you made the error-producing example shorter (but still complete)?

Answer 1

Hi,

The problem seems to be your implementation of the 2 periodic directions. See this link for help. A map that should work in your case is:

# Periodic Boundary Class
class PeriodicBC(SubDomain):
    def __init__(self, tolerance=DOLFIN_EPS, length = 1., length_scaling = 1.):
        SubDomain.__init__(self)
        self.tol = tolerance
        self.length = length
        self.length_scaling = length_scaling

# Left boundary is "target domain" G
def inside(self, x, on_boundary):
    # return True if on left or bottom boundary AND NOT on one of the two corners (0, L) and (L, 0)
    return bool((near(x[0], 0) or near(x[1], 0)) and 
            (not ((near(x[0], 0) and near(x[1], self.length/self.length_scaling)) or 
                    (near(x[0], self.length/self.length_scaling) and near(x[1], 0)))) and on_boundary)

def map(self, x, y):
    L = self.length/self.length_scaling
    if near(x[0], L) and near(x[1], L):
        y[0] = x[0] - L
        y[1] = x[1] - L
    elif near(x[0], L):
        y[0] = x[0] - L
        y[1] = x[1]
    else:   # near(x[1], L)
        y[0] = x[0]
        y[1] = x[1] - L

Jan is correct, the problem is singular because you have a Poisson problem with no boundaries. So if you take your solution phi, then phi +10 will also be a solution. It only looks like you have no problems because your LUsolver is giving you a nice solution that almost averages to 0, probably since you start with initial guess = 0. Try adding this to your code just before you solve it:

phi.vector()[:] = 10.  # initial guess = 10
Solver.parameters['linear_solver'] = 'gmres'
Solver.parameters['preconditioner'] = 'petsc_amg'
Solver.parameters['krylov_solver']['monitor_convergence'] = True
Solver.parameters['krylov_solver']['nonzero_initial_guess'] = True

You will now find the phi + 10 solution. A much used solution to this "problem" is to normalize your solution after solving. Add this after solving and your vector will sum to zero:

normalize(phi.vector())

You should not mix Periodic boundary conditions with Dirichlet conditions. Use normalize, a real space (neuman-poisson demo) or a nullspace (singular demo).

Comments

No, a transient one is working Cahn-Hilliard example. Problematic Poisson example is singular without any Dirichlet condition or Lagrange multiplier - solution is specified up to the constant term.

---

Ah, I was looking at the Cahn-Hilliard stuff. Much too much code, and the CH was not really working...

---

Thank you for your helpful response! You were right, there was a problem with the corners.

---

Thank you so much, Jan and Mikael! Yes, the problem is singular.
By choosing "normalize(phi.vector())" I avoid the undesirable effects at the corners.
Thanks a lot!

Answer 2

Your problem is probably singular. Adding say

bcs = [DirichletBC(W, Expression("k*x[0]*x[0]", k=1e-10), PBC)]

helps. If you want Neumann condition on a whole boundary fix solution in one point

bcs = [DirichletBC(W, 0.0, "x[0]<DOLFIN_EPS & x[1]<DOLFIN_EPS", "pointwise")]

or use Lagrange multiplier approach (see neumann-poisson demo).

Comments

Hello Jan, thank you very much for your helpful response! I add

bcs = [DirichletBC(W, 0.0, "x[0]<DOLFIN_EPS & x[1]<DOLFIN_EPS", "pointwise")]

and it works! But now the solution drops considerably at the corners of the square and thats an effect that I want to avoid. There is something I don't understand: I just wonder why the additional condition is not necessary when solving the poisson problem on the unit square. Here is the code (I'm really sorry, but I really do not know how to reduce it):

from dolfin import *
import numpy as np
import random


# Periodic Boundary Class
class PeriodicBC(SubDomain):
    def __init__(self, tolerance=DOLFIN_EPS, length = 1., length_scaling = 1.):
        SubDomain.__init__(self)
        self.tol = tolerance
        self.length = length
        self.length_scaling = length_scaling
    def inside(self, x, on_boundary):
    # return True if on left or bottom boundary AND NOT on one of the two corners (0, L) and (L, 0)
        return bool((near(x[0], 0) or near(x[1], 0) and on_boundary) and \
                (not ((near(x[0], 0) and near(x[1], self.length/self.length_scaling)) or \
                (near(x[0], self.length/self.length_scaling) and near(x[1], 0)))) and on_boundary)
    def map(self, x, y):
        L = self.length/self.length_scaling
        if near(x[0], L) and near(x[1], L):
            y[0] = x[0] - L
            y[1] = x[1] - L
        elif near(x[0], L):
            y[0] = x[0] - L
            y[1] = x[1]
        else:   # near(x[1], L)
            y[0] = x[0]
            y[1] = x[1] - L

# RHS
class RightHandSide(Expression):
    def __init__(self, center=0.1, radius=0.08):
        self.center = center
        self.radius = radius
        return None
    def eval(self, values, x):
        if ((x[0]-self.center)**2)/((1.*self.radius)**2)+((x[1]-self.center)**2)/((self.radius)**2) < 1.0:
            values[0] = 1.
        else:
            values[0] = 0.                       

# Translation 
def translation(h):
        int_h = h *dx
        average_h = assemble(int_h)
        h_array = h.vector().array() - average_h
        return h_array


# Create Mesh
mesh = RectangleMesh(0.0, 0.0, 1., 1., 64, 64, "crossed")

# Build Function Space and Functions
PBC = PeriodicBC()
W = FunctionSpace(mesh, "CG", 1,  constrained_domain=PBC)
phi = Function(W)
c = Function(W)
v_phi = TrialFunction(W)
w_phi = TestFunction(W)
# Auxiliary Function 
h1 = Function(W)    

# Define Linear Variational Problem
a = inner(grad(v_phi), grad(w_phi))*dx
L = h1*w_phi*dx

# Create Linear Problem and Solver
bcs = []
Problem = LinearVariationalProblem(a, L, phi, bcs) 
Solver = LinearVariationalSolver(Problem)

# Prepare RHS
c.interpolate(RightHandSide())
h1.vector()[:] = translation(c)

# Solve Equation, Plot Solution
Solver.solve()
plot(phi)
interactive()

Could you explain, why the additional condition is not necessary in that case?

---

I'm not sure. It seems to me that the solution is determined only up to the constant term. This would mean that the problem is singular.

---

See updated answer.