NodePartitionedMeshReader.cpp

/*!
  \file
  \author Wenqing Wang
  \date   2014.08
  \brief  Define members of class NodePartitionedMeshReader to read node-wise
  partitioned mesh with MPI functions.

  \copyright
  Copyright (c) 2012-2023, OpenGeoSys Community (http://www.opengeosys.org)
             Distributed under a Modified BSD License.
               See accompanying file LICENSE.txt or
               http://www.opengeosys.org/project/license

 */

#include "NodePartitionedMeshReader.h"

#include <mpi.h>

#include <fstream>
#include <numeric>

#include "BaseLib/FileTools.h"
#include "BaseLib/Logging.h"
#include "BaseLib/RunTime.h"
#include "MeshLib/Elements/Elements.h"
#include "MeshLib/MeshEnums.h"
#include "MeshLib/Properties.h"

// Check if the value can by converted to given type without overflow.
template <typename VALUE, typename TYPE>
bool is_safely_convertable(VALUE const& value)
{
    bool const result = value <= std::numeric_limits<TYPE>::max();
    if (!result)
    {
        ERR("The value {:d} is too large for conversion.", value);
        ERR("Maximum available size is {:d}.",
            std::numeric_limits<TYPE>::max());
    }
    return result;
}

namespace MeshLib
{
namespace IO
{
NodePartitionedMeshReader::NodePartitionedMeshReader(MPI_Comm comm)
    : _mpi_comm(comm)
{
    MPI_Comm_size(_mpi_comm, &_mpi_comm_size);
    MPI_Comm_rank(_mpi_comm, &_mpi_rank);

    registerNodeDataMpiType();
}

NodePartitionedMeshReader::~NodePartitionedMeshReader()
{
    MPI_Type_free(&_mpi_node_type);
}

void NodePartitionedMeshReader::registerNodeDataMpiType()
{
    int const count = 2;
    int blocks[count] = {1, 3};
    MPI_Datatype types[count] = {MPI_UNSIGNED_LONG, MPI_DOUBLE};
    MPI_Aint displacements[count] = {0, sizeof(NodeData::index)};

    MPI_Type_create_struct(count, blocks, displacements, types,
                           &_mpi_node_type);
    MPI_Type_commit(&_mpi_node_type);
}

MeshLib::NodePartitionedMesh* NodePartitionedMeshReader::read(
    const std::string& file_name_base)
{
    BaseLib::RunTime timer;
    timer.start();

    // Always try binary file first
    std::string const fname_new = file_name_base + "_partitioned_msh_cfg" +
                                  std::to_string(_mpi_comm_size) + ".bin";

    if (!BaseLib::IsFileExisting(fname_new))  // binary file does not exist.
    {
        std::string const fname_ascii = file_name_base +
                                        "_partitioned_msh_cfg" +
                                        std::to_string(_mpi_comm_size) + ".msh";
        if (BaseLib::IsFileExisting(fname_ascii))
        {
            ERR("Reading of ASCII mesh file {:s} is not supported since OGS "
                "version 6.3.3.",
                fname_ascii);
        }
        OGS_FATAL("Required binary file {:s} does not exist.\n", fname_new);
    }

    INFO("Reading corresponding part of mesh data from binary file {:s} ...",
         file_name_base);

    MeshLib::NodePartitionedMesh* mesh = readMesh(file_name_base);

    INFO("[time] Reading the mesh took {:f} s.", timer.elapsed());

    MPI_Barrier(_mpi_comm);

    return mesh;
}

template <typename DATA>
bool NodePartitionedMeshReader::readDataFromFile(std::string const& filename,
                                                 MPI_Offset offset,
                                                 MPI_Datatype type,
                                                 DATA& data) const
{
    // Check container size
    if (!is_safely_convertable<std::size_t, int>(data.size()))
    {
        ERR("The container size is too large for MPI_File_read() call.");
        return false;
    }

    // Open file
    MPI_File file;

    char* filename_char = const_cast<char*>(filename.data());
    int const file_status = MPI_File_open(
        _mpi_comm, filename_char, MPI_MODE_RDONLY, MPI_INFO_NULL, &file);

    if (file_status != 0)
    {
        ERR("Error opening file {:s}. MPI error code {:d}", filename,
            file_status);
        return false;
    }

    // Read data
    char file_mode[] = "native";
    MPI_File_set_view(file, offset, type, type, file_mode, MPI_INFO_NULL);
    // The static cast is checked above.
    MPI_File_read(file, data.data(), static_cast<int>(data.size()), type,
                  MPI_STATUS_IGNORE);
    MPI_File_close(&file);

    return true;
}

MeshLib::NodePartitionedMesh* NodePartitionedMeshReader::readMesh(
    const std::string& file_name_base)
{
    //----------------------------------------------------------------------------------
    // Read headers
    const std::string fname_header = file_name_base + "_partitioned_msh_";
    const std::string fname_num_p_ext = std::to_string(_mpi_comm_size) + ".bin";

    // Read the config meta data from *cfg* file into struct PartitionedMeshInfo
    // _mesh_info
    if (!readDataFromFile(
            fname_header + "cfg" + fname_num_p_ext,
            static_cast<MPI_Offset>(static_cast<unsigned>(_mpi_rank) *
                                    sizeof(_mesh_info)),
            MPI_LONG, _mesh_info))
        return nullptr;

    //----------------------------------------------------------------------------------
    // Read Nodes
    std::vector<NodeData> nodes(_mesh_info.number_of_nodes);

    if (!readDataFromFile(fname_header + "nod" + fname_num_p_ext,
                          static_cast<MPI_Offset>(_mesh_info.offset[2]),
                          _mpi_node_type, nodes))
        return nullptr;

    std::vector<MeshLib::Node*> mesh_nodes;
    std::vector<unsigned long> glb_node_ids;
    setNodes(nodes, mesh_nodes, glb_node_ids);

    //----------------------------------------------------------------------------------
    // Read non-ghost elements
    std::vector<unsigned long> elem_data(_mesh_info.number_of_regular_elements +
                                         _mesh_info.offset[0]);
    if (!readDataFromFile(fname_header + "ele" + fname_num_p_ext,
                          static_cast<MPI_Offset>(_mesh_info.offset[3]),
                          MPI_LONG, elem_data))
        return nullptr;

    std::vector<MeshLib::Element*> mesh_elems(
        _mesh_info.number_of_regular_elements +
        _mesh_info.number_of_ghost_elements);
    setElements(mesh_nodes, elem_data, mesh_elems);

    //----------------------------------------------------------------------------------
    // Read ghost element
    std::vector<unsigned long> ghost_elem_data(
        _mesh_info.number_of_ghost_elements + _mesh_info.offset[1]);

    if (!readDataFromFile(fname_header + "ele_g" + fname_num_p_ext,
                          static_cast<MPI_Offset>(_mesh_info.offset[4]),
                          MPI_LONG, ghost_elem_data))
        return nullptr;

    const bool process_ghost = true;
    setElements(mesh_nodes, ghost_elem_data, mesh_elems, process_ghost);

    //----------------------------------------------------------------------------------
    // read the properties
    MeshLib::Properties p(readProperties(file_name_base));

    return newMesh(BaseLib::extractBaseName(file_name_base), mesh_nodes,
                   glb_node_ids, mesh_elems, p);
}

MeshLib::Properties NodePartitionedMeshReader::readProperties(
    const std::string& file_name_base) const
{
    MeshLib::Properties p;
    readProperties(file_name_base, MeshLib::MeshItemType::Node, p);
    readProperties(file_name_base, MeshLib::MeshItemType::Cell, p);
    readProperties(file_name_base, MeshLib::MeshItemType::IntegrationPoint, p);
    return p;
}

void NodePartitionedMeshReader::readProperties(
    const std::string& file_name_base, MeshLib::MeshItemType t,
    MeshLib::Properties& p) const
{
    const std::string item_type = MeshLib::toString(t);

    const std::string fname_cfg = file_name_base + "_partitioned_" + item_type +
                                  "_properties_cfg" +
                                  std::to_string(_mpi_comm_size) + ".bin";
    std::ifstream is(fname_cfg.c_str(), std::ios::binary | std::ios::in);
    if (!is)
    {
        WARN(
            "Could not open file '{:s}'.\n"
            "\tYou can ignore this warning if the mesh does not contain {:s}-"
            "wise property data.",
            fname_cfg, item_type.data());
        return;
    }
    std::size_t number_of_properties = 0;
    is.read(reinterpret_cast<char*>(&number_of_properties),
            sizeof(std::size_t));
    std::vector<std::optional<MeshLib::IO::PropertyVectorMetaData>> vec_pvmd(
        number_of_properties);
    for (std::size_t i(0); i < number_of_properties; ++i)
    {
        vec_pvmd[i] = MeshLib::IO::readPropertyVectorMetaData(is);
        if (!vec_pvmd[i])
        {
            OGS_FATAL(
                "Error in NodePartitionedMeshReader::readProperties: "
                "Could not read the meta data for the PropertyVector {:d}",
                i);
        }
    }
    for (std::size_t i(0); i < number_of_properties; ++i)
    {
        MeshLib::IO::writePropertyVectorMetaData(*(vec_pvmd[i]));
    }
    auto pos = is.tellg();
    auto offset =
        static_cast<long>(pos) +
        static_cast<long>(_mpi_rank *
                          sizeof(MeshLib::IO::PropertyVectorPartitionMetaData));
    is.seekg(offset);
    std::optional<MeshLib::IO::PropertyVectorPartitionMetaData> pvpmd(
        MeshLib::IO::readPropertyVectorPartitionMetaData(is));
    bool pvpmd_read_ok = static_cast<bool>(pvpmd);
    bool all_pvpmd_read_ok;
    MPI_Allreduce(&pvpmd_read_ok, &all_pvpmd_read_ok, 1, MPI_C_BOOL, MPI_LOR,
                  _mpi_comm);
    if (!all_pvpmd_read_ok)
    {
        OGS_FATAL(
            "Error in NodePartitionedMeshReader::readProperties: "
            "Could not read the partition meta data for the mpi process {:d}",
            _mpi_rank);
    }
    DBUG("offset in the PropertyVector: {:d}", pvpmd->offset);
    DBUG("{:d} tuples in partition.", pvpmd->number_of_tuples);
    is.close();

    const std::string fname_val = file_name_base + "_partitioned_" + item_type +
                                  "_properties_val" +
                                  std::to_string(_mpi_comm_size) + ".bin";
    is.open(fname_val.c_str(), std::ios::binary | std::ios::in);
    if (!is)
    {
        ERR("Could not open file '{:s}'\n."
            "\tYou can ignore this warning if the mesh does not contain {:s}-"
            "wise property data.",
            fname_val, item_type.data());
    }

    readDomainSpecificPartOfPropertyVectors(vec_pvmd, *pvpmd, t, is, p);
}

void NodePartitionedMeshReader::readDomainSpecificPartOfPropertyVectors(
    std::vector<std::optional<MeshLib::IO::PropertyVectorMetaData>> const&
        vec_pvmd,
    MeshLib::IO::PropertyVectorPartitionMetaData const& pvpmd,
    MeshLib::MeshItemType t,
    std::istream& is,
    MeshLib::Properties& p) const
{
    unsigned long global_offset = 0;
    std::size_t const number_of_properties = vec_pvmd.size();
    for (std::size_t i(0); i < number_of_properties; ++i)
    {
        DBUG("global offset: {:d}, offset within the PropertyVector: {:d}.",
             global_offset,
             global_offset + pvpmd.offset * vec_pvmd[i]->number_of_components *
                                 vec_pvmd[i]->data_type_size_in_bytes);

        // Special field data such as OGS_VERSION, IntegrationPointMetaData,
        // etc., which are not "real" integration points, are copied "as is"
        // (i.e. fully) for every partition.
        if (vec_pvmd[i]->property_name.find("_ip") == std::string::npos &&
            t == MeshLib::MeshItemType::IntegrationPoint)
        {
            createSpecificPropertyVectorPart<char>(is, *vec_pvmd[i], t,
                                                   global_offset, p);

            global_offset += vec_pvmd[i]->data_type_size_in_bytes *
                             vec_pvmd[i]->number_of_tuples;
            continue;
        }

        if (vec_pvmd[i]->is_int_type)
        {
            if (vec_pvmd[i]->is_data_type_signed)
            {
                if (vec_pvmd[i]->data_type_size_in_bytes == sizeof(char))
                {
                    createPropertyVectorPart<char>(is, *vec_pvmd[i], pvpmd, t,
                                                   global_offset, p);
                }
                else if (vec_pvmd[i]->data_type_size_in_bytes == sizeof(int))
                {
                    createPropertyVectorPart<int>(is, *vec_pvmd[i], pvpmd, t,
                                                  global_offset, p);
                }
                else if (vec_pvmd[i]->data_type_size_in_bytes == sizeof(long))
                    createPropertyVectorPart<long>(is, *vec_pvmd[i], pvpmd, t,
                                                   global_offset, p);
                else
                {
                    WARN(
                        "Implementation for reading signed integer property "
                        "vector '{:s}' is not available.",
                        vec_pvmd[i]->property_name);
                }
            }
            else
            {
                if (vec_pvmd[i]->data_type_size_in_bytes ==
                    sizeof(unsigned char))
                {
                    createPropertyVectorPart<unsigned char>(
                        is, *vec_pvmd[i], pvpmd, t, global_offset, p);
                }
                else if (vec_pvmd[i]->data_type_size_in_bytes ==
                         sizeof(unsigned int))
                    createPropertyVectorPart<unsigned int>(
                        is, *vec_pvmd[i], pvpmd, t, global_offset, p);
                else if (vec_pvmd[i]->data_type_size_in_bytes ==
                         sizeof(unsigned long))
                    createPropertyVectorPart<unsigned long>(
                        is, *vec_pvmd[i], pvpmd, t, global_offset, p);
                else
                {
                    WARN(
                        "Implementation for reading unsigned property vector "
                        "'{:s}' is not available.",
                        vec_pvmd[i]->property_name);
                }
            }
        }
        else
        {
            if (vec_pvmd[i]->data_type_size_in_bytes == sizeof(float))
                createPropertyVectorPart<float>(is, *vec_pvmd[i], pvpmd, t,
                                                global_offset, p);
            else if (vec_pvmd[i]->data_type_size_in_bytes == sizeof(double))
                createPropertyVectorPart<double>(is, *vec_pvmd[i], pvpmd, t,
                                                 global_offset, p);
            else
            {
                WARN(
                    "Implementation for reading floating point property vector "
                    "'{:s}' is not available.",
                    vec_pvmd[i]->property_name);
            }
        }
        global_offset += vec_pvmd[i]->data_type_size_in_bytes *
                         vec_pvmd[i]->number_of_tuples *
                         vec_pvmd[i]->number_of_components;
    }
}

MeshLib::NodePartitionedMesh* NodePartitionedMeshReader::newMesh(
    std::string const& mesh_name,
    std::vector<MeshLib::Node*> const& mesh_nodes,
    std::vector<unsigned long> const& glb_node_ids,
    std::vector<MeshLib::Element*> const& mesh_elems,
    MeshLib::Properties const& properties) const
{
    std::vector<std::size_t> gathered_n_regular_base_nodes(_mpi_comm_size);

    MPI_Allgather(&_mesh_info.number_of_regular_base_nodes,
                  1,
                  MPI_UNSIGNED_LONG,
                  gathered_n_regular_base_nodes.data(),
                  1,
                  MPI_UNSIGNED_LONG,
                  _mpi_comm);

    std::vector<std::size_t> n_regular_base_nodes_at_rank;
    n_regular_base_nodes_at_rank.push_back(0);
    std::partial_sum(begin(gathered_n_regular_base_nodes),
                     end(gathered_n_regular_base_nodes),
                     back_inserter(n_regular_base_nodes_at_rank));

    std::vector<std::size_t> gathered_n_regular_high_order_nodes(
        _mpi_comm_size);
    std::size_t const n_regular_high_order_nodes =
        _mesh_info.number_of_regular_nodes -
        _mesh_info.number_of_regular_base_nodes;
    MPI_Allgather(&n_regular_high_order_nodes,
                  1,
                  MPI_UNSIGNED_LONG,
                  gathered_n_regular_high_order_nodes.data(),
                  1,
                  MPI_UNSIGNED_LONG,
                  _mpi_comm);

    std::vector<std::size_t> n_regular_high_order_nodes_at_rank;
    n_regular_high_order_nodes_at_rank.push_back(0);
    std::partial_sum(begin(gathered_n_regular_high_order_nodes),
                     end(gathered_n_regular_high_order_nodes),
                     back_inserter(n_regular_high_order_nodes_at_rank));

    return new MeshLib::NodePartitionedMesh(
        mesh_name, mesh_nodes, glb_node_ids, mesh_elems, properties,
        _mesh_info.number_of_global_base_nodes,
        _mesh_info.number_of_global_nodes, _mesh_info.number_of_regular_nodes,
        std::move(n_regular_base_nodes_at_rank),
        std::move(n_regular_high_order_nodes_at_rank));
}

void NodePartitionedMeshReader::setNodes(
    const std::vector<NodeData>& node_data,
    std::vector<MeshLib::Node*>& mesh_node,
    std::vector<unsigned long>& glb_node_ids) const
{
    mesh_node.resize(_mesh_info.number_of_nodes);
    glb_node_ids.resize(_mesh_info.number_of_nodes);

    for (std::size_t i = 0; i < mesh_node.size(); i++)
    {
        NodeData const& nd = node_data[i];
        glb_node_ids[i] = nd.index;
        mesh_node[i] = new MeshLib::Node(nd.x, nd.y, nd.z, i);
    }
}

void NodePartitionedMeshReader::setElements(
    const std::vector<MeshLib::Node*>& mesh_nodes,
    const std::vector<unsigned long>& elem_data,
    std::vector<MeshLib::Element*>& mesh_elems, const bool ghost) const
{
    // Number of elements, either ghost or regular
    unsigned long const ne = ghost ? _mesh_info.number_of_ghost_elements
                                   : _mesh_info.number_of_regular_elements;
    unsigned long const id_offset_ghost =
        ghost ? _mesh_info.number_of_regular_elements : 0;

    for (unsigned long i = 0; i < ne; i++)
    {
        unsigned long id_offset_elem = elem_data[i];

        // Unused for now, keep for elem_data documentation purpose here.
        {
            const unsigned mat_idx =
                static_cast<unsigned>(elem_data[id_offset_elem++]);
            (void)mat_idx;
        }
        const unsigned long e_type = elem_data[id_offset_elem++];
        unsigned long const nnodes = elem_data[id_offset_elem++];

        MeshLib::Node** elem_nodes = new MeshLib::Node*[nnodes];
        for (unsigned long k = 0; k < nnodes; k++)
            elem_nodes[k] = mesh_nodes[elem_data[id_offset_elem++]];

        // The element types below are defined by the MeshLib::CellType.
        switch (static_cast<CellType>(e_type))
        {
            case CellType::POINT1:
                mesh_elems[i + id_offset_ghost] =
                    new MeshLib::Point(elem_nodes);
                break;
            case CellType::LINE2:
                mesh_elems[i + id_offset_ghost] = new MeshLib::Line(elem_nodes);
                break;
            case CellType::LINE3:
                mesh_elems[i + id_offset_ghost] =
                    new MeshLib::Line3(elem_nodes);
                break;
            case CellType::QUAD4:
                mesh_elems[i + id_offset_ghost] = new MeshLib::Quad(elem_nodes);
                break;
            case CellType::QUAD8:
                mesh_elems[i + id_offset_ghost] =
                    new MeshLib::Quad8(elem_nodes);
                break;
            case CellType::QUAD9:
                mesh_elems[i + id_offset_ghost] =
                    new MeshLib::Quad9(elem_nodes);
                break;
            case CellType::HEX8:
                mesh_elems[i + id_offset_ghost] = new MeshLib::Hex(elem_nodes);
                break;
            case CellType::HEX20:
                mesh_elems[i + id_offset_ghost] =
                    new MeshLib::Hex20(elem_nodes);
                break;
            case CellType::HEX27:
                OGS_FATAL(
                    "NodePartitionedMeshReader: construction of HEX27 element "
                    "with id {:d} is not implemented.",
                    i);
                break;
            case CellType::TRI3:
                mesh_elems[i + id_offset_ghost] = new MeshLib::Tri(elem_nodes);
                break;
            case CellType::TRI6:
                mesh_elems[i + id_offset_ghost] = new MeshLib::Tri6(elem_nodes);
                break;
            case CellType::TET4:
                mesh_elems[i + id_offset_ghost] = new MeshLib::Tet(elem_nodes);
                break;
            case CellType::TET10:
                mesh_elems[i + id_offset_ghost] =
                    new MeshLib::Tet10(elem_nodes);
                break;
            case CellType::PRISM6:
                mesh_elems[i + id_offset_ghost] =
                    new MeshLib::Prism(elem_nodes);
                break;
            case CellType::PRISM15:
                mesh_elems[i + id_offset_ghost] =
                    new MeshLib::Prism15(elem_nodes);
                break;
            case CellType::PYRAMID5:
                mesh_elems[i + id_offset_ghost] =
                    new MeshLib::Pyramid(elem_nodes);
                break;
            case CellType::PYRAMID13:
                mesh_elems[i + id_offset_ghost] =
                    new MeshLib::Pyramid13(elem_nodes);
                break;
            case CellType::INVALID:
                OGS_FATAL(
                    "NodePartitionedMeshReader: construction of INVALID "
                    "element type with id {:d} is not possible.",
                    i);
                break;
            default:
                OGS_FATAL(
                    "NodePartitionedMeshReader: construction of element type "
                    "{:d} is not implemented.",
                    e_type);
        }
    }
}
}  // namespace IO
}  // namespace MeshLib