Element distribution interpolation

In this example we show how to perform interpolations that keep the same element structure.

Two examples of this are:

  1. Mapping the GLL points of each element into a different distribution.

  2. Perfomring P refinement, i.e., map to points of the same distribution but a different degree.

Import general modules

[1]:

# Import required modules
from mpi4py import MPI #equivalent to the use of MPI_init() in C
import matplotlib.pyplot as plt
import numpy as np

# Get mpi info
comm = MPI.COMM_WORLD

Import modules from pynek

[2]:

from pysemtools.io.ppymech.neksuite import pynekread, pynekwrite
from pysemtools.datatypes.msh import Mesh
from pysemtools.datatypes.field import FieldRegistry

Read the data and build objects

In this instance, we create connectivity for the mesh object, given that we wish to use direct stiffness summation to reduce discontinuities.

[3]:

msh = Mesh(comm, create_connectivity=True)
fld = FieldRegistry(comm)
pynekread('../data/tc_channel0.f00001', comm, data_dtype=np.double, msh=msh, fld=fld)

2024-09-11 21:48:29,145 - Mesh - INFO - Initializing empty Mesh object.
2024-09-11 21:48:29,147 - Field - INFO - Initializing empty Field object
2024-09-11 21:48:29,147 - pynekread - INFO - Reading file: ../data/tc_channel0.f00001
2024-09-11 21:48:29,163 - Mesh - INFO - Initializing Mesh object from x,y,z ndarrays.
2024-09-11 21:48:29,165 - Mesh - INFO - Initializing common attributes.
2024-09-11 21:48:29,167 - Mesh - INFO - Getting vertices
2024-09-11 21:48:29,168 - Mesh - INFO - Getting vertices
2024-09-11 21:48:29,177 - Mesh - INFO - Getting facet centers
2024-09-11 21:48:29,187 - Mesh - INFO - Creating connectivity
2024-09-11 21:48:29,575 - Mesh - INFO - Mesh object initialized.
2024-09-11 21:48:29,576 - Mesh - INFO - Mesh data is of type: float64
2024-09-11 21:48:29,577 - Mesh - INFO - Elapsed time: 0.41460768800000003s
2024-09-11 21:48:29,578 - pynekread - INFO - Reading field data
2024-09-11 21:48:29,586 - pynekread - INFO - File read
2024-09-11 21:48:29,587 - pynekread - INFO - Elapsed time: 0.439429355s

Mapping to a new distribution

We can map to a new distribution of points in any of the dimensions of the element at the same time.

Here we show an example where this is done to an equidistant mesh in the z direction.

For this we use the PMapper class.

[4]:

# Import the module
from pysemtools.interpolation.mesh_to_mesh import PMapper

# Here specify that it should be equal in the 3 direction
mapper = PMapper(n=msh.lx, distribution=['GLL', 'GLL', 'EQ'])

# Create the mesh with the new distribution
eq_msh = mapper.create_mapped_mesh(comm, msh=msh)

# Interpolate the fields. They are passed as a list and returned as a list
mapped_fields = mapper.interpolate_from_field_list(comm, field_list=[fld.registry['u'],fld.registry['v'],fld.registry['w']])

# Create an empty field object for the equidistant fields
eq_fld = FieldRegistry(comm)
eq_fld.add_field(comm, field_name='u', field=mapped_fields[0], dtype = np.double)
eq_fld.add_field(comm, field_name='v', field=mapped_fields[1], dtype = np.double)
eq_fld.add_field(comm, field_name='w', field=mapped_fields[2], dtype = np.double)

# Write the new mesh and fields
fname = "mappedfield0.f00001"
pynekwrite(fname, comm, msh=eq_msh, fld=eq_fld, write_mesh=True, wdsz=4)

2024-09-11 21:48:30,619 - Mesh - INFO - Initializing Mesh object from x,y,z ndarrays.
2024-09-11 21:48:30,619 - Mesh - INFO - Initializing common attributes.
2024-09-11 21:48:30,620 - Mesh - INFO - Getting vertices
2024-09-11 21:48:30,621 - Mesh - INFO - Getting vertices
2024-09-11 21:48:30,630 - Mesh - INFO - Getting facet centers
2024-09-11 21:48:30,637 - Mesh - INFO - Creating connectivity
2024-09-11 21:48:30,986 - Mesh - INFO - Mesh object initialized.
2024-09-11 21:48:30,987 - Mesh - INFO - Mesh data is of type: float64
2024-09-11 21:48:30,987 - Mesh - INFO - Elapsed time: 0.36897678200000006s
2024-09-11 21:48:31,746 - Field - INFO - Initializing empty Field object
2024-09-11 21:48:31,747 - pynekwrite - INFO - Writing file: mappedfield0.f00001
2024-09-11 21:48:31,756 - pynekwrite - INFO - Elapsed time: 0.008410899999999888s

Mapping all fields

In the previous block, we showed how to map specific fields. However if you wish to map all the contents of the field object, you can use the following:

[5]:

eq_fld2 = mapper.create_mapped_field(comm, fld=fld)

for key in eq_fld.registry.keys():
    print(np.allclose(eq_fld.registry[key].data, eq_fld2.registry[key].data))

2024-09-11 21:48:31,766 - Field - INFO - Initializing empty Field object
True
True
True

Performing P refinement

In some instances one wishes to keep the GLL distribution to perform integration or derivation but doing so at higher or lower polynomiar order. For this cases one can employ p refinement/coarsening. We show now how that can be done.

[6]:

# Import the module
from pysemtools.interpolation.mesh_to_mesh import PRefiner

# Here specify that it should be equal in the 3 direction
pr = PRefiner(n_old = msh.lx, n_new = 3)

# Create the mesh with the new distribution
r_msh = pr.create_refined_mesh(comm, msh=msh)

# Interpolate the fields. They are passed as a list and returned as a list
r_fields = pr.interpolate_from_field_list(comm, field_list=[fld.registry['u'],fld.registry['v'],fld.registry['w']])

# Create an empty field object for the equidistant fields
r_fld = FieldRegistry(comm)
r_fld.add_field(comm, field_name='u', field=r_fields[0], dtype = np.double)
r_fld.add_field(comm, field_name='v', field=r_fields[1], dtype = np.double)
r_fld.add_field(comm, field_name='w', field=r_fields[2], dtype = np.double)

# Write the new mesh and fields
fname = "refinedfield0.f00001"
pynekwrite(fname, comm, msh=r_msh, fld=r_fld, write_mesh=True, wdsz=4)

2024-09-11 21:48:33,932 - Mesh - INFO - Initializing Mesh object from x,y,z ndarrays.
2024-09-11 21:48:33,933 - Mesh - INFO - Initializing common attributes.
2024-09-11 21:48:33,933 - Mesh - INFO - Getting vertices
2024-09-11 21:48:33,935 - Mesh - INFO - Getting vertices
2024-09-11 21:48:33,940 - Mesh - INFO - Getting facet centers
2024-09-11 21:48:33,943 - Mesh - INFO - Creating connectivity
2024-09-11 21:48:34,030 - Mesh - INFO - Mesh object initialized.
2024-09-11 21:48:34,031 - Mesh - INFO - Mesh data is of type: float64
2024-09-11 21:48:34,033 - Mesh - INFO - Elapsed time: 0.10139660900000003s
2024-09-11 21:48:34,922 - Field - INFO - Initializing empty Field object
2024-09-11 21:48:34,923 - pynekwrite - INFO - Writing file: refinedfield0.f00001
2024-09-11 21:48:34,926 - pynekwrite - INFO - Elapsed time: 0.002947181000000576s

Refining all fields

In the previous block, we showed how to refine specific fields. However if you wish to refine all the contents of the field object, you can use the following:

[7]:

r_fld2 = pr.create_refined_field(comm, fld=fld)

for key in r_fld.registry.keys():
    print(np.allclose(r_fld.registry[key].data, r_fld2.registry[key].data))

2024-09-11 21:48:34,935 - Field - INFO - Initializing empty Field object
True
True
True