Hello everyone, I used two dolfinx versions #957 and the newest version to solve the general eigenvalue problem. However, they give different results 10.839566377715625 and 2.9713814818807514 respectively, and I’ll past their codes here for comparison. I hope to get your help.
#957
import pygmsh
import meshio
import numpy as np
geom = pygmsh.built_in.Geometry()
side_len = 0.0127
strip_w = 1.27e-3
strip_h = 1.27e-4
die_h = 1.27e-3
point_coords = [[0, 0, 0], [side_len, 0, 0],
[side_len, die_h, 0], [side_len, side_len, 0],
[0, side_len, 0], [0, die_h, 0],
[side_len/2+strip_w/2, die_h, 0], [side_len/2-strip_w/2, die_h, 0],
[side_len/2+strip_w/2, die_h+strip_h, 0], [side_len/2-strip_w/2, die_h+strip_h, 0]]
dense_point_idx = [6, 7, 8, 9]
points = []
for i, point in enumerate(point_coords):
if i in dense_point_idx:
lcar = 5e-5
else:
lcar = 5e-4
points.append(geom.add_point(point, lcar=lcar))
line_tag = [(0, 1), (1, 2), (2, 3), (3, 4), (4, 5),
(5, 0), (7, 5), (6, 7), (2, 6), (6, 8),
(8, 9), (9, 7)]
lines = []
for i, j in line_tag:
lines.append(geom.add_line(points[i], points[j]))
dielectric_loop_tag = [0, 1, 8, 7, 6, 5]
air_loop_tag = [-6, -11, -10, -9, -8, 2, 3, 4]
domain_loop_tag = [0, 1, 2, 3, 4, 5]
strip_loop_tag = [7, 9, 10, 11]
dielectric_loop = geom.add_line_loop([lines[i] for i in dielectric_loop_tag])
air_loop = []
for i in air_loop_tag:
if i >= 0:
air_loop.append(lines[i])
else:
air_loop.append(-lines[-i])
air_loop = geom.add_line_loop(air_loop)
dielectric_domain = geom.add_plane_surface(dielectric_loop)
air_domain = geom.add_plane_surface(air_loop)
geom.add_physical([lines[i] for i in domain_loop_tag + strip_loop_tag], 'PEC')
geom.add_physical(dielectric_domain, 'Diectric')
geom.add_physical(air_domain, 'Air')
geo_file = 'strip.geo'
with open(geo_file, 'w') as f:
f.write(geom.get_code())
mesh = pygmsh.generate_mesh(geom, prune_z_0=True, verbose=True, geo_filename="strip.geo")
points, cells, cell_data, field_data = mesh.points, mesh.cells, mesh.cell_data, mesh.field_data
cells = np.vstack(np.array([cells.data for cells in mesh.cells
if cells.type == "triangle"]))
facet_cells = np.vstack(np.array([cells.data for cells in mesh.cells
if cells.type == "line"]))
cell_data = mesh.cell_data_dict["gmsh:physical"]["triangle"]
facet_data = mesh.cell_data_dict["gmsh:physical"]["line"]
meshio.xdmf.write('tag_triangle.xdmf', meshio.Mesh(
points=points,
cells={'triangle': cells},
cell_data={
'triangle': [cell_data]
},
#field_data=field_data
))
meshio.xdmf.write("tag_all.xdmf", meshio.Mesh(
points=points,
cells={"line": facet_cells},
cell_data={"line": [facet_data]}
))
import numpy as np
import pygmsh
import meshio
from dolfinx.io import XDMFFile
from dolfinx import (MPI, FacetNormal, Function, FunctionSpace,
has_petsc_complex, DirichletBC, cpp)
from dolfinx.fem.assemble import assemble_scalar, assemble_matrix
from dolfinx.fem import locate_dofs_topological
from ufl import (TestFunctions, TrialFunctions, grad, inner, curl, split,
Measure, FiniteElement, dx, Dx, dot)
import dolfinx
from scipy.constants import c
from scipy.sparse import csr_matrix
from scipy.sparse.linalg import eigs
from slepc4py import SLEPc
from petsc4py import PETSc
with XDMFFile(MPI.comm_world, "tag_triangle.xdmf", 'r') as xdmf_infile:
mesh = xdmf_infile.read_mesh(name="Grid")
mesh.create_connectivity(mesh.topology.dim, mesh.topology.dim)
mvc_subdomain = xdmf_infile.read_meshtags(mesh, "Grid")
mesh.create_connectivity(mesh.topology.dim-1, mesh.topology.dim)
with XDMFFile(MPI.comm_world, "tag_all.xdmf", 'r') as xdmf_infile:
mvc_boundaries = xdmf_infile.read_meshtags(mesh, "Grid")
#%%
V = FiniteElement('N1curl', mesh.ufl_cell(), 2)
Q = FiniteElement('P', mesh.ufl_cell(), 3)
W = FunctionSpace(mesh, V*Q)
#%%
S = FunctionSpace(mesh, ('DG', 0))
e_r = Function(S)
x = S.tabulate_dof_coordinates()
for i in range(x.shape[0]):
if mvc_subdomain.values[i] == 2:
e_r.vector.setValueLocal(i, 8.875)
else:
e_r.vector.setValueLocal(i, 1)
#%%
u, p = TrialFunctions(W)
v, q = TestFunctions(W)
#%%
f = 25e9
k0_sqr = (2*np.pi*f/c)**2
electric_wall = Function(W)
with electric_wall.vector.localForm() as bc_local:
bc_local.set(0)
bndry_facets = np.where(mvc_boundaries.values == 1)[0]
bdofs = locate_dofs_topological(W, 1, bndry_facets)
bc= DirichletBC(electric_wall, bdofs)
#%%
theta_sqr = k0_sqr * 8.875
#%%
a = 1/theta_sqr*(inner(curl(u), curl(v)) - k0_sqr*e_r*inner(u, v))*dx
b = inner(grad(p), v)*dx + inner(u, v)*dx + inner(grad(p), grad(q))*dx + inner(u, grad(q))*dx - k0_sqr*e_r*inner(p, q)*dx
c = b + a
#%%
S = assemble_matrix(b, [bc])
T = assemble_matrix(c, [bc])
#%%
S.assemble()
S = S.realPart()
si, sj, sv = S.getValuesCSR()
S_sp = csr_matrix((sv, sj, si))
T.assemble()
T = T.realPart()
ti, tj, tv = T.getValuesCSR()
T_sp = csr_matrix((tv, tj, ti))
#%%
indicators = np.zeros(S_sp.shape[0])
indicators[bc.dof_indices[:, 0]] = 1
free_dofs = np.where(indicators == 0)[0]
S_np = S_sp[free_dofs, :][:, free_dofs]
T_np = T_sp[free_dofs, :][:, free_dofs]
es = SLEPc.EPS(which='LR').create(comm=MPI.comm_world)
es.setDimensions(10)
es.setOperators(S, T)
es.setFromOptions()
es.solve()
#%%
vr, vi = S.getVecs()
lam = es.getEigenpair(0, vr, vi)
((1-1/lam)*theta_sqr)**0.5/k0_sqr**0.5
print(lam)
The newest vestion:
#!/usr/bin/env python
# coding: utf-8
# In[1]:
import pygmsh
import gmsh
import numpy as np
from mpi4py import MPI
from petsc4py import PETSc
from slepc4py import SLEPc
from scipy.constants import c
from ufl import curl, dx, FiniteElement, grad, inner, Measure, TestFunctions, TrialFunctions
from dolfinx import Constant, Function, FunctionSpace, DirichletBC
from dolfinx.fem import assemble_matrix, DofMapRestriction, locate_dofs_topological
from dolfinx.io import XDMFFile
# ### Mesh
# In[2]:
with pygmsh.occ.Geometry() as geom:
side_len = 0.0127
strip_w = 1.27e-3
strip_h = 1.27e-4
die_h = 1.27e-3
point_coords = [[0, 0, 0], [side_len, 0, 0],
[side_len, die_h, 0], [side_len, side_len, 0],
[0, side_len, 0], [0, die_h, 0],
[side_len/2+strip_w/2, die_h, 0], [side_len/2-strip_w/2, die_h, 0],
[side_len/2+strip_w/2, die_h+strip_h, 0], [side_len/2-strip_w/2, die_h+strip_h, 0]]
dense_point_idx = [6, 7, 8, 9]
points = []
for i, point in enumerate(point_coords):
if i in dense_point_idx:
lcar = 5e-5
else:
lcar = 5e-4
points.append(geom.add_point(point, lcar))
line_tag = [(0, 1), (1, 2), (2, 3), (3, 4), (4, 5),
(5, 0), (7, 5), (6, 7), (2, 6), (6, 8),
(8, 9), (9, 7)]
lines = []
for i, j in line_tag:
lines.append(geom.add_line(points[i], points[j]))
dielectric_loop_tag = [0, 1, 8, 7, 6, 5]
air_loop_tag = [-6, -11, -10, -9, -8, 2, 3, 4]
domain_loop_tag = [0, 1, 2, 3, 4, 5]
strip_loop_tag = [7, 9, 10, 11]
dielectric_loop = geom.add_curve_loop([lines[i] for i in dielectric_loop_tag])
air_loop = []
for i in air_loop_tag:
if i >= 0:
air_loop.append(lines[i])
else:
air_loop.append(-lines[-i])
air_loop = geom.add_curve_loop(air_loop)
dielectric_domain = geom.add_plane_surface(dielectric_loop)
air_domain = geom.add_plane_surface(air_loop)
geom.add_physical([lines[i] for i in domain_loop_tag + strip_loop_tag], 'PEC')
geom.add_physical(dielectric_domain, 'Diectric')
geom.add_physical(air_domain, 'Air')
mesh = geom.generate_mesh(2)
from gmsh_helpers import gmsh_model_to_mesh
mesh, cell_tags, facet_tags = gmsh_model_to_mesh(
gmsh.model, cell_data=True, facet_data=True, gdim=2)
V = FiniteElement('N1curl', mesh.ufl_cell(), 2)
Q = FiniteElement('P', mesh.ufl_cell(), 3)
W = FunctionSpace(mesh, V*Q)
dofs_W = np.arange(0, W.dofmap.index_map_bs * (
W.dofmap.index_map.size_local + W.dofmap.index_map.num_ghosts))
fdim = mesh.topology.dim - 1
pec_facets = facet_tags.indices[facet_tags.values == 1]
bdofs_W = locate_dofs_topological(W, fdim, pec_facets)
non_bdofs_W = np.setdiff1d(dofs_W, bdofs_W)
restriction_W_Gamma = DofMapRestriction(W.dofmap, non_bdofs_W) # TODO: remove without fork
electric_wall = Function(W)
with electric_wall.vector.localForm() as bc_local:
bc_local.set(0)
bc = DirichletBC(electric_wall, bdofs_W)
u, p = TrialFunctions(W)
v, q = TestFunctions(W)
# In[7]:
S = FunctionSpace(mesh, ('DG', 0))
e_r = Function(S)
x = S.tabulate_dof_coordinates()
for i in range(x.shape[0]):
if cell_tags.values[i] == 2:
e_r.vector.setValueLocal(i, 8.875)
else:
e_r.vector.setValueLocal(i, 1)
# In[8]:
f = 25e9
k0_sqr = (2*np.pi*f/c)**2
theta_sqr = k0_sqr * 8.875
# In[9]:
print(theta_sqr, k0_sqr)
# In[10]:
a = inner(grad(p), v)*dx + inner(u, v)*dx + inner(grad(p), grad(q))*dx + inner(u, grad(q))*dx - k0_sqr*e_r*inner(p, q)*dx
b = a + 1/theta_sqr*(inner(curl(u), curl(v)) - k0_sqr*e_r*inner(u, v))*dx
# In[11]:
# TODO: requires fork. If you want to use master dolfinx, you should throw away Dirichlet BCs
# using a petsc4py submatrix, see e.g. https://github.com/michalhabera/dolfiny/blob/master/dolfiny/restriction.py#L69
#A = assemble_matrix(a, bcs=[], restriction=(restriction_W_Gamma, restriction_W_Gamma))
A = assemble_matrix(a, bcs=[bc])
A.assemble()
B = assemble_matrix(b, bcs=[bc])
#B = assemble_matrix(b, bcs=[], restriction=(restriction_W_Gamma, restriction_W_Gamma))
B.assemble()
# In[12]:
# Solve
eps = SLEPc.EPS().create()
eps.setOperators(A, B)
eps.setProblemType(SLEPc.EPS.ProblemType.GNHEP)
eps.setDimensions(1, PETSc.DECIDE, PETSc.DECIDE)
eps.setWhichEigenpairs(SLEPc.EPS.Which.LARGEST_REAL)
eps.getST().getKSP().setType("preonly")
eps.getST().getKSP().getPC().setType("lu")
eps.getST().getKSP().getPC().setFactorSolverType("mumps")
eps.solve()
assert eps.getConverged() >= 1
# In[13]:
# Extract leading eigenvalue
eigv = eps.getEigenvalue(0)
r, i = eigv.real, eigv.imag
#assert abs(i) < 1.e-10
print(r, i)
print(((1-1/r)*theta_sqr)**0.5/k0_sqr**0.5)