WriteSystemMatrixHelper.hpp

/*
  Copyright 2019 SINTEF Digital, Mathematics and Cybernetics.

  This file is part of the Open Porous Media project (OPM).

  OPM is free software: you can redistribute it and/or modify
  it under the terms of the GNU General Public License as published by
  the Free Software Foundation, either version 3 of the License, or
  (at your option) any later version.

  OPM is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  GNU General Public License for more details.

  You should have received a copy of the GNU General Public License
  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
*/

#ifndef OPM_WRITESYSTEMMATRIXHELPER_HEADER_INCLUDED
#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 detail {
        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.
            const auto& iterCtx = simulator.problem().iterationContext();
            std::ostringstream oss;
            oss << "prob_"  << simulator.episodeIndex()
                << "_time_" << std::setprecision(15) << std::setw(12) << std::setfill('0') << simulator.time()
                << "_nit_"  << iterCtx.iteration()
                << '_'      << 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

    template <class SimulatorType, class VectorType, class Communicator>
    void writeVector(const SimulatorType& simulator,
                     const VectorType&    rhs,
                     const std::string&   sysName,
                     const Communicator*  comm)
    {
        detail::writeMatrixMarketObject(simulator, rhs, sysName + "_vector", comm);
    }

    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);
    }

    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);
    }

    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);
    }

    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