Find TE and TM resonances of a axisymmetric cavity of arbitrary shape in rz (Minimal working example)

mwe-resonances-TE-TM-vacuum.py

#!/usr/bin/python2
# -*- coding: utf8 -*-
"""
Usage: %(scriptName)s

FEniCS RF cylindrical cavity eigenvalue problem test.
Computes the TE and TM eigenmodes of a cavity using the E field.
No losses are considered. Permittivity is currently not considered.

TODO: Implement induced surface current integral on boundary to estimate Q factor

(c) 2015 by C. Weickhmann <christian.weickhmann@mailbox.org>
Special thanks to Johann Heller
"""

from dolfin import *
import numpy as np
import scipy.constants as co

def convertMesh(geofile, dim=3, noConvert=False):
	"""Convert gmsh geo-file and return strings of newly created files.

	Keyword arguments:
	geofile -- filename of the geo file
	dim     -- enforce dimension (default: 3)

	Returns: Tuple of
	* xml dolfin meshfile
	* xml dolfin pysical regions markers (*)
	* xml dolfin facet markers (*)

	The files marked (*) are only created if physical features are used in the gmsh geo-file.
	The code does not check whether they are actually created.
	"""
	if geofile[-3:] != "geo":
		print("Error: Not a geo file.")
		return -1
	else:
		geofile = geofile[:-4]

	if not noConvert:
		os.system("gmsh %s.geo -%d -smooth 3" % (geofile, dim))
		os.system("dolfin-convert %s.msh %s.xml" % (geofile, geofile))
	return (("%s.xml" % geofile), ("%s_physical_region.xml" % geofile), ("%s_facet_region.xml" % geofile))


parameters["linear_algebra_backend"] = "PETSc" # Force this backend
parameters["form_compiler"]["representation"] = "quadrature"
parameters["form_compiler"]["cpp_optimize"] = True
parameters["form_compiler"]["cpp_optimize_flags"] = '-O3 -funroll-loops'
parameters["form_compiler"]["optimize"] = True
parameters['form_compiler'].add('eliminate_zeros', False)
set_log_level(CRITICAL)

# Shorthand names for physical constants from scipy.constants
mu0  = co.mu_0
eps0 = co.epsilon_0
c0   = co.speed_of_light
eta0 = co.physical_constants['characteristic impedance of vacuum'][0]

# Function definitions for the UFL expressions
def cross_1D(n, A):
	return as_vector([-n[1]*A, n[0]*A])

def cross_2D(x, y):
	return x[1]*y[0]-x[0]*y[1]

def n_x_n_x_A(n, A):
	# 1.) calculate phi = n x A, w/ n and A in rz
	phi = cross_2D(n, A)
	# 2.) calculate rz = n x (n x A), w/ (n x A) in phi and n in rz
	rz_r = -n[1]*phi
	rz_z = n[0]*phi
	return as_vector([rz_r, rz_z])

def inv(eps):
	return 1.0/eps

def curl_p(w_rz):
	c_phi = w_rz[0].dx(1)-w_rz[1].dx(0)
	return c_phi

def curl_rz(w_rz, w_phi, r, k):
	# k is order in phi
	# i is complex number, however for now set i=1
	i = 1
	c_r = i*k/r * w_rz[1] - w_phi.dx(1)
	c_z = 1/r * w_phi + w_phi.dx(0) - i*k/r * w_rz[0]
	return as_vector([c_r, c_z])

# Target frequencies: The mode is obtained by determination of the eigenvalue spectrum.
# Since each mode is identified by its eigenfrequency, the eigenvalue spectrum can be shifted
# to the approximate value and the solution is obtained for a small number of eigenvalues
# close to the target value only. This speeds up the calculation.
target_frequencies = np.array([25.00])*1e9

(geo, regions, borders) = convertMesh("resonator-30ghz.geo", dim=2, noConvert=False)
mesh = Mesh(geo)
conductive = MeshFunction('size_t', mesh, borders)

V1D_phi = FunctionSpace(mesh, "CG", 3)
V2D_rz  = FunctionSpace(mesh, "N2curl", 2)
V = V1D_phi * V2D_rz

(N_phi_i, M_rz_i) = TestFunctions(V)
(N_phi_j, M_rz_j) = TrialFunctions(V)

r = Expression('x[0]')
n = FacetNormal(mesh)
N = 1.0

ds = Measure("ds")[conductive]
cond = 66e6
skin = 1e-7
Zm   = 1.0/(cond*skin)

UV = Function(V)	# full E result

for target_frequency in target_frequencies:
	k_0_sq = ((2*pi*target_frequency/c0)**2)
	
	# REMARK: If k set to k=0, TM modes are calculated with <1e-18 fluctuations in phi space.
	k=0
	
	a_mn =   r*dot( curl_rz(M_rz_i, N_phi_i, r, k), curl_rz(M_rz_j, N_phi_j, r, k) )*dx \
	       + r*dot( curl_p(M_rz_i), curl_p(M_rz_j) )*dx
	b_mn =   r*dot( N_phi_i, N_phi_j )*dx \
	       + r*dot( M_rz_i, M_rz_j )*dx
	g_mn = -eta0/Zm * (   r*dot( cross_2D(n, cross_1D(n, N_phi_i)), N_phi_j )*ds(1)  \
			    + r*dot( n_x_n_x_A(n, M_rz_i),              M_rz_j )*ds(1) )
	
	S_te = PETScMatrix()
	T_te = PETScMatrix()
	
	assemble(a_mn, tensor=S_te)
	assemble(b_mn + g_mn, tensor=T_te) # , keep_diagonal=True)

	te = SLEPcEigenSolver(S_te,T_te)
	te.parameters["tolerance"] = 1e-16
	te.parameters["spectral_shift"] = k_0_sq
	te.parameters["spectral_transform"] = "shift-and-invert"
	te.parameters["spectrum"] = "target real"
	
	te.solve(1)
	
	print("Extracting from %d results." % te.get_number_converged())
	for i in range(te.get_number_converged()):
		k_r, k_i, ev_r, ev_i = te.get_eigenpair(i)
		print k_i/k_r
		if False and k_r <= 0.0: # Remove False to check for negative results
			print("r negative")
			continue
		f = np.sqrt(k_r)*c0/(2*np.pi)
		if False and f <= target_frequency*0.8:	# Remove False to suppress results far off the target frequency
			logging.debug("f smaller than f_target: %e" % f)
			continue
		print("at f=%.2f: %4d: f_calc=%.10f GHz" % (target_frequency/1e9, i, f/1e9))
		
		if True:	# if clause for on/off switching of visualisation ;-)
			UV.vector()[:] = ev_r
			(v_phi, u_rz) = UV.split()
			plot(u_rz, title=("u_rz, f_ana=%f, f_res=%.3f" % (target_frequency/1e9, f/1e9)))
			plot(v_phi, title=("v_phi, f_ana=%f, f_res=%.3f" % (target_frequency/1e9, f/1e9)))
			interactive()

print("Done.")
logging.shutdown()