Maximum value in vector and paraview is different

wwshunan · December 8, 2021, 1:54pm

Dear all,

  I solved a general eigenvalue problem and  got different maximum values of the eigenvector in the code and showed in paraview. The code and the results are:

a_tt = 1/k0_sqr*(inner(curl(u), curl(v)) - k0_sqr*e_r*inner(u, v))*dx
b_tt = inner(u, v)*dx
b_tz = inner(u, grad(q))*dx
b_zt = inner(grad(p), v)*dx
b_zz = inner(grad(p), grad(q))*dx - k0_sqr*e_r*inner(p, q)*dx
s_a = 1 / theta_sqr * k0_sqr * a_tt
s_b = b_tt + b_zz + b_tz + b_zt
s_c = s_b + s_a

B = assemble_matrix(s_b, [bc])
B.assemble()

C = assemble_matrix(s_c, [bc])
C.assemble()
eps = SLEPc.EPS()
eps.create()
eps.setOperators(B, C)
eps.solve()

xr, xi = B.getVecs()
eps.getEigenpair(j, xr, xi)
lam = eps.getEigenpair(j, xr, xi)

x_r = Function(W)
xr.copy(x_r.vector)
x_i = Function(W)
xi.copy(x_i.vector)

xt_r, xz_r = x_r.split()
xt_i, xz_i = x_i.split()
xt_r, xz_r = x_r.sub(0).collapse(), x_r.sub(1).collapse()
print(np.abs(xt_r.vector.array).max())   #3.00074500898824e-05

with XDMFFile(MPI.COMM_WORLD, "Et_r.xdmf", 'w',
        encoding=XDMFFile.Encoding.HDF5) as file:
file.write_mesh(mesh)
file.write_function(xt_r)

The maximum xt_r is about 3e-5, however the maximum value showed in paraview is about 0.12. Why? Thanks for any help.

dokken · December 9, 2021, 10:09am

Please supply a minimal example that can reproduce the error. In the code above, you have excluded import statements and several crucial definitions to make the code executable and reproducable.

wwshunan · December 9, 2021, 10:16am

Dear dokken,
Does Paraview use the raw vector value or the value processed by dolfinx?

wwshunan · December 9, 2021, 10:26am

Dear dokken,
Here is my code

import numpy as np
import gmsh, meshio
from mpi4py import MPI
from slepc4py import SLEPc
from ufl import curl, dx, FiniteElement, grad, inner, Measure, TestFunctions, TrialFunctions
from dolfinx import Constant, Function, FunctionSpace, geometry, cpp, DirichletBC
from dolfinx.fem import locate_dofs_topological, assemble_matrix, Form, assemble_matrix_nest
from scipy.constants import c, epsilon_0
from dolfinx.io import XDMFFile
from petsc4py import PETSc

def create_mesh(mesh, cell_type, prune_z=False):
    cells = mesh.get_cells_type(cell_type)
    cell_data = mesh.get_cell_data("gmsh:physical", cell_type)
    out_mesh = meshio.Mesh(points=mesh.points, cells={cell_type: cells}, cell_data={"name_to_read":[cell_data]})
    if prune_z:
        out_mesh.prune_z_0()
    return out_mesh

with pygmsh.geo.Geometry() as geom:
    core_r = 0.0015
    max_r = 1900
    inner_lcar = 0.00007
    outer_lcar = 40.1
    theta = np.pi / 10
    p0 = [0, 0]
    p1 = [core_r, 0]
    p2 = [core_r * np.cos(theta), core_r * np.sin(theta)]
    p3 = [max_r * np.cos(theta), max_r * np.sin(theta)]
    p4 = [max_r, 0]
    point_coords = [p0, p1, p2, p3, p4]
    lcars = [inner_lcar, inner_lcar, inner_lcar, outer_lcar, outer_lcar]
    points = []

    for p, lcar in zip(point_coords, lcars):
        points.append(geom.add_point([p[0], p[1]], lcar))

    line1 = geom.add_line(points[0], points[1])
    inner_arc = geom.add_circle_arc(points[1], points[0], points[2])
    line2 = geom.add_line(points[2], points[0])
    inner_loop_curves = [line1, inner_arc, line2]
    line3 = geom.add_line(points[2], points[3])
    outer_arc = geom.add_circle_arc(points[3], points[0], points[4])
    line4 = geom.add_line(points[4], points[1])
    outer_loop_curves = [inner_arc, line3, outer_arc, line4]
    inner_curve_loop = geom.add_curve_loop(inner_loop_curves)
    outer_curve_loop = geom.add_curve_loop(outer_loop_curves)
    inner_surface = geom.add_plane_surface(inner_curve_loop)
    outer_surface = geom.add_plane_surface(outer_curve_loop)
    geom.add_physical(outer_arc, 'Dirichlet')
    geom.add_physical(inner_surface, 'Conductor')
    geom.add_physical(outer_surface, 'Air')

    mesh = geom.generate_mesh(2)

    gmsh.write("mesh.msh")
    msh = meshio.read("mesh.msh")

    triangle_mesh = create_mesh(msh, "triangle", prune_z=True)
    line_mesh = create_mesh(msh, "line", prune_z=True)
    meshio.write("mesh.xdmf", triangle_mesh)
    meshio.write("mt.xdmf", line_mesh)


with XDMFFile(MPI.COMM_WORLD, "mesh.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, "mt.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)

num_subs = W.num_sub_spaces()
spaces = []
maps = []
for i in range(num_subs):
    space_i, map_i = W.sub(i).collapse(collapsed_dofs=True)
    spaces.append(space_i)
    maps.append(map_i)


u, p = TrialFunctions(W)
v, q = TestFunctions(W)
omega = 2e6

sigma_core = 8678
sigma_corona = 5e-6
k0_sqr = (omega/c)**2


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, 1-1j*sigma_core/(omega*epsilon_0))
    else:
        e_r.vector.setValueLocal(i, 1)


theta_sqr = 4e-05-0.0002j

a_tt = 1/k0_sqr*(inner(curl(u), curl(v)) - k0_sqr*e_r*inner(u, v))*dx
b_tt = inner(u, v)*dx
b_tz = inner(u, grad(q))*dx
b_zt = inner(grad(p), v)*dx
b_zz = inner(grad(p), grad(q))*dx - k0_sqr*e_r*inner(p, q)*dx
s_a = 1 / theta_sqr * k0_sqr * a_tt
s_b = b_tt + b_zz + b_tz + b_zt
s_c = s_b + s_a

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

bndry_facets = np.where(mvc_boundaries.values == 0)[0]
bdofs = locate_dofs_topological(W, 1, bndry_facets)
bc = DirichletBC(electric_wall, bdofs)
B = assemble_matrix(s_b, [bc])
B.assemble()
C = assemble_matrix(s_c, [bc])
C.assemble()

D = assemble_matrix(b_zt, [bc])
D.assemble()

eps = SLEPc.EPS()
eps.create()
eps.setOperators(B, C)
eps.setDimensions(100)
eps.solve()

n = eps.getConverged()

xr, xi = B.getVecs()
for j in range(n):
    lam = eps.getEigenpair(j, xr, xi)
    if (((1-1/lam)*theta_sqr)**0.5).real > 0.0002 and abs((((1-1/lam)*theta_sqr)**0.5).imag) > 0.0002:
        print(lam)
        print(((1-1/lam)*theta_sqr)**0.5, 'xxx')

        print(np.sum(xr.getArray()))

        x_r = Function(W)
        xr.copy(x_r.vector)
        x_i = Function(W)
        xi.copy(x_i.vector)

        xt_r, xz_r = x_r.split()
        xt_i, xz_i = x_i.split()
        print(np.abs(xt_r.vector.array).max())
        xt_r, xz_r = x_r.sub(0).collapse(), x_r.sub(1).collapse()
        print(np.abs(xt_r.vector.array).max(), xz_r.vector.array.max())

        n = len(maps[1])
        x = PETSc.Vec().createSeq(n)
        b = PETSc.Vec().createSeq(n)
        iset_0 = PETSc.IS().createGeneral(maps[0])
        iset_1 = PETSc.IS().createGeneral(maps[1])
        mat_zt = assemble_matrix(b_zt, [bc])
        mat_zt.assemble()
        mat_zz = assemble_matrix(b_zz, [bc])
        mat_zz.assemble()
        sub_mat_zt = mat_zt.createSubMatrix(iset_1, iset_0)
        sub_mat_zz = mat_zz.createSubMatrix(iset_1, iset_1)
        sub_mat_zt.mult(xt_r.vector, b)
        ksp = PETSc.KSP().create()
        ksp.setOperators(sub_mat_zz)
        ksp.setFromOptions()
        ksp.solve(b, x)

        with XDMFFile(MPI.COMM_WORLD, "Et_r.xdmf", 'w',
                      encoding=XDMFFile.Encoding.HDF5) as file:
            file.write_mesh(mesh)
            file.write_function(xt_r)

        with XDMFFile(MPI.COMM_WORLD, "Ez_r.xdmf", 'w',
                      encoding=XDMFFile.Encoding.HDF5) as file:
            file.write_mesh(mesh)
            file.write_function(xz_r)

        with XDMFFile(MPI.COMM_WORLD, "Et_i.xdmf", 'w',
                      encoding=XDMFFile.Encoding.HDF5) as file:
            file.write_mesh(mesh)
            file.write_function(xt_i)

        with XDMFFile(MPI.COMM_WORLD, "Ez_i.xdmf", 'w',
                      encoding=XDMFFile.Encoding.HDF5) as file:
            file.write_mesh(mesh)
            file.write_function(xz_i)

dokken · December 9, 2021, 10:29am

No, it interpolates the data into CG-1 if you use XDMFFile.

For spaces suchs as N1Curl, you should interpolate the solution into a DG space of the same order, and use VTXWriter: dolfinx/test_adios2.py at main · FEniCS/dolfinx · GitHub

wwshunan · December 9, 2021, 10:56am

Dear dokken,
I cannot import VTXWriter from dolfinx.cpp.io. My dolfinx version is #1212. Could you tell me how to replace the XDMFFile with VTXWriter or VTXFile for ‘N1curl’ space with order 2? xt_r is the sub-vector collapsed with the code x_r.sub(0).collapse() where x_r is the real part of the mixed space function.

with XDMFFile(MPI.COMM_WORLD, "Et_r.xdmf", 'w',
                      encoding=XDMFFile.Encoding.HDF5) as file:
    file.write_mesh(mesh)
    file.write_function(xt_r)

Thank you very much.

dokken · December 9, 2021, 11:33am

You cannot use the version 1212 with VTXWriter as it predates it’s existence: ADIOS2 support (VTX and Fides) by jorgensd · Pull Request #1655 · FEniCS/dolfinx · GitHub

wwshunan · December 9, 2021, 11:48am

Dear dokken,d
Then how to interpolate from N1curl space to DG space? I really encountered the same problem a few years ago. In addition, I want to do the multiply between the submatrix at the left bottom of the whole matrix and the corresponding sub-function. Is the following code quoted from the above code correct?

        n = len(maps[1])
        x = PETSc.Vec().createSeq(n)
        b = PETSc.Vec().createSeq(n)
        iset_0 = PETSc.IS().createGeneral(maps[0])
        iset_1 = PETSc.IS().createGeneral(maps[1])
        mat_zt = assemble_matrix(b_zt, [bc])
        mat_zt.assemble()
        sub_mat_zt = mat_zt.createSubMatrix(iset_1, iset_0)
        sub_mat_zt.mult(xt_r.vector, b)

dokken · December 9, 2021, 12:00pm

Support of interpolation between different spaces was introduced in: Interpolation between functions from different spaces by IgorBaratta · Pull Request #1742 · FEniCS/dolfinx · GitHub

You have already made a separate post about this. I do not know about all petsc functions, so I cannot tell you if this will work or not. I would suggest you try to make a minimal example which can verify your approach:
Like:

making a mixed element
solve poisson in upper left corner
extract the sub matrix
compare it to a problem with a non-mixed function space.

wwshunan · December 9, 2021, 12:10pm

Thank you again. I’d like to ask one more question. Does the list map_i contain the indices for the rows and columns for the ith subspace in the final assembled matrix? Does the x_r.sub(0).collapse() mean projecting to the first subspace?

dokken · December 9, 2021, 12:21pm

See: dolfinx/function.py at main · FEniCS/dolfinx · GitHub

wwshunan · December 9, 2021, 1:38pm

Dear dokken,

b_tt = inner(u, v)*dx
b_tz = inner(u, grad(q))*dx
b_zt = inner(grad(p), v)*dx
b_zz = inner(grad(p), grad(q))*dx - k0_sqr*e_r*inner(p, q)*dx
s_b = b_tt + b_zz + b_tz + b_zt

B = assemble_matrix(s_b, [bc])
B.assemble()

mat_zt = assemble_matrix(b_zt, [bc])
mat_zt.assemble()

num_subs = W.num_sub_spaces()
spaces = []
maps = []
for i in range(num_subs):
    space_i, map_i = W.sub(i).collapse(collapsed_dofs=True)
    spaces.append(space_i)
    maps.append(map_i)

for k in range(1000):
      print(B.getValue(maps[1][k], maps[0][k]), mat_zt.getValue(maps[1][k], maps[0][k]))

The print code gave different results. I used this code to print the element values at the position in the left bottom block of the two matrices. They should give the same values, and it’s beyond my expectation.

wwshunan · December 9, 2021, 2:43pm

I reversed the rows and columns, and can get the same values for the two matrix if I exchange them.

Topic		Replies	Views
Wrong max displacement \| Linear Elasticity General	5	261	November 14, 2023
Problems with paraview I/O	2	571	February 2, 2022
Different solution values for different dolfinx version dolfinx	35	2609	October 19, 2021
Modal Analysis of the beam structure	1	780	September 23, 2021
Wrap by vector in paraview is not consistent I/O paraview	3	220	December 30, 2023

Maximum value in vector and paraview is different

Related topics