Changing the coefficients of a multi-point constraint#
Author Jørgen S. Dokken
A multi-point constraint relates each slave degree of freedom \(u_s\) to a set of masters,
Creating and finalizing a constraint is expensive: the masters are located, communicated and added as ghosts. When only the coefficients \(c_{sj}\) change between solves, for instance in a Floquet-Bloch condition \(u(x + L) = e^{i k L} u(x)\) swept over the wave number \(k\), the constraint can be kept and its coefficients modified instead. This demo sets up a periodic constraint and changes its coefficients in three ways, inspecting the coefficients after each step. Nothing is solved.
A Bloch phase is complex, so we use the complex scalar type matching the mesh precision, whether or not DOLFINx was built with complex PETSc.
dtype = np.result_type(default_real_type, np.complex64).type
mesh = create_unit_square(MPI.COMM_WORLD, 4, 4)
V = fem.functionspace(mesh, ("Lagrange", 1))
Setting up the constraint#
We make the right boundary \(x=1\) periodic with the left one, \(u(1, y) = u(0, y)\). To also show masters that are eliminated by a Dirichlet condition, we fix \(u(0, y) = 1 + y\) on the left boundary only. The slaves on the right boundary then have masters with known values, and these are moved into the offset \(g_s\) of the constraint.
tol = 500 * np.finfo(default_real_type).eps
def left(x):
return np.isclose(x[0], 0, atol=tol)
def right(x):
return np.isclose(x[0], 1, atol=tol)
def periodic_relation(x):
out = x.copy()
out[0] = x[0] - 1
return out
g = fem.Function(V, dtype=dtype)
g.interpolate(lambda x: 1 + x[1])
bc = fem.dirichletbc(g, fem.locate_dofs_geometrical(V, left))
The masters of a constraint are fixed when it is created. By default, a master
whose coefficient is below coefficient_tol times the largest of its slave is
dropped, and could then never be given a coefficient later. Passing
coefficient_tol=0 keeps every master. It comes at a cost,
as each master is a ghost and an entry of the sparsity pattern:
mpc = dolfinx_mpc.MultiPointConstraint(V, dtype=dtype, bcs=[bc])
mpc.create_periodic_constraint_geometrical(V, right, periodic_relation, [bc], scale=dtype(1), coefficient_tol=0)
mpc.finalize()
mpc_cut = dolfinx_mpc.MultiPointConstraint(V, dtype=dtype, bcs=[bc])
mpc_cut.create_periodic_constraint_geometrical(V, right, periodic_relation, [bc], scale=dtype(1))
mpc_cut.finalize()
num_masters = mesh.comm.allreduce(len(mpc.all_masters()), op=MPI.SUM)
num_masters_cut = mesh.comm.allreduce(len(mpc_cut.all_masters()), op=MPI.SUM)
if mesh.comm.rank == 0:
print(f"Masters with coefficient_tol=0: {num_masters}, with the default: {num_masters_cut}")
Masters with coefficient_tol=0: 15, with the default: 5
For a P1 space every slave gets the three vertices of the cell its periodic image lies in, two of them with coefficient zero.
Inspecting the coefficients#
all_coefficients
returns the coefficients of all masters per local degree of freedom, and
all_masters the
corresponding masters. Unlike
coefficients, these
include the masters eliminated by the Dirichlet condition. We mark those by
checking which masters are absent from
masters.
The offset \(g_s\) is given by
constants.
x_dofs = V.tabulate_dof_coordinates()
def print_constraint(mpc, header, num_slaves=3):
"""Print the non-zero terms of the constraint of the first owned slaves of this process."""
coeffs, offsets = mpc.all_coefficients()
masters = mpc.all_masters()
lines = [header]
for s in mpc.slaves[: mpc.num_local_slaves][:num_slaves]:
kept = mpc.masters.links(s)
terms = []
for m, c in zip(masters[offsets[s] : offsets[s + 1]], coeffs[offsets[s] : offsets[s + 1]]):
if abs(c) > tol:
terms.append(f"({c:.3f}) u_{m}" + ("" if m in kept else " [Dirichlet]"))
lines.append(f" y={x_dofs[s, 1]:.2f}: u_{s} = " + " + ".join(terms) + f" g = {mpc.constants[s]:.3f}")
print(f"Rank {mesh.comm.rank}: " + "\n".join(lines), flush=True)
print_constraint(mpc, "Created with scale 1")
Rank 0: Created with scale 1
y=1.00: u_13 = (1.000+0.000j) u_1 [Dirichlet] g = 2.000+0.000j
y=0.75: u_15 = (1.000+0.000j) u_0 [Dirichlet] g = 1.750+0.000j
y=0.50: u_19 = (1.000+0.000j) u_5 [Dirichlet] g = 1.500+0.000j
Each slave has a single master with a non-zero coefficient, the dof at the same height on the left boundary. As that one is fixed by the Dirichlet condition, the constraint reduces to \(u_s = c\, g(0, y) = 1 + y\), which is the offset \(g\).
The figure joins each slave on the right boundary to its masters: the master at the same
height, with coefficient one, by a straight line, and the masters kept by tol=None, with
coefficient zero, by arcs. Each process draws its own cells and slaves, which are gathered
on the first process.
pyvista.global_theme.allow_empty_mesh = True
def link(a, b, bulge, num_points=17):
"""Points along a quadratic arc from `a` to `b`, bulging sideways by `bulge` times its length."""
d = b - a
control = 0.5 * (a + b) + bulge * np.array([-d[1], d[0], 0.0])
t = np.linspace(0, 1, num_points)[:, None]
return (1 - t) ** 2 * a + 2 * t * (1 - t) * control + t**2 * b
x_all = mpc.function_space.tabulate_dof_coordinates()
coeffs, offsets = mpc.all_coefficients()
masters = mpc.all_masters()
links: dict[bool, list[np.ndarray]] = {True: [], False: []}
for s in mpc.slaves[: mpc.num_local_slaves]:
for m, c in zip(masters[offsets[s] : offsets[s + 1]], coeffs[offsets[s] : offsets[s + 1]]):
nonzero = bool(abs(c) > tol)
links[nonzero].append(link(x_all[s], x_all[m], 0.0 if nonzero else 0.15))
slaves = x_all[mpc.slaves[: mpc.num_local_slaves]]
owned_cells = np.arange(mesh.topology.index_map(mesh.topology.dim).size_local, dtype=np.int32)
cells = pyvista.UnstructuredGrid(*plot.vtk_mesh(mesh, entities=owned_cells))
gathered = mesh.comm.gather((cells, links, slaves), root=0)
def polylines(curves):
"""The curves as one PolyData of polylines."""
points = np.vstack(curves)
n = len(curves[0])
lines = np.hstack([[n, *range(i * n, (i + 1) * n)] for i in range(len(curves))])
return pyvista.PolyData(points, lines=lines)
if gathered is not None: # only the first process received the pieces
plotter = pyvista.Plotter(window_size=[600, 600])
for g in gathered:
plotter.add_mesh(g[0], style="wireframe", color="lightgray")
for nonzero, colour, width, label in (
(False, "steelblue", 2, "coefficient 0"),
(True, "orange", 5, "coefficient 1"),
):
curves = [curve for g in gathered for curve in g[1][nonzero]]
plotter.add_mesh(polylines(curves), color=colour, line_width=width, label=label)
plotter.add_points(np.vstack([g[2] for g in gathered]), color="black", point_size=12, render_points_as_spheres=True)
plotter.add_legend(bcolor="white", face="line", size=(0.3, 0.1), loc="upper center")
plotter.view_xy()
# The figure is named after the demo, as the gallery of the documentation expects
figure = Path("demo_scaling_coefficients.py")
if pyvista.OFF_SCREEN:
plotter.screenshot(figure.with_suffix(".png"))
else:
# The interactive scene, for the gallery
plotter.export_vtksz(figure.with_suffix(".vtksz"))
plotter.show(screenshot=figure.with_suffix(".png"))
/tmp/ipykernel_3456/874201897.py:53: PyVistaDeprecationWarning: Plotter.export_vtksz is deprecated. Use plotter.trame.export_vtksz(...) instead. Install via `pip install trame-pyvista`.
plotter.export_vtksz(figure.with_suffix(".vtksz"))
Scaling by a constant#
scale_coefficients
multiplies all coefficients of slave \(s\) by a factor \(f_s\), and recomputes the
offset. Masters eliminated by the Dirichlet condition are scaled too, so the
offset becomes \(f_s c\, g\).
mpc.scale_coefficients(dtype(0.5))
print_constraint(mpc, "Scaled by 0.5")
Rank 0: Scaled by 0.5
y=1.00: u_13 = (0.500+0.000j) u_1 [Dirichlet] g = 1.000+0.000j
y=0.75: u_15 = (0.500+0.000j) u_0 [Dirichlet] g = 0.875+0.000j
y=0.50: u_19 = (0.500+0.000j) u_5 [Dirichlet] g = 0.750+0.000j
Scaling by an expression#
The factor can also vary in space and be given as a UFL expression. It is
interpolated into a function in the space of the constraint, and \(f_s\) is its
value at the slave. For the Bloch phase \(f = e^{ikL}\) we use a
Constant for \(k\), and compile the expression
once, so that it is not recompiled when \(k\) changes.
Repeated calls compound, so the coefficients are now \(0.5 e^{ikL}\).
k = fem.Constant(mesh, dtype(0))
L = 1.0
phase = fem.Expression(ufl.exp(1j * k * L), V.element.interpolation_points, dtype=dtype)
k.value = np.pi / 2
mpc.scale_coefficients(phase)
print_constraint(mpc, "Scaled by exp(i pi/2)")
Rank 0: Scaled by exp(i pi/2)
y=1.00: u_13 = (0.000+0.500j) u_1 [Dirichlet] g = 0.000+1.000j
y=0.75: u_15 = (0.000+0.500j) u_0 [Dirichlet] g = 0.000+0.875j
y=0.50: u_19 = (0.000+0.500j) u_5 [Dirichlet] g = 0.000+0.750j
Replacing the coefficients#
update_coefficients
replaces all coefficients, in the layout of all_coefficients. Here we replace
the coefficient \(c_1\) of the (non-zero) Dirichlet master \(u_k\) of each slave by
\(c_3 = 2\). The slave becomes \(u_s = c_3 u_k = c_3 g_k\), and the offset is
recomputed accordingly. The zero coefficients of the other masters are left
as they are.
coeffs, offsets = mpc.all_coefficients()
masters = mpc.all_masters()
new_coeffs = coeffs.copy()
for s in mpc.slaves:
kept = mpc.masters.links(s)
for j in range(offsets[s], offsets[s + 1]):
if masters[j] not in kept and abs(coeffs[j]) > tol:
new_coeffs[j] = 2
mpc.update_coefficients(new_coeffs)
print_constraint(mpc, "Dirichlet masters set to 2")
Rank 0: Dirichlet masters set to 2
y=1.00: u_13 = (2.000+0.000j) u_1 [Dirichlet] g = 4.000+0.000j
y=0.75: u_15 = (2.000+0.000j) u_0 [Dirichlet] g = 3.500-0.000j
y=0.50: u_19 = (2.000+0.000j) u_5 [Dirichlet] g = 3.000-0.000j
Note that mpc.slaves includes the ghosted slaves: every process must hold the
same coefficients for a slave, so the ghosts are updated as well.
Sweeping a parameter#
As scaling compounds, a sweep over \(k\) stores the original coefficients once,
restores them with update_coefficients, and scales with the current phase.
Only the coefficients change; the masters, index maps and sparsity pattern are
reused for every \(k\).
Rank 0: k = 0.000
y=1.00: u_13 = (1.000+0.000j) u_1 [Dirichlet] g = 2.000+0.000j
Rank 0: k = 0.785
y=1.00: u_13 = (0.707+0.707j) u_1 [Dirichlet] g = 1.414+1.414j
Rank 0: k = 3.142
y=1.00: u_13 = (-1.000+0.000j) u_1 [Dirichlet] g = -2.000+0.000j