Manipulation of real & imag, writing back to vtk

dolfinx
1

Hello again,

sorry for flooding this board, but I’m not able to calculate a new measure based on the real and the imaginary part of a solution ‘uh’ (using dolfinx-complex-mode)

problem = fem.LinearProblem(a, L, u=uh)
problem.solve()

I want to compute, for example,

uh_new = uh.real + 2*uh.imag

This new measure is real-valued only and shall be written back to a vtk-file. But:

  1. I’m failing in getting the real and the imaginary part of ‘uh’
  2. I’m failing in interpolating/projecting uh_new to a (new) FunctionSpace
  3. I’m failing in storing ‘uh_new’ in a vtk-file like
dolfinx.io.VTKFile(MPI.COMM_WORLD, "uh_new.pvd", "w").write(uh_new)

A lot of deficiencies by myself, I know.

Thank’s in advance

2

Please note that you can do this post-processing using the Calculator filter in Paraview.
Press the Filters->Alphabetical->Calculator, and then in the calculator, choose Scalars->"name_of_field_real"
and then add
2 * `Scalars->“name_of_field_complex”.
For instance consider this minimal example:

from dolfinx import UnitSquareMesh, Function, FunctionSpace
from mpi4py import MPI
from dolfinx.io import VTKFile

mesh = UnitSquareMesh(MPI.COMM_WORLD, 16, 16)
u = Function(FunctionSpace(mesh, ("Lagrange", 1)))


def f(x):
    return x[0] + 2j * x[1]


u.interpolate(f)

with VTKFile(mesh.mpi_comm(), "u.pvd", "w") as file:
    file.write_mesh(mesh)
    file.write_function([u._cpp_object])

I open u.pvd in Paraview, then use the Extract Block filter to get the functions, then the Calculator filter with the input: real_f_2 + 2 * imag_f_2 to get output. Note that you could also do the following:

import numpy as np
v = Function(u.function_space)
v.x.array[:] = np.real(u.x.array) + 2 * np.imag(u.x.array)
with VTKFile(mesh.mpi_comm(), "u.pvd", "w") as file:
    file.write_mesh(mesh)
    file.write_function([u._cpp_object, v._cpp_object])

and then you can compare the result of the Paraview calculator with the value in real_f_3

3

@dokken, many thanks. The second option points me in the right direction. VTK-export works quite well but opening these files in paraview may cause problems for 2d-meshes ('Cell type UnknownClass(69) is not a 3D cell).
However, I combined your approach with XDMF-export, and everything works fine. For your background: I want to visualize the complex-valued results in the time domain. My code snippet is:

v = Function(uh.function_space)
v.name = "uh_t"
nums = 20
steps = np.linspace(0., 2*np.pi, num=nums) 
with XDMFFile(MPI.COMM_WORLD, "./uh_time.xdmf", "w",
              encoding=XDMFFile.Encoding.HDF5) as file:
    file.write_mesh(msh)
    for ii in range(nums):
        v.x.array[:] = np.real(uh.x.array)*np.cos(steps[ii])-np.imag(uh.x.array)*np.sin(steps[ii])
        file.write_function(v, steps[ii])
4

The error you get might be due to the version of Paraview you are using (as the example above is a 2D mesh, and I used 5.9.x without getting any errors).

5

You are right. Paraview v5.9 doesn’t show any problems (I used v5.4.1).