eigenvalue problem of newest version

[wwshunan]

Hi everyone, I used the newest version of dolfinx to solve a microstrip eigenvalue problem. It gives a different first eigenvalue from the one calculated by #843 all using the same mesh file. The newest version has made some minor error. Here is the mesh file and my dolfinx code.

mesh file:

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]}
))

dolfinx file in newest version:

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)

dolfinx code for #843:

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") as xdmf_infile:
    mesh = xdmf_infile.read_mesh(cpp.mesh.GhostMode.none)
    mvc_subdomain = xdmf_infile.read_mvc_size_t(mesh)
    # tag_info = xdmf_infile.read_information_int()

with XDMFFile(MPI.comm_world, "tag_all.xdmf") as xdmf_infile:
    mvc_boundaries = xdmf_infile.read_mvc_size_t(mesh)

mf_triangle = cpp.mesh.MeshFunctionSizet(mesh, mvc_subdomain, 0)
mf_line = cpp.mesh.MeshFunctionSizet(mesh, mvc_boundaries, 0)

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 mf_triangle.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(mf_line.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()

T.assemble()
T = T.realPart()

es = SLEPc.EPS(which='LR').create(comm=MPI.comm_world)
es.setDimensions(1)
es.setOperators(S, T)
es.setFromOptions()
es.solve()

vr, vi = S.getVecs()
lam = es.getEigenpair(0, vr, vi)
print(lam)

[francesco-ballarin]

Can you please clarify what you mean by

has made some minor error

?

Do you mean the values of lam are different? If so, please write them in your reply.

[wwshunan]

843 gives "lam" value 6.44970770628454, while the newest version gives 10.839566377715625

[francesco-ballarin]

I suggest that you proceed in this way: 1) using XDMFFile, export the MeshFunction (old version) and MeshTags (new version) that you use in the two files, to confirm with ParaView that the markers are actually the same in both versions. I understand that you expect them to be the same (since they are read from the same file generated by pygmsh): the goal here is to double check that they actually are. 2) confirm that the matrices S and T, which you provide to SLEPc, are actually the same. I would start with comparing some norms and associated condition numbers, and only if they match compare that the actual entries are the same 3) double check that the initialization of the es object is correct, especially for what concerns solver parameters. I cannot comment on the actual physical problem because I am not able to recognize it from the implementation, but I am afraid that there is something wrong in that setup because running the new version on my laptop gives -96.94738627105075+0j, i.e. yet another value.

[wwshunan]

It gives 97.24680980575958-6.808448627744403e-14j in the newest dolfinx. So the mesh changes randomly from one version of dolfinx to another version. If you used the Ubuntu default gmsh rather than the docker gmsh version to generate the mesh, it should give 10.839566377715625 which is not changed with dolfinx version.

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
print(cells, cell_data)
meshio.xdmf.write('tag_triangle.xdmf', meshio.Mesh(
    points=points,
    cells={'triangle': cells[1].data},
    cell_data={
        'triangle': [cell_data['gmsh:physical'][1]]
    },
    #field_data=field_data
))
line_tags = cell_data['gmsh:physical'][0]
meshio.xdmf.write("tag_all.xdmf", meshio.Mesh(
    points=points,
    cells={"line": cells[0].data},
    cell_data={"line": [line_tags]}
))

[francesco-ballarin]

Please try and carry out items 1, 2 and 3 in my post above, because we had four different results in the last three posts, and without breaking the checks into smaller parts it is impossible to understand where the issue is.

[wwshunan]

I'm not familiar with many methods in dolfinx and don't know how to write meshtags to xdmf as I tried with failure. It's not an easy task for me to check them. I cannot finish them if you not give detailed indications.

[francesco-ballarin]

1) since XDMFFile changed between your "old" installation and the current one, you might want to have a look at https://github.com/FEniCS/dolfinx/blob/e2a96c4211bb0e4fda9256657ce23afc756c05a5/python/dolfinx/io.py#L147 for what concerns the old installation, and https://github.com/FEniCS/dolfinx/blob/master/cpp/dolfinx/io/XDMFFile.h#L56 for the current one. You will easily find on your own what is the name of the appropriate write method for the task.

2) petsc4py has methods for each of the operation I asked, see https://www.mcs.anl.gov/petsc/petsc4py-current/docs/apiref/petsc4py.PETSc.Mat-class.html

3) you might find helpful to have a look at how I setup an EPS object in one of my demos https://gitlab.com/multiphenics/multiphenics/-/blob/fenicsx/tutorials/04_infsup_stokes/tutorial_infsup_stokes.ipynb

[wwshunan]

Hi, I have already finished the first step. The meshtag file imported to paraview with the option "Xdmf3ReaderS" only has the default white background color, while it shows the read and blue color representing different tag values with the option "Xdmf reader" for the newest dolfinx version. It behaves normal for the #843 version. I don't know if this difference causes problem.

[francesco-ballarin]

You should probably use Xdmf3ReaderT, so that you can select/deselect blocks. In any case, if the marked entities look the same in both cases, then step 1 is probably ok.

[wwshunan]

The option "Xdmf3ReaderS" did not show the red and blue color for the meshtag file for the newest dolfinx version. Is it OK?

[wwshunan]

The 2nd step norm result comparison:

843:

1169.9355192266012 393301.4261144103

newest dolfinx: 1166.2513958086834 391283.29239720985

[wwshunan]

Hi, is this problem solved?

[francesco-ballarin]

I think there is nothing wrong with current dolfinx master. I attach a notebook in which I am able to obtain the desired value of 6.443994747149748.

To my understanding, your main issue is that you are not throwing away rows/cols associated to Dirichlet boundary conditions. Due to time constraints, in the notebook I provide you a code that works with my fork, which has additional functionalities that are helpful in carrying out this task. Still, I marked with TODO the two things that you need to change to have a code that is compatible with dolfinx master (by avoiding using my additional functionalities), assuming that you can work on these two TODOs on your own.

I am closing this, because IMHO it is not a programming issue inside dolfinx. If you think instead that there still is some dolfinx issue associated to your problem, feel free to reopen this issue, but if possible try and make an effort to provide a more minimal example that highlights a specific problem.

issue.zip

[wwshunan]

Now the newest branch maintained by francesco-ballarin still has the problem to solve the same eigenvalue problem correctly. I cannot get the desired value 6.443994747149748 with the same code issue.py. Please see if there is anything wrong with the new code.