group_solution_handler.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/distributed/solution_transfer.h>
#include <deal.II/lac/affine_constraints.h>
#include <deal.II/lac/la_parallel_block_vector.h>
#include <deal.II/lac/la_parallel_vector.h>
#include <deal.II/matrix_free/matrix_free.h>
#include <deal.II/multigrid/mg_transfer_global_coarsening.h>
 
#include <prismspf/core/conditional_ostreams.h>
#include <prismspf/core/constraint_manager.h>
#include <prismspf/core/dependency_extents.h>
#include <prismspf/core/dof_manager.h>
#include <prismspf/core/field_attributes.h>
#include <prismspf/core/matrix_free_manager.h>
#include <prismspf/core/solve_block.h>
#include <prismspf/core/timer.h>
#include <prismspf/core/type_enums.h>
#include <prismspf/core/types.h>
 
#include <prismspf/config.h>
 
#include <vector>
 
PRISMS_PF_BEGIN_NAMESPACE

template <unsigned int dim, typename number>
struct SolutionLevel
{
  BlockVector<number>              solutions;
  std::vector<BlockVector<number>> old_solutions;
};

template <unsigned int dim, typename number>
class GroupSolutionHandler
{
public:
#if DEAL_II_VERSION_MAJOR >= 9 && DEAL_II_VERSION_MINOR >= 7
  using SolutionTransfer = dealii::SolutionTransfer<dim, SolutionVector<number>>;
#else
  using SolutionTransfer =
    dealii::parallel::distributed::SolutionTransfer<dim, SolutionVector<number>>;
#endif

  GroupSolutionHandler(SolveBlock                            _solve_block,
                       const std::vector<FieldAttributes>   &_attributes_list,
                       const MatrixFreeManager<dim, number> &_matrix_free_manager);

  [[nodiscard]] BlockVector<number> &
  get_solution_full_vector();

  [[nodiscard]] const BlockVector<number> &
  get_solution_full_vector() const;

  [[nodiscard]] SolutionVector<number> &
  get_solution_vector(unsigned int global_index);

  [[nodiscard]] const SolutionVector<number> &
  get_solution_vector(unsigned int global_index) const;

  [[nodiscard]] BlockVector<number> &
  get_old_solution_full_vector(unsigned int age);

  [[nodiscard]] const BlockVector<number> &
  get_old_solution_full_vector(unsigned int age) const;

  [[nodiscard]] SolutionVector<number> &
  get_old_solution_vector(unsigned int age, unsigned int global_index);

  [[nodiscard]] const SolutionVector<number> &
  get_old_solution_vector(unsigned int age, unsigned int global_index) const;

  [[nodiscard]] SolutionLevel<dim, number> &
  get_primary_solutions();

  [[nodiscard]] const SolutionLevel<dim, number> &
  get_primary_solutions() const;

  [[nodiscard]] SolutionLevel<dim, number> &
  get_solution_level(unsigned int relative_level = -1);

  [[nodiscard]] const SolutionLevel<dim, number> &
  get_solution_level(unsigned int relative_level = -1) const;

  const std::vector<MatrixFree<dim, number>> &
  get_matrix_free_levels() const;

  [[nodiscard]] unsigned int
  get_block_index(unsigned int global_index) const;

  [[nodiscard]] const SolveBlock &
  get_solve_block() const;

  [[nodiscard]] const std::vector<unsigned int> &
  get_global_to_block_index() const;

  [[nodiscard]] const std::vector<unsigned int> &
  get_block_to_global_index() const;

  void
  init(const NewDependencyExtents &extents);

  void
  reinit();

  void
  update_ghosts() const;

  void
  zero_out_ghosts() const;

  void
  apply_constraints();

  void
  apply_constraints(BlockVector<number> &solution_vectors);

  void
  apply_constraints_to_all();

  void
  apply_initial_condition_for_old_fields();

  void
  update();

  void
  init_solution_transfer();

  void
  prepare_for_solution_transfer();

  void
  execute_solution_transfer();

  template <unsigned int degree>
  void
  reinit_mg_transfer(const DoFManager<dim, degree> &dof_manager);

  template <unsigned int degree>
  void
  mg_transfer_down(const DoFManager<dim, degree> &dof_manager,
                   unsigned int                   finest_level,
                   bool                           transfer_old_solutions = false);

  [[nodiscard]] unsigned int
  num_levels();

  void
  print_solution_full_vector(std::ostream &out) const;
 
private:
  SolveBlock solve_block;

  std::vector<unsigned int> block_to_global_index;

  std::vector<unsigned int> global_to_block_index;

  SolutionLevel<dim, number> primary_solutions;

  std::vector<SolutionLevel<dim, number>> solution_levels;

  const MatrixFreeManager<dim, number> *matrix_free_manager = nullptr;

  std::vector<SolutionTransfer> block_solution_transfer;

  std::vector<dealii::MGConstrainedDoFs> mg_constraints;

  std::vector<dealii::MGTransferMatrixFree<dim, number>> mg_transfer;
};
 
template <unsigned int dim, typename number>
template <unsigned int degree>
inline void
GroupSolutionHandler<dim, number>::reinit_mg_transfer(
  const DoFManager<dim, degree> &dof_manager)
{
  // 1. Initialize constraints
  const unsigned int num_blocks = solve_block.field_indices.size();
  mg_constraints                = std::vector<dealii::MGConstrainedDoFs>(num_blocks);
  mg_transfer = std::vector<dealii::MGTransferMatrixFree<dim, number>>(num_blocks);
  for (unsigned int block_index = 0; block_index < num_blocks; block_index++)
    {
      unsigned int field_index = block_to_global_index[block_index];
      mg_constraints[block_index].initialize(
        dof_manager.get_field_dof_handler(field_index));
      // 2. Initialize MG Transfer
      mg_transfer[block_index].initialize_constraints(mg_constraints[block_index]);
      mg_transfer[block_index].build(dof_manager.get_field_dof_handler(field_index));
    }
}
 
template <unsigned int dim, typename number>
template <unsigned int degree>
inline void
GroupSolutionHandler<dim, number>::mg_transfer_down(
  const DoFManager<dim, degree> &dof_manager,
  unsigned int                   finest_level,
  bool                           transfer_old_solutions)
{
  Timer::start_section("MG Transfer LHS Dependencies");
  const unsigned int num_blocks = solve_block.field_indices.size();
  for (unsigned int block_index = 0; block_index < num_blocks; block_index++)
    {
      unsigned int field_index = block_to_global_index[block_index];
      dealii::MGLevelObject<SolutionVector<number>> temp_mg_solutions(
        1 + finest_level - solution_levels.size(),
        finest_level);
 
      const auto &dof_handler = dof_manager.get_field_dof_handler(field_index);
 
      // transfer regular solutions to mg levels
      mg_transfer[block_index].interpolate_to_mg(dof_handler,
                                                 temp_mg_solutions,
                                                 primary_solutions.solutions.block(
                                                   block_index));
      // swap to actual mg solution vectors
      for (unsigned int relative_level = 0; relative_level < solution_levels.size();
           ++relative_level)
        {
          unsigned int level = finest_level - relative_level;
          solution_levels[relative_level]
            .solutions.block(block_index)
            .swap(temp_mg_solutions[level]);
          Timer::start_section("MG LHS Dependencies Update Ghosts");
          solution_levels[relative_level].solutions.update_ghost_values();
          Timer::end_section("MG LHS Dependencies Update Ghosts");
        }
 
      if (!transfer_old_solutions)
        {
          Timer::end_section("MG Transfer LHS Dependencies");
          return;
        }
      // Transfer old solutions
      for (unsigned int age_index = 0;
           age_index < solution_levels[0].old_solutions.size();
           ++age_index)
        {
          mg_transfer[block_index].interpolate_to_mg(
            dof_handler,
            temp_mg_solutions,
            solution_levels[0].old_solutions[age_index].block(block_index));
          for (unsigned int relative_level = 0; relative_level < solution_levels.size();
               ++relative_level)
            {
              unsigned int level = finest_level - relative_level;
              if (age_index < solution_levels[relative_level].old_solutions.size())
                {
                  solution_levels[relative_level]
                    .old_solutions[age_index]
                    .block(block_index)
                    .swap(temp_mg_solutions[level]);
                  Timer::start_section("MG LHS Dependencies Update Ghosts");
                  solution_levels[relative_level]
                    .old_solutions[age_index]
                    .update_ghost_values();
                  Timer::end_section("MG LHS Dependencies Update Ghosts");
                }
            }
        }
    }
  Timer::end_section("MG Transfer LHS Dependencies");
}
 
PRISMS_PF_END_NAMESPACE