Different solution values for different dolfinx version

Since early last year we removed eigen as a dependency of dolfinx (we use xtensor instead in the newer releases). Seems like you Need to install eigen.

I find it strange that some of the outputs you have presented finds slepc, others doesn’t.

Dear dokken,
Thank you. I’ll try to install eigen. Is the slepc necessary to dolfinx?

Dear dokken,
How to install eigen? Slepc can be found in commits later than #1202 (include), but cannot be found in commits #1201, #1200, …

Dear dokken,
I have installed the old version #980 and got the desired eigenvalue now. The problem encountered in the installation is caused by the residue of the uninstallation of the newest version. After I removed them with the commands, I can install the older version. Additionally, why does anyone maintaining dolfinx check the results between the #980 and the newest dolfinx?

sudo rm -rf /usr/local/include/dolfinx*
sudo rm /usr/local/lib/cmake/dolfinx/DOLFINXTargets.cmake
sudo rm -rf /usr/local/cmake/dolfinx
sudo rm -rf /usr/local/lib/cmake/dolfinx

You are free to make a test that checks this, and then add a pull request to dolfinx with it. The point of the tests are that they are simple, or solve known problems with known solutions for verification. See for instance all the tests in: dolfinx/test_fem_pipeline.py at main · FEniCS/dolfinx · GitHub
which checks assembly and solve for a large variety of finite elements and degrees.

As I previously said, I’m happy to try to help you if you make the codes look similar between versions of dolfinx.

I pointed out a some major differences in Different solution values for different dolfinx version - #4 by dokken
which you still have not explained to me.

Dear dokken,
I need to install the branch francesco-ballarin/dolfinx at migration (github.com) first to repeat my result. Because francesco-ballarin answer my question about 2020.6.5 eigenvalue problem of newest version · Issue #957 · FEniCS/dolfinx (github.com), so I installed his branch of that time and used his code, but an error reports “cannot import DofMapRestriction from dolfinx.fem”.

You need to submit an issue on his fork then.

There is no such class in dolfinx. It appears in multiphenicsx, not in dolfinx.

Dear dokken,
Francesco gave me a code working with his branch last year 2020.6.5 GitHub - francesco-ballarin/dolfinx at migration. I tried a lot, but can’t find his branch code at that time. Could you tell me how to get it? Thank you in advance.

Dear dokken,
I have compared the same code results executed by the two versions of dolfinx #1042 and #1756. They give different results 6.4439947471528845 and 4740.365909597999.
The code

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, locate_dofs_topological
from dolfinx.io import XDMFFile

with XDMFFile(MPI.COMM_WORLD, "tag_triangle.xdmf", 'r') as xdmf_infile:
    mesh = xdmf_infile.read_mesh(name="Grid")
    mesh.topology.create_connectivity(mesh.topology.dim, mesh.topology.dim)
    mvc_subdomain = xdmf_infile.read_meshtags(mesh, "Grid")

mesh.topology.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)

bdofs_W = locate_dofs_topological(W, mvc_boundaries.dim, mvc_boundaries.indices[mvc_boundaries.values == 1] )

u, p = TrialFunctions(W)
v, q = TestFunctions(W)

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)

f = 25e9
k0_sqr = (2*np.pi*f/c)**2
theta_sqr = k0_sqr * 8.875

electric_wall = Function(W)

with electric_wall.vector.localForm() as bc_local:
    bc_local.set(0)

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

boundary_dofs = locate_dofs_topological(W, mvc_boundaries.dim,
                                        mvc_boundaries.indices[mvc_boundaries.values == 1])
bc = DirichletBC(electric_wall, boundary_dofs)

A = assemble_matrix(a, bcs=[bc])
A.assemble()
B = assemble_matrix(b, bcs=[bc])
B.assemble()

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

eigv = eps.getEigenvalue(0)
r, i = eigv.real, eigv.imag
#assert abs(i) < 1.e-10
print(r)

The xdmf files can be downloaded from here different eigenvalues in different dolfinx versions · Issue #1761 · FEniCS/dolfinx (github.com)

Could you give me the output of


import numpy as np
from mpi4py import MPI
from ufl import dx, FiniteElement, grad, inner, TestFunctions, TrialFunctions, MixedElement, TrialFunction, TestFunction
from dolfinx import FunctionSpace, UnitSquareMesh
from dolfinx.fem import assemble_matrix


N = 1
mesh = UnitSquareMesh(MPI.COMM_WORLD, N, N)


V = FiniteElement('N1curl', mesh.ufl_cell(), 2)
Q = FiniteElement('P', mesh.ufl_cell(), 3)
VQ = MixedElement([V, Q])

def mixed():
    W = FunctionSpace(mesh, VQ)

    u, p = TrialFunctions(W)
    v, q = TestFunctions(W)
    a = inner(u, v) * dx
    A = assemble_matrix(a, bcs=[])
    A.assemble()
    Anp =  A[:, :]
    np.set_printoptions(precision=3)
    print(A.norm(0))


def single():
    W = FunctionSpace(mesh, V)

    u = TrialFunction(W)
    v = TestFunction(W)
    a = inner(u, v) * dx
    A = assemble_matrix(a, bcs=[])
    A.assemble()
    Anp =  A[:, :]
    np.set_printoptions(precision=3)
    print(A.norm(0))


mixed()
single()

on your #1042 installation?

On the latest version of dolfinx, you get:

root@de701b1b0fa3:~/shared# python3 main.py 
8.622222222222236
8.622222222222199

One key difference from older versions of dolfin-x and the later versions, is that we no longer use point evaluations for Nedelec spaces (they are defined through integral moments). This could result in different matrices, and thus different eigenvalues.

Dear dokken,
I have changed to #1212 now as it can also give me the desired result. Your script gives the results:
8.044444444444462
8.04444444444443

As mentioned previously, how Nedelec finite element functions are defined (through edge integration, not point evaluation) might change the structure of your matrix (This change happened in january/februar 2021). @mscroggs might be able to shed some more light on this.

If you can pinpoint a commit where the eigenvalues in your problem go from what you mean is correct, to wrong, it would be helpful. The commit you are at now is still almost a year back (fix typo introduced to FindSLEPc.cmake by drew-parsons · Pull Request #1212 · FEniCS/dolfinx · GitHub)

Dear dokken,
6.443994747 is the correct one.

My point is that you need to pinpoint at what commit this change is introduced. Giving me a span on 1 year of commits is not something I have time to sit down and debug.

Dear dokken,
I have something to do this days and will try this when I’m free.