Defining internal boundaries as PEC

General
1

Hi,

I am solving vector Helmholtz equation in 3D and need perfect electric conductors (PEC) at both external and internal boundaries in my geometry. For external boundaries, I am able to realize it with Nitche’s method. However, I am unsure about internal boundaries due to field discontinuities. For now, I just add Nitche’s term on both sides of the edges to the variational form.

Nonetheless, I check my approach with a simple toy model. I divide a vacuum box into two halves separated by a PEC and launch a plane wave in only one. If the internal PEC works reliably, half of the volume should be dark at noise levels. But I get observable amplitudes in the supposedly dark box too. Below is the paraview output:

image
image1370×376 41.6 KB

I will appreciate if knowledgeable users could help me understand how the Nitche’s approach is modified for internal facets or point me to literature where I might find a solution. I have the book of Jian-Ming Jin but it seems silent on this.

Below is the simple MWE which produced the result shown above. To avoid mesh tags from changing (hopefully), the gmsh mesh file can be downloaded from here.

import dolfinx, petsc4py, mpi4py, gmsh, ufl
from dolfinx.io import gmshio
import numpy as np
from dolfinx.fem.petsc import LinearProblem
import time

#%% Load mesh
gmsh.initialize()
gmsh.option.setNumber("General.Terminal", 0) # disable output messages
gmsh.clear()
gmsh.merge('mwe.msh')
model_rank = 0
gdim = 3
mesh, ct, ft = gmshio.model_to_mesh(gmsh.model, mpi4py.MPI.COMM_WORLD, model_rank, gdim)

#%% input paramters
port_bnd_idx = (11, 10) # incident, transmitted
pec_int_bnd_idx = (2, )
pec_ext_bnd_idx = (9, 4)
k0 = 2*np.pi/0.5
fe_order = 2

#%% Variational form
Omega = dolfinx.fem.functionspace(mesh, ("N1curl", fe_order))
U = ufl.TrialFunction(Omega)
Ut = ufl.TestFunction(Omega)

F = (ufl.inner(ufl.curl(U), ufl.curl(Ut)) - k0**2*ufl.inner(U, Ut))*ufl.dx
F += -ufl.inner(dolfinx.fem.Constant(mesh, petsc4py.PETSc.ScalarType((0.0, 0.0,0.0))), Ut)*ufl.dx # source term

#%% Port boundary at the top (excitation) and bottom surfaces 
ds = ufl.Measure("ds", domain=mesh, subdomain_data=ft) # boundary differential element

n_3d = ufl.FacetNormal(mesh)

ki_3d = ufl.as_vector((0, 0, -k0)) # incident wave
kr_3d = ufl.as_vector((0, 0, k0)) # reflected wave, if any
kt_3d = ufl.as_vector((0, 0, -k0)) # transmitted wave

# port boundary at the incident port
# reflected wave
F_bnd_term = ufl.inner(ufl.cross(n_3d, ufl.cross(1j*kr_3d, U)), Ut)*ds(port_bnd_idx[0])
F += F_bnd_term
# excitation
E_inc = ufl.as_vector((0, 1, 0)) # y-polarized
F_bnd_term_exc = ufl.inner(ufl.cross(n_3d, ufl.cross(1j*(ki_3d-kr_3d), E_inc)), Ut)*ds(port_bnd_idx[0])
F += F_bnd_term_exc

# port boundary at the transmitted port
F_bc_tr = ufl.inner(ufl.cross(n_3d, ufl.cross(1j*kt_3d, U)), Ut)*ds(port_bnd_idx[1])
F += F_bc_tr

#%% PEC boundaries using Nitche's method
# hE = 2 * ufl.Circumradius(mesh)
hE = ufl.CellDiameter(mesh)
alpha = dolfinx.fem.Constant(mesh, petsc4py.PETSc.ScalarType(10))
U_pec = ufl.as_vector((U[0],0,U[2])) # Non-zero components will be forced zero

F_pec_ext = ufl.inner(ufl.cross(n_3d, ufl.curl(U)), Ut)*ds(pec_ext_bnd_idx) + alpha/hE*ufl.inner(U_pec, Ut)*ds(pec_ext_bnd_idx) - \
ufl.inner(U_pec, ufl.cross(n_3d, ufl.curl(Ut)))*ds(pec_ext_bnd_idx)

dS = ufl.Measure("dS", domain=mesh, subdomain_data=ft) # boundary differential element

F_pec_int_pos = ufl.inner(ufl.cross(n_3d, ufl.curl(U))('+'), Ut('+'))*dS(pec_int_bnd_idx) + alpha/hE('+')*ufl.inner(U_pec('+'), Ut('+'))*dS(pec_int_bnd_idx) - \
ufl.inner(U_pec('+'), ufl.cross(n_3d, ufl.curl(Ut))('+'))*dS(pec_int_bnd_idx)

F_pec_int_neg = ufl.inner(ufl.cross(n_3d, ufl.curl(U))('-'), Ut('-'))*dS(pec_int_bnd_idx) + alpha/hE('-')*ufl.inner(U_pec('-'), Ut('-'))*dS(pec_int_bnd_idx) - \
ufl.inner(U_pec('-'), ufl.cross(n_3d, ufl.curl(Ut))('-'))*dS(pec_int_bnd_idx)

F += F_pec_ext
F += F_pec_int_pos
F += F_pec_int_neg

a = ufl.lhs(F)
L = ufl.rhs(F)

#%% Solve

U_sol = dolfinx.fem.Function(Omega)

petsc_options = {"ksp_type": "preonly", "pc_type": "lu", "pc_factor_mat_solver_type":"mumps"}

Lin_system = LinearProblem(a, L, bcs = [], u=U_sol, petsc_options=petsc_options)

solver = Lin_system.solver
t1 = time.perf_counter()
Lin_system.solve()
t2 = time.perf_counter()
print("Solved in {}s".format(t2-t1))

Omega_sol = dolfinx.fem.functionspace(mesh, ("DG", fe_order, (mesh.topology.dim,)))
U_sol_dg = dolfinx.fem.Function(Omega_sol)
U_sol_dg.interpolate(U_sol)

with dolfinx.io.VTXWriter(mpi4py.MPI.COMM_WORLD, "E_field.bp", [U_sol_dg], engine="BP4") as vtx:
    vtx.write(0.0)
2

I am not an expert on Nitsche’s Method, but having used weak-form boundary condition implementations in the past, I know that they are not “perfect”. Since you are applying a boundary condition in a least-squares sense (or using some other norm), the boundary condition will exhibit some leakage. Using the traditional Dirichlet approach (that removes DoFs on the boundary), explicitly constrains the FE space to obey the boundary conditions (in your case, zero) and you would not see any leakage.

1 Like
3

I understand the note about perfection but I still had better hopes from the imperfect.

As far as I know, the implementation of the edge elements in the current dolfinx only allows setting the entire field to zero through Dirichlet condition. Can we do it to specific components too? Thanks this to this impression (right or wrong), I did not even try it before.

4

Turns out I was indeed wrong in my understanding.

Strongly enforcing all field components to zero only causes \hat{n}\times\mathbf{E}=0 at the boundary. Now the model works to my satisfaction.

Thanks a lot for your comment, it triggered my thoughts in ways perhaps neither of us would have anticipated.