solution_output.h

Go to the documentation of this file.

// SPDX-FileCopyrightText: © 2025 PRISMS Center at the University of Michigan
// SPDX-License-Identifier: GNU Lesser General Public Version 2.1
 
#pragma once
 
#include <deal.II/base/data_out_base.h>
#include <deal.II/base/exceptions.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/lac/la_parallel_vector.h>
#include <deal.II/numerics/data_component_interpretation.h>
#include <deal.II/numerics/data_out.h>
 
#include <prismspf/core/dof_manager.h>
#include <prismspf/core/field_attributes.h>
#include <prismspf/core/simulation_timer.h>
#include <prismspf/core/solution_indexer.h>
#include <prismspf/core/solution_output.h>
#include <prismspf/core/type_enums.h>
 
#include <prismspf/solvers/solve_context.h>
 
#include <prismspf/user_inputs/output_parameters.h>
#include <prismspf/user_inputs/user_input_parameters.h>
 
#include <prismspf/config.h>
 
#include <filesystem>
#include <fstream>
#include <iomanip>
#include <mpi.h>
#include <sstream>
#include <string>
#include <vector>
 
PRISMS_PF_BEGIN_NAMESPACE
 
template <unsigned int dim>
class UserInputParameters;

template <unsigned int dim, unsigned int degree, typename number>
class SolutionOutput
{
public:
  using VectorType = SolutionVector<number>;

  SolutionOutput(const std::vector<FieldAttributes> &field_attributes,
                 const SolutionIndexer<dim, number> &solution_indexer,
                 const SimulationTimer              &sim_timer,
                 const DoFManager<dim, degree>      &dof_manager,
                 const std::string                  &file_prefix,
                 const UserInputParameters<dim>     &user_inputs)
  {
    const OutputParameters &output_parameters = user_inputs.output_parameters;
    // Some stuff to determine the actual name of the output file.
    const auto n_trailing_digits = static_cast<unsigned int>(
      std::floor(std::log10(user_inputs.temporal_discretization.num_increments)) + 1);
 
    // Init data out
    dealii::DataOut<dim> data_out;
 
    // Add data vectors
    for (unsigned int index = 0; index < field_attributes.size(); ++index)
      {
        const FieldAttributes &field    = field_attributes[index];
        auto                  &solution = solution_indexer.get_solution_vector(index);
        solution.update_ghost_values();
 
        // Mark field as Scalar/Vector
        const unsigned int n_components =
          (field.field_type == TensorRank::Scalar) ? 1 : dim;
 
        const std::vector<
          dealii::DataComponentInterpretation::DataComponentInterpretation>
          data_type(n_components,
                    (n_components == 1)
                      ? dealii::DataComponentInterpretation::component_is_scalar
                      : dealii::DataComponentInterpretation::component_is_part_of_vector);
 
        const std::vector<std::string> names(n_components, field.name);
 
        data_out.add_data_vector(dof_manager.get_field_dof_handler(index),
                                 solution,
                                 names,
                                 data_type);
 
        solution.zero_out_ghost_values();
      }
 
    // Build patches to linearly interpolate from higher order element degrees. Note that
    // this essentially converts the element to an equal amount of subdivisions in the
    // output. This does not make subdivisions and element degree equivalent in the
    // simulation!
    const unsigned int n_divisions = output_parameters.patch_subdivisions == 0
                                       ? degree
                                       : output_parameters.patch_subdivisions;
    data_out.build_patches(n_divisions);
 
    // Set some flags for data output
    dealii::DataOutBase::VtkFlags flags;
    flags.time                = sim_timer.get_time();
    flags.cycle               = sim_timer.get_increment();
    flags.print_date_and_time = true;
#ifdef PRISMS_PF_WITH_ZLIB
    flags.compression_level = dealii::DataOutBase::CompressionLevel::best_speed;
#endif
    data_out.set_flags(flags);
 
    // Write to file based on the user input.
    const unsigned int increment = sim_timer.get_increment();
    const std::string &file_type = output_parameters.file_type;
    if (file_type == "vtu")
      {
        std::ostringstream increment_stream;
        increment_stream << std::setw(static_cast<int>(n_trailing_digits))
                         << std::setfill('0') << increment;
        const std::string filename = file_prefix + "_" + increment_stream.str() + ".vtu";
        data_out.write_vtu_in_parallel(filename, MPI_COMM_WORLD);
      }
    else if (file_type == "pvtu")
      {
        std::filesystem::path output_path = file_prefix;
        std::string           filename    = output_path.filename();
        std::string           directory   = output_path.remove_filename();
        data_out.write_vtu_with_pvtu_record(directory,
                                            filename,
                                            increment,
                                            MPI_COMM_WORLD,
                                            n_trailing_digits);
      }
    else if (file_type == "vtk")
      {
        std::ostringstream increment_stream;
        increment_stream << std::setw(static_cast<int>(n_trailing_digits))
                         << std::setfill('0') << increment;
        const std::string filename = file_prefix + "_" + increment_stream.str() + ".vtk";
        std::ofstream     vtk_output(filename);
        data_out.write_vtk(vtk_output);
      }
    else if (file_type == "xdmf")
      {
#ifdef DEAL_II_WITH_HDF5
        std::ostringstream increment_stream;
        increment_stream << std::setw(static_cast<int>(n_trailing_digits))
                         << std::setfill('0') << increment;
 
        const std::string h5_filename =
          file_prefix + "_" + increment_stream.str() + ".h5";
        const std::string xdmf_filename =
          file_prefix + "_" + increment_stream.str() + ".xdmf";
 
        // Prepare the data filter
        dealii::DataOutBase::DataOutFilter data_filter(
          dealii::DataOutBase::DataOutFilterFlags(true, true));
        data_out.write_filtered_data(data_filter);
 
        // Write binary HDF5
        data_out.write_hdf5_parallel(data_filter, h5_filename, MPI_COMM_WORLD);
 
        // Create the XDMF wrapper for this timestep
        std::vector<dealii::XDMFEntry> xdmf_entries;
        xdmf_entries.push_back(data_out.create_xdmf_entry(data_filter,
                                                          h5_filename,
                                                          increment,
                                                          MPI_COMM_WORLD));
 
        data_out.write_xdmf_file(xdmf_entries, xdmf_filename, MPI_COMM_WORLD);
#else
        AssertThrow(
          false,
          dealii::ExcMessage(
            "You are trying to write an XDMF file as an output; however, deal.II "
            "was not built with HDF5. Please reconfig deal.II with HDF5."));
#endif
      }
    else
      {
        AssertThrow(false, UnreachableCode());
      }
 
    // Update the ghost values again to allow for read access
    for (unsigned int index = 0; index < field_attributes.size(); ++index)
      {
        solution_indexer.get_solution_vector(index).update_ghost_values();
      }
  }
};
 
PRISMS_PF_END_NAMESPACE