In using the “allreduce” statement to generate an integral over a boundary surface that computes the surface area of a boundary, something seems to have changed in the way this is done. (It works in DolfinX 0.9.0.) The boundaries in question are sub-boundaries labeled as such in the GMSH code (labels 2 and 3).
The problematic lines are at the end of this almost-running example
import sys
import numpy as np
from mpi4py import MPI as M
from petsc4py import PETSc as P
import gmsh
import ufl as U
import basix.ufl
import dolfinx
from dolfinx.io import gmsh as gmshio
from dolfinx.fem import petsc
from dolfinx import mesh as msh
from dolfinx import fem
from dolfinx import io
from dolfinx import plot
import scifem
import time
try:
import pyvista
have_pyvista = True
except ModuleNotFoundError:
print("Pyvista required to visualize solution.")
have_pyvista = False
t_start = time.perf_counter()
eps = 1.0e-4
a = 22.86 #WG width
b = 10.16 #WG height
l = 50.8 # WG length
posx = a/3
posz = l/2
rad = 1.0
f = 16.5 #frequency of operation GHz
c = 300.0 #speed of light
#lm = 0.565
lm = 1.2
mode = 1
model = None
comm = M.COMM_WORLD
mpiRank = comm.rank
mpiSize = comm.size
print("MPI rank = {0}, MPI size = {1}.".format(mpiRank, mpiSize))
modelRank = 0 # where GMSH runs
k0 = 2.0 * np.pi * f / c #wavenumber wrt frequency
if mpiRank == modelRank:
WithObstruction = 1
gmsh.initialize()
gmsh.option.setNumber('General.Terminal', 1)
gmsh.model.add('Waveguide')
if(WithObstruction != 0):
gmsh.model.occ.addBox(0, 0, 0, a, b, l, 1)
gmsh.model.occ.addCylinder(posx, 0, posz, 0, b, 0, rad, 2, np.pi * 2.0)
gmsh.model.occ.addCylinder(posx/2, 0, posz, 0, b, 0, rad, 20, np.pi * 2.0)
gmsh.model.occ.cut([(3,1)], [(3,2),(3,20)], 3, removeObject=True, removeTool=True)
else:
gmsh.model.occ.addBox(0, 0, 0, a, b, l, 3)
gmsh.model.occ.synchronize()
pt = gmsh.model.getEntities(dim=0)
gmsh.model.mesh.setSize(pt, lm)
gmsh.option.setNumber('Mesh.MeshSizeMin', 0.01)
gmsh.option.setNumber('Mesh.MeshSizeMax', 5.2)
gmsh.option.setNumber('Mesh.Algorithm', 6) #1=Mesh Adapt, 2=Auto, 3=Initial mesh only, 5=Delaunay, 6=Frontal-Delaunay
gmsh.option.setNumber('Mesh.MinimumCirclePoints', 30)
gmsh.option.setNumber('Mesh.CharacteristicLengthFromCurvature', 1)
gmsh.option.setNumber('Mesh.Tetrahedra', 0)
OUT = []
IN = []
PEC = []
bb = gmsh.model.getEntities(dim=3)
pt = gmsh.model.getBoundary(bb, combined=True, oriented=False, recursive=False)
for bnd in pt:
CoM = gmsh.model.occ.getCenterOfMass(bnd[0], bnd[1])
if np.abs(CoM[2]) < eps:
IN.append(bnd[1])
elif np.abs(CoM[2] - l) < eps:
OUT.append(bnd[1])
else:
PEC.append(bnd[1])
gmsh.model.addPhysicalGroup(3, [3], 1) #WG
gmsh.model.setPhysicalName(3, 1, "Waveguide")
gmsh.model.addPhysicalGroup(2, PEC, 1)
gmsh.model.setPhysicalName(2, 1, "PEC")
gmsh.model.addPhysicalGroup(2, IN, 2)
gmsh.model.setPhysicalName(2, 2, "INPORT")
gmsh.model.addPhysicalGroup(2, OUT, 3)
gmsh.model.setPhysicalName(2, 3, "OUTPORT")
gmsh.model.mesh.generate(3)
model = gmsh.model
model = comm.bcast(model, root=modelRank)
mesh_data = gmshio.model_to_mesh(model, comm, modelRank, gdim=3)
if mpiRank == modelRank:
gmsh.finalize()
comm.barrier()
mesh_data.mesh.topology.create_connectivity(mesh_data.mesh.topology.dim-1, mesh_data.mesh.topology.dim)
exterior_facets=msh.exterior_facet_indices(mesh_data.mesh.topology)
fm = mesh_data.facet_tags
ct = mesh_data.cell_tags
#Real function space to be used on surface tags 2 and 3 for mode 1
Z1 = scifem.create_real_functionspace(mesh_data.mesh, value_shape=(4,))
zu = U.TrialFunction(Z1)
zv = U.TestFunction(Z1)
#Internal E field solution space
elem = basix.ufl.element('N1curl', mesh_data.mesh.basix_cell(), degree=2)
V = fem.functionspace(mesh_data.mesh, elem)
u = U.TrialFunction(V)
v = U.TestFunction(V)
n = U.FacetNormal(mesh_data.mesh)
# The integration domain over the input and output ports (boundary tags 2 and 3)
ds2 = U.Measure("ds", domain=mesh_data.mesh, subdomain_data=fm, subdomain_id=2)
ds3 = U.Measure("ds", domain=mesh_data.mesh, subdomain_data=fm, subdomain_id=3)
# Weak form. Need areas over the input and output ports in order to normalize mode coefficients
# This is where the program crashes.
Area1 = mesh_data.mesh.comm.allreduce(fem.assemble_scalar(fem.form(U.inner(n, n) * ds2)), op=M.SUM)
Area2 = mesh_data.mesh.comm.allreduce(fem.assemble_scalar(fem.form(U.inner(n, n) * ds3)), op=M.SUM)
print("Area 1 = {0}, Area 2 = {1}".format(Area1, Area2))
sys.exit(0)
I get an “incompatible function arguments” error:
Traceback (most recent call last):
File "/home/bill/WaveguideBC/Test.py", line 120, in <module>
Area1 = mesh_data.mesh.comm.allreduce(fem.assemble_scalar(fem.form(U.inner(n, n) * ds2)), op=M.SUM)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
TypeError: __init__(): incompatible function arguments. The following argument types are supported:
1. __init__(self, spaces: collections.abc.Sequence[dolfinx.cpp.fem.FunctionSpace_float64], integrals: collections.abc.Mapping[dolfinx.cpp.fem._IntegralType, collections.abc.Sequence[tuple[int, int, ndarray[dtype=int32, shape=(*), order='C', writable=False], ndarray[dtype=int8, shape=(*), order='C', writable=False]]]], coefficients: collections.abc.Sequence[dolfinx.cpp.fem.Function_complex128], constants: collections.abc.Sequence[dolfinx.cpp.fem.Constant_complex128], need_permutation_data: bool, entity_maps: collections.abc.Sequence[dolfinx.cpp.mesh.EntityMap], mesh: dolfinx.cpp.mesh.Mesh_float64) -> None
2. __init__(self, form: collections.abc.Sequence[int], spaces: collections.abc.Sequence[dolfinx.cpp.fem.FunctionSpace_float64], coefficients: collections.abc.Sequence[dolfinx.cpp.fem.Function_complex128], constants: collections.abc.Sequence[dolfinx.cpp.fem.Constant_complex128], subdomains: collections.abc.Mapping[dolfinx.cpp.fem._IntegralType, collections.abc.Sequence[tuple[int, ndarray[dtype=int32, shape=(*), order='C', writable=False]]]], entity_maps: collections.abc.Sequence[dolfinx.cpp.mesh.EntityMap], mesh: dolfinx.cpp.mesh.Mesh_float64 | None) -> None
I looked at some tutorial examples that use this form integration method and I cannot see what I am doing incorrectly. What am I missing?