Periodic conditions for vector valued spaces

cristian · October 10, 2024, 10:04am

Hi, I’m trying to solve a periodic problem with curl elements. I get the following error: “RuntimeError: Periodic conditions for vector valued spaces are not implemented”.
I’m using dolfinx_mpc v0.8.1 through the docker image provided by Jørgen.
I would be happy if it worked just for conformal mesh.

Here below is my MWE. I’m using dolfinx-complex-mode (i.e., source /usr/local/bin/dolfinx-complex-mode):

import dolfinx
import ufl
import dolfinx_mpc
import gmsh
from dolfinx.fem.petsc import LinearProblem
from dolfinx.io.gmshio import model_to_mesh
from petsc4py import PETSc
from mpi4py import MPI
import numpy as np
from dolfinx import fem, mesh, common
from basix.ufl import element


Lx = 0.1
Ly = 1


# Define region tags
air_tag = 1
port_1_tag = 2 #input port
port_2_tag = 3 #output port
left_tag = 4
right_tag = 5

# Generate mesh
gmsh.initialize()
gmsh.model.add("geometry")
r1 = gmsh.model.occ.addRectangle(-Lx/2, -Ly/2, 0, Lx, Ly)
gmsh.model.occ.synchronize()
gmsh.model.addPhysicalGroup(2, [1], tag=1, name="air")  
gmsh.model.addPhysicalGroup(1, [1], tag=2, name="bottom")  
gmsh.model.addPhysicalGroup(1, [2], tag=3, name="left")  
gmsh.model.addPhysicalGroup(1, [3], tag=4, name="right")    
gmsh.model.addPhysicalGroup(1, [4], tag=5, name="top")  

gmsh.option.setNumber("Mesh.MeshSizeMax", 0.01)

gmsh.model.mesh.generate(2)
gmsh.write("mesh.msh")
domain, cell_tags, facet_tags = model_to_mesh(gmsh.model, MPI.COMM_WORLD, 0, gdim=3)
gmsh.finalize()


k0 = dolfinx.fem.Constant(domain, PETSc.ScalarType(2 * np.pi / 0.8))
epsr = dolfinx.fem.Constant(domain, PETSc.ScalarType(1))


E0=ufl.as_vector((1, 0, 0))

#Define elements, spaces and measures
degree = 2
curl_el = element('N1curl', domain.basix_cell(), degree)
V = dolfinx.fem.functionspace(domain, curl_el)
dx = ufl.Measure("dx", domain, subdomain_data=cell_tags, metadata={'quadrature_degree': 4})
ds = ufl.Measure("exterior_facet", domain, subdomain_data=facet_tags,metadata={'quadrature_degree': 4})

E = ufl.TrialFunction(V)
testE = ufl.TestFunction(V)

n = ufl.FacetNormal(domain)

# Define problem
F = (- ufl.inner(ufl.curl(E), ufl.curl(testE)) * dx
    + epsr * k0**2 * ufl.inner(E, testE)  * dx(air_tag)
    + 1j*k0 * ufl.inner(ufl.cross(E, n), ufl.cross(testE, n)) * ds((port_1_tag,port_2_tag))
    - 2j*k0 * ufl.inner(ufl.cross(E0, n), ufl.cross(testE, n)) * ds(port_1_tag)
)
a, L = ufl.lhs(F), ufl.rhs(F)


# Periodic Boundary condition
bcs=[]
def periodic_relation(x):
    out_x = np.zeros(x.shape)
    out_x[0] = Lx-x[0]
    out_x[1] = x[1]
    out_x[2] = x[2]
    return out_x

facets = facet_tags.find(right_tag)
arg_sort = np.argsort(facets)
mt = mesh.meshtags(domain, domain.topology.dim - 1, facets[arg_sort], np.full(len(facets), 2, dtype=np.int32))
with common.Timer("~~Periodic: Compute mpc condition"):
    mpc = dolfinx_mpc.MultiPointConstraint(V)
    mpc.create_periodic_constraint_topological(V, mt, 2, periodic_relation, bcs)
    mpc.finalize()

# Solve problem
problem = dolfinx_mpc.LinearProblem(a, L,  mpc, bcs, petsc_options={"ksp_type": "preonly", "pc_type": "lu","pc_factor_mat_solver_type": "mumps"})#
E=problem.solve()

dokken · October 10, 2024, 1:44pm

For a conforming grid this could probably be done manually (with create general constraint). Some work is still required to extend it in general.

I have been thinking about an alternative approach where we create a «periodic» mesh. I will keep you posted on this development.

cristian · October 17, 2024, 1:30pm

Thank you Jorgen. I’m trying to implement this with create_general_constraint but I’m not sure how to define the dof’s location as Nedelec elements don’t have a pointwise evaluation (tabulate_dof_coordinates gives a runtime error). Any suggestion on how could I implement it?

dokken · October 17, 2024, 3:27pm

If you know what facets should be connected together, one should be able to find the entity closure dof of each facet, and their index (maybe with a -1) coefficient due to the normal vector).

cristian · October 18, 2024, 12:36pm

Not sure I follow, I still need to pass to create_general_constraint a dictionary the dof’s coordinates, not the dof’s index.

dokken · October 18, 2024, 3:10pm

You can use the lowermost interface of MPC for dof indices: dolfinx_mpc — DOLFINx-MPC: An extension to DOLFINx for multi point constraints

cristian · October 24, 2024, 7:15pm

Hi Jorgen. I think I was able to figure it out (see code). However there seems to be a random issue where the solutions is not correct for a given meshing, for example if I set gmsh.option.setNumber(“Mesh.MeshSizeMax”, 0.021) it will work other if I use gmsh.option.setNumber(“Mesh.MeshSizeMax”, 0.02) it doesn’t (see figures). Any idea what the problem could be?

import ufl
import dolfinx_mpc
import gmsh
from dolfinx.io.gmshio import model_to_mesh
from petsc4py import PETSc
from mpi4py import MPI
import numpy as np
from dolfinx import fem
from basix.ufl import element
from mpi4py import MPI
from dolfinx.io import VTXWriter

comm = MPI.COMM_WORLD

#routine to save fields
def save(E,filename):
    gdim=domain.geometry.dim
    W = fem.functionspace(domain,("Discontinuous Lagrange", 2, (3,)))
    E_dg = fem.Function(W)
    E_dg.interpolate(fem.Expression(E, W.element.interpolation_points()))

    with VTXWriter(domain.comm, "sols/E.bp", E_dg) as f:
        f.write(0.0)
    if MPI.COMM_WORLD.rank == 0:
        print("Info    : Saved solution as "+filename)


# Define region tags
air_tag = 1
port_1_tag = 2 #input port
port_2_tag = 4 #output port
left_tag = 3
right_tag = 5
pec_tag = (3,5)

# Generate mesh
Lx = 0.2
Ly = 1
gmsh.initialize()
gmsh.model.add("geometry")
r1 = gmsh.model.occ.addRectangle(-Lx/2, -Ly/2, 0, Lx, Ly)
gmsh.model.occ.synchronize()
gmsh.model.addPhysicalGroup(2, [1], tag=1, name="air")  
gmsh.model.addPhysicalGroup(1, [1], tag=2, name="bottom")  
gmsh.model.addPhysicalGroup(1, [2], tag=3, name="left")  
gmsh.model.addPhysicalGroup(1, [3], tag=4, name="right")    
gmsh.model.addPhysicalGroup(1, [4], tag=5, name="top")  

translation = [1, 0, 0, Lx,   # X-axis translation component
               0, 1, 0, 0,   # Y-axis translation component
               0, 0, 1, 0,   # Z-axis translation component
               0, 0, 0, 1]
gmsh.model.mesh.setPeriodic(1, [3], [2], translation)

gmsh.option.setNumber("Mesh.MeshSizeMax", 0.02)

gmsh.model.mesh.generate(2)
gmsh.write("mesh.msh")
domain, cell_tags, facet_tags = model_to_mesh(gmsh.model, MPI.COMM_WORLD, 0, gdim=3)
gmsh.finalize()

#Define elements, spaces and measures
degree = 2
curl_el = element('N1curl', domain.basix_cell(), degree)
V = dolfinx.fem.functionspace(domain, curl_el)
dx = ufl.Measure("dx", domain, subdomain_data=cell_tags, metadata={'quadrature_degree': 4})
ds = ufl.Measure("exterior_facet", domain, subdomain_data=facet_tags,metadata={'quadrature_degree': 4})

E = ufl.TrialFunction(V)
testE = ufl.TestFunction(V)

n = ufl.FacetNormal(domain)
k0 = dolfinx.fem.Constant(domain, PETSc.ScalarType(2 * np.pi / 0.4))

# Define problem
F = (- ufl.inner(ufl.curl(E), ufl.curl(testE)) * dx
    + k0**2 * ufl.inner(E, testE)  * dx
    + 1j*k0 * ufl.inner(ufl.cross(E, n), ufl.cross(testE, n)) * ds((port_1_tag,port_2_tag))
    + 2j*k0 * ufl.inner(ufl.cross(ufl.as_vector((1, 0,0)), n), ufl.cross(testE, n)) * ds(port_2_tag)
)
a, L = ufl.lhs(F), ufl.rhs(F)

# mpc constraints
mpc = dolfinx_mpc.MultiPointConstraint(V)
master_dofs = dolfinx.fem.locate_dofs_topological(V, domain.topology.dim - 1, facet_tags.find(left_tag))
owners = np.full_like(master_dofs, comm.rank)
slave_dofs = dolfinx.fem.locate_dofs_topological(V, domain.topology.dim - 1, facet_tags.find(right_tag))
coeffs = np.ones_like(master_dofs, dtype=PETSc.ScalarType)
offsets = np.arange(master_dofs.size+1, dtype=np.int64)
mpc.add_constraint(V, slave_dofs, master_dofs, coeffs, owners, offsets)
mpc.finalize()

# Solve problem
problem = dolfinx_mpc.LinearProblem(a, L, mpc, bcs=[], petsc_options={"ksp_type": "preonly", "pc_type": "lu","pc_factor_mat_solver_type": "mumps"})#
E=problem.solve()

# Save solution
save(E,'E.bp')

cristian · October 25, 2024, 9:07pm

I’m trying to make the code work with mph. I’ve modified the mpc section as follow but I keep getting a “malloc(): invalid size (unsorted)” error.

# Find the master dof 
master_dofs = fem.locate_dofs_topological(V, domain.topology.dim - 1, facet_tags.find(left_tag))
owners = np.full_like(master_dofs, domain.comm.rank)
imap = V.dofmap.index_map
master_dofs_global=imap.local_to_global(master_dofs)

# Find the slave dofs
slave_dofs = fem.locate_dofs_topological(V, domain.topology.dim - 1, facet_tags.find(right_tag))

offsets = np.arange(master_dofs_global.size+1, dtype=np.int64)
coeffs = np.ones_like(master_dofs, dtype=PETSc.ScalarType)

# Create constraint and finalize mpc
mpc.add_constraint(V, slave_dofs, master_dofs_global, coeffs, owners, offsets)
mpc.finalize()

dokken · November 4, 2024, 12:53pm

@cristian

I tried to get the MPC strategy to work with Nedelec elements, but got into similar issues as you did (it is probably related to the facet orientation and dof orientation).

Therefore, I decided to go down another route that I have wanted to try for a long time, which is to make the mesh topology periodic, as one would no longer need specialized assembly routines.

Following is a prototype that works in serial:

# Create a periodic mesh in serial
# SPDX-License-Identifier: MIT
# Author: Jørgen S. Dokken

from mpi4py import MPI
import numpy as np
import dolfinx
import ufl


def create_periodic_mesh(mesh, indicator, mapping_function):
    """
    Create a periodic mesh that takes all facets that satisfy the `indicator` function,
    and map the vertices of these facets to the vertices that satisfies the mapping function.

    Example:

        .. code-block:: python
        .. highlight:: python
    
              mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 7, 19)
              def indicator(x):
                  return numpy.isclose(x[1], 1)
      
              def map(x):
                  values = x.copy()
                  values[1] -= 1
                  return values
      
              periodic_mesh = create_periodic_mesh(mesh, indicator, map)        
    """    

    geometry = mesh.geometry._cpp_object
    topology = mesh.topology

    mesh.topology.create_connectivity(mesh.topology.dim, mesh.topology.dim-1)
    mesh.topology.create_connectivity(mesh.topology.dim-1, mesh.topology.dim)

    # Map left side to right side
    left_facets = dolfinx.mesh.locate_entities_boundary(mesh, mesh.topology.dim - 1, indicator)
    left_vertices = dolfinx.mesh.compute_incident_entities(mesh.topology, left_facets, mesh.topology.dim -1, 0)
    left_midpoints = dolfinx.mesh.compute_midpoints(mesh, 0, left_vertices)

    right_midpoints = mapping_function(left_midpoints.T).T
    right_vertices = dolfinx.mesh.locate_entities_boundary(mesh, 0, lambda x: x[0] > 1.0 - 1e-10)

    # Find closest vertex on right side
    bb_tree = dolfinx.geometry.bb_tree(mesh,0, right_vertices)
    mid_tree = dolfinx.geometry.create_midpoint_tree(mesh, 0, right_vertices)
    closest_vertex = dolfinx.geometry.compute_closest_entity(bb_tree, mid_tree, mesh, right_midpoints)

    # Keep only left side vertices
    num_vertices_local = mesh.topology.index_map(0).size_local + mesh.topology.index_map(0).num_ghosts
    keep_vertices = np.ones(num_vertices_local, dtype=np.bool_)
    keep_vertices[right_vertices] = False
    # Create submap
    new_vertices = np.flatnonzero(keep_vertices)
    new_vertex_map, sub_to_parent = dolfinx.cpp.common.create_sub_index_map(mesh.topology.index_map(0), new_vertices, allow_owner_change=False)

    # Invert map
    num_vertices_local = mesh.topology.index_map(0).size_local + mesh.topology.index_map(0).num_ghosts
    parent_to_sub = np.full(num_vertices_local, -1, dtype=np.int32)
    parent_to_sub[sub_to_parent] = np.arange(sub_to_parent.size, dtype=np.int32)

    # Create map from right to left vertices for replacement
    replace_map = np.arange(num_vertices_local, dtype=np.int32)
    replace_map[closest_vertex] = left_vertices

    # First map vertices from right to left vertices
    c_to_v = mesh.topology.connectivity(mesh.topology.dim, 0)
    new_c = replace_map[c_to_v.array]
    # Then map vertices to new indices
    new_c = parent_to_sub[new_c]
    new_o = c_to_v.offsets.copy()
    new_c_to_v = dolfinx.graph.adjacencylist(new_c, new_o)
    new_v_to_v = dolfinx.graph.adjacencylist(np.arange(len(sub_to_parent), dtype=np.int32))

    topology = dolfinx.cpp.mesh.Topology(MPI.COMM_WORLD, mesh.topology.cell_type)
    topology.set_index_map(0, new_vertex_map)
    topology.set_index_map(mesh.topology.dim, mesh.topology.index_map(mesh.topology.dim))
    topology.set_connectivity(new_v_to_v, 0,0)
    topology.set_connectivity(new_c_to_v, mesh.topology.dim, 0)
    c_el = dolfinx.fem.coordinate_element(mesh._ufl_domain.ufl_coordinate_element().basix_element)
    geometry = dolfinx.mesh.create_geometry(mesh.geometry.index_map(), mesh.geometry.dofmap, c_el._cpp_object, mesh.geometry.x[:, :mesh.geometry.dim].copy(),  mesh.geometry.input_global_indices)
    if mesh.geometry.x.dtype == np.float64:
        cpp_mesh = dolfinx.cpp.mesh.Mesh_float64(mesh.comm, topology, geometry._cpp_object)
    elif mesh.geometry.x.dtype == np.float32:
        cpp_mesh = dolfinx.cpp.mesh.Mesh_float32(mesh.comm, topology, geometry._cpp_object)
    else:
        raise RuntimeError(f"Unsupported dtype for mesh {mesh.geometry.x.dtype}")

    new_mesh = dolfinx.mesh.Mesh(cpp_mesh, domain = ufl.Mesh(mesh._ufl_domain.ufl_coordinate_element()))

    return new_mesh



mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 19, 13)



def indicator(x):
    return np.isclose(x[0], 0.0)

def mapping(x):
    values = x.copy()
    values[0] += 1
    return values


new_mesh = create_periodic_mesh(mesh, indicator, mapping)




with dolfinx.io.XDMFFile(MPI.COMM_WORLD, "periodic_mesh.xdmf", "w") as xdmf:
    xdmf.write_mesh(new_mesh)


V = dolfinx.fem.functionspace(new_mesh, ("Lagrange", 2, (new_mesh.geometry.dim, )))
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
a = ufl.inner(u, v) * ufl.dx
x = ufl.SpatialCoordinate(new_mesh)
f = ufl.as_vector([5*x[0]*ufl.sin(3*np.pi * x[0])+2*x[1], x[0]*ufl.sin(3*np.pi * x[1])])
L = ufl.inner(f, v) * ufl.dx
import dolfinx.fem.petsc
problem = dolfinx.fem.petsc.LinearProblem(a, L, bcs=[], petsc_options = {"ksp_type": "preonly", "pc_type": "lu", "pc_factor_mat_solver_type": "mumps"})
uh = problem.solve()


V_out = dolfinx.fem.functionspace(new_mesh, ("DG", 2, (new_mesh.geometry.dim, )))
u_out = dolfinx.fem.Function(V_out)
u_out.interpolate(uh)
with dolfinx.io.VTXWriter(new_mesh.comm, "u_periodic.bp", [u_out]) as writer:
    writer.write(0.0)

yielding:

for second order Lagrange, and with N1curl (considering tangential continuity)

Topic		Replies	Views
Periodic Boundary Conditons for Mixed Functionspace with dolfinx_mpc General	10	1118	January 4, 2023
Periodic boundary constraints for vector-valued functionspaces in dolfinx_mpc, when V.sub(i) is insufficient dolfinx	0	227	August 26, 2023
Periodic boundary condition with DG and MPC General	6	361	December 6, 2023
Periodic boundary condition with non-periodic mesh using dolfinx_mpc General	1	446	July 21, 2023
Simple periodic boundary condition without dolfinx_mpc dolfinx boundary-conditions	1	36	November 4, 2024

Periodic conditions for vector valued spaces

Related Topics