Add Helper Functions for Outputting Vectors/Matrices Separately

While here, also add wrapper functions for outputting systems of
linear equations arising in geo-mechanical applications.
This commit is contained in:
Halvor M Nilsen
2025-01-31 18:09:10 +01:00
committed by Bård Skaflestad
parent 381182dedd
commit 08e39f7479

View File

@@ -21,68 +21,268 @@
#define OPM_WRITESYSTEMMATRIXHELPER_HEADER_INCLUDED
#include <dune/istl/matrixmarket.hh>
#include <opm/simulators/linalg/MatrixMarketSpecializations.hpp>
#include <filesystem>
#include <iomanip>
#include <sstream>
#include <string>
namespace Opm::Helper {
namespace Opm
{
namespace Helper
{
template <class SimulatorType, class MatrixType, class VectorType, class Communicator>
void writeSystem(const SimulatorType& simulator,
const MatrixType& matrix,
const VectorType& rhs,
[[maybe_unused]] const Communicator* comm)
namespace detail {
/// Output a linear algebra object to file in MatrixMarket format.
///
/// Common implementation function.
///
/// \tparam SimulatorType Packaged simulator type. Expected to
/// provide, among other services, an episode index, a simulation
/// time, and a physical model.
///
/// \tparam LinalgObjectType Object type for a linear algebra
/// object. Expected to be a coefficient matrix type or a linear
/// system vector type from the Dune ISTL.
///
/// \tparam Communicator ISTL communication type. Typically
/// \code Dune::OwnerOverlapCopyCommunication<> \endcode or similar.
///
/// \param[in] simulator Simulator object.
///
/// \param[in] linalgObject Linear algebra object. Expected to be a
/// Dune ISTL matrix or vector.
///
/// \param[in] objName Name of linear algebra object. Incorporated
/// into the name of the MatrixMarket output file.
///
/// \param[in] comm Dune ISTL communication object. Unused in
/// builds without MPI support. You may use a nullptr value in a
/// sequential run of an MPI-enabled build to signify that parallel
/// output is unneeded.
template <class SimulatorType, class LinalgObjectType, class Communicator>
void writeMatrixMarketObject(const SimulatorType& simulator,
const LinalgObjectType& linalgObject,
const std::string& objName,
[[maybe_unused]] const Communicator* comm)
{
namespace fs = std::filesystem;
const auto output_dir = fs::path {simulator.problem().outputDir()} / "reports";
if (! fs::exists(output_dir)) {
fs::create_directory(output_dir);
}
// Combine and return.
std::ostringstream oss;
oss << "prob_" << simulator.episodeIndex()
<< "_time_" << std::setprecision(15) << std::setw(12) << std::setfill('0') << simulator.time()
<< "_nit_" << simulator.model().newtonMethod().numIterations()
<< '_' << objName << "_istl";
const auto filename = (output_dir / oss.str()).generic_string();
#if HAVE_MPI
if (comm != nullptr) {
Dune::storeMatrixMarket(linalgObject, filename, *comm);
}
else
#endif // HAVE_MPI
{
Dune::storeMatrixMarket(linalgObject, filename + ".mm");
}
}
} // namespace detail
/// Output a Dune ISTL vector to file in MatrixMarket format.
///
/// \tparam SimulatorType Packaged simulator type. Expected to
/// provide, among other services, an episode index, a simulation
/// time, and a physical model.
///
/// \tparam VectorType Object type for a linear system vector object.
/// Typically a linear system vector type from the Dune ISTL.
///
/// \tparam Communicator ISTL communication type. Typically
/// \code Dune::OwnerOverlapCopyCommunication<> \endcode or similar.
///
/// \param[in] simulator Simulator object.
///
/// \param[in] rhs Linear system right-hand side vector. Expected to be
/// a Dune ISTL vector.
///
/// \param[in] sysName Name of linear system/physical model.
/// Incorporated into the name of the MatrixMarket output file.
///
/// \param[in] comm Dune ISTL communication object. Unused in builds
/// without MPI support. You may use a nullptr value in a sequential
/// run of an MPI-enabled build to signify that parallel output is
/// unneeded.
template <class SimulatorType, class VectorType, class Communicator>
void writeVector(const SimulatorType& simulator,
const VectorType& rhs,
const std::string& sysName,
const Communicator* comm)
{
std::string dir = simulator.problem().outputDir();
if (dir == ".") {
dir = "";
} else if (!dir.empty() && dir.back() != '/') {
dir += "/";
}
namespace fs = ::std::filesystem;
fs::path output_dir(dir);
fs::path subdir("reports");
output_dir = output_dir / subdir;
if (!(fs::exists(output_dir))) {
fs::create_directory(output_dir);
}
// Combine and return.
std::ostringstream oss;
oss << "prob_" << simulator.episodeIndex() << "_time_";
oss << std::setprecision(15) << std::setw(12) << std::setfill('0') << simulator.time() << "_";
int nit = simulator.model().newtonMethod().numIterations();
oss << "_nit_" << nit << "_";
std::string output_file(oss.str());
fs::path full_path = output_dir / output_file;
std::string prefix = full_path.string();
{
std::string filename = prefix + "matrix_istl";
#if HAVE_MPI
if (comm != nullptr) { // comm is not set in serial runs
Dune::storeMatrixMarket(matrix, filename, *comm, true);
} else
#endif
{
Dune::storeMatrixMarket(matrix, filename + ".mm");
}
}
{
std::string filename = prefix + "rhs_istl";
#if HAVE_MPI
if (comm != nullptr) { // comm is not set in serial runs
Dune::storeMatrixMarket(rhs, filename, *comm, true);
} else
#endif
{
Dune::storeMatrixMarket(rhs, filename + ".mm");
}
}
detail::writeMatrixMarketObject(simulator, rhs, sysName + "_vector", comm);
}
/// Output a Dune ISTL matrix to file in MatrixMarket format.
///
/// \tparam SimulatorType Packaged simulator type. Expected to
/// provide, among other services, an episode index, a simulation
/// time, and a physical model.
///
/// \tparam MatrixType Object type for a linear system matrix object.
/// Typically a linear system matrix type from the Dune ISTL.
///
/// \tparam Communicator ISTL communication type. Typically
/// \code Dune::OwnerOverlapCopyCommunication<> \endcode or similar.
///
/// \param[in] simulator Simulator object.
///
/// \param[in] matrix Linear system coefficient matrix. Expected to be
/// a Dune ISTL (sparse) matrix.
///
/// \param[in] sysName Name of linear system/physical model.
/// Incorporated into the name of the MatrixMarket output file.
///
/// \param[in] comm Dune ISTL communication object. Unused in builds
/// without MPI support. You may use a nullptr value in a sequential
/// run of an MPI-enabled build to signify that parallel output is
/// unneeded.
template <class SimulatorType, class MatrixType, class Communicator>
void writeMatrix(const SimulatorType& simulator,
const MatrixType& matrix,
const std::string& sysName,
const Communicator* comm)
{
detail::writeMatrixMarketObject(simulator, matrix, sysName + "_matrix", comm);
}
} // namespace Helper
} // namespace Opm
#endif
/// Output a Dune ISTL linear system to files in MatrixMarket format.
///
/// This function will create one output file for the coefficient matrix
/// and another output file for the system right-hand side vector.
///
/// \tparam SimulatorType Packaged simulator type. Expected to
/// provide, among other services, an episode index, a simulation
/// time, and a physical model.
///
/// \tparam MatrixType Object type for a linear system matrix object.
/// Typically a linear system matrix type from the Dune ISTL.
///
/// \tparam VectorType Object type for a linear system vector object.
/// Typically a linear system vector type from the Dune ISTL.
///
/// \tparam Communicator ISTL communication type. Typically
/// \code Dune::OwnerOverlapCopyCommunication<> \endcode or similar.
///
/// \param[in] simulator Simulator object.
///
/// \param[in] matrix Linear system coefficient matrix. Expected to be
/// a Dune ISTL (sparse) matrix.
///
/// \param[in] rhs Linear system right-hand side vector. Expected to be
/// a Dune ISTL vector.
///
/// \param[in] sysName Name of linear system/physical model.
/// Incorporated into the name of the MatrixMarket output file.
///
/// \param[in] comm Dune ISTL communication object. Unused in builds
/// without MPI support. You may use a nullptr value in a sequential
/// run of an MPI-enabled build to signify that parallel output is
/// unneeded.
template <class SimulatorType, class MatrixType, class VectorType, class Communicator>
void writeSystem(const SimulatorType& simulator,
const MatrixType& matrix,
const VectorType& rhs,
const std::string& sysName,
const Communicator* comm)
{
writeMatrix(simulator, matrix, sysName, comm);
writeVector(simulator, rhs, sysName, comm);
}
/// Output a Dune ISTL linear system from linearising a set of flow
/// equations to files in MatrixMarket format.
///
/// This function will create one output file for the coefficient matrix
/// and another output file for the system right-hand side vector.
///
/// \tparam SimulatorType Packaged simulator type. Expected to
/// provide, among other services, an episode index, a simulation
/// time, and a physical model.
///
/// \tparam MatrixType Object type for a linear system matrix object.
/// Typically a linear system matrix type from the Dune ISTL.
///
/// \tparam VectorType Object type for a linear system vector object.
/// Typically a linear system vector type from the Dune ISTL.
///
/// \tparam Communicator ISTL communication type. Typically
/// \code Dune::OwnerOverlapCopyCommunication<> \endcode or similar.
///
/// \param[in] simulator Simulator object.
///
/// \param[in] matrix Linear system coefficient matrix. Expected to be
/// a Dune ISTL (sparse) matrix.
///
/// \param[in] rhs Linear system right-hand side vector. Expected to be
/// a Dune ISTL vector.
///
/// \param[in] comm Dune ISTL communication object. Unused in builds
/// without MPI support. You may use a nullptr value in a sequential
/// run of an MPI-enabled build to signify that parallel output is
/// unneeded.
template <class SimulatorType, class MatrixType, class VectorType, class Communicator>
void writeSystem(const SimulatorType& simulator,
const MatrixType& matrix,
const VectorType& rhs,
const Communicator* comm)
{
writeSystem(simulator, matrix, rhs, "flow_", comm);
}
/// Output a Dune ISTL linear system from linearising a set of
/// geo-mechanical equations to files in MatrixMarket format.
///
/// This function will create one output file for the coefficient matrix
/// and another output file for the system right-hand side vector.
///
/// \tparam SimulatorType Packaged simulator type. Expected to
/// provide, among other services, an episode index, a simulation
/// time, and a physical model.
///
/// \tparam MatrixType Object type for a linear system matrix object.
/// Typically a linear system matrix type from the Dune ISTL.
///
/// \tparam VectorType Object type for a linear system vector object.
/// Typically a linear system vector type from the Dune ISTL.
///
/// \tparam Communicator ISTL communication type. Typically
/// \code Dune::OwnerOverlapCopyCommunication<> \endcode or similar.
///
/// \param[in] simulator Simulator object.
///
/// \param[in] matrix Linear system coefficient matrix. Expected to be
/// a Dune ISTL (sparse) matrix.
///
/// \param[in] rhs Linear system right-hand side vector. Expected to be
/// a Dune ISTL vector.
///
/// \param[in] comm Dune ISTL communication object. Unused in builds
/// without MPI support. You may use a nullptr value in a sequential
/// run of an MPI-enabled build to signify that parallel output is
/// unneeded.
template <class SimulatorType, class MatrixType, class VectorType, class Communicator>
void writeMechSystem(const SimulatorType& simulator,
const MatrixType& matrix,
const VectorType& rhs,
const Communicator* comm)
{
writeSystem(simulator, matrix, rhs, "mech_", comm);
}
} // namespace Opm::Helper
#endif // OPM_WRITESYSTEMMATRIXHELPER_HEADER_INCLUDED