Good day everyone,
I am trying to build a right hand side vector where I need to use a conditional statement to exclude a possible “divide by zero” condition. The key code elements are the Hinc class at the beginning and the “Build RHS” segment at the end of the code snippet. Code follows:
import sys
import numpy as np
from scipy import optimize as opt
from scipy import special as sp
from mpi4py import MPI as nMPI
from petsc4py import PETSc
import gmsh
import meshio
import ufl
import basix.ufl
from dolfinx.mesh import locate_entities_boundary, meshtags
from dolfinx.io import gmshio
from dolfinx import fem
from dolfinx.fem.petsc import LinearProblem
from dolfinx import io
eps = 1.0e-4 # global geometric tolerance
# Incident field
class Hinc:
def __init__(self, b, kc, H):
self.beta = b
self.kc = kc
self.H = H
def eval(self, x):
rr = np.sqrt(x[0] * x[0] + x[1] * x[1])
s = ufl.conditional(ufl.gt(rr, eps), x[1]/rr, 0.0)
c = ufl.conditional(ufl.gt(rr, eps), x[0]/rr, 1.0)
hphi = ufl.conditional(ufl.gt(rr, eps), 1j * self.beta * self.H * sp.jvp(1, rr, 0) * s / (self.kc * self.kc * rr), 0.0)
hr = -1j * self.beta * self.H * sp.jvp(1, rr, 1) * c / self.kc
hz = np.full_like(hr, 0.0+0j, dtype=np.complex128)
hx = hr * c - hphi * s
hy = hr * s + hphi * c
return(hx, hy, hz)
def Model(x):
comm = nMPI.COMM_WORLD
mpiRank = comm.rank
modelRank = 0 # where GMSH runs
a = x[0] # Waveguide radius
h = x[1] # horn length
b = x[2] # horn mouth radius
wl1 = 1.0
wl2 = 1.0 # waveguide lengths
r = 6.0 # rad sphere radius
t = 0.1 # WG thickness
Plot = 1
k = 2.0 * np.pi * x[3] / 3.0e10 # k0 as a function of frequency
eta0 = 377.0
kc = 1.84 / a
beta = np.sqrt(k * k - kc * kc)
lm = 0.3
ls = 0.025
if mpiRank == modelRank:
print("Waveguide radius = {0:<f}, Horn length = {1:<f}, Horn radius = {2:<f}, Freq = {3:<f}".format(a, h, b, x[3]))
gmsh.initialize()
gmsh.option.setNumber('General.Terminal', 1)
gmsh.model.add('Feed Horn')
#Rad hemisphere
gmsh.model.occ.addSphere(0, 0, 0, r, 1)
gmsh.model.occ.addBox(-r, -r, 0, 2*r, 2*r, r, 2)
gmsh.model.occ.intersect([(3,1)], [(3,2)], 3, removeObject=True, removeTool=True)
# Stick short waveguide on bottom
gmsh.model.occ.addCylinder(0, 0, -wl2, 0, 0, wl2, a, 2, 2.0*np.pi)
gmsh.model.occ.fuse([(3,3)], [(3,2)], 1, removeObject=True, removeTool=True)
# Make cutout for horn geometry
gmsh.model.occ.addCylinder(0, 0, 0, 0, 0, wl1+h, b+t, 2, 2.0*np.pi)
gmsh.model.occ.addCylinder(0, 0, 0, 0, 0, wl1, a, 3, 2.0*np.pi)
gmsh.model.occ.addCone(0, 0, wl1, 0, 0, h, a, b, 4, 2.0*np.pi)
gmsh.model.occ.cut([(3,2)], [(3,3), (3,4)], 5, removeObject=True, removeTool=True)
gmsh.model.occ.cut([(3,1)], [(3,5)], 6, removeObject=True, removeTool=True)
# Let's only model 1/4 of the geometry to reduce compute needs.
gmsh.model.occ.addBox(0, 0, -wl2, r, r, r+wl2, 2)
gmsh.model.occ.intersect([(3,6)], [(3,2)], 1, removeObject=True, removeTool=True)
gmsh.model.occ.synchronize()
pt = gmsh.model.getEntities(dim=0)
gmsh.model.mesh.setSize(pt, lm)
pt = gmsh.model.getEntitiesInBoundingBox(-eps, -eps, wl1+h-eps, b+t+eps, b+t+eps, wl1+h+eps)
gmsh.model.mesh.setSize(pt, ls)
pt
gmsh.option.setNumber('Mesh.MeshSizeMin', 0.01)
gmsh.option.setNumber('Mesh.MeshSizeMax', 1.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.MshFileVersion', 2.2)
gmsh.option.setNumber('Mesh.Format', 1)
gmsh.option.setNumber('Mesh.MinimumCirclePoints', 30)
gmsh.option.setNumber('Mesh.CharacteristicLengthFromCurvature', 1)
gmsh.option.setNumber('Mesh.Tetrahedra', 0)
bb = gmsh.model.getEntities(dim=3)
pt = gmsh.model.getBoundary(bb, combined=True, oriented=False, recursive=False)
CoaxPort = []
RAD = []
PMC = []
PEC = []
for bnd in pt:
CoM = gmsh.model.occ.getCenterOfMass(bnd[0], bnd[1])
print(CoM)
if np.abs(CoM[2] + wl2) < eps:
CoaxPort.append(bnd[1])
elif np.abs(CoM[0]) < eps:
PMC.append(bnd[1])
elif (np.abs(CoM[2]) > wl1 + h + eps) and (np.abs(CoM[0]) > eps) and (np.abs(CoM[1]) > eps):
RAD.append(bnd[1])
else:
PEC.append(bnd[1])
gmsh.model.addPhysicalGroup(3, [1], 1) # rad zone
gmsh.model.setPhysicalName(3, 1, "RadZone")
gmsh.model.addPhysicalGroup(2, PEC, 1)
gmsh.model.setPhysicalName(2, 1, "PEC")
gmsh.model.addPhysicalGroup(2, PMC, 2)
gmsh.model.setPhysicalName(2, 2, "PMC")
gmsh.model.addPhysicalGroup(2, CoaxPort, 3)
gmsh.model.setPhysicalName(2, 3, "In Port")
gmsh.model.addPhysicalGroup(2, RAD, 4)
gmsh.model.setPhysicalName(2, 4, "Rad Surface")
gmsh.model.mesh.generate(3)
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 = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
n = ufl.FacetNormal(mesh)
ds = ufl.Measure("ds", domain=mesh, subdomain_data=fm)
# Dirichlet BCs on PECs
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)
f = Hinc(beta, kc, 1.0+0j)
uh.interpolate(f.eval)
finc = 2.0j * k0 * eta0 * ufl.inner(ufl.cross(n, uh), v) * ds(3)
uh.name = "Hinc"
return 0
x=np.array([1.115, 1.0, 2.23, 1.0e10])
e = Model(x)
print(e)
sys.exit(0)
I know that I am doing something wrong in the UFL conditional statements. How do I do the testing correctly? I get the following error.
Traceback (most recent call last):
File "/usr/lib/petscdir/petsc3.15/x86_64-linux-gnu-complex/lib/python3/dist-packages/dolfinx/fem/function.py", line 459, in interpolate
_interpolate(u0)
File "/usr/lib/python3.10/functools.py", line 889, in wrapper
return dispatch(args[0].__class__)(*args, **kw)
File "/usr/lib/petscdir/petsc3.15/x86_64-linux-gnu-complex/lib/python3/dist-packages/dolfinx/fem/function.py", line 440, in _interpolate
self._cpp_object.interpolate(u0, cells0, cells1) # type: ignore
TypeError: interpolate(): incompatible function arguments. The following argument types are supported:
1. interpolate(self, f: ndarray[dtype=complex128, shape=(*), order='C', writable=False], cells: ndarray[dtype=int32, shape=(*), order='C', writable=False]) -> None
2. interpolate(self, f: ndarray[dtype=complex128, shape=(*, *), order='C', writable=False], cells: ndarray[dtype=int32, shape=(*), order='C', writable=False]) -> None
3. interpolate(self, u: dolfinx.cpp.fem.Function_complex128, cells0: ndarray[dtype=int32, shape=(*), order='C', writable=False], cells1: ndarray[dtype=int32, shape=(*), order='C', writable=False]) -> None
4. interpolate(self, u: dolfinx.cpp.fem.Function_complex128, cells: ndarray[dtype=int32, shape=(*), order='C', writable=False], interpolation_data: dolfinx.cpp.geometry.PointOwnershipData_float64) -> None
5. interpolate(self, e0: dolfinx.cpp.fem.Expression_complex128, cells0: ndarray[dtype=int32, order='C', writable=False], cells1: ndarray[dtype=int32, order='C', writable=False]) -> None
Invoked with types: dolfinx.cpp.fem.Function_complex128, method, ndarray, ndarray
During handling of the above exception, another exception occurred:
Traceback (most recent call last):
File "/home/bill/Cluster/Fenics2020/FeedHorn/SimpleCircHorn.py", line 194, in <module>
Traceback (most recent call last):
File "/usr/lib/petscdir/petsc3.15/x86_64-linux-gnu-complex/lib/python3/dist-packages/dolfinx/fem/function.py", line 459, in interpolate
_interpolate(u0)
File "/usr/lib/python3.10/functools.py", line 889, in wrapper
return dispatch(args[0].__class__)(*args, **kw)
File "/usr/lib/petscdir/petsc3.15/x86_64-linux-gnu-complex/lib/python3/dist-packages/dolfinx/fem/function.py", line 440, in _interpolate
self._cpp_object.interpolate(u0, cells0, cells1) # type: ignore
TypeError: interpolate(): incompatible function arguments. The following argument types are supported:
1. interpolate(self, f: ndarray[dtype=complex128, shape=(*), order='C', writable=False], cells: ndarray[dtype=int32, shape=(*), order='C', writable=False]) -> None
2. interpolate(self, f: ndarray[dtype=complex128, shape=(*, *), order='C', writable=False], cells: ndarray[dtype=int32, shape=(*), order='C', writable=False]) -> None
3. interpolate(self, u: dolfinx.cpp.fem.Function_complex128, cells0: ndarray[dtype=int32, shape=(*), order='C', writable=False], cells1: ndarray[dtype=int32, shape=(*), order='C', writable=False]) -> None
4. interpolate(self, u: dolfinx.cpp.fem.Function_complex128, cells: ndarray[dtype=int32, shape=(*), order='C', writable=False], interpolation_data: dolfinx.cpp.geometry.PointOwnershipData_float64) -> None
5. interpolate(self, e0: dolfinx.cpp.fem.Expression_complex128, cells0: ndarray[dtype=int32, order='C', writable=False], cells1: ndarray[dtype=int32, order='C', writable=False]) -> None
Invoked with types: dolfinx.cpp.fem.Function_complex128, method, ndarray, ndarray
During handling of the above exception, another exception occurred:
Traceback (most recent call last):
File "/home/bill/Cluster/Fenics2020/FeedHorn/SimpleCircHorn.py", line 194, in <module>
e = Model(x)
File "/home/bill/Cluster/Fenics2020/FeedHorn/SimpleCircHorn.py", line 154, in Model
Traceback (most recent call last):
File "/usr/lib/petscdir/petsc3.15/x86_64-linux-gnu-complex/lib/python3/dist-packages/dolfinx/fem/function.py", line 459, in interpolate
_interpolate(u0)
File "/usr/lib/python3.10/functools.py", line 889, in wrapper
return dispatch(args[0].__class__)(*args, **kw)
File "/usr/lib/petscdir/petsc3.15/x86_64-linux-gnu-complex/lib/python3/dist-packages/dolfinx/fem/function.py", line 440, in _interpolate
self._cpp_object.interpolate(u0, cells0, cells1) # type: ignore
TypeError: interpolate(): incompatible function arguments. The following argument types are supported:
1. interpolate(self, f: ndarray[dtype=complex128, shape=(*), order='C', writable=False], cells: ndarray[dtype=int32, shape=(*), order='C', writable=False]) -> None
2. interpolate(self, f: ndarray[dtype=complex128, shape=(*, *), order='C', writable=False], cells: ndarray[dtype=int32, shape=(*), order='C', writable=False]) -> None
3. interpolate(self, u: dolfinx.cpp.fem.Function_complex128, cells0: ndarray[dtype=int32, shape=(*), order='C', writable=False], cells1: ndarray[dtype=int32, shape=(*), order='C', writable=False]) -> None
4. interpolate(self, u: dolfinx.cpp.fem.Function_complex128, cells: ndarray[dtype=int32, shape=(*), order='C', writable=False], interpolation_data: dolfinx.cpp.geometry.PointOwnershipData_float64) -> None
5. interpolate(self, e0: dolfinx.cpp.fem.Expression_complex128, cells0: ndarray[dtype=int32, order='C', writable=False], cells1: ndarray[dtype=int32, order='C', writable=False]) -> None
Invoked with types: dolfinx.cpp.fem.Function_complex128, method, ndarray, ndarray
During handling of the above exception, another exception occurred:
Traceback (most recent call last):
File "/home/bill/Cluster/Fenics2020/FeedHorn/SimpleCircHorn.py", line 194, in <module>
e = Model(x)
EDIT:
I found the error. I fell into the rabbit hole of mixing Numpy and UFL! It turns out I needed to use the generalized conditional “numpy.where” for Numpy arrays in the Hinc function like so:
class Hinc:
def __init__(self, b, kc, H):
self.beta = b
self.kc = kc
self.H = H
def eval(self, x):
rr = np.sqrt(x[0] * x[0] + x[1] * x[1])
s = np.where(rr > eps, x[1]/rr, 0.0)
c = np.where(rr > eps, x[0]/rr, 1.0)
hphi = np.where(rr > eps, 1j * self.beta * self.H * sp.jvp(1, rr, 0) * s / (self.kc * self.kc * rr), 0.0)
hr = -1j * self.beta * self.H * sp.jvp(1, rr, 1) * c / self.kc
hz = np.full_like(hr, 0.0+0j, dtype=np.complex128)
hx = hr * c - hphi * s
hy = hr * s + hphi * c
return(hx, hy, hz)
The code now runs fine! Cheers!