Hello again
After some much apprciated assistance on this discussion, I have managed to get a Microstrip line eigenvalue solver.
Ignoring all the hard stuff ahead of me (calculating H-fields, impedances, etc.), I rushed to see if I can also solve for a coaxial mode, as I had singularity issues with it.
After running the following code:
import sys
from mpi4py import MPI
import numpy as np
import dolfinx
print(f"DOLFINx version: {dolfinx.__version__} based on GIT commit: {dolfinx.git_commit_hash} of https://github.com/FEniCS/dolfinx/")
import scipy
import ufl
from basix.ufl import element, mixed_element
from dolfinx import default_scalar_type, fem, io, plot
from dolfinx.fem.petsc import assemble_matrix
from dolfinx.io import gmshio
try:
import pyvista
have_pyvista = True
except ModuleNotFoundError:
print("pyvista and pyvistaqt are required to visualise the solution")
have_pyvista = False
try:
from slepc4py import SLEPc
except ModuleNotFoundError:
print("slepc4py is required for this demo")
sys.exit(0)
f0 = 2e9
er_teflon = 2.5
er_air = 1.0
numSolutions = 5
eigenMode_num = 0
vector_basisFunction_order = 2
scalar_basisFunction_order = 3
msh, cell_tags, facet_tags = gmshio.read_from_msh( "COAX.msh", MPI.COMM_WORLD, gdim=2)
print("Import Gud!! Number of dimensions: {}".format(msh.topology.dim))
# Scalar and vector function spaces
FS_V = element("N1curl", msh.basix_cell(), vector_basisFunction_order)
FS_S = element("Lagrange", msh.basix_cell(), scalar_basisFunction_order)
combined_space = fem.functionspace(msh, mixed_element([FS_V, FS_S]))
# Define trial and test functions
Et_bf, Ez_bf = ufl.TrialFunctions(combined_space)
Et_tf, Ez_tf = ufl.TestFunctions(combined_space)
m0 = scipy.constants.mu_0
e0 = scipy.constants.epsilon_0
c0 = scipy.constants.speed_of_light
k0 = 2.0*np.pi*f0/c0
# Find tags for relevant facets and cells
facets_gnd = facet_tags.find(1)
facets_line = facet_tags.find(2)
cells_teflon = facet_tags.find(3)
# Concatenate metals
facet_metals = np.concatenate((facets_gnd, facets_line))
# Set a function space for the epsilon values
DG0 = fem.FunctionSpace(msh, ("DG", 0))
# Function for relative epsilon
e_r = fem.Function(DG0)
m_r = fem.Function(DG0)
# The 'x' here does not mean only x. This is the class that returns
# the "Function degree-of-freedom vector".
e_r.x.array[cells_teflon] = np.full_like(cells_teflon, er_teflon, dtype=default_scalar_type)
m_r.x.array[cells_teflon] = np.full_like(cells_teflon, 1.0, dtype=default_scalar_type)
# Formulate FEM problem
a_tt = ((1/m_r)*ufl.inner(ufl.curl(Et_bf), ufl.curl(Et_tf)) - (k0**2) * e_r * ufl.inner(Et_bf, Et_tf)) * ufl.dx
b_tt = (1/m_r)*ufl.inner(Et_bf, Et_tf) * ufl.dx
b_tz = (1/m_r)*ufl.inner(Et_bf, ufl.grad(Ez_tf)) * ufl.dx
b_zt = ufl.inner(ufl.grad(Ez_bf), Et_tf) * ufl.dx
b_zz = ((1/m_r)*ufl.inner(ufl.grad(Ez_bf), ufl.grad(Ez_tf)) - (k0**2) * e_r * ufl.inner(Ez_bf, Ez_tf)) * ufl.dx
Aform = fem.form(a_tt)
Bform = fem.form(b_tt + b_tz + b_zt + b_zz)
# Add boundary conditions
# Original:
# bc_facets = exterior_facet_indices(msh.topology)
# bc_dofs = fem.locate_dofs_topological(V, msh.topology.dim - 1, bc_facets)
# u_bc = fem.Function(V)
# with u_bc.vector.localForm() as loc:
# loc.set(0)
# bc = fem.dirichletbc(u_bc, bc_dofs)
bc_dofs = fem.locate_dofs_topological(combined_space, msh.topology.dim - 1, facet_metals)
PEC_bc = fem.Function(combined_space)
with PEC_bc.vector.localForm() as loc:
loc.set(0)
bc = fem.dirichletbc(PEC_bc, bc_dofs)
# Now we can solve the problem with SLEPc. First of all, we need to
# assemble our and matrices with PETSc in this way:
A = assemble_matrix(Aform, bcs=[bc])
A.assemble()
B = assemble_matrix(Bform, bcs=[bc])
B.assemble()
# Now, we need to create the eigenvalue problem in SLEPc. Our problem is a
# linear eigenvalue problem, that in SLEPc is solved with the EPS module.
# We can initialize this solver in the following way:
eps = SLEPc.EPS().create(msh.comm)
# We can pass to EPS our matrices by using the setOperators routine:
eps.setOperators(A, B)
# Next, we need to specify a tolerance for the iterative solver,
# so that it knows when to stop:
tol = 1e-9
eps.setTolerances(tol=tol)
# Set solver type
eps.setType(SLEPc.EPS.Type.KRYLOVSCHUR)
# Get ST context from eps
st = eps.getST()
# Set shift-and-invert transformation
st.setType(SLEPc.ST.Type.SINVERT)
# Target only real values
# eps.setWhichEigenpairs(SLEPc.EPS.Which.TARGET_REAL)
# Target value:
bt2 = k0*k0*er_teflon
eps.setTarget(-bt2)
# Number of eigenvalues
eps.setDimensions(nev = numSolutions)
eps.solve()
eps.view()
eps.errorView()
################################
################################
bt2_vals = [(i, -eps.getEigenvalue(i)) for i in range(eps.getConverged())]
kz2_vect = []
for eigIdx in range(len(bt2_vals)):
kz2_vect.append(bt2_vals[eigIdx][1])
print(bt2_vals[eigIdx])
# Start visualizing the solution
eh = fem.Function(combined_space)
# Save eigenvector in eh
eps.getEigenpair(eigenMode_num, eh.vector)
kz = np.sqrt(kz2_vect[eigenMode_num])
print(f"eigenvalue: {-(kz**2)}")
print(f"kz: {kz}")
print(f"kz/k0: {kz / k0}")
eh.x.scatter_forward()
eth, ezh = eh.split()
eth = eh.sub(0).collapse()
ez = eh.sub(1).collapse()
# Transform eth, ezh into Et and Ez
eth.x.array[:] = eth.x.array[:] / kz
ezh.x.array[:] = ezh.x.array[:] * 1j
# Prepare container to write Et
gdim = msh.geometry.dim
# "Discontinuous Lagrange" or "DG", representing scalar discontinuous Lagrange finite elements (discontinuous piecewise polynomial functions);
OD_dg = fem.functionspace(msh, ("DQ", vector_basisFunction_order, (gdim,)))
Et_dg = fem.Function(OD_dg)
Et_dg.interpolate(eth)
# Save solutions
# with io.VTXWriter(msh.comm, f"sols/Et_{eigenMode_num}.bp", Et_dg) as fHdl:
# io.VTKFile(msh.comm, f"sols/Et_{eigenMode_num}.pos", "w").write(Et_dg)
with io.VTKFile(msh.comm, f"sols/Et_{eigenMode_num}.pvd", "w") as file:
file.write_mesh(msh)
file.write_function([Et_dg._cpp_object])
I got this error:
petsc4py.PETSc.Error: error code 95
[0] EPSSolve() at /usr/local/slepc/src/eps/interface/epssolve.c:147
[0] EPSSolve_KrylovSchur_Default() at /usr/local/slepc/src/eps/impls/krylov/krylovschur/krylovschur.c:259
[0] BVMatArnoldi() at /usr/local/slepc/src/sys/classes/bv/interface/bvkrylov.c:91
[0] BVOrthonormalizeColumn() at /usr/local/slepc/src/sys/classes/bv/interface/bvorthog.c:402
[0] BVOrthogonalizeGS() at /usr/local/slepc/src/sys/classes/bv/interface/bvorthog.c:177
[0] BVOrthogonalizeCGS1() at /usr/local/slepc/src/sys/classes/bv/interface/bvorthog.c:101
[0] BV_SquareRoot_Default() at /usr/local/slepc/include/slepc/private/bvimpl.h:383
[0] BV_SafeSqrt() at /usr/local/slepc/include/slepc/private/bvimpl.h:132
[0] Missing or incorrect user input
[0] The inner product is not well defined: indefinite matrix -nan.
[0]PETSC ERROR: --------------------- Error Message --------------------------------------------------------------
[0]PETSC ERROR: Operation done in wrong order
[0]PETSC ERROR: Need to call MatDenseRestoreSubMatrix() first
[0]PETSC ERROR: See https://petsc.org/release/faq/ for trouble shooting.
[0]PETSC ERROR: Petsc Release Version 3.20.0, unknown
[0]PETSC ERROR: Unknown Name on a linux-gnu-real64-32 named 20e406027a47 by Unknown Fri Apr 5 20:50:50 2024
[0]PETSC ERROR: Configure options PETSC_ARCH=linux-gnu-real64-32 --COPTFLAGS=-O2 --CXXOPTFLAGS=-O2 --FOPTFLAGS=-O2 --with-64-bit-indices=no --with-debugging=no --with-fortran-bindings=no --with-shared-libraries --download-hypre --download-metis --download-mumps --download-ptscotch --download-scalapack --download-spai --download-suitesparse --download-superlu --download-superlu_dist --with-scalar-type=real --with-precision=double
[0]PETSC ERROR: #1 MatDestroy_SeqDense() at /usr/local/petsc/src/mat/impls/dense/seq/dense.c:1695
[0]PETSC ERROR: #2 MatDestroy() at /usr/local/petsc/src/mat/interface/matrix.c:1401
[0]PETSC ERROR: #3 DSReset() at /usr/local/slepc/src/sys/classes/ds/interface/dsbasic.c:887
[0]PETSC ERROR: #4 DSDestroy() at /usr/local/slepc/src/sys/classes/ds/interface/dsbasic.c:910
[0]PETSC ERROR: #5 EPSDestroy() at /usr/local/slepc/src/eps/interface/epsbasic.c:343
[0]PETSC ERROR: #6 PetscObjectDestroy() at /usr/local/petsc/src/sys/objects/destroy.c:51
[0]PETSC ERROR: #7 PetscObjectDelayedDestroy() at /usr/local/petsc/src/sys/objects/garbage.c:78
Traceback (most recent call last):
File "petsc4py/PETSc/PETSc.pyx", line 79, in petsc4py.PETSc.CHKERR
petsc4py.PETSc.Error: error code 58
[0] EPSSolve() at /usr/local/slepc/src/eps/interface/epssolve.c:147
[0] EPSSolve_KrylovSchur_Default() at /usr/local/slepc/src/eps/impls/krylov/krylovschur/krylovschur.c:259
[0] BVMatArnoldi() at /usr/local/slepc/src/sys/classes/bv/interface/bvkrylov.c:91
[0] BVOrthonormalizeColumn() at /usr/local/slepc/src/sys/classes/bv/interface/bvorthog.c:402
[0] BVOrthogonalizeGS() at /usr/local/slepc/src/sys/classes/bv/interface/bvorthog.c:177
[0] BVOrthogonalizeCGS1() at /usr/local/slepc/src/sys/classes/bv/interface/bvorthog.c:101
[0] BV_SquareRoot_Default() at /usr/local/slepc/include/slepc/private/bvimpl.h:383
[0] BV_SafeSqrt() at /usr/local/slepc/include/slepc/private/bvimpl.h:132
[0] Missing or incorrect user input
[0] The inner product is not well defined: indefinite matrix -nan.
Exception ignored in: 'petsc4py.PETSc.Object.__dealloc__'
Traceback (most recent call last):
File "petsc4py/PETSc/PETSc.pyx", line 79, in petsc4py.PETSc.CHKERR
petsc4py.PETSc.Error: error code 58
[0] EPSSolve() at /usr/local/slepc/src/eps/interface/epssolve.c:147
[0] EPSSolve_KrylovSchur_Default() at /usr/local/slepc/src/eps/impls/krylov/krylovschur/krylovschur.c:259
[0] BVMatArnoldi() at /usr/local/slepc/src/sys/classes/bv/interface/bvkrylov.c:91
[0] BVOrthonormalizeColumn() at /usr/local/slepc/src/sys/classes/bv/interface/bvorthog.c:402
[0] BVOrthogonalizeGS() at /usr/local/slepc/src/sys/classes/bv/interface/bvorthog.c:177
[0] BVOrthogonalizeCGS1() at /usr/local/slepc/src/sys/classes/bv/interface/bvorthog.c:101
[0] BV_SquareRoot_Default() at /usr/local/slepc/include/slepc/private/bvimpl.h:383
[0] BV_SafeSqrt() at /usr/local/slepc/include/slepc/private/bvimpl.h:132
[0] Missing or incorrect user input
[0] The inner product is not well defined: indefinite matrix -nan.
I have read a bit in the SLEPc guides, and I am pretty sure I have this has to do with the spectral transformation method I used, but I have a lot more reading ahead of me, to solve this.
If you can advise regarding how to approach this, I would appreciate it, a lot.
The mesh file can be found here.