Output.cpp

Go to the documentation of this file.

 
#include "Output.h"
 
#include <cassert>
#include <exception>
#include <fstream>
#include <vector>
 
#include "Applications/InSituLib/Adaptor.h"
#include "BaseLib/FileTools.h"
#include "BaseLib/Logging.h"
#include "BaseLib/RunTime.h"
#include "ProcessLib/Process.h"
 
namespace
{
int convertVtkDataMode(std::string const& data_mode)
{
    if (data_mode == "Ascii")
    {
        return 0;
    }
    if (data_mode == "Binary")
    {
        return 1;
    }
    if (data_mode == "Appended")
    {
        return 2;
    }
    OGS_FATAL(
        "Unsupported vtk output file data mode '{:s}'. Expected Ascii, Binary, "
        "or Appended.",
        data_mode);
}
 
std::string constructPVDName(std::string const& output_directory,
                             std::string const& output_file_prefix,
                             std::string const& mesh_name)
{
    return BaseLib::joinPaths(output_directory,
                              BaseLib::constructFormattedFileName(
                                  output_file_prefix, mesh_name, 0, 0, 0) +
                                  ".pvd");
}
}  // namespace
 
namespace ProcessLib
{
bool Output::shallDoOutput(int timestep, double const t)
{
    auto const fixed_output_time = std::lower_bound(
        cbegin(_fixed_output_times), cend(_fixed_output_times), t);
    if ((fixed_output_time != cend(_fixed_output_times)) &&
        (std::abs(*fixed_output_time - t) <
         std::numeric_limits<double>::epsilon()))
    {
        return true;
    }
 
    int each_steps = 1;
 
    for (auto const& pair : _repeats_each_steps)
    {
        each_steps = pair.each_steps;
 
        if (timestep > pair.repeat * each_steps)
        {
            timestep -= pair.repeat * each_steps;
        }
        else
        {
            break;
        }
    }
 
    if (timestep % each_steps == 0)
    {
        return true;
    }
 
    return false;
}
 
Output::Output(std::string directory, OutputType file_type,
               std::string file_prefix, std::string file_suffix,
               bool const compress_output, unsigned int const n_files,
               std::string const& data_mode,
               bool const output_nonlinear_iteration_results,
               std::vector<PairRepeatEachSteps> repeats_each_steps,
               std::vector<double>&& fixed_output_times,
               OutputDataSpecification&& output_data_specification,
               std::vector<std::string>&& mesh_names_for_output,
               std::vector<std::unique_ptr<MeshLib::Mesh>> const& meshes)
    : _output_directory(std::move(directory)),
      _output_file_type(file_type),
      _output_file_prefix(std::move(file_prefix)),
      _output_file_suffix(std::move(file_suffix)),
      _output_file_compression(compress_output),
      _n_files(n_files),
      _output_file_data_mode(convertVtkDataMode(data_mode)),
      _output_nonlinear_iteration_results(output_nonlinear_iteration_results),
      _repeats_each_steps(std::move(repeats_each_steps)),
      _fixed_output_times(std::move(fixed_output_times)),
      _output_data_specification(std::move(output_data_specification)),
      _mesh_names_for_output(mesh_names_for_output),
      _meshes(meshes)
{
    if (!std::is_sorted(cbegin(_fixed_output_times), cend(_fixed_output_times)))
    {
        OGS_FATAL(
            "Vector of fixed output time steps passed to the Output "
            "constructor must be sorted");
    }
}
 
void Output::addProcess(ProcessLib::Process const& process)
{
    if (_mesh_names_for_output.empty())
    {
        _mesh_names_for_output.push_back(process.getMesh().getName());
    }
 
    for (auto const& mesh_output_name : _mesh_names_for_output)
    {
        auto const filename = constructPVDName(
            _output_directory, _output_file_prefix, mesh_output_name);
        _process_to_pvd_file.emplace(std::piecewise_construct,
                                     std::forward_as_tuple(&process),
                                     std::forward_as_tuple(filename));
    }
}
 
// TODO return a reference.
MeshLib::IO::PVDFile* Output::findPVDFile(
    Process const& process,
    const int process_id,
    std::string const& mesh_name_for_output)
{
    auto const filename = constructPVDName(
        _output_directory, _output_file_prefix, mesh_name_for_output);
    auto range = _process_to_pvd_file.equal_range(&process);
    int counter = 0;
    MeshLib::IO::PVDFile* pvd_file = nullptr;
    for (auto spd_it = range.first; spd_it != range.second; ++spd_it)
    {
        if (spd_it->second.pvd_filename == filename)
        {
            if (counter == process_id)
            {
                pvd_file = &spd_it->second;
                break;
            }
            counter++;
        }
    }
    if (pvd_file == nullptr)
    {
        OGS_FATAL(
            "The given process is not contained in the output configuration. "
            "Aborting.");
    }
 
    return pvd_file;
}
 
struct Output::OutputFile
{
    OutputFile(std::string const& directory, OutputType const type,
               std::string const& prefix, std::string const& suffix,
               std::string const& mesh_name, int const timestep, double const t,
               int const iteration, int const data_mode_,
               bool const compression_,
               std::set<std::string> const& outputnames,
               unsigned int const n_files)
        : name(constructFilename(type, prefix, suffix, mesh_name, timestep, t,
                                 iteration)),
          path(BaseLib::joinPaths(directory, name)),
          type(type),
          data_mode(data_mode_),
          compression(compression_),
          outputnames(outputnames),
          n_files(n_files)
    {
    }
 
    std::string const name;
    std::string const path;
    OutputType const type;
    int const data_mode;
    bool const compression;
    std::set<std::string> outputnames;
    unsigned int n_files;
    static std::string constructFilename(OutputType const type,
                                         std::string prefix, std::string suffix,
                                         std::string mesh_name,
                                         int const timestep, double const t,
                                         int const iteration)
    {
        std::map<OutputType, std::string> filetype_to_extension = {
            {OutputType::vtk, "vtu"}, {OutputType::xdmf, "xdmf"}};
 
        try
        {
            std::string extension = filetype_to_extension.at(type);
            return BaseLib::constructFormattedFileName(prefix, mesh_name,
                                                       timestep, t, iteration) +
                   BaseLib::constructFormattedFileName(suffix, mesh_name,
                                                       timestep, t, iteration) +
                   "." + extension;
        }
        catch (std::out_of_range&)
        {
            OGS_FATAL(
                "No supported file type provided. Read `{:s}' from <output><type> \
                    in prj file. Supported: VTK, XDMF.",
                type);
        }
    }
};
 
void Output::outputMeshXdmf(
    OutputFile const& output_file,
    std::vector<std::reference_wrapper<const MeshLib::Mesh>>
        meshes,
    int const timestep,
    double const t)
{
    // \TODO (tm) Refactor to a dedicated VTKOutput and XdmfOutput
    // The XdmfOutput will create on construction the XdmfHdfWriter
    if (!_mesh_xdmf_hdf_writer)
    {
        std::filesystem::path path(output_file.path);
        _mesh_xdmf_hdf_writer = std::make_unique<MeshLib::IO::XdmfHdfWriter>(
            std::move(meshes), path, timestep, t,
            _output_data_specification.output_variables,
            output_file.compression, output_file.n_files);
    }
    else
    {
        _mesh_xdmf_hdf_writer->writeStep(t);
    };
}
 
void Output::outputMesh(OutputFile const& output_file,
                        MeshLib::IO::PVDFile* const pvd_file,
                        MeshLib::Mesh const& mesh,
                        double const t)
{
    DBUG("output to {:s}", output_file.path);
 
    if (output_file.type == OutputType::vtk)
    {
        pvd_file->addVTUFile(output_file.name, t);
    }
 
    makeOutput(output_file.path, mesh, output_file.compression,
               output_file.data_mode);
}
 
void Output::doOutputAlways(Process const& process,
                            const int process_id,
                            int const timestep,
                            const double t,
                            int const iteration,
                            std::vector<GlobalVector*> const& x)
{
    BaseLib::RunTime time_output;
    time_output.start();
 
    std::vector<NumLib::LocalToGlobalIndexMap const*> dof_tables;
    dof_tables.reserve(x.size());
    for (std::size_t i = 0; i < x.size(); ++i)
    {
        dof_tables.push_back(&process.getDOFTable(i));
    }
 
    bool output_secondary_variable = true;
    // Need to add variables of process to vtu even no output takes place.
    addProcessDataToMesh(t, x, process_id, process.getMesh(), dof_tables,
                         dof_tables, process.getProcessVariables(process_id),
                         process.getSecondaryVariables(),
                         output_secondary_variable,
                         process.getIntegrationPointWriter(process.getMesh()),
                         _output_data_specification);
 
    // For the staggered scheme for the coupling, only the last process, which
    // gives the latest solution within a coupling loop, is allowed to make
    // output.
    if (!(process_id == static_cast<int>(_process_to_pvd_file.size() /
                                         _mesh_names_for_output.size()) -
                            1 ||
          process.isMonolithicSchemeUsed()))
    {
        return;
    }
 
    auto output_bulk_mesh =
        [&](std::vector<std::reference_wrapper<const MeshLib::Mesh>> meshes)
    {
        MeshLib::IO::PVDFile* pvd_file = nullptr;
        if (_output_file_type == ProcessLib::OutputType::vtk)
        {
            for (auto const& mesh : meshes)
            {
                OutputFile const file(
                    _output_directory, _output_file_type, _output_file_prefix,
                    _output_file_suffix, mesh.get().getName(), timestep, t,
                    iteration, _output_file_data_mode, _output_file_compression,
                    _output_data_specification.output_variables, 1);
 
                pvd_file =
                    findPVDFile(process, process_id, mesh.get().getName());
                outputMesh(file, pvd_file, mesh, t);
            }
        }
        else if (_output_file_type == ProcessLib::OutputType::xdmf)
        {
            std::string name = meshes[0].get().getName();
            OutputFile const file(
                _output_directory, _output_file_type, _output_file_prefix, "",
                name, timestep, t, iteration, _output_file_data_mode,
                _output_file_compression,
                _output_data_specification.output_variables, _n_files);
 
            outputMeshXdmf(std::move(file), std::move(meshes), timestep, t);
        }
    };
 
    std::vector<std::reference_wrapper<const MeshLib::Mesh>> output_meshes;
    for (auto const& mesh_output_name : _mesh_names_for_output)
    {
        // process related output
        if (process.getMesh().getName() == mesh_output_name)
        {
            output_meshes.push_back(process.getMesh());
            continue;
        }
        // mesh related output
        auto& non_bulk_mesh = *BaseLib::findElementOrError(
            begin(_meshes), end(_meshes),
            [&mesh_output_name](auto const& m)
            { return m->getName() == mesh_output_name; },
            "Need mesh '" + mesh_output_name + "' for the output.");
 
        std::vector<MeshLib::Node*> const& nodes = non_bulk_mesh.getNodes();
        DBUG("Found {:d} nodes for output at mesh '{:s}'.", nodes.size(),
             non_bulk_mesh.getName());
 
        std::vector<std::unique_ptr<NumLib::LocalToGlobalIndexMap>>
            mesh_dof_tables;
        mesh_dof_tables.reserve(x.size());
        for (std::size_t i = 0; i < x.size(); ++i)
        {
            mesh_dof_tables.push_back(
                process.getDOFTable(i).deriveBoundaryConstrainedMap(
                    MeshLib::MeshSubset{non_bulk_mesh, nodes}));
        }
        std::vector<NumLib::LocalToGlobalIndexMap const*>
            mesh_dof_table_pointers;
        mesh_dof_table_pointers.reserve(mesh_dof_tables.size());
        transform(cbegin(mesh_dof_tables), cend(mesh_dof_tables),
                  back_inserter(mesh_dof_table_pointers),
                  [](std::unique_ptr<NumLib::LocalToGlobalIndexMap> const& p)
                  { return p.get(); });
 
        output_secondary_variable = false;
        addProcessDataToMesh(
            t, x, process_id, non_bulk_mesh, dof_tables,
            mesh_dof_table_pointers, process.getProcessVariables(process_id),
            process.getSecondaryVariables(), output_secondary_variable,
            process.getIntegrationPointWriter(non_bulk_mesh),
            _output_data_specification);
 
        output_meshes.push_back(non_bulk_mesh);
    }
 
    output_bulk_mesh(std::move(output_meshes));
    INFO("[time] Output of timestep {:d} took {:g} s.", timestep,
         time_output.elapsed());
}
 
void Output::doOutput(Process const& process,
                      const int process_id,
                      int const timestep,
                      const double t,
                      int const iteration,
                      std::vector<GlobalVector*> const& x)
{
    if (shallDoOutput(timestep, t))
    {
        doOutputAlways(process, process_id, timestep, t, iteration, x);
    }
#ifdef USE_INSITU
    // Note: last time step may be output twice: here and in
    // doOutputLastTimestep() which throws a warning.
    InSituLib::CoProcess(process.getMesh(), t, timestep, false,
                         _output_directory);
#endif
}
 
void Output::doOutputLastTimestep(Process const& process,
                                  const int process_id,
                                  int const timestep,
                                  const double t,
                                  int const iteration,
                                  std::vector<GlobalVector*> const& x)
{
    if (!shallDoOutput(timestep, t))
    {
        doOutputAlways(process, process_id, timestep, t, iteration, x);
    }
#ifdef USE_INSITU
    InSituLib::CoProcess(process.getMesh(), t, timestep, true,
                         _output_directory);
#endif
}
 
void Output::doOutputNonlinearIteration(Process const& process,
                                        const int process_id,
                                        int const timestep, const double t,
                                        int const iteration,
                                        std::vector<GlobalVector*> const& x)
{
    if (!_output_nonlinear_iteration_results)
    {
        return;
    }
 
    BaseLib::RunTime time_output;
    time_output.start();
 
    std::vector<NumLib::LocalToGlobalIndexMap const*> dof_tables;
    for (std::size_t i = 0; i < x.size(); ++i)
    {
        dof_tables.push_back(&process.getDOFTable(i));
    }
 
    bool const output_secondary_variable = true;
    addProcessDataToMesh(t, x, process_id, process.getMesh(), dof_tables,
                         dof_tables, process.getProcessVariables(process_id),
                         process.getSecondaryVariables(),
                         output_secondary_variable,
                         process.getIntegrationPointWriter(process.getMesh()),
                         _output_data_specification);
 
    // For the staggered scheme for the coupling, only the last process, which
    // gives the latest solution within a coupling loop, is allowed to make
    // output.
    if (!(process_id == static_cast<int>(_process_to_pvd_file.size()) - 1 ||
          process.isMonolithicSchemeUsed()))
    {
        return;
    }
 
    // Only check whether a process data is available for output.
    findPVDFile(process, process_id, process.getMesh().getName());
 
    std::string const output_file_name = OutputFile::constructFilename(
        _output_file_type, _output_file_prefix, _output_file_suffix,
        process.getMesh().getName(), timestep, t, iteration);
 
    std::string const output_file_path =
        BaseLib::joinPaths(_output_directory, output_file_name);
 
    DBUG("output iteration results to {:s}", output_file_path);
    makeOutput(output_file_path, process.getMesh(), _output_file_compression,
               _output_file_data_mode);
    INFO("[time] Output took {:g} s.", time_output.elapsed());
}
}  // namespace ProcessLib