# Maximum value in vector and paraview is different

**URL:** <https://fenicsproject.discourse.group/t/maximum-value-in-vector-and-paraview-is-different/7281>\
**Category:** dolfinx\
**Created:** [December 8, 2021, 1:54pm UTC](https://fenicsproject.discourse.group/t/maximum-value-in-vector-and-paraview-is-different/7281 "2021-12-08T13:54:36Z")\
**Posts on this page:** 13\
**Page:** 1

<div class="post-metadata">

**Author:** ![wwshunan](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/wwshunan/32/852_2.png) [@wwshunan](https://fenicsproject.discourse.group/u/wwshunan)\
**Post date:** [December 8, 2021, 1:54pm UTC](https://fenicsproject.discourse.group/t/maximum-value-in-vector-and-paraview-is-different/7281/1 "2021-12-08T13:54:37Z")

</div>

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:

```

```auto
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.

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [December 9, 2021, 10:09am UTC](https://fenicsproject.discourse.group/t/maximum-value-in-vector-and-paraview-is-different/7281/2 "2021-12-09T10:09:14Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![wwshunan](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/wwshunan/32/852_2.png) [@wwshunan](https://fenicsproject.discourse.group/u/wwshunan)\
**Post date:** [December 9, 2021, 10:16am UTC](https://fenicsproject.discourse.group/t/maximum-value-in-vector-and-paraview-is-different/7281/3 "2021-12-09T10:16:04Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![wwshunan](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/wwshunan/32/852_2.png) [@wwshunan](https://fenicsproject.discourse.group/u/wwshunan)\
**Post date:** [December 9, 2021, 10:26am UTC](https://fenicsproject.discourse.group/t/maximum-value-in-vector-and-paraview-is-different/7281/4 "2021-12-09T10:26:31Z")

</div>

Dear dokken,  
Here is my code

```auto
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)

```

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [December 9, 2021, 10:29am UTC](https://fenicsproject.discourse.group/t/maximum-value-in-vector-and-paraview-is-different/7281/5 "2021-12-09T10:29:27Z")

</div>

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](https://github.com/FEniCS/dolfinx/blob/main/python/test/unit/io/test_adios2.py#L224)

---

<div class="post-metadata">

**Author:** ![wwshunan](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/wwshunan/32/852_2.png) [@wwshunan](https://fenicsproject.discourse.group/u/wwshunan)\
**Post date:** [December 9, 2021, 10:56am UTC](https://fenicsproject.discourse.group/t/maximum-value-in-vector-and-paraview-is-different/7281/6 "2021-12-09T10:56:57Z")

</div>

Dear dokken,  
I cannot import VTXWriter from [dolfinx.cpp.io](http://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.

```auto
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.

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [December 9, 2021, 11:33am UTC](https://fenicsproject.discourse.group/t/maximum-value-in-vector-and-paraview-is-different/7281/7 "2021-12-09T11:33:19Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![wwshunan](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/wwshunan/32/852_2.png) [@wwshunan](https://fenicsproject.discourse.group/u/wwshunan)\
**Post date:** [December 9, 2021, 11:48am UTC](https://fenicsproject.discourse.group/t/maximum-value-in-vector-and-paraview-is-different/7281/8 "2021-12-09T11:48:27Z")

</div>

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?

```auto
        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)

```

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [December 9, 2021, 12:00pm UTC](https://fenicsproject.discourse.group/t/maximum-value-in-vector-and-paraview-is-different/7281/9 "2021-12-09T12:00:39Z")

</div>

Support of interpolation between different spaces was introduced in: [Interpolation between functions from different spaces by IgorBaratta · Pull Request #1742 · FEniCS/dolfinx · GitHub](https://github.com/FEniCS/dolfinx/pull/1742)

> [@wwshunan](#):
>
> , 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?

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:

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

---

<div class="post-metadata">

**Author:** ![wwshunan](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/wwshunan/32/852_2.png) [@wwshunan](https://fenicsproject.discourse.group/u/wwshunan)\
**Post date:** [December 9, 2021, 12:10pm UTC](https://fenicsproject.discourse.group/t/maximum-value-in-vector-and-paraview-is-different/7281/10 "2021-12-09T12:10:45Z")

</div>

> [@wwshunan](#):
>
> `space_i, map_i = W.sub(i).collapse(collapsed_dofs=True)`

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?

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [December 9, 2021, 12:21pm UTC](https://fenicsproject.discourse.group/t/maximum-value-in-vector-and-paraview-is-different/7281/11 "2021-12-09T12:21:12Z")

</div>

See: [dolfinx/function.py at main · FEniCS/dolfinx · GitHub](https://github.com/FEniCS/dolfinx/blob/main/python/dolfinx/fem/function.py#L519-L530)

---

<div class="post-metadata">

**Author:** ![wwshunan](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/wwshunan/32/852_2.png) [@wwshunan](https://fenicsproject.discourse.group/u/wwshunan)\
**Post date:** [December 9, 2021, 1:38pm UTC](https://fenicsproject.discourse.group/t/maximum-value-in-vector-and-paraview-is-different/7281/12 "2021-12-09T13:38:09Z")

</div>

Dear dokken,

```auto
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.

---

<div class="post-metadata">

**Author:** ![wwshunan](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/wwshunan/32/852_2.png) [@wwshunan](https://fenicsproject.discourse.group/u/wwshunan)\
**Post date:** [December 9, 2021, 2:43pm UTC](https://fenicsproject.discourse.group/t/maximum-value-in-vector-and-paraview-is-different/7281/13 "2021-12-09T14:43:30Z")

</div>

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