Sem subdomains

In this example we will see how one can choose a subdomain of the SEM mesh.

We see to utilities of this.

  1. One wishes to write a reduced set of data to reduce storage

  2. One wishes to reduce the scope of an analysis, in turn speeding it up.

Import general modules

mpi4py is always required when using these tools. Numpy is always good to have if any manipulation is to be done.

[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 pysemtools

In this case we will import all the data types that we currently support, as well as io functions that are required to populate them.

[2]:

# Data types
from pysemtools.datatypes.msh import Mesh
from pysemtools.datatypes.coef import Coef
from pysemtools.datatypes.field import FieldRegistry

# Readers
from pysemtools.io.ppymech.neksuite import pynekread

# Writers
from pysemtools.io.ppymech.neksuite import pynekwrite

fname = '../data/rbc0.f00001'

Writing subdomains with a utility function

We have written a utility function that allows to write only a subdomain of a field.

The function works by identidying which element in each rank satisfy a condition and then write those elements. Some aspects to note:

  1. Each rank will check independantly and write its own data. Therefore if the subdomain is all owned by one rank, then only that rank will write.

  2. The subdomains will not be available in memory. If one wants to have the subdomain partitioned across ranks, then one must write and read the subdomain.

In the following example we will write out a subdomain of the input file that is contained between the values x = [-1,1], y=[-1,1], z=[0,0.5].

[3]:

# Utils
from pysemtools.datatypes.utils import write_fld_subdomain_from_list

# Instance the empty objects
msh = Mesh(comm, create_connectivity=False)
fld = FieldRegistry(comm)

# Read the data
pynekread(fname, comm, data_dtype=np.single, msh=msh, fld = fld)

# Write the data in a subdomain and with a different order than what was read
fout = 'subdomains0.f00001'

write_fld_subdomain_from_list(fout, comm, msh, field_list=[fld.registry['u'],fld.registry['v'],fld.registry['w']], subdomain=[[-1, 1], [-1, 1], [0, 0.5]])

2024-08-25 19:39:32,840 - Mesh - INFO - Initializing empty Mesh object.
2024-08-25 19:39:32,840 - Field - INFO - Initializing empty Field object
2024-08-25 19:39:32,841 - pynekread - INFO - Reading file: ../data/rbc0.f00001
2024-08-25 19:39:32,852 - Mesh - INFO - Initializing Mesh object from x,y,z ndarrays.
2024-08-25 19:39:32,853 - Mesh - INFO - Initializing common attributes.
2024-08-25 19:39:32,853 - Mesh - INFO - Mesh object initialized.
2024-08-25 19:39:32,854 - Mesh - INFO - Mesh data is of type: float32
2024-08-25 19:39:32,855 - Mesh - INFO - Elapsed time: 0.00260147s
2024-08-25 19:39:32,855 - pynekread - INFO - Reading field data
2024-08-25 19:39:32,862 - pynekread - INFO - File read
2024-08-25 19:39:32,862 - pynekread - INFO - Elapsed time: 0.021230949s
2024-08-25 19:39:32,872 - Mesh - INFO - Initializing Mesh object from x,y,z ndarrays.
2024-08-25 19:39:32,873 - Mesh - INFO - Initializing common attributes.
2024-08-25 19:39:32,874 - Mesh - INFO - Creating connectivity
2024-08-25 19:39:33,011 - Mesh - INFO - Mesh object initialized.
2024-08-25 19:39:33,012 - Mesh - INFO - Mesh data is of type: float32
2024-08-25 19:39:33,012 - Mesh - INFO - Elapsed time: 0.14000945999999997s
2024-08-25 19:39:33,013 - Field - INFO - Initializing empty Field object
Writing file: subdomains0.f00001

If you visualize the field, you will observe that the subdomain should have been written. Please note that this function might not be very memory efficient.

Obtaining subdomains and partitioning the resuling elements

The naive approach

While the previous method works perfectly fine, there are instances in whihc one might wish to do this in memory.

Choosing the elements that satisfy certain conditions in each rank is very easy and can be seen in the following snippet, where we choose a subdomain that contains values where z is smaller than 0.1 and larger than 0.9 for the rbc data set.

[4]:

# Instance the empty objects
msh = Mesh(comm, create_connectivity=False)
fld = FieldRegistry(comm)

# Read the data
pynekread(fname, comm, data_dtype=np.single, msh=msh, fld = fld)

# Choose a condition that you want the subdomain to satisfy
condition1 = msh.z < 0.1
conidtion2 = msh.z > 0.9
cond = condition1 | conidtion2 # Logical OR
# Get a list of elements that satisfy the condition
condition = np.all([cond], axis=0)
ce = np.unique(np.where(condition)[0])

# Create new object
msh_sub = Mesh(comm, x=msh.x[ce], y=msh.y[ce], z=msh.z[ce], create_connectivity=False)

print(f'Elements in rank {comm.rank}: {msh_sub.nelv}')

2024-08-25 19:39:33,044 - Mesh - INFO - Initializing empty Mesh object.
2024-08-25 19:39:33,045 - Field - INFO - Initializing empty Field object
2024-08-25 19:39:33,046 - pynekread - INFO - Reading file: ../data/rbc0.f00001
2024-08-25 19:39:33,050 - Mesh - INFO - Initializing Mesh object from x,y,z ndarrays.
2024-08-25 19:39:33,051 - Mesh - INFO - Initializing common attributes.
2024-08-25 19:39:33,052 - Mesh - INFO - Mesh object initialized.
2024-08-25 19:39:33,052 - Mesh - INFO - Mesh data is of type: float32
2024-08-25 19:39:33,053 - Mesh - INFO - Elapsed time: 0.0029262059999999313s
2024-08-25 19:39:33,053 - pynekread - INFO - Reading field data
2024-08-25 19:39:33,056 - pynekread - INFO - File read
2024-08-25 19:39:33,057 - pynekread - INFO - Elapsed time: 0.011478025000000058s
2024-08-25 19:39:33,060 - Mesh - INFO - Initializing Mesh object from x,y,z ndarrays.
2024-08-25 19:39:33,060 - Mesh - INFO - Initializing common attributes.
2024-08-25 19:39:33,061 - Mesh - INFO - Mesh object initialized.
2024-08-25 19:39:33,061 - Mesh - INFO - Mesh data is of type: float32
2024-08-25 19:39:33,062 - Mesh - INFO - Elapsed time: 0.002286441999999944s
Elements in rank 0: 120

If you run this code in parallel, you might observe that some ranks will have all the elements while other will not have any. This means that the load is unbalanced and if you wish to operate on the data, many cores might be idling.

Partitoning the elements to keep the load balanced

To fight this we provide a mesh partitoning object that will take charge of selecting the elements that satisfy the condition and then redistribute them to all ranks such tha the load becomes balanced again. The procedure is the following:

[5]:

# Import the mesh partitioner
from pysemtools.datatypes.msh_partitioning import MeshPartitioner

# Instance the empty objects
msh = Mesh(comm, create_connectivity=False)
fld = FieldRegistry(comm)

# Read the data
pynekread(fname, comm, data_dtype=np.single, msh=msh, fld = fld)

# Choose a condition that you want the subdomain to satisfy
condition1 = msh.z < 0.1
conidtion2 = msh.z > 0.9
cond = condition1 | conidtion2 # Logical OR

# Initialize the mesh partitioner with the given condition
mp = MeshPartitioner(comm, msh=msh, conditions=[cond])

# Create the properly partitioned sub mesh and field
partitioned_mesh = mp.create_partitioned_mesh(msh, partitioning_algorithm="load_balanced_linear", create_conectivity=False)
partitioned_field = mp.create_partitioned_field(fld, partitioning_algorithm="load_balanced_linear")

fname = "partitioned_field0.f00001"
pynekwrite(fname, comm, msh=partitioned_mesh, fld=partitioned_field, write_mesh=True)

print(f'Elements in rank {comm.rank}: {partitioned_mesh.nelv}')

2024-08-25 19:39:33,072 - Mesh - INFO - Initializing empty Mesh object.
2024-08-25 19:39:33,072 - Field - INFO - Initializing empty Field object
2024-08-25 19:39:33,073 - pynekread - INFO - Reading file: ../data/rbc0.f00001
2024-08-25 19:39:33,076 - Mesh - INFO - Initializing Mesh object from x,y,z ndarrays.
2024-08-25 19:39:33,077 - Mesh - INFO - Initializing common attributes.
2024-08-25 19:39:33,078 - Mesh - INFO - Mesh object initialized.
2024-08-25 19:39:33,078 - Mesh - INFO - Mesh data is of type: float32
2024-08-25 19:39:33,079 - Mesh - INFO - Elapsed time: 0.002756633000000064s
2024-08-25 19:39:33,079 - pynekread - INFO - Reading field data
2024-08-25 19:39:33,083 - pynekread - INFO - File read
2024-08-25 19:39:33,087 - pynekread - INFO - Elapsed time: 0.014093703000000013s
2024-08-25 19:39:33,089 - Mesh Partitioner - INFO - Initializing Mesh Partitioner
2024-08-25 19:39:33,091 - Mesh - INFO - Initializing Mesh object from x,y,z ndarrays.
2024-08-25 19:39:33,092 - Mesh - INFO - Initializing common attributes.
2024-08-25 19:39:33,093 - Mesh - INFO - Mesh object initialized.
2024-08-25 19:39:33,094 - Mesh - INFO - Mesh data is of type: float32
2024-08-25 19:39:33,094 - Mesh - INFO - Elapsed time: 0.0030915799999999827s
2024-08-25 19:39:33,096 - Mesh Partitioner - INFO - Partitioning the mesh coordinates with load_balanced_linear algorithm
2024-08-25 19:39:33,098 - Mesh Partitioner - INFO - Creating mesh object
2024-08-25 19:39:33,099 - Mesh - INFO - Initializing Mesh object from x,y,z ndarrays.
2024-08-25 19:39:33,099 - Mesh - INFO - Initializing common attributes.
2024-08-25 19:39:33,100 - Mesh - INFO - Mesh object initialized.
2024-08-25 19:39:33,101 - Mesh - INFO - Mesh data is of type: float32
2024-08-25 19:39:33,101 - Mesh - INFO - Elapsed time: 0.002486525999999989s
2024-08-25 19:39:33,102 - Mesh Partitioner - INFO - Partitioning the field object with load_balanced_linear algorithm
2024-08-25 19:39:33,102 - Field - INFO - Initializing empty Field object
2024-08-25 19:39:33,105 - Mesh Partitioner - INFO - done
2024-08-25 19:39:33,107 - pynekwrite - INFO - Writing file: partitioned_field0.f00001
2024-08-25 19:39:33,115 - pynekwrite - INFO - Elapsed time: 0.007852660000000067s
Elements in rank 0: 120

Here we decided to write the partitioned mesh, but one can easily operate on the data. For example, the coef object can be created with the submesh

[6]:

coef_sub = Coef(partitioned_mesh, comm)

2024-08-25 19:39:33,126 - Coef - INFO - Initializing Coef object
2024-08-25 19:39:33,127 - Coef - INFO - Getting derivative matrices
2024-08-25 19:39:33,130 - Coef - INFO - Calculating the components of the jacobian
2024-08-25 19:39:33,145 - Coef - INFO - Calculating the jacobian determinant and inverse of the jacobian matrix
2024-08-25 19:39:33,147 - Coef - INFO - Calculating the mass matrix
2024-08-25 19:39:33,148 - Coef - INFO - Coef object initialized
2024-08-25 19:39:33,148 - Coef - INFO - Coef data is of type: float32
2024-08-25 19:39:33,149 - Coef - INFO - Elapsed time: 0.02330738499999996s

From here on you can do any operation on the partitioned mesh and fields