Generate an inner product between RHS vector and solution vector in Dolfinx 0.9.0

Linear Algebra
1

Hello all,

A question about doing an inner product on RHS and solution vectors: How do I do this, taking into account the scattering of vector components across processes? Code is below. Pertinent code block is at the bottom.

# Example to test multi-mode waveguide absorbing boundary condition
# for a simple rectangular waveguide
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
from dolfinx.mesh import locate_entities_boundary, meshtags
from dolfinx.io import gmshio
from dolfinx.fem import petsc
from dolfinx import fem
from dolfinx import io
import time

eta0 = 377.0 # free space wave impedance

class ModesTE:
    def __init__(self, beta, k, a, b, m, n):
        self.beta = beta
        self.k = k
        self.a = a
        self.b = b
        self.m = m
        self.n = n
    def eval_H(self, x):
        ZTE = eta0 * self.k / self.beta
        if((self.m != 0) and (self.n != 0)):
            B = np.sqrt(self.n * self.m * 4.0 / (ZTE * self.a * self.b))
        elif self.m == 0:
            B = np.sqrt(self.n * 2.0 / (ZTE * self.a * self.b))
        else:
            B = np.sqrt(self.m * 2.0 / (ZTE * self.a * self.b))
        hx = B * np.sin((np.pi * self.m / self.a) * x[0]) * np.cos((np.pi * self.n / self.b) * x[1])
        hy = B * np.cos((np.pi * self.m / self.a) * x[0]) * np.sin((np.pi * self.n / self.b) * x[1])
        hz = np.full_like(hx, 0.0+0j, dtype=np.complex128)
        return(hx, hy, hz)
    def eval_E(self, x):
        if((self.m != 0) and (self.n != 0)):
            B = np.sqrt(self.n * self.m * 4.0 / (ZTE * self.a * self.b))
        elif self.m == 0:
            B = np.sqrt(self.n * 2.0 / (ZTE * self.a * self.b))
        else:
            B = np.sqrt(self.m * 2.0 / (ZTE * self.a * self.b))
        ZTE = eta0 * self.k / self.beta
        ey = -ZTE * B * np.sin((np.pi * self.m / self.a) * x[0]) * np.cos((np.pi * self.n / self.b) * x[1])
        ex = ZTE * B * np.cos((np.pi * self.m / self.a) * x[0]) * np.sin((np.pi * self.n / self.b) * x[1])
        ez = np.full_like(hx, 0.0+0j, dtype=np.complex128)
        return(ex, ey, ez)

eps = 1.0e-4
a = 22.86 #WG width
b = 10.16 #WG height
l = 50.8 # WG length
f = 15.5 #frequency of operation GHz
c = 300.0 #speed of light
lm = 1.00

comm = M.COMM_WORLD
mpiRank = comm.rank
modelRank = 0 # where GMSH runs
k0 = 2.0 * np.pi * f / c #wavenumber wrt frequency

if mpiRank == modelRank:
    gmsh.initialize()
    gmsh.option.setNumber('General.Terminal', 1)
    gmsh.model.add('Waveguide')
    gmsh.model.occ.addBox(0, 0, 0, a, b, l, 1)
    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', 2.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, [1], 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)
    gmsh.fltk.run()


mesh, ct, fm = gmshio.model_to_mesh(gmsh.model, comm, modelRank, gdim=3)
if mpiRank == modelRank:
    gmsh.finalize()

mesh.topology.create_connectivity(mesh.topology.dim-1, mesh.topology.dim)
elem = basix.ufl.element('Nedelec 1st kind H(curl)', mesh.basix_cell(), degree=2)
V = fem.functionspace(mesh, elem)
u = U.TrialFunction(V)
v = U.TestFunction(V)
normal = U.FacetNormal(mesh)
ds = U.Measure("ds", domain=mesh, subdomain_data = fm)
# Dirichlet BCs
facets = fm.find(1)
ubc = fem.Function(V)
ubc.interpolate(lambda x: np.array([0*x[0], 0*x[1], 0*x[2]], dtype=np.complex128))
dofs = fem.locate_dofs_topological(V, mesh.topology.dim-1, facets)
bc = fem.dirichletbc(ubc, dofs)
# Build rhs
uh = fem.Function(V, dtype=np.complex128)
kz = np.sqrt(k0 * k0 - (np.pi / a)**2)
f = ModesTE(kz, k0, a, b, 1, 0)
uh.interpolate(f.eval_H)
finc = 2.0j * k0 * eta0 * U.inner(U.cross(normal, uh), v) * ds(2)
uh.name = "H_inc"
print(k0, kz)
# Generate mode inner product for projection operator
ua = fem.Function(V, dtype=np.complex128)
kz = np.sqrt(k0 * k0 - (np.pi / a)**2)
Zm = eta0 * k0 / kz
f = ModesTE(kz, k0, a, b, 1, 0) #Mode 1
ua.interpolate(f.eval_H)
um1a = fem.form(U.inner(U.cross(ua, normal), v) * ds(2)) # mode on input plane
um1b = fem.form(U.inner(U.cross(ua, normal), v) * ds(3)) #mode on output plane
t_start = time.perf_counter()

# Weak form
L = (U.inner(U.curl(u), U.curl(v)) - k0 * k0 * U.inner(u, v)) * U.dx
Df = fem.form(L)
Rf = fem.form(finc)

num_dofs_local = V.dofmap.index_map.size_local * V.dofmap.index_map_bs
num_dofs_global = V.dofmap.index_map.size_global * V.dofmap.index_map_bs
print(f"Number of dofs (owned) by rank {mpiRank}: {num_dofs_local}")
if mpiRank == modelRank:
    print(f"Number of dofs global: {num_dofs_global}")

# Build RHSs for mode vector solves
def VecCreate(linear_form, bilin_form, BCs):
    vec = petsc.create_vector(linear_form)
    petsc.assemble_vector(vec, linear_form)
    petsc.apply_lifting(vec, [bilin_form], [[BCs]])
    vec.ghostUpdate(addv=P.InsertMode.ADD_VALUES, mode=P.ScatterMode.REVERSE)
    petsc.set_bc(vec, [BCs])
    vec.ghostUpdate(addv=P.InsertMode.INSERT, mode=P.ScatterMode.FORWARD)
    return vec

#A = petsc.create_matrix(Df)
#A.zeroEntries()
A = petsc.assemble_matrix(Df, bcs=[bc])
A.assemble()
b = VecCreate(Rf, Df, bc)
u1 = VecCreate(um1a, Df, bc)
u2 = VecCreate(um1b, Df, bc)

petsc_options = {
        "ksp_type": "preonly",
        "pc_type": "lu",
        "pc_factor_mat_solver": "mumps",
        "pc_factor_nonzeros_along_diagonal": 1e-8,
        "ksp_error_if_not_converged": False,
        "ksp_monitor": None }
opts = P.Options()
for key, value in petsc_options.items():
    opts[key] = value

solver = P.KSP().create(mesh.comm)
solver.setOperators(A)
solver.setFromOptions()
solver.setUp()
t_stop = time.perf_counter()

print("Solver setup time = {0}".format(t_stop - t_start))
sol = fem.Function(V)
sol.name="Solution"
sol1 = fem.Function(V)
sol1.name="BC1"
sol2 = fem.Function(V)
sol2.name="BC2"

solver.solve(b, sol.x.petsc_vec)
sol.x.petsc_vec.ghostUpdate(addv=P.InsertMode.INSERT, mode=P.ScatterMode.FORWARD)
t_solve = time.perf_counter()

print("Solve time = {0}".format(t_solve - t_stop))

solver.solve(u1, sol1.x.petsc_vec)
sol1.x.petsc_vec.ghostUpdate(addv=P.InsertMode.INSERT, mode=P.ScatterMode.FORWARD)
t_solve1 = time.perf_counter()
P1 = np.vdot(sol.x.array, sol.x.array) #This works OK
print(mpiRank)
print(u1.array)
print(sol.x.array)

#XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
#u1 is a RHS vector constructed with dolfinx.petsc.assemble_vector
#sol is a solution vector.  How do I compute the following inner product?
print(np.vdot(u1.array, sol.x.array)) # Would like to calculate this vector inner product over all processes.
#XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX

print("Solve1 time = {0}".format(t_solve1 - t_solve))

solver.solve(u2, sol2.x.petsc_vec)
sol2.x.petsc_vec.ghostUpdate(addv=P.InsertMode.INSERT, mode=P.ScatterMode.FORWARD)
t_solve2 = time.perf_counter()

print("Solve1 time = {0}".format(t_solve2 - t_solve1))
sys.exit(0)
```

Cheers,

Bill

2

Since you use the vectors wrapped as PETSc.Vec, you can use dot (petsc4py.PETSc.Vec — Python 3.24.0 documentation) to compute the “l2” inner product.

3

I tried using

L2norm = u1.array.dot(sol.x.array)

and it gives the error

ValueError:shapes (45132,) and (46814,) not aligned 45132 (dim 0) != 47692 (dim=0).

(and a similar error for each process)

The vectors should have the same length. Is there something with “ghosted” values that I need to account for? (Full code for EZ running is below.) I need to compute the inner product of the full vector and redistribute it to all processes.


# Example to test multi-mode waveguide absorbing boundary condition
# for a simple rectangular waveguide
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
from dolfinx.mesh import locate_entities_boundary, meshtags
from dolfinx.io import gmshio
from dolfinx.fem import petsc
from dolfinx import fem
from dolfinx import io
import time

eta0 = 377.0 # free space wave impedance

class ModesTE:
    def __init__(self, beta, k, a, b, m, n):
        self.beta = beta
        self.k = k
        self.a = a
        self.b = b
        self.m = m
        self.n = n
    def eval_H(self, x):
        ZTE = eta0 * self.k / self.beta
        if((self.m != 0) and (self.n != 0)):
            B = np.sqrt(self.n * self.m * 4.0 / (ZTE * self.a * self.b))
        elif self.m == 0:
            B = np.sqrt(self.n * 2.0 / (ZTE * self.a * self.b))
        else:
            B = np.sqrt(self.m * 2.0 / (ZTE * self.a * self.b))
        hx = B * np.sin((np.pi * self.m / self.a) * x[0]) * np.cos((np.pi * self.n / self.b) * x[1])
        hy = B * np.cos((np.pi * self.m / self.a) * x[0]) * np.sin((np.pi * self.n / self.b) * x[1])
        hz = np.full_like(hx, 0.0+0j, dtype=np.complex128)
        return(hx, hy, hz)
    def eval_E(self, x):
        if((self.m != 0) and (self.n != 0)):
            B = np.sqrt(self.n * self.m * 4.0 / (ZTE * self.a * self.b))
        elif self.m == 0:
            B = np.sqrt(self.n * 2.0 / (ZTE * self.a * self.b))
        else:
            B = np.sqrt(self.m * 2.0 / (ZTE * self.a * self.b))
        ZTE = eta0 * self.k / self.beta
        ey = -ZTE * B * np.sin((np.pi * self.m / self.a) * x[0]) * np.cos((np.pi * self.n / self.b) * x[1])
        ex = ZTE * B * np.cos((np.pi * self.m / self.a) * x[0]) * np.sin((np.pi * self.n / self.b) * x[1])
        ez = np.full_like(hx, 0.0+0j, dtype=np.complex128)
        return(ex, ey, ez)

eps = 1.0e-4
a = 22.86 #WG width
b = 10.16 #WG height
l = 50.8 # WG length
f = 15.5 #frequency of operation GHz
c = 300.0 #speed of light
lm = 1.00

comm = M.COMM_WORLD
mpiRank = comm.rank
modelRank = 0 # where GMSH runs
k0 = 2.0 * np.pi * f / c #wavenumber wrt frequency

if mpiRank == modelRank:
    gmsh.initialize()
    gmsh.option.setNumber('General.Terminal', 1)
    gmsh.model.add('Waveguide')
    gmsh.model.occ.addBox(0, 0, 0, a, b, l, 1)
    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', 2.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, [1], 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)
    gmsh.fltk.run()


mesh, ct, fm = gmshio.model_to_mesh(gmsh.model, comm, modelRank, gdim=3)
if mpiRank == modelRank:
    gmsh.finalize()

mesh.topology.create_connectivity(mesh.topology.dim-1, mesh.topology.dim)
elem = basix.ufl.element('Nedelec 1st kind H(curl)', mesh.basix_cell(), degree=2)
V = fem.functionspace(mesh, elem)
u = U.TrialFunction(V)
v = U.TestFunction(V)
normal = U.FacetNormal(mesh)
ds = U.Measure("ds", domain=mesh, subdomain_data = fm)
# Dirichlet BCs
facets = fm.find(1)
ubc = fem.Function(V)
ubc.interpolate(lambda x: np.array([0*x[0], 0*x[1], 0*x[2]], dtype=np.complex128))
dofs = fem.locate_dofs_topological(V, mesh.topology.dim-1, facets)
bc = fem.dirichletbc(ubc, dofs)
# Build rhs
uh = fem.Function(V, dtype=np.complex128)
kz = np.sqrt(k0 * k0 - (np.pi / a)**2)
f = ModesTE(kz, k0, a, b, 1, 0)
uh.interpolate(f.eval_H)
finc = 2.0j * k0 * eta0 * U.inner(U.cross(normal, uh), v) * ds(2)
uh.name = "H_inc"
print(k0, kz)
# Generate mode inner product for projection operator
ua = fem.Function(V, dtype=np.complex128)
kz = np.sqrt(k0 * k0 - (np.pi / a)**2)
Zm = eta0 * k0 / kz
f = ModesTE(kz, k0, a, b, 1, 0) #Mode 1
ua.interpolate(f.eval_H)
um1a = fem.form(U.inner(U.cross(ua, normal), v) * ds(2)) # mode on input plane
um1b = fem.form(U.inner(U.cross(ua, normal), v) * ds(3)) #mode on output plane
t_start = time.perf_counter()

# Weak form
L = (U.inner(U.curl(u), U.curl(v)) - k0 * k0 * U.inner(u, v)) * U.dx
Df = fem.form(L)
Rf = fem.form(finc)

num_dofs_local = V.dofmap.index_map.size_local * V.dofmap.index_map_bs
num_dofs_global = V.dofmap.index_map.size_global * V.dofmap.index_map_bs
print(f"Number of dofs (owned) by rank {mpiRank}: {num_dofs_local}")
if mpiRank == modelRank:
    print(f"Number of dofs global: {num_dofs_global}")

# Build RHSs for mode vector solves
def VecCreate(linear_form, bilin_form, BCs):
    vec = petsc.create_vector(linear_form)
    petsc.assemble_vector(vec, linear_form)
    petsc.apply_lifting(vec, [bilin_form], [[BCs]])
    vec.ghostUpdate(addv=P.InsertMode.ADD_VALUES, mode=P.ScatterMode.REVERSE)
    petsc.set_bc(vec, [BCs])
    vec.ghostUpdate(addv=P.InsertMode.INSERT, mode=P.ScatterMode.FORWARD)
    return vec

#A = petsc.create_matrix(Df)
#A.zeroEntries()
A = petsc.assemble_matrix(Df, bcs=[bc])
A.assemble()
b = VecCreate(Rf, Df, bc)
u1 = VecCreate(um1a, Df, bc)
u2 = VecCreate(um1b, Df, bc)

petsc_options = {
        "ksp_type": "preonly",
        "pc_type": "lu",
        "pc_factor_mat_solver": "mumps",
        "pc_factor_nonzeros_along_diagonal": 1e-8,
        "ksp_error_if_not_converged": False,
        "ksp_monitor": None }
opts = P.Options()
for key, value in petsc_options.items():
    opts[key] = value

solver = P.KSP().create(mesh.comm)
solver.setOperators(A)
solver.setFromOptions()
solver.setUp()
t_stop = time.perf_counter()

print("Solver setup time = {0}".format(t_stop - t_start))
sol = fem.Function(V)
sol.name="Solution"
sol1 = fem.Function(V)
sol1.name="BC1"
sol2 = fem.Function(V)
sol2.name="BC2"

solver.solve(b, sol.x.petsc_vec)
sol.x.petsc_vec.ghostUpdate(addv=P.InsertMode.INSERT, mode=P.ScatterMode.FORWARD)
t_solve = time.perf_counter()

print("Solve time = {0}".format(t_solve - t_stop))

solver.solve(u1, sol1.x.petsc_vec)
sol1.x.petsc_vec.ghostUpdate(addv=P.InsertMode.INSERT, mode=P.ScatterMode.FORWARD)
t_solve1 = time.perf_counter()
P1 = np.vdot(sol.x.array, sol.x.array)
print(mpiRank)
print(u1.array)
print(sol.x.array)
#u1 is a RHS vector constructed with dolfinx.petsc.assemble_vector
#sol is a solution vector.  How do I compute the following inner product?
print(u1.array.dot(sol.x.array)) # Would like to calculate this vector inner product over all processes.
print("Solve1 time = {0}".format(t_solve1 - t_solve))

solver.solve(u2, sol2.x.petsc_vec)
sol2.x.petsc_vec.ghostUpdate(addv=P.InsertMode.INSERT, mode=P.ScatterMode.FORWARD)
t_solve2 = time.perf_counter()

print("Solve1 time = {0}".format(t_solve2 - t_solve1))
sys.exit(0)
```
4

Try using u1.dot(sol.x.petsc_vec), as you want it to work on the PETSc Vectors, not the underlying numpy arrays.

5

Yes! That seems to work. Thanks!

cheers!