text:
/* Copyright (c) 2016 - 2021, the adamantine authors.
 *
 * This file is subject to the Modified BSD License and may not be distributed
 * without copyright and license information. Please refer to the file LICENSE
 * for the text and further information on this license.
 */

#ifndef THERMAL_PHYSICS_TEMPLATES_HH
#define THERMAL_PHYSICS_TEMPLATES_HH

#include <ThermalOperator.hh>
#if defined(ADAMANTINE_HAVE_CUDA) && defined(__CUDACC__)
#include <ThermalOperatorDevice.hh>
#endif
#include <CubeHeatSource.hh>
#include <ElectronBeamHeatSource.hh>
#include <GoldakHeatSource.hh>
#include <ThermalPhysics.hh>

#include <deal.II/base/geometry_info.h>
#include <deal.II/distributed/cell_data_transfer.templates.h>
#include <deal.II/dofs/dof_tools.h>
#include <deal.II/fe/fe_nothing.h>
#include <deal.II/fe/fe_q.h>
#include <deal.II/grid/filtered_iterator.h>
#include <deal.II/hp/fe_values.h>
#include <deal.II/hp/q_collection.h>
#include <deal.II/lac/precondition.h>
#include <deal.II/lac/solver_gmres.h>

#ifdef ADAMANTINE_WITH_CALIPER
#include <caliper/cali.h>
#endif

#include <algorithm>
#include <execution>

namespace adamantine
{
namespace
{
template <int dim, int fe_degree, typename MemorySpaceType,
          std::enable_if_t<
              std::is_same<MemorySpaceType, dealii::MemorySpace::Host>::value,
              int> = 0>
dealii::LA::distributed::Vector<double, MemorySpaceType>
evaluate_thermal_physics_impl(
    std::shared_ptr<ThermalOperatorBase<dim, MemorySpaceType>> thermal_operator,
    double const t, double const current_source_height,
    dealii::LA::distributed::Vector<double, MemorySpaceType> const &y,
    std::vector<Timer> &timers)
{
  timers[evol_time_eval_mat_prop].start();
  thermal_operator->evaluate_material_properties(y);
  timers[evol_time_eval_mat_prop].stop();

  timers[evol_time_eval_th_ph].start();
  thermal_operator->set_time_and_source_height(t, current_source_height);

  dealii::LA::distributed::Vector<double, dealii::MemorySpace::Host> value(
      y.get_partitioner());
  value = 0.;
  // Apply the Thermal Operator.
  thermal_operator->vmult_add(value, y);

  // Multiply by the inverse of the mass matrix.
  value.scale(*thermal_operator->get_inverse_mass_matrix());

  timers[evol_time_eval_th_ph].stop();

  return value;
}

template <int dim, int fe_degree, typename MemorySpaceType,
          std::enable_if_t<
              std::is_same<MemorySpaceType, dealii::MemorySpace::Host>::value,
              int> = 0>
void init_dof_vector(
    double const value,
    dealii::LinearAlgebra::distributed::Vector<double, MemorySpaceType> &vector)
{
  unsigned int const local_size = vector.locally_owned_size();
  for (unsigned int i = 0; i < local_size; ++i)
    vector.local_element(i) = value;
}

#ifdef ADAMANTINE_HAVE_CUDA
template <int dim, int fe_degree, typename MemorySpaceType,
          std::enable_if_t<
              std::is_same<MemorySpaceType, dealii::MemorySpace::CUDA>::value,
              int> = 0>
dealii::LA::distributed::Vector<double, MemorySpaceType> vmult_and_scale(
    std::shared_ptr<ThermalOperatorBase<dim, MemorySpaceType>> const
        &thermal_operator,
    dealii::LA::distributed::Vector<double, MemorySpaceType> const &y,
    dealii::LA::distributed::Vector<double, dealii::MemorySpace::Host> &value,
    std::vector<Timer> &timers)
{
  dealii::LA::distributed::Vector<double, MemorySpaceType> value_dev(
      value.get_partitioner());
  value_dev.import(value, dealii::VectorOperation::insert);

  // Apply the Thermal Operator.
  thermal_operator->vmult_add(value_dev, y);

  // Multiply by the inverse of the mass matrix.
  value_dev.scale(*thermal_operator->get_inverse_mass_matrix());

  timers[evol_time_eval_th_ph].stop();

  return value_dev;
}

template <int dim, int fe_degree, typename MemorySpaceType,
          std::enable_if_t<
              std::is_same<MemorySpaceType, dealii::MemorySpace::CUDA>::value,
              int> = 0>
dealii::LA::distributed::Vector<double, MemorySpaceType>
evaluate_thermal_physics_impl(
    std::shared_ptr<ThermalOperatorBase<dim, MemorySpaceType>> const
        &thermal_operator,
    dealii::hp::FECollection<dim> const &fe_collection, double const t,
    dealii::DoFHandler<dim> const &dof_handler,
    std::vector<std::shared_ptr<HeatSource<dim>>> const &heat_sources,
    double current_source_height, BoundaryType boundary_type,
    std::shared_ptr<MaterialProperty<dim>> const &material_properties,
    dealii::AffineConstraints<double> const &affine_constraints,
    dealii::LA::distributed::Vector<double, MemorySpaceType> const &y,
    std::vector<Timer> &timers)
{
  timers[evol_time_eval_mat_prop].start();
  // TODO do this on the GPU
  dealii::LA::distributed::Vector<double, dealii::MemorySpace::Host> y_host(
      y.get_partitioner());
  y_host.import(y, dealii::VectorOperation::insert);
  thermal_operator->evaluate_material_properties(y_host);
  timers[evol_time_eval_mat_prop].stop();

  timers[evol_time_eval_th_ph].start();
  dealii::LA::distributed::Vector<double, dealii::MemorySpace::Host> source(
      y.get_partitioner());
  source = 0.;

  // Compute the source term.
  dealii::hp::QCollection<dim> source_q_collection;
  source_q_collection.push_back(dealii::QGauss<dim>(fe_degree + 1));
  source_q_collection.push_back(dealii::QGauss<dim>(1));
  dealii::hp::FEValues<dim> hp_fe_values(fe_collection, source_q_collection,
                                         dealii::update_quadrature_points |
                                             dealii::update_values |
                                             dealii::update_JxW_values);
  unsigned int const dofs_per_cell = fe_collection.max_dofs_per_cell();
  unsigned int const n_q_points = source_q_collection.max_n_quadrature_points();
  std::vector<dealii::types::global_dof_index> local_dof_indices(dofs_per_cell);
  dealii::QGauss<dim - 1> face_quadrature(fe_degree + 1);
  dealii::FEFaceValues<dim> fe_face_values(
      fe_collection[0], face_quadrature,
      dealii::update_values | dealii::update_quadrature_points |
          dealii::update_JxW_values);
  unsigned int const n_face_q_points = face_quadrature.size();
  dealii::Vector<double> cell_source(dofs_per_cell);

  // Loop over the locally owned cells with an active FE index of zero
  for (auto const &cell : dealii::filter_iterators(
           dof_handler.active_cell_iterators(),
           dealii::IteratorFilters::LocallyOwnedCell(),
           dealii::IteratorFilters::ActiveFEIndexEqualTo(0)))
  {
    cell_source = 0.;
    hp_fe_values.reinit(cell);
    dealii::FEValues<dim> const &fe_values =
        hp_fe_values.get_present_fe_values();

    for (unsigned int i = 0; i < dofs_per_cell; ++i)
    {
      for (unsigned int q = 0; q < n_q_points; ++q)
      {
        double const inv_rho_cp = thermal_operator->get_inv_rho_cp(cell, q);
        double quad_pt_source = 0.;
        dealii::Point<dim> const &q_point = fe_values.quadrature_point(q);
        for (auto &beam : heat_sources)
          quad_pt_source += beam->value(q_point, t, current_source_height);

        cell_source[i] += inv_rho_cp * quad_pt_source *
                          fe_values.shape_value(i, q) * fe_values.JxW(q);
      }
    }

    // If we don't have a adiabatic boundary conditions, we need to add boundary
    // conditions
    if (!(boundary_type & BoundaryType::adiabatic))
    {
      for (unsigned int f = 0; f < dealii::GeometryInfo<dim>::faces_per_cell;
           ++f)
      {
        // We need to add the boundary conditions on the faces on the boundary
        // but also on the faces at the interface with FE_Nothing
        auto const &face = cell->face(f);
        if ((face->at_boundary()) &&
            ((!face->at_boundary()) &&
             (cell->neighbor(f)->active_fe_index() != 0)))
        {
          double conv_temperature_infty = 0.;
          double conv_heat_transfer_coef = 0.;
          double rad_temperature_infty = 0.;
          double rad_heat_transfer_coef = 0.;
          if (boundary_type & BoundaryType::convective)
          {
            conv_temperature_infty = material_properties->get_cell_value(
                cell, Property::convection_temperature_infty);
            conv_heat_transfer_coef = material_properties->get_cell_value(
                cell, StateProperty::convection_heat_transfer_coef);
          }
          if (boundary_type & BoundaryType::radiative)
          {
            rad_temperature_infty = material_properties->get_cell_value(
                cell, Property::radiation_temperature_infty);
            rad_heat_transfer_coef = material_properties->get_cell_value(
                cell, StateProperty::radiation_heat_transfer_coef);
          }

          fe_face_values.reinit(cell, face);
          for (unsigned int i = 0; i < dofs_per_cell; ++i)
          {
            for (unsigned int q = 0; q < n_face_q_points; ++q)
            {
              double const inv_rho_cp =
                  thermal_operator->get_inv_rho_cp(cell, q);
              cell_source[i] +=
                  inv_rho_cp *
                  (conv_heat_transfer_coef * conv_temperature_infty +
                   rad_heat_transfer_coef * rad_temperature_infty) *
                  fe_face_values.shape_value(i, q) * fe_face_values.JxW(q);
            }
          }
        }
      }
    }
    cell->get_dof_indices(local_dof_indices);
    affine_constraints.distribute_local_to_global(cell_source,
                                                  local_dof_indices, source);
  }

  source.compress(dealii::VectorOperation::add);

  return vmult_and_scale<dim, fe_degree, MemorySpaceType>(thermal_operator, y,
                                                          source, timers);
}

template <int dim, int fe_degree, typename MemorySpaceType,
          std::enable_if_t<
              std::is_same<MemorySpaceType, dealii::MemorySpace::CUDA>::value,
              int> = 0>
void init_dof_vector(
    double const value,
    dealii::LinearAlgebra::distributed::Vector<double, MemorySpaceType> &vector)
{
  dealii::LA::distributed::Vector<double, dealii::MemorySpace::Host>
      vector_host(vector.get_partitioner());
  unsigned int const local_size = vector_host.local_size();
  for (unsigned int i = 0; i < local_size; ++i)
    vector_host.local_element(i) = value;

  vector.import(vector_host, dealii::VectorOperation::insert);
}
#endif
} // namespace

template <int dim, int fe_degree, typename MemorySpaceType,
          typename QuadratureType>
ThermalPhysics<dim, fe_degree, MemorySpaceType, QuadratureType>::ThermalPhysics(
    MPI_Comm const &communicator, boost::property_tree::ptree const &database,
    Geometry<dim> &geometry)
    : _boundary_type(BoundaryType::invalid), _geometry(geometry),
      _dof_handler(_geometry.get_triangulation())
{
  // Create the FECollection
  _fe_collection.push_back(dealii::FE_Q<dim>(fe_degree));
  _fe_collection.push_back(dealii::FE_Nothing<dim>());

  // Create the QCollection
  _q_collection.push_back(QuadratureType(fe_degree + 1));
  _q_collection.push_back(QuadratureType(1));

  // Create the material properties
  boost::property_tree::ptree const &material_database =
      database.get_child("materials");
  _material_properties.reset(new MaterialProperty<dim>(
      communicator, _geometry.get_triangulation(), material_database));

  // Create the heat sources
  boost::property_tree::ptree const &source_database =
      database.get_child("sources");
  // PropertyTreeInput sources.n_beams
  unsigned int const n_beams = source_database.get<unsigned int>("n_beams");
  _heat_sources.resize(n_beams);
  for (unsigned int i = 0; i < n_beams; ++i)
  {
    // PropertyTreeInput sources.beam_X.type
    boost::property_tree::ptree const &beam_database =
        source_database.get_child("beam_" + std::to_string(i));
    std::string type = beam_database.get<std::string>("type");
    if (type == "goldak")
    {
      _heat_sources[i] = std::make_shared<GoldakHeatSource<dim>>(beam_database);
    }
    else if (type == "electron_beam")
    {
      _heat_sources[i] =
          std::make_shared<ElectronBeamHeatSource<dim>>(beam_database);
    }
    else if (type == "cube")
    {
      _heat_sources[i] = std::make_shared<CubeHeatSource<dim>>(beam_database);
    }
    else
    {
      ASSERT_THROW(false, "Error: Beam type '" +
                              beam_database.get<std::string>("type") +
                              "' not recognized.");
    }
  }

  // Create the boundary condition type
  // PropertyTreeInput boundary.type
  std::string boundary_type_str = database.get<std::string>("boundary.type");
  // Parse the string
  size_t pos_str = 0;
  std::string boundary;
  std::string delimiter = ",";
  auto parse_boundary_type = [&](std::string const &boundary) {
    if (boundary == "adiabatic")
    {
      ASSERT_THROW(_boundary_type == BoundaryType::invalid,
                   "Adiabatic condition cannot be combined with another type.");
      _boundary_type = BoundaryType::adiabatic;
    }
    else
    {
      ASSERT_THROW(_boundary_type != BoundaryType::adiabatic,
                   "Adiabatic condition cannot be combined with another type.");

      if (boundary == "radiative")
      {
        _boundary_type |= BoundaryType::radiative;
      }
      else if (boundary == "convective")
      {
        _boundary_type |= BoundaryType::convective;
      }
      else
      {
        ASSERT_THROW(false, "Unknown boundary type.");
      }
    }
  };
  while ((pos_str = boundary_type_str.find(delimiter)) != std::string::npos)
  {
    boundary = boundary_type_str.substr(0, pos_str);
    parse_boundary_type(boundary);
    boundary_type_str.erase(0, pos_str + delimiter.length());
  }
  parse_boundary_type(boundary_type_str);

  // Create the thermal operator
  if (std::is_same<MemorySpaceType, dealii::MemorySpace::Host>::value)
    _thermal_operator =
        std::make_shared<ThermalOperator<dim, fe_degree, MemorySpaceType>>(
            communicator, _boundary_type, _material_properties, _heat_sources);
#if defined(ADAMANTINE_HAVE_CUDA) && defined(__CUDACC__)
  else
    _thermal_operator = std::make_shared<
        ThermalOperatorDevice<dim, fe_degree, MemorySpaceType>>(
        communicator, _material_properties);
#endif

  // Create the time stepping scheme
  boost::property_tree::ptree const &time_stepping_database =
      database.get_child("time_stepping");
  // PropertyTreeInput time_stepping.method
  std::string method = time_stepping_database.get<std::string>("method");
  std::transform(method.begin(), method.end(), method.begin(),
                 [](unsigned char c) { return std::tolower(c); });
  if (method.compare("forward_euler") == 0)
    _time_stepping =
        std::make_unique<dealii::TimeStepping::ExplicitRungeKutta<LA_Vector>>(
            dealii::TimeStepping::FORWARD_EULER);
  else if (method.compare("rk_third_order") == 0)
    _time_stepping =
        std::make_unique<dealii::TimeStepping::ExplicitRungeKutta<LA_Vector>>(
            dealii::TimeStepping::RK_THIRD_ORDER);
  else if (method.compare("rk_fourth_order") == 0)
    _time_stepping =
        std::make_unique<dealii::TimeStepping::ExplicitRungeKutta<LA_Vector>>(
            dealii::TimeStepping::RK_CLASSIC_FOURTH_ORDER);
  else if (method.compare("heun_euler") == 0)
  {
    _time_stepping = std::make_unique<
        dealii::TimeStepping::EmbeddedExplicitRungeKutta<LA_Vector>>(
        dealii::TimeStepping::HEUN_EULER);
    _embedded_method = true;
  }
  else if (method.compare("bogacki_shampine") == 0)
  {
    _time_stepping = std::make_unique<
        dealii::TimeStepping::EmbeddedExplicitRungeKutta<LA_Vector>>(
        dealii::TimeStepping::BOGACKI_SHAMPINE);
    _embedded_method = true;
  }
  else if (method.compare("dopri") == 0)
  {
    _time_stepping = std::make_unique<
        dealii::TimeStepping::EmbeddedExplicitRungeKutta<LA_Vector>>(
        dealii::TimeStepping::DOPRI);
    _embedded_method = true;
  }
  else if (method.compare("fehlberg") == 0)
  {
    _time_stepping = std::make_unique<
        dealii::TimeStepping::EmbeddedExplicitRungeKutta<LA_Vector>>(
        dealii::TimeStepping::FEHLBERG);
    _embedded_method = true;
  }
  else if (method.compare("cash_karp") == 0)
  {
    _time_stepping = std::make_unique<
        dealii::TimeStepping::EmbeddedExplicitRungeKutta<LA_Vector>>(
        dealii::TimeStepping::CASH_KARP);
    _embedded_method = true;
  }
  else if (method.compare("backward_euler") == 0)
  {
    _time_stepping =
        std::make_unique<dealii::TimeStepping::ImplicitRungeKutta<LA_Vector>>(
            dealii::TimeStepping::BACKWARD_EULER);
    _implicit_method = true;
  }
  else if (method.compare("implicit_midpoint") == 0)
  {
    _time_stepping =
        std::make_unique<dealii::TimeStepping::ImplicitRungeKutta<LA_Vector>>(
            dealii::TimeStepping::IMPLICIT_MIDPOINT);
    _implicit_method = true;
  }
  else if (method.compare("crank_nicolson") == 0)
  {
    _time_stepping =
        std::make_unique<dealii::TimeStepping::ImplicitRungeKutta<LA_Vector>>(
            dealii::TimeStepping::CRANK_NICOLSON);
    _implicit_method = true;
  }
  else if (method.compare("sdirk2") == 0)
  {
    _time_stepping =
        std::make_unique<dealii::TimeStepping::ImplicitRungeKutta<LA_Vector>>(
            dealii::TimeStepping::SDIRK_TWO_STAGES);
    _implicit_method = true;
  }

  if (_embedded_method == true)
  {
    // PropertyTreeInput time_steppping.coarsening_parameter
    double coarsen_param =
        time_stepping_database.get("coarsening_parameter", 1.2);
    // PropertyTreeInput time_steppping.refining_parameter
    double refine_param = time_stepping_database.get("refining_parameter", 0.8);
    // PropertyTreeInput time_stepping.min_time_step
    double min_delta = time_stepping_database.get("min_time_step", 1e-14);
    // PropertyTreeInput time_stepping.max_time_step
    double max_delta = time_stepping_database.get("max_time_step", 1e100);
    // PropertyTreeInput time_stepping.refining_tolerance
    double refine_tol = time_stepping_database.get("refining_tolerance", 1e-8);
    // PropertyTreeInput time_stepping.coarsening_tolerance
    double coarsen_tol =
        time_stepping_database.get("coarsening_tolerance", 1e-12);
    dealii::TimeStepping::EmbeddedExplicitRungeKutta<LA_Vector> *embedded_rk =
        static_cast<
            dealii::TimeStepping::EmbeddedExplicitRungeKutta<LA_Vector> *>(
            _time_stepping.get());
    embedded_rk->set_time_adaptation_parameters(coarsen_param, refine_param,
                                                min_delta, max_delta,
                                                refine_tol, coarsen_tol);
  }

  // If the time stepping scheme is implicit, set the parameters for the solver
  // and create the implicit operator.
  if (_implicit_method == true)
  {
    // PropertyTreeInput time_stepping.max_iteration
    _max_iter = time_stepping_database.get("max_iteration", 1000);
    // PropertyTreeInput time_stepping.tolerance
    _tolerance = time_stepping_database.get("tolerance", 1e-12);
    // PropertyTreeInput time_stepping.right_preconditioning
    _right_preconditioning =
        time_stepping_database.get("right_preconditioning", false);
    // PropertyTreeInput time_stepping.n_tmp_vectors
    _max_n_tmp_vectors = time_stepping_database.get("n_tmp_vectors", 30);
    // PropertyTreeInput time_stepping.newton_max_iteration
    unsigned int newton_max_iter =
        time_stepping_database.get("newton_max_iteration", 100);
    // PropertyTreeInput time_stepping.newton_tolerance
    double newton_tolerance =
        time_stepping_database.get("newton_tolerance", 1e-6);
    dealii::TimeStepping::ImplicitRungeKutta<LA_Vector> *implicit_rk =
        static_cast<dealii::TimeStepping::ImplicitRungeKutta<LA_Vector> *>(
            _time_stepping.get());
    implicit_rk->set_newton_solver_parameters(newton_max_iter,
                                              newton_tolerance);

    // PropertyTreeInput time_stepping.jfnk
    bool jfnk = time_stepping_database.get("jfnk", false);
    _implicit_operator = std::make_unique<ImplicitOperator<MemorySpaceType>>(
        _thermal_operator, jfnk);
  }

  // Set material on part of the domain
  // PropertyTreeInput geometry.material_height
  double const material_height = database.get("geometry.material_height", 1e9);
  for (auto const &cell :
       dealii::filter_iterators(_dof_handler.active_cell_iterators(),
                                dealii::IteratorFilters::LocallyOwnedCell()))
  {
    // If the center of the cell is below material_height, it contains material
    // otherwise it does not.
    if (cell->center()[axis<dim>::z] < material_height)
      cell->set_active_fe_index(0);
    else
      cell->set_active_fe_index(1);
  }

  // Set the initial height of the heat source. Right now this is just the
  // maximum heat source height, which can lead to unexpected behavior for
  // different sources with different heights.
  double temp_height = std::numeric_limits<double>::lowest();
  for (auto const &source : _heat_sources)
  {
    temp_height = std::max(temp_height, source->get_current_height(0.0));
  }
  _current_source_height = temp_height;
}

template <int dim, int fe_degree, typename MemorySpaceType,
          typename QuadratureType>
void ThermalPhysics<dim, fe_degree, MemorySpaceType,
                    QuadratureType>::setup_dofs()
{
  _dof_handler.distribute_dofs(_fe_collection);
  dealii::IndexSet locally_relevant_dofs;
  dealii::DoFTools::extract_locally_relevant_dofs(_dof_handler,
                                                  locally_relevant_dofs);
  _affine_constraints.clear();
  _affine_constraints.reinit(locally_relevant_dofs);
  dealii::DoFTools::make_hanging_node_constraints(_dof_handler,
                                                  _affine_constraints);
  _affine_constraints.close();

  _thermal_operator->reinit(_dof_handler, _affine_constraints, _q_collection);
}

template <int dim, int fe_degree, typename MemorySpaceType,
          typename QuadratureType>
void ThermalPhysics<dim, fe_degree, MemorySpaceType,
                    QuadratureType>::compute_inverse_mass_matrix()
{
  _thermal_operator->compute_inverse_mass_matrix(
      _dof_handler, _affine_constraints, _fe_collection);
  if (_implicit_method == true)
    _implicit_operator->set_inverse_mass_matrix(
        _thermal_operator->get_inverse_mass_matrix());
}

template <int dim, int fe_degree, typename MemorySpaceType,
          typename QuadratureType>
void ThermalPhysics<dim, fe_degree, MemorySpaceType, QuadratureType>::
    add_material(
        std::vector<
            typename dealii::DoFHandler<dim>::active_cell_iterator> const
            &elements_to_activate,
        double new_material_temperature,
        dealii::LA::distributed::Vector<double, MemorySpaceType> &solution)
{
#ifdef ADAMANTINE_WITH_CALIPER
  CALI_CXX_MARK_FUNCTION;
#endif
  std::vector<dealii::Vector<double>> data_to_transfer;
  unsigned int const dofs_per_cell = _dof_handler.get_fe().n_dofs_per_cell();
  dealii::Vector<double> cell_solution(dofs_per_cell);
  dealii::Vector<double> dummy_cell_solution(dofs_per_cell);
  for (auto &val : dummy_cell_solution)
  {
    val = std::numeric_limits<double>::infinity();
  }
  solution.update_ghost_values();
  std::vector<dealii::types::global_dof_index> dof_indices(dofs_per_cell);
  for (auto const &cell : _dof_handler.active_cell_iterators())
  {
    if ((cell->is_locally_owned()) && (cell->active_fe_index() == 0))
    {
      cell->get_dof_values(solution, cell_solution);
      data_to_transfer.push_back(cell_solution);
    }
    else
    {
      data_to_transfer.push_back(dummy_cell_solution);
    }
  }

  // Activate elements by updating the fe_index
  for (auto const &cell : elements_to_activate)
  {
    cell->set_active_fe_index(0);
  }

  dealii::parallel::distributed::Triangulation<dim> &triangulation =
      dynamic_cast<dealii::parallel::distributed::Triangulation<dim> &>(
          const_cast<dealii::Triangulation<dim> &>(
              _dof_handler.get_triangulation()));
  triangulation.prepare_coarsening_and_refinement();
  dealii::parallel::distributed::CellDataTransfer<
      dim, dim, std::vector<dealii::Vector<double>>>
      cell_data_trans(triangulation);
  cell_data_trans.prepare_for_coarsening_and_refinement(data_to_transfer);
  triangulation.execute_coarsening_and_refinement();

  setup_dofs();

  // Update MaterialProperty DoFHandler and resize the state vectors
  _material_properties->reinit_dofs();

  // Recompute the inverse of the mass matrix
  compute_inverse_mass_matrix();

  initialize_dof_vector(new_material_temperature, solution);
  std::vector<dealii::Vector<double>> transferred_data(
      triangulation.n_active_cells(), dealii::Vector<double>(dofs_per_cell));
  cell_data_trans.unpack(transferred_data);

  unsigned int cell_i = 0;
  for (auto const &cell : _dof_handler.active_cell_iterators())
  {
    if ((cell->is_locally_owned()) && (transferred_data[cell_i][0] !=
                                       std::numeric_limits<double>::infinity()))
    {
      cell->set_dof_values(transferred_data[cell_i], solution);
    }
    ++cell_i;
  }

  // Communicate the results.
  solution.compress(dealii::VectorOperation::min);

  // Set the value to the newly create DoFs. Here we need to be careful with the
  // hanging nodes. When there is a hanging node, the dofs at the vertices are
  // "doubled": there is a dof associated to the coarse cell and a dof
  // associated to the fine cell. The final value is decided by
  // AffineConstraints. Thus, we need to make sure that the newly activated
  // cells are at the same level than their neighbors.
  std::for_each(std::execution::par_unseq, solution.begin(), solution.end(),
                [&](double &val) {
                  if (val == std::numeric_limits<double>::infinity())
                  {
                    val = new_material_temperature;
                  }
                });
}

template <int dim, int fe_degree, typename MemorySpaceType,
          typename QuadratureType>
double ThermalPhysics<dim, fe_degree, MemorySpaceType, QuadratureType>::
    evolve_one_time_step(
        double t, double delta_t,
        dealii::LA::distributed::Vector<double, MemorySpaceType> &solution,
        std::vector<Timer> &timers)
{

  // Update the height of the heat source. Right now this is just the
  // maximum heat source height, which can lead to unexpected behavior for
  // different sources with different heights.
  double temp_height = std::numeric_limits<double>::lowest();
  for (auto const &source : _heat_sources)
  {
    temp_height = std::max(temp_height, source->get_current_height(t));
  }
  _current_source_height = temp_height;

  auto eval = [&](double const t, LA_Vector const &y) {
    return evaluate_thermal_physics(t, y, timers);
  };
  auto id_m_Jinv = [&](double const t, double const tau, LA_Vector const &y) {
    return id_minus_tau_J_inverse(t, tau, y, timers);
  };

  double time = _time_stepping->evolve_one_time_step(eval, id_m_Jinv, t,
                                                     delta_t, solution);

  // If the method is embedded, get the next time step. Otherwise, just use the
  // current time step.
  if (_embedded_method == false)
    _delta_t_guess = delta_t;
  else
  {
    dealii::TimeStepping::EmbeddedExplicitRungeKutta<LA_Vector> *embedded_rk =
        static_cast<
            dealii::TimeStepping::EmbeddedExplicitRungeKutta<LA_Vector> *>(
            _time_stepping.get());
    _delta_t_guess = embedded_rk->get_status().delta_t_guess;
  }

  // Return the time at the end of the time step. This may be different than
  // t+delta_t for embedded methods.
  return time;
}

template <int dim, int fe_degree, typename MemorySpaceType,
          typename QuadratureType>
void ThermalPhysics<dim, fe_degree, MemorySpaceType, QuadratureType>::
    initialize_dof_vector(
        dealii::LA::distributed::Vector<double, MemorySpaceType> &vector) const
{
  _thermal_operator->initialize_dof_vector(vector);
}

template <int dim, int fe_degree, typename MemorySpaceType,
          typename QuadratureType>
void ThermalPhysics<dim, fe_degree, MemorySpaceType, QuadratureType>::
    initialize_dof_vector(
        double const value,
        dealii::LA::distributed::Vector<double, MemorySpaceType> &vector) const
{
  // Resize the vector
  _thermal_operator->initialize_dof_vector(vector);

  init_dof_vector<dim, fe_degree, MemorySpaceType>(value, vector);
}

template <int dim, int fe_degree, typename MemorySpaceType,
          typename QuadratureType>
void ThermalPhysics<dim, fe_degree, MemorySpaceType,
                    QuadratureType>::get_state_from_material_properties()
{
  _thermal_operator->get_state_from_material_properties();
}

template <int dim, int fe_degree, typename MemorySpaceType,
          typename QuadratureType>
void ThermalPhysics<dim, fe_degree, MemorySpaceType,
                    QuadratureType>::set_state_to_material_properties()
{
  _thermal_operator->set_state_to_material_properties();
}

template <int dim, int fe_degree, typename MemorySpaceType,
          typename QuadratureType>
dealii::LA::distributed::Vector<double, MemorySpaceType>
ThermalPhysics<dim, fe_degree, MemorySpaceType, QuadratureType>::
    evaluate_thermal_physics(
        double const t,
        dealii::LA::distributed::Vector<double, MemorySpaceType> const &y,
        std::vector<Timer> &timers) const
{
#ifdef ADAMANTINE_WITH_CALIPER
  CALI_CXX_MARK_FUNCTION;
#endif
  if constexpr (std::is_same<MemorySpaceType, dealii::MemorySpace::Host>::value)
  {
    _thermal_operator->evaluate_material_properties(y);
    return evaluate_thermal_physics_impl<dim, fe_degree, MemorySpaceType>(
        _thermal_operator, t, _current_source_height, y, timers);
  }
  else
  {
    // TODO do this on the GPU
    dealii::LA::distributed::Vector<double, dealii::MemorySpace::Host> y_host(
        y.get_partitioner());
    y_host.import(y, dealii::VectorOperation::insert);
    _thermal_operator->evaluate_material_properties(y_host);
    return evaluate_thermal_physics_impl<dim, fe_degree, MemorySpaceType>(
        _thermal_operator, _fe_collection, t, _dof_handler, _heat_sources,
        _current_source_height, _boundary_type, _material_properties,
        _affine_constraints, y, timers);
  }

  // Dummy to silence warning
  return dealii::LA::distributed::Vector<double, MemorySpaceType>();
}

template <int dim, int fe_degree, typename MemorySpaceType,
          typename QuadratureType>
dealii::LA::distributed::Vector<double, MemorySpaceType>
ThermalPhysics<dim, fe_degree, MemorySpaceType, QuadratureType>::
    id_minus_tau_J_inverse(
        double const /*t*/, double const tau,
        dealii::LA::distributed::Vector<double, MemorySpaceType> const &y,
        std::vector<Timer> &timers) const
{
  timers[evol_time_J_inv].start();
  _implicit_operator->set_tau(tau);
  dealii::LA::distributed::Vector<double, MemorySpaceType> solution(
      y.get_partitioner());

  // TODO Add a geometric multigrid preconditioner.
  dealii::PreconditionIdentity preconditioner;

  dealii::SolverControl solver_control(_max_iter, _tolerance * y.l2_norm());
  // We need to inverse (I - tau M^{-1} J). While M^{-1} and J are SPD,
  // (I - tau M^{-1} J) is symmetric indefinite in the general case.
  typename dealii::SolverGMRES<
      dealii::LA::distributed::Vector<double, MemorySpaceType>>::AdditionalData
      additional_data(_max_n_tmp_vectors, _right_preconditioning);
  dealii::SolverGMRES<dealii::LA::distributed::Vector<double, MemorySpaceType>>
      solver(solver_control, additional_data);
  solver.solve(*_implicit_operator, solution, y, preconditioner);

  timers[evol_time_J_inv].stop();

  return solution;
}
} // namespace adamantine

#endif

ext: hh

max_stars_repo_name: Rombur/adamantine

max_issues_repo_name: Rombur/adamantine

max_forks_repo_name: Rombur/adamantine


text:
/*
 * This file is a part of the TChecker project.
 *
 * See files AUTHORS and LICENSE for copyright details.
 *
 */

#ifndef TCHECKER_PARSING_DECLARATION_HH
#define TCHECKER_PARSING_DECLARATION_HH

#include <cassert>
#include <iostream>
#include <string>
#include <unordered_map>
#include <vector>

#include <boost/core/noncopyable.hpp>

#include "tchecker/basictypes.hh"
#include "tchecker/utils/iterator.hh"

/*!
 \file declaration.hh
 \brief System declarations
 */

namespace tchecker {

  namespace parsing {

    /* Forward declaration */

    class declaration_visitor_t;

    /*!
     \class declaration_t
     \brief Declaration from input file
     \note Abstract base class
     */
    class declaration_t : private boost::noncopyable {
    public:
      /*!
       \brief Constructor
       */
      declaration_t() = default;

      /*!
       \brief Destructor
       */
      virtual ~declaration_t() = default;

      /*!
       \brief Assignment operator (DELETED)
       */
      tchecker::parsing::declaration_t & operator= (tchecker::parsing::declaration_t const &) = delete;

      /*!
       \brief Move assignment operator (DELETED)
       */
      tchecker::parsing::declaration_t & operator= (tchecker::parsing::declaration_t &&) = delete;

      /*!
       \brief Clone
       \return Returns a clone of this
       */
      tchecker::parsing::declaration_t * clone() const;

      /*!
       \brief Visit
       \param v : visitor
       \post this has been visited by v
       */
      void visit(tchecker::parsing::declaration_visitor_t & v) const;
    private:
      /*!
       \brief Clone
       */
      virtual tchecker::parsing::declaration_t * do_clone() const = 0;

      /*!
       \brief Visit
       */
      virtual void do_visit(tchecker::parsing::declaration_visitor_t & v) const = 0;

      /*!
       \brief Output the declaration
       \param os : output stream
       */
      virtual std::ostream & do_output(std::ostream & os) const = 0;

      friend std::ostream & operator<< (std::ostream & os, tchecker::parsing::declaration_t const & decl);
    };

    /*!
     \brief Output declaration
     \param os : output stream
     \param decl : declaration
     \post decl has been output to os
     \return os after outputing decl has been output
     */
    inline std::ostream & operator<< (std::ostream & os, tchecker::parsing::declaration_t const & decl)
    {
      return decl.do_output(os);
    }

    /*!
     \class inner_declaration_t
     */
    class inner_declaration_t : public tchecker::parsing::declaration_t {
    public:
      /*!
       \brief Destructor
       */
      virtual ~inner_declaration_t() = default;
    };

    /*!
     \class clock_declaration_t
     \brief Declaration of an (array of) clock variable
     */
    class clock_declaration_t final : public tchecker::parsing::inner_declaration_t {
    public:
      /*!
       \brief Constructor
       \param name : clock name
       \param size : size of array
       \pre name is not empty and size >= 1
       \throw std::invalid_argument : if name is empty or size < 1
       */
      clock_declaration_t(std::string const & name, unsigned int size);

      /*!
       \brief Destructor
       */
      virtual ~clock_declaration_t() = default;

      /*!
       \brief Accessor
       \return Name
       */
      inline std::string const & name() const
      {
        return _name;
      }

      /*!
       \brief Accessor
       \return Size
       */
      inline unsigned int size() const
      {
        return _size;
      }
    private:
      /*!
       \brief Clone
       \return A clone of this
       */
      virtual tchecker::parsing::declaration_t * do_clone() const;

      /*!
       \brief Visit
       \post this has been visited by v
       */
      virtual void do_visit(tchecker::parsing::declaration_visitor_t & v) const;

      /*!
       \brief Output the declaration
       \param os : output stream
       \post this declaration has been output to os
       \return os after this declaration has been output
       */
      virtual std::ostream & do_output(std::ostream & os) const;

      std::string const _name;    /*!< Name */
      unsigned int const _size;   /*!< Size */
    };

    /*!
     \class int_declaration_t
     \brief Declaration of an (array of) int variable
     */
    class int_declaration_t final : public tchecker::parsing::inner_declaration_t {
    public:
      /*!
       \brief Constructor
       \param name : int name
       \param size : size of array
       \param min : minimal value
       \param max : maximal value
       \param init : initial value
       \pre name is not empty, size >= 1, min <= init <= max
       \throws std::invalid_argument : name is empty, or if size < 1, or
       not (min <= init <= max)
       */
      int_declaration_t(std::string const & name,
                        unsigned int size,
                        tchecker::integer_t min,
                        tchecker::integer_t max,
                        tchecker::integer_t init);

      /*!
       \brief Destructor
       */
      virtual ~int_declaration_t() = default;

      /*!
       \brief Accessor
       \return Name
       */
      inline std::string const & name() const
      {
        return _name;
      }

      /*!
       \brief Accessor
       \return Size
       */
      inline unsigned int size() const
      {
        return _size;
      }

      /*!
       \brief Accessor
       \return Min value
       */
      inline tchecker::integer_t min() const
      {
        return _min;
      }

      /*!
       \brief Accessor
       \return Max value
       */
      inline tchecker::integer_t max() const
      {
        return _max;
      }

      /*!
       \brief Accessor
       \return Initial value
       */
      inline tchecker::integer_t init() const
      {
        return _init;
      }
    private:
      /*!
       \brief Clone
       \returns A clone of this
       */
      virtual tchecker::parsing::declaration_t * do_clone() const;

      /*!
       \brief Visit
       \post this has been visited by v
       */
      virtual void do_visit(tchecker::parsing::declaration_visitor_t & v) const;

      /*!
       \brief Output the declaration
       \param os : output stream
       \post this declaration has been output to os
       \return os after this declaration has been output
       */
      virtual std::ostream & do_output(std::ostream & os) const;

      std::string const _name;           /*!< Name */
      unsigned int const _size;          /*!< Size */
      tchecker::integer_t const _min;    /*!< Min value */
      tchecker::integer_t const _max;    /*!< Max value */
      tchecker::integer_t const _init;   /*!< Initial value */
    };

    /*
     \class process_declaration_t
     \brief Declaration of a process
     */
    class process_declaration_t final : public tchecker::parsing::inner_declaration_t {
    public:
      /*!
       \brief Constructor
       \param name : process name
       \pre name is not empty
       \throws std::invalid_argument : if name is empty
       */
      process_declaration_t(std::string const & name);

      /*!
       \brief Destructor
       */
      virtual ~process_declaration_t() = default;

      /*!
       \brief Accessor
       \return Name
       */
      inline std::string const & name() const
      {
        return _name;
      }
    private:
      /*!
       \brief Clone
       \return A clone of this
       */
      virtual tchecker::parsing::declaration_t * do_clone() const;

      /*!
       \brief Visit
       \post this has been visited by v
       */
      virtual void do_visit(tchecker::parsing::declaration_visitor_t & v) const;

      /*!
       \brief Output the declaration
       \param os : output stream
       \post this declaration has been output to os
       \return os after this declaration has been output
       */
      virtual std::ostream & do_output(std::ostream & os) const;

      std::string const _name;  /*!< Process name */
    };

    /*!
     \class event_declaration_t
     \brief Declaration of an event
     */
    class event_declaration_t final : public tchecker::parsing::inner_declaration_t {
    public:
      /*!
       \brief Constructor
       \param name : event name
       \pre name is not empty
       \throws std::invalid_argument : if name is empty
       */
      event_declaration_t(std::string const & name);

      /*!
       \brief Destructor
       */
      virtual ~event_declaration_t() = default;

      /*!
       \brief Accessor
       \return Name
       */
      inline std::string const & name() const
      {
        return _name;
      }
    private:
      /*!
       \brief Clone
       \return A clone of this
       */
      virtual tchecker::parsing::declaration_t * do_clone() const;

      /*!
       \brief Visit
       \post this has been visited by v
       */
      virtual void do_visit(tchecker::parsing::declaration_visitor_t & v) const;

      /*!
       \brief Output the declaration
       \param os : output stream
       \post this declaration has been output to os
       \return os after this declaration has been output
       */
      virtual std::ostream & do_output(std::ostream & os) const;

      std::string const _name;   /*!< Event name */
    };

    /*!
     \class attr_t
     \brief Attribute for locations and edges
     */
    class attr_t {
    public:
      /*!
       \brief Constructor
       \param key : attribute key
       \param value : attribute value
       \param key_context : contextual information for key (position in input
       file, etc)
       \param value_context : contextual information for value (position in
       input file, etc)
       */
      attr_t(std::string const & key,
             std::string const & value,
             std::string const & key_context,
             std::string const & value_context);

      /*!
       \brief Copy constructor
       \param attr : attribute
       \post this is a copy of attr
       */
      attr_t(tchecker::parsing::attr_t const & attr) = default;

      /*!
       \brief Move constructor
       \param attr : attribute
       \post attr has been moved to this
       */
      attr_t(tchecker::parsing::attr_t && attr) = default;

      /*!
       \brief Destructor
       */
      ~attr_t() = default;

      /*!
       \brief Assignment operator
       \param attr : attribute
       \post this is a copy of attr
       */
      tchecker::parsing::attr_t & operator= (tchecker::parsing::attr_t const & attr) = default;

      /*!
       \brief Move assignment operator
       \param attr : attribute
       \post attr has been moved to this
       */
      tchecker::parsing::attr_t & operator= (tchecker::parsing::attr_t && attr) = default;

      /*!
       \brief Accessor
       \return Attribute key
       */
      inline std::string const & key() const
      {
        return _key;
      }

      /*!
       \brief Accessor
       \return Attribute value
       */
      inline std::string const & value() const
      {
        return _value;
      }

      /*!
       \brief Accessor
       \return Contextual information for key
       */
      inline std::string const & key_context() const
      {
        return _key_context;
      }

      /*!
       \brief Accessor
       \return Contextual information for value
       */
      inline std::string const & value_context() const
      {
        return _value_context;
      }
    private:
      std::string _key;            /*!< Key */
      std::string _value;          /*!< Value */
      std::string _key_context;    /*!< Contextual information for key */
      std::string _value_context;  /*!< Contextual information for value */
    };

    /*!
     \brief Output operator
     \param os : output stream
     \param attr : attribute
     \post attr has been output to os
     \return os after attr has been output
     */
    std::ostream & operator<< (std::ostream & os, tchecker::parsing::attr_t const & attr);

    /*!
     \class attributes_t
     \brief Attributes map for locations and edges
     */
    class attributes_t {
    public:
      /*!
       \brief Type of attributes map
       */
      using map_t = std::unordered_map<std::string, tchecker::parsing::attr_t const *>;

      /*!
       \brief Constructor
       \post this is an empty attributes map
       */
      attributes_t() = default;

      /*!
       \brief Copy constructor
       \param attr : attributes to copy from
       \post this is a copy of attr
       */
      attributes_t(attributes_t const & attr);

      /*!
       \brief Move constructor
       \param attr : attributes to move
       \post attr has been moved to this
       */
      attributes_t(attributes_t && attr);

      /*!
       \brief Destructor
       \post All attributes in the map have been deleted
       */
      ~attributes_t();

      /*!
       \brief Assignement operator
       \param attr : attributes to assign from
       \post this is a copy of attr
       \return this
       */
      attributes_t & operator= (attributes_t const & attr);

      /*!
       \brief Move assignement operator
       \param attr : attributes to move
       \post attr has been moved to this
       \return this
       */
      attributes_t & operator= (attributes_t && attr);

      /*!
       \brief Empties the map
       \post this map is empty
       */
      void clear();

      /*!
       \brief Emptiness check
       \return true if this map is empty, false otherwise
       */
      bool empty() const;

      /*!
       \brief Accessor
       \return size of the map
       */
      std::size_t size() const;

      /*!
       \brief Insert
       \param attr : attribute to insert in the map
       \return true if attr has been inserted in the map, false
       otherwise (when an attribute with the same key is already in the map)
       \note this takes ownership on attr
       */
      bool insert(tchecker::parsing::attr_t const * attr);

      /*!
       \class const_iterator_t
       \brief Const iterator on the map
       */
      class const_iterator_t : public map_t::const_iterator {
      public:
        /*!
         \brief Constructor
         \param it : map iterator
         \post this iterator point to it
         */
        const_iterator_t(map_t::const_iterator const & it);

        /*!
         \brief Accessor
         \return attribute at pointed map entry
         */
        inline tchecker::parsing::attr_t const & operator*() const
        {
          return *(*this)->second;
        }
      };

      /*!
       \brief Accessor
       \param key : key of attribute to access
       \return attribute at key in the map if any
       \throw std::out_of_range if there is no attibute with
       key in map
       */
      tchecker::parsing::attr_t const & at(std::string const & key) const;

      /*!
       \brief Accessor
       \param key : key of attribute to access
       \return iterator on atribute at key if any, end()
       otherwise
       */
      const_iterator_t find(std::string const & key) const;

      /*!
       \brief Accessor
       \return iterator to the first entry of the map (if any)
       */
      const_iterator_t begin() const;

      /*!
       \brief Accessor
       \return iteratot to the past-the-end "entry" of the map
       */
      const_iterator_t end() const;
    private:
      map_t _attr;  /*!< Attributes */
    };

    /*!
     \brief Output operator
     \param os : output stream
     \param attr : attributes map
     \post attr has been output to os
     \return os after attr has been output
     */
    std::ostream & operator<< (std::ostream & os, tchecker::parsing::attributes_t const & attr);

    /*!
     \class location_declaration_t
     \brief Declaration of a location
     */
    class location_declaration_t final : public tchecker::parsing::inner_declaration_t {
    public:
      /*!
       \brief Constructor
       \param name : location name
       \param process : location process
       \param attr : location attributes
       \pre name is not empty
       \throws std::invalid_argument : if name is empty
       \note attr is moved to this.
       */
      location_declaration_t(std::string const & name,
                             tchecker::parsing::process_declaration_t const & process,
                             tchecker::parsing::attributes_t && attr);

      /*!
       \brief Destructor
       */
      virtual ~location_declaration_t() = default;

      /*!
       \brief Accessor
       \return Name
       */
      inline std::string const & name() const
      {
        return _name;
      }

      /*!
       \brief Accessor
       \return Location process
       */
      inline tchecker::parsing::process_declaration_t const & process() const
      {
        return _process;
      }

      /*!
       \brief Accessor
       \return Attributes
       */
      inline tchecker::parsing::attributes_t const & attributes() const
      {
        return _attr;
      }
    private:
      /*!
       \brief Clone
       \return A clone of this
       \note The clone shares pointers with this
       */
      virtual tchecker::parsing::declaration_t * do_clone() const;

      /*!
       \brief Visit
       \post this has been visited by v
       */
      virtual void do_visit(tchecker::parsing::declaration_visitor_t & v) const;

      /*!
       \brief Output the declaration
       \param os : output stream
       \post this declaration has been output to os
       \return os after this declaration has been output
       */
      virtual std::ostream & do_output(std::ostream & os) const;

      std::string const _name;                                    /*!< Name */
      tchecker::parsing::process_declaration_t const & _process;  /*!< Process */
      tchecker::parsing::attributes_t const _attr;                /*!< Attributes */
    };

    /*!
     \class edge_declaration_t
     \brief Declaration of an edge
     */
    class edge_declaration_t final : public tchecker::parsing::inner_declaration_t {
    public:
      /*!
       \brief Constructor
       \param process : process
       \param src : source location
       \param tgt : target location
       \param event : edge event
       \param attr : edge attributes
       \pre src and tgt belong to process
       \throws std::invalid_argument : if src and tgt do not belong to process
       \note attr is moved to this
       */
      edge_declaration_t(tchecker::parsing::process_declaration_t const & process,
                         tchecker::parsing::location_declaration_t const & src,
                         tchecker::parsing::location_declaration_t const & tgt,
                         tchecker::parsing::event_declaration_t const & event,
                         tchecker::parsing::attributes_t && attr);

      /*!
       \brief Destructor
       */
      virtual ~edge_declaration_t() = default;

      /*!
       \brief Accessor
       \return Process
       */
      inline tchecker::parsing::process_declaration_t const & process() const
      {
        return _process;
      }

      /*!
       \brief Accessor
       \return Source location
       */
      inline tchecker::parsing::location_declaration_t const & src() const
      {
        return _src;
      }

      /*!
       \brief Accessor
       \return Target location
       */
      inline tchecker::parsing::location_declaration_t const & tgt() const
      {
        return _tgt;
      }

      /*!
       \brief Accessor
       \return Event
       */
      inline tchecker::parsing::event_declaration_t const & event() const
      {
        return _event;
      }

      /*!
       \brief Accessor
       \return Attributes
       */
      inline tchecker::parsing::attributes_t const & attributes() const
      {
        return _attr;
      }
    private:
      /*!
       \brief Clone
       \return A clone of this
       \note The clone shares pointers with this
       */
      virtual tchecker::parsing::declaration_t * do_clone() const;

      /*!
       \brief Visit
       \post this has been visited by v
       */
      virtual void do_visit(tchecker::parsing::declaration_visitor_t & v) const;

      /*!
       \brief Output the declaration
       \param os : output stream
       \post this declaration has been output to os
       \return os after this declaration has been output
       */
      virtual std::ostream & do_output(std::ostream & os) const;

      tchecker::parsing::process_declaration_t const & _process;  /*!< Process */
      tchecker::parsing::location_declaration_t const & _src;     /*!< Source loc */
      tchecker::parsing::location_declaration_t const & _tgt;     /*!< Target loc */
      tchecker::parsing::event_declaration_t const & _event;      /*!< Event */
      tchecker::parsing::attributes_t const _attr;                /*!< Attributes */
    };

    /*!
     \class sync_constraint_t
     \brief Synchronization constraint
     */
    class sync_constraint_t {
    public:
      /*!
       \brief Constructor
       \param process : synchronized process
       \param event : synchronized event
       \param strength : strength of synchronization
       */
      sync_constraint_t(tchecker::parsing::process_declaration_t const & process,
                        tchecker::parsing::event_declaration_t const & event,
                        enum tchecker::sync_strength_t strength);

      /*!
       \brief Copy constructor
       \param sync : synchronization constraint
       \post this is copy of sync
       */
      sync_constraint_t(tchecker::parsing::sync_constraint_t const & sync) = default;

      /*!
       \brief Move constructor
       \param sync : synchronization constraint
       \post sync has been moved to this
       */
      sync_constraint_t(tchecker::parsing::sync_constraint_t && sync) = default;

      /*!
       \brief Destructor
       */
      ~sync_constraint_t() = default;

      /*!
       \brief Assignment operator (DELETED)
       */
      tchecker::parsing::sync_constraint_t & operator= (tchecker::parsing::sync_constraint_t const &) = delete;

      /*!
       \brief Move assignment operator (DELETED)
       */
      tchecker::parsing::sync_constraint_t & operator= (tchecker::parsing::sync_constraint_t &&) = delete;

      /*!
       \brief Clone
       \return A clone of this
       */
      sync_constraint_t * clone() const;

      /*!
       \brief Accessor
       \return Process
       */
      inline tchecker::parsing::process_declaration_t const & process() const
      {
        return _process;
      }

      /*!
       \brief Accessor
       \return Event
       */
      inline tchecker::parsing::event_declaration_t const & event() const
      {
        return _event;
      }

      /*!
       \brief Accessor
       \return Strength
       */
      inline enum tchecker::sync_strength_t strength() const
      {
        return _strength;
      }
    private:
      tchecker::parsing::process_declaration_t const & _process;   /*!< Process */
      tchecker::parsing::event_declaration_t const & _event;       /*!< Event */
      enum tchecker::sync_strength_t const _strength;              /*!< Strength */

      friend std::ostream & operator<< (std::ostream & os, tchecker::parsing::sync_constraint_t const & c);
    };

    /*!
     \brief Output operator
     \param os : output stream
     \param c : synchronization constraint
     \return os after c has been output
     */
    std::ostream & operator<< (std::ostream & os, tchecker::parsing::sync_constraint_t const & c);

    /*!
     \class sync_declaration_t
     \brief Declaration of synchronized process events
     */
    class sync_declaration_t final : public tchecker::parsing::inner_declaration_t {
    public:
      /*!
       \brief Constructor
       \param syncs : synchronization constraints
       \pre syncs has size > 0, syncs does not contain nullptr, and there
       is no two events with the same process in syncs
       \throws std::invalid_argument : if syncs size == 0, or if syncs
       contains nullptr, or two events in syncs have the same process
       \note syncs has been moved to this
       */
      sync_declaration_t(std::vector<tchecker::parsing::sync_constraint_t const *> && syncs);

      /*!
       \brief Destructor
       */
      virtual ~sync_declaration_t();

      /*!
       \brief Iterator over attributes
       */
      using const_iterator_t = std::vector<tchecker::parsing::sync_constraint_t const *>::const_iterator;

      /*!
       \brief Accessor
       \return Range [begin;end) of process events
       */
      inline tchecker::range_t<const_iterator_t> sync_constraints() const
      {
        return tchecker::make_range(_syncs.begin(), _syncs.end());
      }
    private:
      /*!
       \brief Clone
       \return A clone of this
       \note The clone shares pointers with this
       */
      virtual tchecker::parsing::declaration_t * do_clone() const;

      /*!
       \brief Visit
       \post this has been visited by v
       */
      virtual void do_visit(tchecker::parsing::declaration_visitor_t & v) const;

      /*!
       \brief Output the declaration
       \param os : output stream
       \post this declaration has been output to os
       \return os after this declaration has been output
       */
      virtual std::ostream & do_output(std::ostream & os) const;

      std::vector<tchecker::parsing::sync_constraint_t const *> const _syncs;  /*!< Synchronization vector */
    };

    /*!
     \class system_declaration_t
     \brief System declaration
     */
    class system_declaration_t final : public tchecker::parsing::declaration_t {
    public:
      /*!
       \brief Constructor
       \param name : system name
       \pre name is not empty
       \throw std::invalid_argument : if name is empty
       */
      system_declaration_t(std::string const & name);

      /*!
       \brief Destructor
       \post
       */
      virtual ~system_declaration_t();

      /*!
       \brief Accessor
       \return Name
       */
      inline std::string const & name() const
      {
        return _name;
      }

      /*!
       \brief Type of iterator on declarations
       */
      using const_iterator_t = std::vector<tchecker::parsing::inner_declaration_t const *>::const_iterator;

      /*!
       \brief Accessor
       \return Range [begin;end) of declarations
       */
      inline tchecker::range_t<const_iterator_t> declarations() const
      {
        return tchecker::make_range(_decl.begin(), _decl.end());
      }

      /*!
       \brief Accessor
       \param name : name of declaration
       \return int declaration with name if any, nullptr otherwise
       */
      inline tchecker::parsing::int_declaration_t const * get_int_declaration(std::string const & name) const
      {
        return get_declaration(name, _ints);
      }

      /*!
       \brief Insert an int declaration
       \param d : int declaration
       \pre d != nullptr (checked by assertion)
       \return true if d has been inserted, false otherwise
       \note this takes owenship on d if inserted
       */
      inline bool insert_int_declaration(tchecker::parsing::int_declaration_t const * d)
      {
        assert(d != nullptr);
        return insert_declaration(d->name(), d, _ints);
      }

      /*!
       \brief Accessor
       \param name : name of declaration
       \return clock declaration with name if any, nullptr otherwise
       */
      inline tchecker::parsing::clock_declaration_t const * get_clock_declaration(std::string const & name) const
      {
        return get_declaration(name, _clocks);
      }

      /*!
       \brief Insert a clock declaration
       \param d : clock declaration
       \pre d != nullptr (checked by assertion)
       \return true if d has been inserted, false otherwise
       \note this takes ownership on d if inserted
       */
      inline bool insert_clock_declaration(tchecker::parsing::clock_declaration_t const * d)
      {
        assert(d != nullptr);
        return insert_declaration(d->name(), d, _clocks);
      }

      /*!
       \brief Accessor
       \param name : name of declaration
       \return process declaration with name if any, nullptr otherwise
       */
      inline tchecker::parsing::process_declaration_t const * get_process_declaration(std::string const & name) const
      {
        return get_declaration(name, _procs);
      }

      /*!
       \brief Insert a process declaration
       \param d : process declaration
       \pre d != nullptr (checked by assertion)
       \return true if d has been inserted, false otherwise
       \note this takes iwnership on d if inserted
       */
      inline bool insert_process_declaration(tchecker::parsing::process_declaration_t const * d)
      {
        assert(d != nullptr);
        return insert_declaration(d->name(), d, _procs);
      }

      /*!
       \brief Accessor
       \param name : name of declaration
       \return event declaration with name if any, nullptr otherwise
       */
      inline tchecker::parsing::event_declaration_t const * get_event_declaration(std::string const & name) const
      {
        return get_declaration(name, _events);
      }

      /*!
       \brief Insert an event declaration
       \param d : event declaration
       \pre d != nullptr (checked by assertion)
       \return true if d has been inserted, false otherwise
       \note this takes ownership on d if inserted
       */
      inline bool insert_event_declaration(tchecker::parsing::event_declaration_t const * d)
      {
        assert(d != nullptr);
        return insert_declaration(d->name(), d, _events);
      }

      /*!
       \brief Accessor
       \param proc : process name
       \param name : location name
       \return location declaration with process proc and name if any, nullptr
       otherwise
       */
      inline tchecker::parsing::location_declaration_t const * get_location_declaration(std::string const & proc,
                                                                                std::string const & name) const
      {
        return get_declaration(location_map_key(proc, name), _locs);
      }

      /*!
       \brief Insert a location declaration
       \param d : location declaration
       \pre d != nullptr (checked by assertion)
       \return true if d has been inserted, false otherwise
       \note this takes ownership on d if inserted
       */
      inline bool insert_location_declaration(tchecker::parsing::location_declaration_t const * d)
      {
        assert(d != nullptr);
        auto key = location_map_key(d->process().name(), d->name());
        return insert_declaration(key, d, _locs);
      }

      /*!
       \brief Insert an edge declaration
       \param d : edge declaration
       \pre d != nullptr (checked by assertion)
       \return true if d has been inserted, false otherwise
       \note this takes ownership on d if inserted
       */
      inline bool insert_edge_declaration(tchecker::parsing::edge_declaration_t const * d)
      {
        assert(d != nullptr);
        _decl.push_back(d);
        return true;
      }

      /*!
       \brief Insert a synchronization declaration
       \param d : synchronization declaration
       \pre d != nullptr (checked by assertion)
       \return true if d has been inserted, false otherwise
       \note this takes ownership on d if inserted
       */
      inline bool insert_sync_declaration(tchecker::parsing::sync_declaration_t const * d)
      {
        assert(d != nullptr);
        _decl.push_back(d);
        return true;
      }
    private:
      /*!
       \brief Clone
       \return A clone of this
       */
      virtual tchecker::parsing::declaration_t * do_clone() const;

      /*!
       \brief Visit
       \post this has been visited by v
       */
      virtual void do_visit(tchecker::parsing::declaration_visitor_t & v) const;

      /*!
       \brief Output the declaration
       \param os : output stream
       \post this declaration has been output to os
       \return os after this declaration has been output
       */
      virtual std::ostream & do_output(std::ostream & os) const;

      /*!
       \brief Type of declarations index
       */
      template <class T>
      using declaration_map_t = std::unordered_map<std::string, T const *>;

      /*!
       \brief Accessor
       \tparam T : type of declaration
       \param name : name of declaration
       \param m : declaration map
       \return declaration of type T with key name in map m, nullptr if no such
       declaration in m
       */
      template <class T>
      T const * get_declaration(std::string const & name, declaration_map_t<T> const & m) const
      {
        auto it = m.find(name);
        if (it == m.end())
          return nullptr;
        return it->second;
      }

      /*!
       \brief Insert a declaration
       \tparam T : type of declaration
       \param name : name of declaration
       \param d : declaration
       \param m : declaration map
       \pre d is not nullptr (checked by assertion), m does not already contain
       a declaration with key name
       \return true if d has been inserted in m with key name, false otherwise
       (m already contains a declaration with key name)
       \post d has been inserted in m with key name and in the list of
       declarations if returned value is true
       */
      template <class T>
      bool insert_declaration(std::string const & name, T const * d, declaration_map_t<T> & m)
      {
        assert( d != nullptr );
        if ( ! m.insert({name, d}).second )
          return false;
        _decl.push_back(d);
        return true;
      }

      /*!
       \brief Accessor
       \param process_name : process name
       \param name : location name
       \return key to location map for location name in process process_name
       */
      inline static std::string location_map_key(std::string const & process_name, std::string const & name)
      {
        return (process_name + ":" + name);
      }

      std::string const _name;                                             /*!< Name */
      std::vector<inner_declaration_t const *> _decl;                      /*!< Declarations */
      declaration_map_t<tchecker::parsing::int_declaration_t> _ints;       /*!< int declarations index */
      declaration_map_t<tchecker::parsing::clock_declaration_t> _clocks;   /*!< clock declaration index */
      declaration_map_t<tchecker::parsing::process_declaration_t> _procs;  /*!< process declaration index */
      declaration_map_t<tchecker::parsing::event_declaration_t> _events;   /*!< event declaration index */
      declaration_map_t<tchecker::parsing::location_declaration_t> _locs;  /*!< location declaration index */
    };

    /*!
     \class declaration_visitor_t
     \brief Visitor for declarations
     */
    class declaration_visitor_t {
    public:
      /*!
       \brief Constructor
       */
      declaration_visitor_t() = default;

      /*!
       \brief Copy constructor
       */
      declaration_visitor_t(tchecker::parsing::declaration_visitor_t const &) = default;

      /*!
       \brief Move constructor
       */
      declaration_visitor_t(tchecker::parsing::declaration_visitor_t &&) = default;

      /*!
       \brief Destructor
       */
      virtual ~declaration_visitor_t() = default;

      /*!
       \brief Assignment operator
       */
      tchecker::parsing::declaration_visitor_t & operator= (tchecker::parsing::declaration_visitor_t const &) = default;

      /*!
       \brief Move assignment operator
       */
      tchecker::parsing::declaration_visitor_t & operator= (tchecker::parsing::declaration_visitor_t &&) = default;

      /*!
       \brief Visitors
       */
      virtual void visit(tchecker::parsing::system_declaration_t const & d) = 0;
      virtual void visit(tchecker::parsing::clock_declaration_t const & d) = 0;
      virtual void visit(tchecker::parsing::int_declaration_t const & d) = 0;
      virtual void visit(tchecker::parsing::process_declaration_t const & d) = 0;
      virtual void visit(tchecker::parsing::event_declaration_t const & d) = 0;
      virtual void visit(tchecker::parsing::location_declaration_t const & d) = 0;
      virtual void visit(tchecker::parsing::edge_declaration_t const & d) = 0;
      virtual void visit(tchecker::parsing::sync_declaration_t const & d) = 0;
    };

  } // end of namespace parsing

} // end of namespace tchecker

#endif // TCHECKER_PARSING_DECLARATION_HH

max_stars_repo_name: pictavien/tchecker

max_issues_repo_name: pictavien/tchecker

max_forks_repo_name: pictavien/tchecker


text:
#include "IncrementalModuleSet.h"
#include "TreeInfixPrinter.h"
#include "../common/Timer.h"
#include "../common/Logger.h"

#include <iostream>
#include <algorithm>
#include <boost/lambda/lambda.hpp>

using namespace boost::lambda;

namespace hydla {
namespace hierarchy {

IncrementalModuleSet::IncrementalModuleSet()
{}

IncrementalModuleSet::IncrementalModuleSet(ModuleSet ms):
  ModuleSetContainer(ms)
{}

IncrementalModuleSet::IncrementalModuleSet(const IncrementalModuleSet& im)
{
  stronger_modules_ = im.stronger_modules_;
  weaker_modules_ = im.weaker_modules_;
  maximal_module_set_ = im.maximal_module_set_;
}

IncrementalModuleSet::~IncrementalModuleSet()
{}

/**
 * 対象のモジュール集合から除去可能なモジュールの集合を返す
 * @param current_ms 対象のモジュール集合
 * @param ms 矛盾の原因となったモジュール集合
 */
std::set<ModuleSet> IncrementalModuleSet::get_removable_module_sets(ModuleSet &current_ms, const ModuleSet &ms)
{
  HYDLA_LOGGER_DEBUG("%% candidate modules : ", ms.get_name(), "\n");
  HYDLA_LOGGER_DEBUG("%% current modules : ", current_ms.get_name(), "\n");
  // msは矛盾集合
  std::set<ModuleSet> removable;
  std::set<ModuleSet> ret;

  for ( auto it : ms )
  {
    // Check if "it" is required module
    if (!stronger_modules_.count(it)) continue;
    // Check if "it" is included in current module set
    if (current_ms.find(it) != current_ms.end()) continue;
    ModuleSet removable_modules;
    // insert "it" to empty set
    removable_modules.add_module(it);
    // insert weaker modules than "it" to "removable_modules"
    removable_modules.insert(weaker_modules_[it]);
    // "removable_modules" = { M | M = "it" \/ M << "it" }
    removable.insert(removable_modules);
  }

  // select minimal module sets
  for ( auto it : removable )
  {
    bool is_minimal = true;
    for ( auto it2 : removable )
    {
      // If it == it2 then continue roop
      if (it.including(it2) && it2.including(it)) continue;
      // Check if it include other module sets in removable
      if (it.including(it2))
      {
        is_minimal = false;
        break;
      }
    }
    // "ret" are minimal module sets of "removable"
    if (is_minimal) ret.insert(it);
  }

  // string for debug
  std::string str = "";
  for (auto it : ret)
  {
    str += it.get_name();
    str += " , ";
  }
  HYDLA_LOGGER_DEBUG("%% removable modules : ", str, "\n");
  return ret;
}

ModuleSet IncrementalModuleSet::unadopted_module_set()
{
  return *ms_to_visit_.begin();
}

void IncrementalModuleSet::add_parallel(IncrementalModuleSet& parallel_module_set)
{
  add_order_data(parallel_module_set);

  // 今まで出現したすべてのモジュールの集合を保持
  maximal_module_set_.insert(parallel_module_set.maximal_module_set_);
}

void IncrementalModuleSet::add_weak(IncrementalModuleSet& weak_module_set_list)
{
  add_order_data(weak_module_set_list);

  // thisに含まれるすべてのモジュールが
  // weak_module_setに含まれるすべてのモジュールよりも強いという情報を保持
  for (auto this_module : maximal_module_set_)
  {
    for (auto weaker_module : weak_module_set_list.maximal_module_set_)
    {
      stronger_modules_[weaker_module].add_module(this_module);
      weaker_modules_[this_module].add_module(weaker_module);
    }
  }

  // 今まで出現したすべてのモジュールの集合を保持
  maximal_module_set_.insert(weak_module_set_list.maximal_module_set_);
}

void IncrementalModuleSet::add_order_data(IncrementalModuleSet& ims)
{
  for (auto sm : ims.stronger_modules_)
  {
    stronger_modules_[sm.first].insert(sm.second);
  }
  for (auto wm : ims.weaker_modules_)
  {
    weaker_modules_[wm.first].insert(wm.second);
  }
}

std::ostream& IncrementalModuleSet::dump_module_sets_for_graphviz(std::ostream& s)
{
  reset();
  s << "digraph module_set { " << std::endl;
  s << "  edge [dir=back];" << std::endl;
  module_set_set_t mss;
  s << "  \"" << get_max_module_set().get_name() << "\"" << std::endl;
  while(has_next())
  {
    ModuleSet tmp = get_max_module_set();
    tmp.erase(unadopted_module_set());
    generate_new_ms(mss, tmp);
    for (auto ms : ms_to_visit_)
    {
      ModuleSet child = get_max_module_set();
      child.erase(ms);
      if (tmp.including(child))
      {
        s << "  \"" << tmp.get_name() << "\" -> \"" << child.get_name() << "\"" << std::endl;
      }
    }
  }
  s << "}" << std::endl;
  reset();
  return s;
}

std::ostream& IncrementalModuleSet::dump_priority_data_for_graphviz(std::ostream& s) const
{
  s << "digraph priority_data { " << std::endl;
  s << "  edge [dir=back];" << std::endl;
  for (auto m : weaker_modules_)
  {
    for (auto wm : m.second)
    {
      s << "  \"" << m.first.first << "\" -> \"" << wm.first << "\";" << std::endl;
    }
  }
  s << "}" << std::endl;
  return s;
}
std::ostream& IncrementalModuleSet::dump(std::ostream& s) const
{
  return s;
}

std::ostream& IncrementalModuleSet::dump_node_names(std::ostream& s) const
{
  return s;
}

std::ostream& IncrementalModuleSet::dump_node_trees(std::ostream& s) const
{
  return s;
}

/**
 * この mark_nodes は調べているモジュール集合が
 * 無矛盾だった場合に呼ばれ、
 * そのモジュール集合が包含しているモジュール集合を
 * 探索対象から外す
 */

bool IncrementalModuleSet::check_same_ms_generated(module_set_set_t &new_mss, ModuleSet &ms)
{
  return new_mss.count(ms);
}

void IncrementalModuleSet::update_by_new_mss(module_set_set_t &new_mss)
{
  ms_to_visit_.clear();
  for (auto new_ms : new_mss)
  {
    ms_to_visit_.insert(new_ms);
  }
}

void IncrementalModuleSet::init()
{
}

void IncrementalModuleSet::remove_included_ms_by_current_ms()
{
  /// current は現在のモジュール集合
  ModuleSet current = unadopted_module_set();
  module_set_set_t::iterator lit = ms_to_visit_.begin();
  while(lit != ms_to_visit_.end())
  {
    /**
     * ms_to_visit_内のモジュール集合で
     * currentが包含するモジュール集合を削除
     */
    if (lit->including(current))
    {
      lit = ms_to_visit_.erase(lit);
    }
    else lit++;
  }
}

/**
 * この mark_nodes は調べているモジュール集合が
 * 矛盾した場合に呼ばれ、新たなモジュール集合を生成する
 * @param mms そのフェーズにおおける極大無矛盾集合
 * @param ms 矛盾の原因となったモジュール集合
 */
void IncrementalModuleSet::generate_new_ms(const module_set_set_t& mcss, const ModuleSet& ms)
{
  HYDLA_LOGGER_DEBUG("%% inconsistent module set : ", ms.get_name());
  std::string for_debug = "";
  for (auto mit : mcss)
  {
    for_debug += mit.get_name() + " , ";
  }
  HYDLA_LOGGER_DEBUG("%% maximal consistent module set : ", for_debug);
  // 制約の優先度を崩さない範囲で除去可能なモジュールを要素として持つ集合を得る
  // 結果用のモジュール集合の集合
  module_set_set_t new_mss;
  ModuleSet inconsistent_ms = ms;
  for (auto lit : ms_to_visit_)
  {
    // 探索対象のモジュール集合がmsを含んでいる場合
    // そのモジュール集合は矛盾するため
    // 新たなモジュール集合を生成する
    if (lit.disjoint(inconsistent_ms))
    {
      // get vector of removable module set
      std::set<ModuleSet> rm = get_removable_module_sets(lit,inconsistent_ms);
      for (auto removable_module_set : rm)
      {
        // remove removable module set of lit.
        // make new module set
        // The set has no module which is in removable module set
        ModuleSet new_ms = lit;
        new_ms.insert(removable_module_set);
        // check weather the new_ms is included by maximal consistent module sets.
        bool checked = false;
        for (auto mcs : mcss)
        {
         // 生成されたモジュール集合が極大無矛盾集合に包含されている場合checkedをtrueにしておく
          if (new_ms.including(mcs))
          {
            checked = true;
            break;
          }
        }
       // checkedがtrueでない場合、生成したモジュール集合を探索対象に追加
        if (!checked)
        {
          new_mss.insert(new_ms);
          HYDLA_LOGGER_DEBUG("%% new ms : ", new_ms.get_name());
        }
      }
    }else{
      // 探索対象がmsを含んでいない場合そのまま残しておく
      if (!check_same_ms_generated(new_mss, lit))
      {
        new_mss.insert(lit);
      }
    }
  }
  // update ms_to_visit_ by generated module sets
  update_by_new_mss(new_mss);
}

/// 最も要素数の多いモジュール集合を返す
  ModuleSet IncrementalModuleSet::get_max_module_set() const{
    ModuleSet ret = maximal_module_set_;
    return ret;
  }

// 探索対象を引数のmssが示す状態にする
  void IncrementalModuleSet::reset(const module_set_set_t &mss)
  {
    ms_to_visit_ = mss;
  }

// 探索対象を初期状態に戻す
  void IncrementalModuleSet::reset()
  {
    ms_to_visit_.clear();
    ms_to_visit_.insert(ModuleSet());
  }

} // namespace hierarchy
} // namespace hydla

ext: cpp

max_stars_repo_name: takafumihoriuchi/HyLaGI

max_issues_repo_name: takafumihoriuchi/HyLaGI

max_forks_repo_name: takafumihoriuchi/HyLaGI