Skip to content
Snippets Groups Projects
Select Git revision
  • e464abba248926cb97e826ab9ad32a3760adf77f
  • develop default protected
  • feature/variational-hydro
  • origin/stage/bouguettaia
  • feature/gmsh-reader
  • feature/reconstruction
  • save_clemence
  • feature/kinetic-schemes
  • feature/local-dt-fsi
  • feature/composite-scheme-sources
  • feature/composite-scheme-other-fluxes
  • feature/serraille
  • feature/composite-scheme
  • hyperplastic
  • feature/polynomials
  • feature/gks
  • feature/implicit-solver-o2
  • feature/coupling_module
  • feature/implicit-solver
  • feature/merge-local-dt-fsi
  • master protected
  • v0.5.0 protected
  • v0.4.1 protected
  • v0.4.0 protected
  • v0.3.0 protected
  • v0.2.0 protected
  • v0.1.0 protected
  • Kidder
  • v0.0.4 protected
  • v0.0.3 protected
  • v0.0.2 protected
  • v0 protected
  • v0.0.1 protected
33 results

CheckpointUtils.hpp

Blame
  • GnuplotWriter1D.cpp 11.40 KiB
    #include <output/GnuplotWriter1D.hpp>
    
    #include <mesh/Connectivity.hpp>
    #include <mesh/ItemValue.hpp>
    #include <mesh/Mesh.hpp>
    #include <mesh/MeshData.hpp>
    #include <mesh/MeshDataManager.hpp>
    #include <utils/Messenger.hpp>
    #include <utils/PugsTraits.hpp>
    #include <utils/RevisionInfo.hpp>
    #include <utils/Stringify.hpp>
    
    #include <utils/Demangle.hpp>
    
    #include <fstream>
    #include <iomanip>
    
    std::string
    GnuplotWriter1D::_getDateAndVersionComment() const
    {
      std::ostringstream os;
    
      std::time_t now = std::time(nullptr);
      os << "#  Generated by pugs: " << std::ctime(&now);
      os << "#  version: " << RevisionInfo::version() << '\n';
      os << "#  tag:  " << RevisionInfo::gitTag() << '\n';
      os << "#  HEAD: " << RevisionInfo::gitHead() << '\n';
      os << "#  hash: " << RevisionInfo::gitHash() << " (" << ((RevisionInfo::gitIsClean()) ? "clean" : "dirty") << ")\n";
      os << '\n';
    
      return os.str();
    }
    
    std::string
    GnuplotWriter1D::_getFilename() const
    {
      std::ostringstream sout;
      sout << m_base_filename;
      if (m_period_manager.has_value()) {
        sout << '.' << std::setfill('0') << std::setw(4) << m_period_manager->nbSavedTimes();
      }
      sout << ".gnu";
      return sout.str();
    }
    
    template <typename ItemDataType>
    bool
    GnuplotWriter1D::_is_cell_data(const ItemDataType&) const
    {
      return ItemDataType::item_t == ItemType::cell;
    }
    
    template <typename ItemDataType>
    bool
    GnuplotWriter1D::_is_node_data(const ItemDataType&) const
    {
      return ItemDataType::item_t == ItemType::node;
    }
    
    void
    GnuplotWriter1D::_writePreamble(const OutputNamedItemDataSet& output_named_item_data_set, std::ostream& fout) const
    {
      fout << "# list of data\n";
      fout << "# 1:x";
      uint64_t i = 2;
      for (const auto& i_named_item_data : output_named_item_data_set) {
        const std::string name        = i_named_item_data.first;
        const auto& item_data_variant = i_named_item_data.second;
        std::visit(
          [&](auto&& item_data) {
            using ItemDataType = std::decay_t<decltype(item_data)>;
            using DataType     = std::decay_t<typename ItemDataType::data_type>;
            if constexpr (is_item_value_v<ItemDataType>) {
              if constexpr (std::is_arithmetic_v<DataType>) {
                fout << ' ' << i++ << ':' << name;
              } else if constexpr (is_tiny_vector_v<DataType>) {
                for (size_t j = 0; j < DataType{}.dimension(); ++j) {
                  fout << ' ' << i++ << ':' << name << '[' << j << ']';
                }
              } else if constexpr (is_tiny_matrix_v<DataType>) {
                for (size_t j = 0; j < DataType{}.numberOfRows(); ++j) {
                  for (size_t k = 0; k < DataType{}.numberOfColumns(); ++k) {
                    fout << ' ' << i++ << ':' << name << '(' << j << ',' << k << ')';
                  }
                }
              } else {
                throw UnexpectedError("invalid data type");
              }
            } else if constexpr (is_item_array_v<ItemDataType>) {
              if constexpr (std::is_arithmetic_v<DataType>) {
                for (size_t j = 0; j < item_data.sizeOfArrays(); ++j) {
                  fout << ' ' << i++ << ':' << name << '[' << j << ']';
                }
              } else {
                throw UnexpectedError("invalid data type");
              }
            } else {
              throw UnexpectedError("invalid ItemData type");
            }
          },
          item_data_variant);
      }
      fout << "\n\n";
    }
    
    template <typename DataType, ItemType item_type>
    size_t
    GnuplotWriter1D::_itemDataNbRow(const ItemValue<DataType, item_type>&) const
    {
      if constexpr (std::is_arithmetic_v<DataType>) {
        return 1;
      } else if constexpr (is_tiny_vector_v<std::decay_t<DataType>>) {
        return DataType::Dimension;
      } else if constexpr (is_tiny_matrix_v<std::decay_t<DataType>>) {
        return DataType{}.dimension();
      } else {
        throw UnexpectedError("invalid data type for cell value output: " + demangle<DataType>());
      }
    }
    
    template <typename DataType, ItemType item_type>
    size_t
    GnuplotWriter1D::_itemDataNbRow(const ItemArray<DataType, item_type>& item_array) const
    {
      return item_array.sizeOfArrays();
    }
    
    template <typename MeshType, ItemType item_type>
    void
    GnuplotWriter1D::_writeItemDatas(const std::shared_ptr<const MeshType>& mesh,
                                     const OutputNamedItemDataSet& output_named_item_data_set,
                                     std::ostream& fout) const
    {
      using ItemId = ItemIdT<item_type>;
    
      const size_t& number_of_columns = [&] {
        size_t number_of_columns = 1;
        for (auto [name, item_data] : output_named_item_data_set) {
          std::visit([&](auto&& value) { number_of_columns += _itemDataNbRow(value); }, item_data);
        }
        return number_of_columns;
      }();
    
      auto is_owned = mesh->connectivity().template isOwned<item_type>();
    
      const size_t& number_of_owned_lines = [&]() {
        if (parallel::size() > 1) {
          size_t number_of_owned_items = 0;
          for (ItemId item_id = 0; item_id < mesh->template numberOf<item_type>(); ++item_id) {
            if (is_owned[item_id]) {
              ++number_of_owned_items;
            }
          }
    
          return number_of_owned_items;
        } else {
          return mesh->template numberOf<item_type>();
        }
      }();
    
      Array<double> values{number_of_columns * number_of_owned_lines};
    
      if constexpr (item_type == ItemType::cell) {
        auto& mesh_data         = MeshDataManager::instance().getMeshData(*mesh);
        const auto& cell_center = mesh_data.xj();
    
        size_t index = 0;
        for (ItemId item_id = 0; item_id < mesh->template numberOf<item_type>(); ++item_id) {
          if (is_owned[item_id]) {
            values[number_of_columns * index++] = cell_center[item_id][0];
          }
        }
      } else if constexpr (item_type == ItemType::node) {
        const auto& node_position = mesh->xr();
    
        size_t index = 0;
        for (ItemId item_id = 0; item_id < mesh->template numberOf<item_type>(); ++item_id) {
          if (is_owned[item_id]) {
            values[number_of_columns * index++] = node_position[item_id][0];
          }
        }
      } else {
        throw UnexpectedError("invalid item type");
      }
    
      size_t column_number = 1;
      for (auto [name, output_item_data] : output_named_item_data_set) {
        std::visit(
          [&](auto&& item_data) {
            using ItemDataT = std::decay_t<decltype(item_data)>;
            if constexpr (ItemDataT::item_t == item_type) {
              if constexpr (is_item_value_v<ItemDataT>) {
                using DataT  = std::decay_t<typename ItemDataT::data_type>;
                size_t index = 0;
                for (ItemId item_id = 0; item_id < item_data.numberOfItems(); ++item_id) {
                  if (is_owned[item_id]) {
                    if constexpr (std::is_arithmetic_v<DataT>) {
                      values[number_of_columns * index + column_number] = item_data[item_id];
                    } else if constexpr (is_tiny_vector_v<DataT>) {
                      const size_t k = number_of_columns * index + column_number;
                      for (size_t j = 0; j < DataT::Dimension; ++j) {
                        values[k + j] = item_data[item_id][j];
                      }
                    } else if constexpr (is_tiny_matrix_v<DataT>) {
                      size_t k = number_of_columns * index + column_number;
                      for (size_t i = 0; i < DataT{}.numberOfRows(); ++i) {
                        for (size_t j = 0; j < DataT{}.numberOfColumns(); ++j) {
                          values[k++] = item_data[item_id](i, j);
                        }
                      }
                    }
                    ++index;
                  }
                }
              } else {
                using DataT  = std::decay_t<typename ItemDataT::data_type>;
                size_t index = 0;
                for (ItemId item_id = 0; item_id < item_data.numberOfItems(); ++item_id) {
                  if (is_owned[item_id]) {
                    if constexpr (std::is_arithmetic_v<DataT>) {
                      const size_t k = number_of_columns * index + column_number;
                      for (size_t j = 0; j < item_data.sizeOfArrays(); ++j) {
                        values[k + j] = item_data[item_id][j];
                      }
                    }
                    ++index;
                  }
                }
              }
            }
            column_number += _itemDataNbRow(item_data);
          },
          output_item_data);
      }
    
      if (parallel::size() > 1) {
        values = parallel::gatherVariable(values, 0);
      }
    
      if (parallel::rank() == 0) {
        Assert(values.size() % number_of_columns == 0);
    
        std::vector<size_t> line_numbers(values.size() / number_of_columns);
        for (size_t i = 0; i < line_numbers.size(); ++i) {
          line_numbers[i] = i;
        }
    
        std::sort(line_numbers.begin(), line_numbers.end(),
                  [&](size_t i, size_t j) { return values[i * number_of_columns] < values[j * number_of_columns]; });
    
        for (auto i_line : line_numbers) {
          fout << values[i_line * number_of_columns];
          for (size_t j = 1; j < number_of_columns; ++j) {
            fout << ' ' << values[i_line * number_of_columns + j];
          }
          fout << '\n';
        }
      }
    }
    
    template <typename MeshType>
    void
    GnuplotWriter1D::_writeAtTime(const std::shared_ptr<const MeshType>& mesh,
                                  const OutputNamedItemDataSet& output_named_item_data_set,
                                  double time) const
    {
      bool has_cell_data = false;
      for (const auto& [name, item_data_variant] : output_named_item_data_set) {
        has_cell_data |= std::visit([&](auto&& item_data) { return this->_is_cell_data(item_data); }, item_data_variant);
      }
    
      bool has_node_data = false;
      for (const auto& [name, item_data_variant] : output_named_item_data_set) {
        has_node_data |= std::visit([&](auto&& item_data) { return this->_is_node_data(item_data); }, item_data_variant);
      }
    
      if (has_cell_data and has_node_data) {
        throw NormalError("cannot store both node and cell data in the same gnuplot file");
      }
    
      std::ofstream fout;
    
      if (parallel::rank() == 0) {
        fout.open(_getFilename());
        fout.precision(15);
        fout.setf(std::ios_base::scientific);
        fout << _getDateAndVersionComment();
    
        fout << "# time = " << time << "\n\n";
    
        _writePreamble(output_named_item_data_set, fout);
      }
    
      if (has_cell_data) {
        this->_writeItemDatas<MeshType, ItemType::cell>(mesh, output_named_item_data_set, fout);
      } else {   // has_node_value
        this->_writeItemDatas<MeshType, ItemType::node>(mesh, output_named_item_data_set, fout);
      }
    }
    
    void
    GnuplotWriter1D::writeMesh(const std::shared_ptr<const IMesh>&) const
    {
      std::ostringstream errorMsg;
      errorMsg << "gnuplot_1d_writer does not write meshes\n"
               << rang::style::bold << "note:" << rang::style::reset << " one can use " << rang::fgB::blue
               << "gnuplot_writer" << rang::fg::reset << "  instead";
      throw NormalError(errorMsg.str());
    }
    
    void
    GnuplotWriter1D::_writeAtTime(
      const std::vector<std::shared_ptr<const NamedDiscreteFunction>>& named_discrete_function_list,
      double time) const
    {
      std::shared_ptr mesh = this->_getMesh(named_discrete_function_list);
    
      OutputNamedItemDataSet output_named_item_data_set = this->_getOutputNamedItemDataSet(named_discrete_function_list);
    
      switch (mesh->dimension()) {
      case 1: {
        this->_writeAtTime(std::dynamic_pointer_cast<const Mesh<Connectivity<1>>>(mesh), output_named_item_data_set, time);
        break;
      }
      case 2: {
        std::ostringstream errorMsg;
        errorMsg << "gnuplot_1d_writer is not available in dimension " << stringify(mesh->dimension()) << '\n'
                 << rang::style::bold << "note:" << rang::style::reset << " one can use " << rang::fgB::blue
                 << "gnuplot_writer" << rang::fg::reset << " in dimension 2";
        throw NormalError(errorMsg.str());
      }
      default: {
        throw NormalError("gnuplot format is not available in dimension " + stringify(mesh->dimension()));
      }
      }
    }