Are there any ready-made packages in fenicsx to help imply PML on the boundary of domain when doing simulations of electromagnetic field？
See for instance Electromagnetic scattering from a wire with perfectly matched layer condition — DOLFINx 0.9.0 documentation
I believe this is @bay_swiss speciality:)
2 Likes
Few years ago I implemented the Locally-Conformal PML, based on the work of Ozgun et al, that allows to start from an existing geometry and get an automatically extruded PML together with the stretch functions and tensor necessary to modify the Helmholtz equation. You will find the code here:
github.com/bayswiss/autogen_LC_PML
And a detailed explanation of how it works (the same that I gave at the FENICS23 presentation) here:
undabit.com/locally-conformal-perfectly-matched-layer-implementation-in-fenicsx
In the link above, you will find also some references, where the method was taken from.
While I used it in acoustic applications the \Lambda_{PML} tensor and the \alpha_{PML} functions do not depend on the field type. So you should be able to adapt my scripts to the Electromagnetic case.
2 Likes
Thanks for your help every much. I tried to download air.step,autogen_LC_PML.py,helmholtz_with_PML.py and ran helmholtz_with_PML.py. It seems that the code is not adoptable in fenicsx 0.9, for example VectorFunctionSpace and TensorFunctionSpace seem not available now.
However, after I changed them into a available form, python gave a feedback like this, which seems that complex number is not supported. 
Traceback (most recent call last):
File "/mnt/d/pypro/THz/helmholtz_with_PML.py", line 96, in <module>
problem = LinearProblem(a, L, u=uh,
File "/usr/lib/petsc/lib/python3/dist-packages/dolfinx/fem/petsc.py", line 782, in __init__
self._a = _create_form(
File "/usr/lib/petsc/lib/python3/dist-packages/dolfinx/fem/forms.py", line 337, in form
return _create_form(form)
File "/usr/lib/petsc/lib/python3/dist-packages/dolfinx/fem/forms.py", line 331, in _create_form
return _form(form)
File "/usr/lib/petsc/lib/python3/dist-packages/dolfinx/fem/forms.py", line 254, in _form
ufcx_form, module, code = jit.ffcx_jit(
File "/usr/lib/petsc/lib/python3/dist-packages/dolfinx/jit.py", line 62, in mpi_jit
return local_jit(*args, **kwargs)
File "/usr/lib/petsc/lib/python3/dist-packages/dolfinx/jit.py", line 212, in ffcx_jit
r = ffcx.codegeneration.jit.compile_forms([ufl_object], options=p_ffcx, **p_jit)
File "/usr/lib/python3/dist-packages/ffcx/codegeneration/jit.py", line 225, in compile_forms
raise e
File "/usr/lib/python3/dist-packages/ffcx/codegeneration/jit.py", line 205, in compile_forms
impl = _compile_objects(
File "/usr/lib/python3/dist-packages/ffcx/codegeneration/jit.py", line 330, in _compile_objects
_, code_body = ffcx.compiler.compile_ufl_objects(
File "/usr/lib/python3/dist-packages/ffcx/compiler.py", line 108, in compile_ufl_objects
analysis = analyze_ufl_objects(ufl_objects, options["scalar_type"]) # type: ignore
File "/usr/lib/python3/dist-packages/ffcx/analysis.py", line 94, in analyze_ufl_objects
form_data = tuple(_analyze_form(form, scalar_type) for form in forms)
File "/usr/lib/python3/dist-packages/ffcx/analysis.py", line 94, in <genexpr>
form_data = tuple(_analyze_form(form, scalar_type) for form in forms)
File "/usr/lib/python3/dist-packages/ffcx/analysis.py", line 180, in _analyze_form
form_data = ufl.algorithms.compute_form_data(
File "/usr/lib/python3/dist-packages/ufl/algorithms/compute_form_data.py", line 278, in compute_form_data
form = preprocess_form(form, complex_mode)
File "/usr/lib/python3/dist-packages/ufl/algorithms/compute_form_data.py", line 238, in preprocess_form
form = remove_complex_nodes(form)
File "/usr/lib/python3/dist-packages/ufl/algorithms/remove_complex_nodes.py", line 40, in remove_complex_nodes
return map_integrand_dags(ComplexNodeRemoval(), expr)
File "/usr/lib/python3/dist-packages/ufl/algorithms/map_integrands.py", line 86, in map_integrand_dags
return map_integrands(
File "/usr/lib/python3/dist-packages/ufl/algorithms/map_integrands.py", line 29, in map_integrands
mapped_integrals = [
File "/usr/lib/python3/dist-packages/ufl/algorithms/map_integrands.py", line 30, in <listcomp>
map_integrands(function, itg, only_integral_type) for itg in form.integrals()
File "/usr/lib/python3/dist-packages/ufl/algorithms/map_integrands.py", line 39, in map_integrands
return itg.reconstruct(function(itg.integrand()))
File "/usr/lib/python3/dist-packages/ufl/algorithms/map_integrands.py", line 87, in <lambda>
lambda expr: map_expr_dag(function, expr, compress), form, only_integral_type
File "/usr/lib/python3/dist-packages/ufl/corealg/map_dag.py", line 35, in map_expr_dag
(result,) = map_expr_dags(
File "/usr/lib/python3/dist-packages/ufl/corealg/map_dag.py", line 103, in map_expr_dags
r = handlers[v._ufl_typecode_](v, *[vcache[u] for u in v.ufl_operands])
File "/usr/lib/python3/dist-packages/ufl/algorithms/remove_complex_nodes.py", line 28, in terminal
raise ValueError("Unexpected complex value in real expression.")
ValueError: Unexpected complex value in real expression.
Exception ignored in: <function LinearProblem.__del__ at 0x7fcf406c9cf0>
Traceback (most recent call last):
File "/usr/lib/petsc/lib/python3/dist-packages/dolfinx/fem/petsc.py", line 831, in __del__
self._solver.destroy()
AttributeError: 'LinearProblem' object has no attribute '_solver'
This only requires minor adaptations (using dolfinx.fem.functionspace, and the appropriate element based on the syntax at: Introduction to the Unified Form Language — FEniCS Workshop)
This depends on how you install dolfinx and PETSC. As PETSC can’t support different floating types simultaneously one has to decide this at the time of installation. How did you install dolfinx?
1 Like
I use Ubuntu 22.04.5 LTS to install fenicsx. Are there any solutions?
There maybe some problems with PETSC, since I find that PETSC is a none project file in fenicsx 0.9, and PETSc.ScalarType also does not exist.
some code in autogen_LC_PML.py
import numpy as np
from dolfinx.fem import Function, functionspace, Constant
from dolfinx.io import XDMFFile, gmshio
from ufl import as_matrix, ln
from mpi4py import MPI
from petsc4py import PETSc
msh, cell_tags, facet_tags = gmshio.read_from_msh(mesh_name_prefix + "TOT.msh", MPI.COMM_SELF, 0, gdim=3)
V = functionspace(msh, ("CG", elem_degree))
VV = functionspace(msh, ("CG", elem_degree, (3, )))
indexes_dolfinx = msh.geometry.input_global_indices
omega = Constant(msh, PETSc.ScalarType(1))
k0 = Constant(msh, PETSc.ScalarType(1))
See the complex numbers paragraph at Download | FEniCS Project
Thanks for your kind help.
When doing this
/mnt/d/pypro/THz$ PETSC_DIR=/usr/lib/petscdir/petsc-complex python3 helmholtz_with_PML.py,
the page stop here and WSL occupy a lot of my CPU resources, is it ok for this condition?
Info : No ill-shaped tets in the mesh :-)
Info : 0.00 < quality < 0.10 : 0 elements
Info : 0.10 < quality < 0.20 : 0 elements
Info : 0.20 < quality < 0.30 : 1 elements
Info : 0.30 < quality < 0.40 : 698 elements
Info : 0.40 < quality < 0.50 : 978 elements
Info : 0.50 < quality < 0.60 : 1722 elements
Info : 0.60 < quality < 0.70 : 3986 elements
Info : 0.70 < quality < 0.80 : 8735 elements
Info : 0.80 < quality < 0.90 : 13295 elements
Info : 0.90 < quality < 1.00 : 6561 elements
Info : Done optimizing mesh (Wall 0.0632079s, CPU 0.063203s)
Info : 29305 nodes 188441 elements
Info : Writing 'airTOT.msh'...
Info : Done writing 'airTOT.msh'
Info : Reading 'airTOT.msh'...
Info : 327 entities
Info : 29281 nodes
Info : 173877 elements
Info : Done reading 'airTOT.msh'
Computing...
Computing frequency: 100 Hz...
Computing frequency: 300 Hz...
Computing frequency: 500 Hz...
Computing frequency: 700 Hz...
Computing frequency: 900 Hz...
Computing frequency: 1100 Hz...
Dear bay_swiss,
Thanks for your work!
Since I want to use PML to give a boundary condition to EM field and it seems too complex to understand all the details, I just want to know whether this code could fit my code well.
I have a domain and I could calculate E and H at any point in the domain, and in the model, E and H will propagate freely outside of domain. I wonder is there any easy way to attach the PML package into this model. Since I am a totally green hand, so I may ask someeasy questions. 
Can you tell something about the dimension (2D, 3D) and PML shape? I might have a working model that you could directly use.
Thanks for your help!
The dimension is 3D, and maybe finally I will need to shape of the core domain to be a cross shape(3D), but now I just set core domain as a unit cube. So I guess a square PML layer is ok.
I am sorry that I may not reply you quickly, since I am about to sleep.
Below is the model for a patch antenna which uses PMLs of cartesian geometry. I hope in near future I will get time to publish my exemplary models after necessary refinements and explanations. For now, I hope it will be useful as is.
I encourage you to do some sanity checks to ensure the PML is optimally tuned. It has been a while and I am unable to recall details.
#!/usr/bin/env python
# coding: utf-8
# # Patch antenna model
#
# This notebook reproduces in FEniCSx the simple patch antenna model from CST [tutorial](https://www.youtube.com/watch?v=EDkQV4JHR_o). The PML formulation is identical to the one used in FEniCSx [tutorial](https://docs.fenicsproject.org/dolfinx/v0.8.0/python/demos/demo_axis.html).
#
# Author: Shakeeb Bin Hasan
#
# Version: Fenicsx 0.9
import dolfinx, ufl, mpi4py, gmsh, petsc4py
from dolfinx.io import gmshio
import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import griddata
# ## CAD model
# In[2]:
def create_cad_model(fov_R=100, dpml=10, symmetry=False, refine=0, max_size=None):
'''Create the CAd_sub geometry using Gmsh. Default unit is mm'''
gmsh.initialize()
gmsh.option.setNumber("General.Terminal", 0) # disable output messages
gmsh.clear()
gdim = 3
occ = gmsh.model.occ
gmsh.merge('cad/patch_antenna_cart_pml.step')
gmsh.model.occ.synchronize()
all_doms = occ.getEntities(gdim)
frag, _ = occ.fragment([all_doms[0]], all_doms[1:])
occ.synchronize()
int_box = occ.addBox(-fov_R, 0, -fov_R, 2*fov_R, fov_R, 2*fov_R)
inters, _ = occ.intersect(frag, [(gdim, int_box)])
occ.synchronize()
# add_sub physical tags
all_doms = occ.getEntities(gdim)
for j, dom in enumerate(all_doms):
gmsh.model.addPhysicalGroup(dom[0], [dom[1]], j + 1) # create the main group/node
# number all boundaries
all_edges = gmsh.model.getEntities(gdim - 1)
for j, edge in enumerate(all_edges):
gmsh.model.addPhysicalGroup(edge[0], [edge[1]], edge[1]) # create the main group/node
# create mesh field_sub to have higher resolution near magnets
mpi_rank = 0
if max_size is not None:
gmsh.model.mesh.setSize(gmsh.model.getEntities(), max_size)
gmsh.option.setNumber("Mesh.ElementOrder", 2)
gmsh.model.mesh.generate(gdim)
# gmsh.model.mesh.optimize("Netgen")
for j in range(refine):
gmsh.model.mesh.refine()
mesh, ct, ft = gmshio.model_to_mesh(gmsh.model, mpi4py.MPI.COMM_WORLD, mpi_rank, gdim)
gmsh.write('cad/patch_antenna_cart_pml.msh')
return mesh, ct, ft
# ## Input parameters
# In[3]:
dom_idx = [(14,15,16,17), (8,13), # substrate, air,
(2,19,), (11,), (7,9), # pml_x, pml_y, pml_z,
(5,22), (10, 12,), (1,3,18,20), # pml_xy, pml_yz, pml_xz
(4,6,21,23)] # , pml_xyz
pml_dom_idx = dom_idx[2:]
sigma_vec = np.zeros(len(dom_idx))*0.0j
epsr_vec = np.ones(len(dom_idx))+0.0j
epsr_vec[0] = 4.3
mur_vec = np.ones(len(dom_idx))+0.0j
geom_pec_bnd_idx = (85,176,177,179,180,183)
# scat_bnd_idx = 26
port_bnd_idx = 178 # lumped_sub port excitation
aux_bnd_idx = (65,66,68,70)
ff_dom_idx = 13
ff_bnd_idx = aux_bnd_idx
freq = 2.4e9
fem_order = 3
use_pml = True
fov_R = 100
dpml = 10
mesh, ct, ft = create_cad_model(refine=0, max_size=3e8/freq/10*1e3)
# ## Automatic detection of PML domains
# In[4]:
fov = 200
dpml = 10
gdim = 3
eps = 1e-3
all_doms = gmsh.model.getEntities(gdim)
doms_in_bb = gmsh.model.getEntitiesInBoundingBox(-fov/2+eps, -fov/2+eps, -fov/2+eps, fov/2-eps, fov/2-eps, fov/2-eps, gdim)
pml_doms = tuple(set(all_doms) - set(doms_in_bb))
# pml_doms = np.sort(pml_doms, axis=1)
pml_doms = np.sort([(j[1]) for j in pml_doms])
pml_doms = [(gdim, j) for j in pml_doms]
# print("PML domains are: ", pml_doms)
# pml along x
pml_nx = gmsh.model.getEntitiesInBoundingBox(-fov/2-eps, -fov/2+dpml-eps, -fov/2+dpml-eps, -fov/2+dpml+eps, fov/2-dpml+eps, fov/2-dpml+eps, gdim)
pml_px = gmsh.model.getEntitiesInBoundingBox(fov/2-dpml-eps, -fov/2+dpml-eps, -fov/2+dpml-eps, fov/2+eps, fov/2-dpml+eps, fov/2-dpml+eps, gdim)
print(pml_nx, pml_px)
# pml along y
pml_ny = gmsh.model.getEntitiesInBoundingBox(-fov/2+dpml-eps, -fov/2-eps, -fov/2+dpml-eps, fov/2-dpml+eps, -fov/2+dpml+eps, fov/2-dpml+eps, gdim)
pml_py = gmsh.model.getEntitiesInBoundingBox(-fov/2+dpml-eps, fov/2-dpml-eps, -fov/2+dpml-eps, fov/2-dpml+eps, fov/2+eps, fov/2-dpml+eps, gdim)
print(pml_ny, pml_py)
# pml along z
pml_nz = gmsh.model.getEntitiesInBoundingBox(-fov/2+dpml-eps, -fov/2+dpml-eps, -fov/2-eps, fov/2-dpml+eps, fov/2-dpml+eps, -fov/2+dpml+eps, gdim)
pml_pz = gmsh.model.getEntitiesInBoundingBox(-fov/2+dpml-eps, -fov/2+dpml-eps, fov/2-dpml-eps, fov/2-dpml+eps, fov/2-dpml+eps, fov/2+eps, gdim)
print(pml_nz, pml_pz)
# ## Domain variables
# In[5]:
Omega_dg = dolfinx.fem.functionspace(mesh, ("DG", 0))
epsr_param = dolfinx.fem.Function(Omega_dg)
mur_param = dolfinx.fem.Function(Omega_dg)
sigma_param = dolfinx.fem.Function(Omega_dg)
for epsr, mur, sigma, dom in zip(epsr_vec, mur_vec, sigma_vec, dom_idx):
for j in dom:
cells = ct.find(j)
epsr_param.x.array[cells] = epsr
mur_param.x.array[cells] = mur
sigma_param.x.array[cells] = sigma
# ## Weak form
#
# The strong form of the complex vector Helmholtz equation is
#
# $$
# \nabla\times\mu_\mathrm{r}^{-1}\nabla\times\mathbf{E} - k_0^2\left(\frac{\mathrm{i}\sigma}{\omega\varepsilon_0}+\varepsilon_\mathrm{r}\right)\mathbf{E} = \mathrm{i}k_0 Z_0\mathbf{J}.
# $$
#
# Mutliplying by the test function $\mathbf{V}$ and integrating gives
#
# $$
# \intop_\Omega \mathbf{dr} \mu_\mathrm{r}^{-1} \left(\nabla\times\mathbf{E}\right) . \left(\nabla\times\mathbf{V}\right) - k_0^2\intop_\Omega \mathbf{dr}\left(\frac{\mathrm{i}\sigma}{\omega\varepsilon_0}+\varepsilon_\mathrm{r}\right)\mathbf{E}.\mathbf{V} + \intop_{\partial\Omega}\mathrm{ds}\hat{n}.\left[\left(\mu_\mathrm{r}^{-1}\nabla\times\mathbf{E}\right)\times\mathbf{V}\right] = \mathrm{i}k_0 Z_0\intop_\Omega \mathbf{dr}\mathbf{J}.\mathbf{V}.
# $$
#
# For convenience, let us eliminate the frequency by writing in terms of impedance $Z_0$ and wavenumber $k_0$ as
#
# $$
# \intop_\Omega \mathbf{dr} \mu_\mathrm{r}^{-1} \left(\nabla\times\mathbf{E}\right) . \left(\nabla\times\mathbf{V}\right) - k_0^2\intop_\Omega \mathbf{dr}\left(\frac{\mathrm{i}Z_0\sigma}{k_0}+\varepsilon_\mathrm{r}\right)\mathbf{E}.\mathbf{V} + \intop_{\partial\Omega}\mathrm{ds}\hat{n}.\left[\left(\mu_\mathrm{r}^{-1}\nabla\times\mathbf{E}\right)\times\mathbf{V}\right] = \mathrm{i}k_0 Z_0\intop_\Omega \mathbf{dr}\mathbf{J}.\mathbf{V},
# $$
#
# $$
# \intop_\Omega \mathbf{dr} \mu_\mathrm{r}^{-1} \left(\nabla\times\mathbf{E}\right) . \left(\nabla\times\mathbf{V}\right) - k_0^2\intop_\Omega \mathbf{dr}\left(\frac{\mathrm{i}Z_0\sigma}{k_0}+\varepsilon_\mathrm{r}\right)\mathbf{E}.\mathbf{V} + \intop_{\partial\Omega}\mathrm{ds}\left[\hat{n}\times\left(\mu_\mathrm{r}^{-1}\nabla\times\mathbf{E}\right)\right].\mathbf{V} = \mathrm{i}k_0 Z_0\intop_\Omega \mathbf{dr}\mathbf{J}.\mathbf{V}.
# $$
# In[6]:
Omega = dolfinx.fem.functionspace(mesh, ("N1curl", fem_order))
c0_const = 2.9979e8*1e3 # m/s
eps0_const = 8.845e-12*1e3 # in F/m
mu0_const = 4*np.pi*1e-7*1e3 # in H/m
Z0_const = 376.73 # free-space impedance, Ohm
lbda = c0_const/freq
k0_param = dolfinx.fem.Constant(mesh, 2*np.pi/lbda + 0.0j)
w_freq = k0_param*c0_const
omega_param = k0_param*c0_const # angular frequency
U = ufl.TrialFunction(Omega)
Ut = ufl.TestFunction(Omega)
fs = dolfinx.fem.Constant(mesh, mesh.geometry.dim*(0.0+0.0j,)) # dummy source
dx = ufl.Measure('dx', domain=mesh, subdomain_data=ct, metadata={'quadrature_degree':8})
ds = ufl.Measure('ds', domain=mesh, subdomain_data=ft, metadata={'quadrature_degree':8})
dS = ufl.Measure('dS', domain=mesh, subdomain_data=ft, metadata={'quadrature_degree':8})
F = ufl.inner(1/mur_param * ufl.curl(U), ufl.curl(Ut))*dx(tuple(np.hstack(dom_idx[:2]))) - \
k0_param**2*(1j*sigma_param*Z0_const/(k0_param) + epsr_param)*ufl.inner(U, Ut)*dx(tuple(np.hstack(dom_idx[:2])))
F += -ufl.inner(fs, Ut)*dx
# ## Uniaxial PML
#
# We incorporate the PML by modifying the petmirrivity and permeability as follows
#
# $$
# \overline{\overline{S}}=\left(\begin{array}{ccc}
# \frac{s_{y}s_{z}}{s_{x}} & 0 & 0\\
# 0 & \frac{s_{x}s_{z}}{s_{y}} & 0\\
# 0 & 0 & \frac{s_{x}s_{y}}{s_{z}}
# \end{array}\right),
# $$
#
# $$
# \overline{\overline{\varepsilon}}_{\mathrm{pml}}=\overline{\overline{\varepsilon}}_{\mathrm{r}}\overline{\overline{S}},
# $$
#
# $$
# \overline{\overline{\mu}}_{\mathrm{pml}}=\overline{\overline{\mu}}_{\mathrm{r}}\overline{\overline{S}}.
# $$
# In[7]:
def upml_epsilon_mu(epsr, mur, pml_x=True, pml_y=True, pml_z=True, n_pml=2, alpha=1, R_pml=1e-3):
'''create material parameters for uniaxial PML'''
x_ufl, y_ufl, z_ufl = ufl.SpatialCoordinate(mesh)
if pml_x:
sx = 1 - 1j*alpha/k0_param*(n_pml+1)*np.log(R_pml)/(2*dpml)*((ufl.algebra.Abs(x_ufl)-x0)/dpml)**n_pml
else:
sx = 1
if pml_y:
sy = 1 - 1j*alpha/k0_param*(n_pml+1)*np.log(R_pml)/(2*dpml)*((ufl.algebra.Abs(y_ufl)-y0)/dpml)**n_pml
else:
sy = 1
if pml_z:
sz = 1 - 1j*alpha/k0_param*(n_pml+1)*np.log(R_pml)/(2*dpml)*((ufl.algebra.Abs(z_ufl)-z0)/dpml)**n_pml
else:
sz = 1
S_pml = ufl.as_matrix([[sy*sz/sx, 0, 0], [0, sx*sz/sy, 0], [0, 0, sx*sy/sz]])
epsr_pml = epsr*S_pml
mur_pml = mur*S_pml
return epsr_pml, mur_pml
# In[8]:
if use_pml:
n_pml = 2
alpha_pml = 1
x0 = fov_R - dpml
y0 = fov_R - dpml
z0 = fov_R - dpml
pml_dir = (
(True, False, False), (False, True, False), (False, False, True), # x, y, z
(True, True, False), (False, True, True), (True, False, True), # xy, yz, xz
(True, True, True), # xyz
)
for idx, dir in enumerate(pml_dir):
dx_ = dx(pml_dom_idx[idx])
epsr_pml, mur_pml = upml_epsilon_mu(epsr_param, mur_param, *dir, n_pml=n_pml, alpha=alpha_pml)
F_pml = ufl.inner(ufl.inv(mur_pml)*(ufl.curl(U)), ufl.curl(Ut))*dx_ - \
k0_param**2*ufl.inner(epsr_pml*U, Ut)*dx_
F += F_pml
else:
F_pml = (ufl.inner(ufl.inv(mur_param)*(ufl.curl(U)), ufl.curl(Ut)) - \
k0_param**2*epsr_param*ufl.inner(U, Ut))*dx(tuple(np.hstack(pml_dom_idx)))
F += F_pml
# ## PEC boundaries
#
# The geometry is made of PEC which is imposed strongly.
# In[ ]:
U_dbc = dolfinx.fem.Function(Omega)
U_dbc.x.array[:]= 0.0j
dbc = []
pec_facets = np.hstack([ft.find(j) for j in geom_pec_bnd_idx])
pec_dofs = dolfinx.fem.locate_dofs_topological(Omega, mesh.geometry.dim-1, pec_facets)
dbc.append(dolfinx.fem.dirichletbc(U_dbc, pec_dofs))
# ## Current source for the lumped port
#
# $$
# \mathbf{J}_\mathrm{s} = \frac{V Y}{h}(E_z \hat{a}_z - \mathbf{E}_\mathrm{p}),
# $$
#
# where $\mathbf{E}_\mathrm{p}$ is the electric field sustained by the port.
#
#
# In[10]:
import dolfinx.fem.function
Z0 = 376.73 # impedance of free space, ohms
Z_ref = 50 # circuit element's impedance
V0 = 1 # applied voltage
w_line = 0.74
d_sub = 1.5
eta_ref = Z_ref*w_line/d_sub # impedance of the guided mode
E0_phase = 0
E0 = ufl.as_vector((0, 0, -V0/d_sub))
J_lp = (U - 2*E0)/eta_ref
F_lp = 1j*k0_param*Z0*ufl.avg(ufl.inner(J_lp, Ut))*dS(port_bnd_idx)
F += -F_lp
# ## Solver
# In[12]:
U_sol = dolfinx.fem.Function(Omega)
print("number of dofs:",U_sol.x.array.size)
# In[13]:
# rtol = 1e-3
petsc_options = {"ksp_type": "preonly", "pc_type": "lu", "pc_factor_mat_solver_type":"mumps"}
solver_idx = 2
from dolfinx.fem.petsc import LinearProblem
Lin_system = LinearProblem(ufl.lhs(F), ufl.rhs(F), bcs=dbc, u=U_sol, petsc_options=petsc_options)
solver = Lin_system.solver
import time
t1 = time.perf_counter()
Lin_system.solve()
t2 = time.perf_counter()
print("Solved in {}s".format(t2-t1))
it = solver.getIterationNumber()
print("Constrained_sub solver iterations {0:d}".format(it))
# solver.view()
# ## Visualization
# In[14]:
Omega_sol = dolfinx.fem.functionspace(mesh, ("DG", fem_order, (mesh.geometry.dim,)))
U_sol_dg = dolfinx.fem.Function(Omega_sol)
U_sol_dg.interpolate(U_sol)
with dolfinx.io.VTXWriter(mpi4py.MPI.COMM_WORLD, "results/E_antenna.bp", [U_sol_dg], engine="BP4") as vtx:
vtx.write(0.0)
I have already set the environment variable PETSC_DIR to point at the required PETSc version located under /usr/lib/petscdir
export PETSC_DIR=/usr/lib/petscdir/petsc-complex
echo $PETSC_DIR
/usr/lib/petscdir/petsc-complex
However, after I run this in Ubuntu
/mnt/d/pypro/THz$ PETSC_DIR=/usr/lib/petscdir/petsc-complex python3 helmholtz_with_PML.py
It seems that the program starts to run directly in Ubuntu. And pycharm’s return still show that fenicsx doesn’t support complex(I use fenicsx as a package of python). Could you please give some advices?
You need to set PYTHONPATH and PETSC_DIR within pycharm.
Thanks for your help!
Now I set PETSC_DIR=/usr/lib/petscdir/petsc-complex，
but what should I set PYTHONPATH for?
Try with just the PETSC dir first, see if that helps.
It still couldn’t recognize dolfinx
Traceback (most recent call last):
File "/mnt/d/pypro/THz/THz.py", line 3, in <module>
import dolfinx.fem.petsc
ModuleNotFoundError: No module named 'dolfinx'
When you run python through terminal, where dolfinx is working, could you list the following:
which python3
python3 -c "import dolfinx; print(dolfinx)"
Yeah, I could run it.
fenics@LAPTOP-IBUP1ARC:~$ which python3
python3 -c "import dolfinx; print(dolfinx)"
/usr/bin/python3
<module 'dolfinx' from '/usr/lib/petscdir/petsc-complex/lib/python3/dist-packages/dolfinx/__init__.py'>