FEM Failures, Finale. Featuring Dolphins And Whales!

The final chapter is upon us. But will a solution be found? Stay tuned, to find out more.

We left off with a seemingly singular solution. Who knew, that the coaxial cable, of all… It should be the simplest one, no? I mean, you can solve it analytically…

So how do we proceed from here. I have pretty much exhausted the options using Sparselizards python interface. Maybe reverting to the C++ library will allow more options. Another idea I had, was to try a library that is still maintained. Maybe earn the support of the community there. For that purpose, I went with DOLFINx, a newer, python based derivative of FENicS

So this was great, for two reason. The first would be that it’s better to move in any direction when you are stuck on a problem. The second is that it was to learn something new. I thought it was only two things, but then I was exposed to the third. I found out about Docker.

Now I know all you experienced developers are hissing at your desks right now. Docker isn’t exactly news. But us, copper welders from the hardware realm, we don’t deal with this stuff usually. But if you are like myself, let me introduce this to you:

I’m sure most of you were acquainted to a virtual machine once in your life. Some resources, somewhere, are allocated to run an OS of your choice. This is like a computer, for all extents and pourposes, as it keeps running as long as it’s (virtual) power is on.
Docker is a slightly different beast. The docker container is kind of equivalent to your virtual machine, in the sense that it defines the OS that will operate. However that’s not all it has. It also has all the necessary libraries, files, drivers, etc., that you will need to start working, immediately.

Forget this subject, for a second. Consider how easy it would be if instead of bringing up for hours\days on a new project, you could just start a docker container and let it rip. Well, this idea also sounded awesome.

So pretty darn fast, I was able to run my first FENicS (or should I say DOLFINx) project. I still couldn’t get any display going (it goes through linking an X-server to windows, or something), so for now I will print the files. Has to be said, though. Although I had most of the template written down, at this point,

Ok, now it’s time to run the coaxial 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)

try:
    from petsc4py import PETSc
except ModuleNotFoundError:
    print("petsc4py is required for this demo")
    sys.exit(0)

####################
# Input parameters #
####################

f0 = 2e9

er_teflon = 2.5
er_air = 1.0

numSolutions = 2
eigenMode_num = 0

vector_basisFunction_order = 2
scalar_basisFunction_order = 3

# Import mesh file
msh, cell_tags, facet_tags = gmshio.read_from_msh( "COAX.msh", MPI.COMM_WORLD, gdim=2)

# My hello world
print("Import Good!! 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)

# Tags for dielectric fill
cells_teflon = facet_tags.find(3)

# Concatenate metals
facet_metals = np.concatenate((facets_gnd, facets_line))

# Set a function space for the epsilon and mu values (for interpolation)
DG0 = fem.FunctionSpace(msh, ("DG", 0))

# Function for relative epsilon
e_r  = fem.Function(DG0)
m_r = fem.Function(DG0)

# The 'x' here doesn't 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 dirichlet boundary conditions
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)

############################
# Regular old shift invert #
############################
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_MAGNITUDE)
eps.setDimensions(nev = numSolutions)

# Target value:
bt2 = k0*k0*er_teflon

# Target the magnitude of the eigenvalues
eps.setTarget(-bt2);
eps.setWhichEigenpairs(SLEPc.EPS.Which.TARGET_MAGNITUDE);

eps.solve()
eps.view()
eps.errorView()

Which yields, of course

  File "slepc4py/SLEPc/EPS.pyx", line 1266, in slepc4py.SLEPc.EPS.solve
  File "petsc4py/PETSc/PETSc.pyx", line 79, in petsc4py.PETSc.CHKERR
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 Thu Jul 25 07:24:14 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.

Well, an error. This is new! In Sparselizard there wasn’t any runtime error… An erroneous solution isn’t the same as an actual error. So I’m going to do something unfriendly. I’ll try to replicate Sparselizards eigensolve configuration into this code.

// Tell slepc how many eigs we want:
EPSSetDimensions(eps, numeigenvaluestocompute, PETSC_DECIDE, PETSC_DECIDE);
// Set tolerance and max num of iterations allowed:
EPSSetTolerances(eps, 1e-6, 100);
// Set the eigenvalue solver:
EPSSetType(eps, EPSKRYLOVSCHUR);

EPSSetFromOptions(eps);

// We use a shift and invert transform:
ST st;
EPSGetST(eps, &st);
STSetType(st, STSINVERT);

// We target the eigenvalues with a given magnitude:
EPSSetTarget(eps, targeteigenvaluemagnitude);
EPSSetWhichEigenpairs(eps, EPS_TARGET_MAGNITUDE);

// MUMPS petsc solver context:
KSP ksp;
STGetKSP(st, &ksp);
KSPSetType(ksp, "preonly");
PC pc;
KSPGetPC(ksp, &pc);
PCSetType(pc, PCLU);
PCFactorSetMatSolverType(pc, universe::solvertype);

And here is the python equivalent.

tol = 1e-3
bt2 = k0*k0*er_teflon  # Target value

eps.setTolerances(tol=tol, max_it=100);

# Set solver type
eps.setType(SLEPc.EPS.Type.KRYLOVSCHUR)
# eps.setType(SLEPc.EPS.Type.LAPACK)

eps.setFromOptions()

eps.setProblemType(SLEPc.EPS.ProblemType.GHIEP)
# eps.setProblemType(SLEPc.EPS.ProblemType.GNHEP)

# Get ST context from eps
st = eps.getST()

# Set shift-and-invert transformation
st.setType(SLEPc.ST.Type.SINVERT)

# Setup eigen-problem solver
eps.setTarget(-bt2);
eps.setWhichEigenpairs(SLEPc.EPS.Which.TARGET_MAGNITUDE);

# Directly edit the PETSc linear system solver
ksp = st.getKSP()
ksp.setType(PETSc.KSP.Type.PREONLY)

# Pre-conditioner setup
pc = ksp.getPC()
pc.setType(PETSc.PC.Type.LU)
pc.setFactorSolverType(PETSc.Mat.SolverType.MUMPS)

# Number of eigenvalues
eps.setDimensions(nev = numSolutions,ncv = PETSc.DECIDE, mpd = PETSc.DECIDE)

Now let’s run it again.

  File "slepc4py/SLEPc/EPS.pyx", line 1266, in slepc4py.SLEPc.EPS.solve
  File "petsc4py/PETSc/PETSc.pyx", line 79, in petsc4py.PETSc.CHKERR
petsc4py.PETSc.Error: error code 76
[0] EPSSolve() at /usr/local/slepc/src/eps/interface/epssolve.c:134
[0] EPSSetUp() at /usr/local/slepc/src/eps/interface/epssetup.c:381
[0] STSetUp() at /usr/local/slepc/src/sys/classes/st/interface/stsolve.c:524
[0] STSetUp_Sinvert() at /usr/local/slepc/src/sys/classes/st/impls/sinvert/sinvert.c:116
[0] KSPSetUp() at /usr/local/petsc/src/ksp/ksp/interface/itfunc.c:415
[0] PCSetUp() at /usr/local/petsc/src/ksp/pc/interface/precon.c:1068
[0] PCSetUp_LU() at /usr/local/petsc/src/ksp/pc/impls/factor/lu/lu.c:121
[0] MatLUFactorNumeric() at /usr/local/petsc/src/mat/interface/matrix.c:3252
[0] MatFactorNumeric_MUMPS() at /usr/local/petsc/src/mat/impls/aij/mpi/mumps/mumps.c:1964
[0] Error in external library
[0] MUMPS error in numerical factorization: INFOG(1)=-9, INFO(2)=196415 (see users manual https://mumps-solver.org/index.php?page=doc "Error and warning diagnostics")

Shorter answer, always better. There it is. This means that matrix is completely singular. Meaning, that even if some dirty trick was taken to stabilize the solution, it would probably result in a spurious solution. Starting to make sense now?

No happy endings. Not for me, not today. Not for a while, now, in this specific aspect, at least. I couldn’t solve this problem yet.
I have narrowed it down to a few options:

  • The mesh is repetative. This creates a matrix with similar rows. I tried making the mesh non-uniform, didn’t work.
  • This mode has no normal (Z) component at all. Maybe this can be resolved with the alternative formulation for this problem, for example: \mathbf{A}_{tt}\mathbf{e}_t = k_z^2 \mathbf{B}_{tt}^' \mathbf{e}_t, where $latex \mathbf{B}^’_{tt}=\mathbf{B}_{tz}\mathbf{B}_{zz}^{-1}\mathbf{B}_{zt}-\mathbf{B}_{tt}.
  • Different solver. Well, I guess a way to start would be to export the matrices and run some tests on Octave, but for some reason, I doubt this is the issue. Two different solvers, no solution.

however I was not succesful in implementing it, yet, in neither software.

Well, this is it, for now. I am looking for the help of experts much wiser than myself, or at least people lucky enough to have encountered the solution before.

I hope you enjoyed the journey, so far. I sure as hell learned quite a bit from this. As usual, the codes are in my GitHub.

See you all on the next adventure!

Leave a Reply

Your email address will not be published. Required fields are marked *