Periodic boundary condition with DG and MPC

General
1

Hi everyone,

I recently installed dolfinx_mpc with the hope to be able to impose periodic boundary conditions together with the DG discretization to solve an eigenvalue problem.
I had a look at the different demos (in particular demo_periodic_gep.py) and none of them concerns the DG discretization. Is is actually possible to impose periodic boundary condition when the space of interest is FunctionSpace(msh, ufl.VectorElement("DG", msh.ufl_cell(), 2, 3))? I did several tries with both create_periodic_constraint_geometrical and create_periodic_constraint_topological but the computations did not run successfully (i.e. I did not obtain the expected eigenvalues).

So my questions are:

  1. it is possible to use dolfinx_mpc with the DG discretization (e.g. with a space FunctionSpace(msh, ufl.VectorElement("DG", msh.ufl_cell(), 2, 3))) to solve an eigenproblem?
  2. if so, is there somewhere such an example?

Best,

Lucas

2

Could you please provide a minimal reproducible example? As it is hard to give any guidance based of just your text.

One thing to note is that you would have to use create_periodic_constraint_geometrical, and not topological as all dofs of DG spaces are associated with the cell, and not facets.

3

So at the moment I don’t have such a MWE (but I would like to ^^!).

However I modified the script of the post https://fenicsproject.discourse.group/t/discontinuous-1d-eigenvector-using-dg-elements/10491 to use the function assemble_matrix provided by dolfinx_mpc together with an empty MultiPointConstraint instance:

# From https://fenicsproject.discourse.group/t/discontinuous-1d-eigenvector-using-dg-elements/10491
from dolfinx.fem import FunctionSpace, form, Function,dirichletbc,locate_dofs_topological
from dolfinx.mesh import meshtags,locate_entities, create_interval
from ufl import TestFunction,TrialFunction,Measure,inner,grad,jump,avg,Circumradius, FacetNormal
from dolfinx_mpc import MultiPointConstraint, assemble_matrix
import numpy as np
from slepc4py import SLEPc
from mpi4py import MPI

def OneDimensionalSetup(n_elem, x_MP = 0.2, r_outer = 0.25):
    # 1 ------------------------------------ 2
    #                 x=r_outer
    """ This function builds one dimensional setup.
        For boundaries, Tag 1 specifies left end and Tag 2 specifies right end. 
    Returns:
        mesh, subdomains, facet_tags
    """
    mesh = create_interval(MPI.COMM_WORLD, n_elem, [0.0, r_outer])
    def subdomain_func(x, x_MP=x_MP, eps=1e-16):
        x = x[0]
        return np.logical_and(x_MP - eps <= x, x <= x_MP + eps)
    tdim = mesh.topology.dim
    marked_cells = locate_entities(mesh, tdim, subdomain_func)
    fl = 0
    subdomains = meshtags(mesh, tdim, marked_cells, np.full(len(marked_cells), fl, dtype=np.int32))
    boundaries = [  (0, lambda x: np.isclose(x[0], x_MP)),
                    (1, lambda x: np.isclose(x[0], 0)),
                    (2, lambda x: np.isclose(x[0], r_outer))]
    facet_indices, facet_markers = [], []
    fdim = mesh.topology.dim - 1
    for (marker, locator) in boundaries:
        facets = locate_entities(mesh, fdim, locator)
        facet_indices.append(facets)
        facet_markers.append(np.full(len(facets), marker))
    facet_indices = np.array(np.hstack(facet_indices), dtype=np.int32)
    facet_markers = np.array(np.hstack(facet_markers), dtype=np.int32)
    sorted_facets = np.argsort(facet_indices)
    facet_tags = meshtags(mesh, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets])
    return mesh, subdomains, facet_tags

n_elem = 50
mesh, subdomains, facet_tags = OneDimensionalSetup(n_elem)
dx = Measure("dx", domain=mesh, subdomain_data=subdomains)
ds = Measure("ds", domain=mesh, subdomain_data=facet_tags)
dS = Measure("dS", domain=mesh, subdomain_data=facet_tags)
V = FunctionSpace(mesh, ("DG", 1))

# Dirichlet boundary conditions for left and right ends (equals to zero)
fdim = mesh.topology.dim-1
u_bc_left = Function(V)
u_bc_right = Function(V)
facets_left = np.array(facet_tags.indices[facet_tags.values == 1])
facets_right = np.array(facet_tags.indices[facet_tags.values == 2])
dofs_left = locate_dofs_topological(V, fdim, facets_left)
dofs_right = locate_dofs_topological(V, fdim, facets_right)
bc_left = dirichletbc(u_bc_left, dofs_left)
bc_right = dirichletbc(u_bc_right, dofs_right)
bcs=[bc_left,bc_right]

u = TrialFunction(V)
v = TestFunction(V)
n = FacetNormal(mesh)
h = 2 * Circumradius(mesh)
c = 347
alpha = 10
gamma = 10
h_avg = avg(h)

mass  = -c**2 * (inner(grad(u), grad(v))*dx + inner(inner(grad(u), n), v) * ds)
# Add DG/IP terms
mass +=  c**2 * (inner(jump(u, n), avg(grad(v)) )*dS + inner(avg(grad(u)), jump(v, n))*dS)
mass += -c**2*gamma/h_avg*inner(jump(u, n), jump(v, n))*dS
# Add Nitsche terms
mass +=  c**2*(inner(u, inner(grad(v), n)) * ds - alpha / h * inner(u, v) * ds)

mass_form = form(mass)

mpc = MultiPointConstraint(V)
mpc.finalize()

A = assemble_matrix(mass_form, mpc, bcs=bcs)
A.assemble()

stiffness = inner(u , v) * dx
stiffnes_form = form(stiffness)

C = assemble_matrix(stiffnes_form, mpc)
C.assemble()

target = 694*2*np.pi

nev=6

# This function prints the results of the eps_solver1 function below
def results(E):
    print("******************************")
    print("*** SLEPc Solution Results ***")
    print("******************************")
    print("Number of iterations of the method: %d" % E.getIterationNumber())
    print("Solution method: %s" % E.getType())
    nev, ncv, mpd = E.getDimensions()
    print("Number of requested eigenvalues: %d" % nev)
    print("Stopping condition: tol=%.4g, maxit=%d" % (E.getTolerances()))
    print("Number of converged eigenpairs %d" % E.getConverged())
    if E.getConverged() > 0:
        print()
    for i in range(E.getConverged()):
        k = E.getEigenvalue(i)
        w = np.sqrt(k)
        f = w/np.pi/2
        print("%6f, %6f" % (f.real, f.imag))
    print()

# This function solves eigensystem
def eps_solver1(A, C, target, nev, two_sided=False, print_results=False):
    E = SLEPc.EPS().create(MPI.COMM_WORLD)
    C = - C
    E.setOperators(A, C)
    st = E.getST()
    st.setType('sinvert')
    E.setTarget(target)
    E.setWhichEigenpairs(SLEPc.EPS.Which.TARGET_MAGNITUDE)  # TARGET_REAL or TARGET_IMAGINARY
    E.setDimensions(nev, SLEPc.DECIDE)
    E.setTolerances(1e-15)
    E.setFromOptions()
    E.solve()
    if print_results:
        results(E)
    return E

E = eps_solver1(A, C, target**2, nev=nev, print_results= True)

vr = Function(V)
vi = Function(V)

ind = 0
E.getEigenvector(ind, vr.vector, vi.vector)

import matplotlib.pyplot as plt
fig, ax = plt.subplots(2, figsize=(12, 6))
ax[0].plot(V.tabulate_dof_coordinates()[:, 0], vr.x.array.real)
ax[1].plot(V.tabulate_dof_coordinates()[:, 0], vr.x.array.imag)
plt.show()

Since mpc is empty, I expect to get the same eigenvalues and eigenvectors as in the original post. However, the computed eigenvector looks to be quite spurious. However if do the same thing with a “CG” FunctionSpace I obtain the same eigenvectors with and without an empty mpc.

PS: I use doflinx in complex mode.

4

You cannot use:

with DG, for the same reason you cannot use it for creating MPCs, as explained above.
You need to use locate_dofs_geometrical when working with DG, as the degrees of freedom are not associated with vertices, edges or facets, but cells.

5

Sorry my bad, I just copied the code from there. Below is the corrected code where we still obtain a spurious eigenvector (while if I use the assemble_matrix of dolfinx.fem.petsc, a “physical” eigenvector is plotted):

# From https://fenicsproject.discourse.group/t/discontinuous-1d-eigenvector-using-dg-elements/10491
from dolfinx.fem import FunctionSpace, form, Function,dirichletbc,locate_dofs_geometrical
from dolfinx.mesh import meshtags,locate_entities, create_interval
from ufl import TestFunction,TrialFunction,Measure,inner,grad,jump,avg,Circumradius, FacetNormal
from dolfinx_mpc import MultiPointConstraint, assemble_matrix
import numpy as np
from slepc4py import SLEPc
from mpi4py import MPI

def OneDimensionalSetup(n_elem, x_MP = 0.2, r_outer = 0.25):
    # 1 ------------------------------------ 2
    #                 x=r_outer
    """ This function builds one dimensional setup.
        For boundaries, Tag 1 specifies left end and Tag 2 specifies right end. 
    Returns:
        mesh, subdomains, facet_tags
    """
    mesh = create_interval(MPI.COMM_WORLD, n_elem, [0.0, r_outer])
    def subdomain_func(x, x_MP=x_MP, eps=1e-16):
        x = x[0]
        return np.logical_and(x_MP - eps <= x, x <= x_MP + eps)
    tdim = mesh.topology.dim
    marked_cells = locate_entities(mesh, tdim, subdomain_func)
    fl = 0
    subdomains = meshtags(mesh, tdim, marked_cells, np.full(len(marked_cells), fl, dtype=np.int32))
    boundaries = [  (0, lambda x: np.isclose(x[0], x_MP)),
                    (1, lambda x: np.isclose(x[0], 0)),
                    (2, lambda x: np.isclose(x[0], r_outer))]
    facet_indices, facet_markers = [], []
    fdim = mesh.topology.dim - 1
    for (marker, locator) in boundaries:
        facets = locate_entities(mesh, fdim, locator)
        facet_indices.append(facets)
        facet_markers.append(np.full(len(facets), marker))
    facet_indices = np.array(np.hstack(facet_indices), dtype=np.int32)
    facet_markers = np.array(np.hstack(facet_markers), dtype=np.int32)
    sorted_facets = np.argsort(facet_indices)
    facet_tags = meshtags(mesh, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets])
    return mesh, subdomains, facet_tags

n_elem = 50
r_outer = 0.25
mesh, subdomains, facet_tags = OneDimensionalSetup(n_elem, r_outer)
dx = Measure("dx", domain=mesh, subdomain_data=subdomains)
ds = Measure("ds", domain=mesh, subdomain_data=facet_tags)
dS = Measure("dS", domain=mesh, subdomain_data=facet_tags)
V = FunctionSpace(mesh, ("DG", 1))

# Dirichlet boundary conditions for left and right ends (equals to zero)
fdim = mesh.topology.dim-1
u_bc_left = Function(V)
u_bc_right = Function(V)
dofs_left = locate_dofs_geometrical(V, lambda x: np.isclose(x[0], 0.0))
dofs_right = locate_dofs_geometrical(V, lambda x: np.isclose(x[0], r_outer))
assert(dofs_left.size==dofs_right.size)
bc_left = dirichletbc(u_bc_left, dofs_left)
bc_right = dirichletbc(u_bc_right, dofs_right)
bcs=[bc_left,bc_right]

u = TrialFunction(V)
v = TestFunction(V)
n = FacetNormal(mesh)
h = 2 * Circumradius(mesh)
c = 347
alpha = 10
gamma = 10
h_avg = avg(h)

mass  = -c**2 * (inner(grad(u), grad(v))*dx + inner(inner(grad(u), n), v) * ds)
# Add DG/IP terms
mass +=  c**2 * (inner(jump(u, n), avg(grad(v)) )*dS + inner(avg(grad(u)), jump(v, n))*dS)
mass += -c**2*gamma/h_avg*inner(jump(u, n), jump(v, n))*dS
# Add Nitsche terms
mass +=  c**2*(inner(u, inner(grad(v), n)) * ds - alpha / h * inner(u, v) * ds)

mass_form = form(mass)

mpc = MultiPointConstraint(V)
mpc.finalize()

A = assemble_matrix(mass_form, mpc, bcs=bcs)
A.assemble()

stiffness = inner(u , v) * dx
stiffnes_form = form(stiffness)

C = assemble_matrix(stiffnes_form, mpc)
C.assemble()

target = 694*2*np.pi

nev=6

# This function prints the results of the eps_solver1 function below
def results(E):
    print("******************************")
    print("*** SLEPc Solution Results ***")
    print("******************************")
    print("Number of iterations of the method: %d" % E.getIterationNumber())
    print("Solution method: %s" % E.getType())
    nev, ncv, mpd = E.getDimensions()
    print("Number of requested eigenvalues: %d" % nev)
    print("Stopping condition: tol=%.4g, maxit=%d" % (E.getTolerances()))
    print("Number of converged eigenpairs %d" % E.getConverged())
    if E.getConverged() > 0:
        print()
    for i in range(E.getConverged()):
        k = E.getEigenvalue(i)
        w = np.sqrt(k)
        f = w/np.pi/2
        print("%6f, %6f" % (f.real, f.imag))
    print()

# This function solves eigensystem
def eps_solver1(A, C, target, nev, two_sided=False, print_results=False):
    E = SLEPc.EPS().create(MPI.COMM_WORLD)
    C = - C
    E.setOperators(A, C)
    st = E.getST()
    st.setType('sinvert')
    E.setTarget(target)
    E.setWhichEigenpairs(SLEPc.EPS.Which.TARGET_MAGNITUDE)  # TARGET_REAL or TARGET_IMAGINARY
    E.setDimensions(nev, SLEPc.DECIDE)
    E.setTolerances(1e-15)
    E.setFromOptions()
    E.solve()
    if print_results:
        results(E)
    return E

E = eps_solver1(A, C, target**2, nev=nev, print_results= True)

vr = Function(V)
vi = Function(V)

ind = 0
E.getEigenvector(ind, vr.vector, vi.vector)

import matplotlib.pyplot as plt
fig, ax = plt.subplots(2, figsize=(12, 6))
ax[0].plot(V.tabulate_dof_coordinates()[:, 0], vr.x.array.real)
ax[1].plot(V.tabulate_dof_coordinates()[:, 0], vr.x.array.imag)
plt.show()
6

First of all, the code is not valid when you set X_MP to r_outer.
You are passing it in wrongly to your function:

as r_outer is interpreted as X_MP.
This should be:
mesh, subdomains, facet_tags = OneDimensionalSetup(n_elem, r_outer=r_outer)

However, irregardless, you are right the the integrals aren’t correct. This is due to me not implementing dS integrals as part of the integration kernels: dolfinx_mpc/cpp/assemble_matrix.cpp at main · jorgensd/dolfinx_mpc · GitHub
Feel free to open an issue in the Github repo, and I’ll see if I have time to implement it over christmas.

1 Like
7

Thanks for the answer. I created an issue.