[BUG] - writing/reading of discontinuous mesh variables

[jcgraciosa]

Describe the bug When reading a discontinuous mesh variable(in this case, dP1 or dP2) using the read_timestep() function, the dimensions of the coordinates and data obtained from the h5 file are incorrect.

What is the bug? When the write() function in discretisation.py is executed (line 1833-1894), it seems that the dimensionality of the discontinuous mesh variable is tied to that of the number of mesh elements (Number of 2-cells per rank). This causes and error to be thrown during read_timestep.

What version code? This issue is found in the development and uw_constants branches

Steps to reproduce: The following can be run the reproduce the issue:

import os
import petsc4py
import underworld3 as uw
from petsc4py import PETSc

res = 4
width = 1.
height = 1.

Pdeg = 1
Pcont = False

outfile = f"test_out"
outdir = f"./test_outdir"

if uw.mpi.rank == 0:
    os.makedirs(f"{outdir}", exist_ok=True)

meshbox = uw.meshing.UnstructuredSimplexBox(
                                                minCoords=(0.0, 0.0), 
                                                maxCoords=(width, height), 
                                                cellSize=1.0 / res, 
                                                regular=False, 
                                                qdegree = 3
                                        )
p_soln = uw.discretisation.MeshVariable("P", meshbox, 1, degree=Pdeg, continuous = Pcont)

meshbox.dm.view()

# check the dimensionality of the mesh variable
with meshbox.access():
    print(p_soln.data.shape)
    print(p_soln.coords.shape)

# saving of discontinous mesh variable
meshbox.write_timestep(
            outfile,
            meshUpdates = True,
            meshVars = [p_soln],
            outputPath = outdir,
            index = 0,
        )

p_soln2 = uw.discretisation.MeshVariable("P2", meshbox, 1, degree=Pdeg, continuous = Pcont) 

# this will result to the error
p_soln2.read_timestep(data_filename = outfile, data_name = "P", index = 0, outputPath = outdir)

This is a portion of the read_timestep function that shows that the data stored in the h5 files now have incorrect dimensions.

import h5py 

output_base_name = os.path.join(outdir, outfile)
data_file = output_base_name + f".mesh.P.{0:05}.h5"
print(data_file)

data_name = "P"

h5f = h5py.File(data_file)

D = h5f["fields"][data_name][()]
X = h5f["fields"]["coordinates"][()]
h5f.close()

print(D.shape) # three columns instead of 1
print(X.shape) # six columns instead of 2

[lmoresi]

I assumed that the petsc writer and h5 reader would match up automatically. You can try the fix I just checked in and see if it works more generally (it seems to fix this code).

[jcgraciosa]

Thanks for the fix! I thought of fixing in the write side, but doing it in the read side makes sense. I've tested the fix and it works. Closing the issue.