solver_base.h

Go to the documentation of this file.

// SPDX-FileCopyrightText: © 2025 PRISMS Center at the University of Michigan
// SPDX-License-Identifier: GNU Lesser General Public Version 2.1
 
#pragma once
 
#include <deal.II/base/exceptions.h>
#include <deal.II/numerics/vector_tools.h>
 
#include <boost/geometry/core/cs.hpp>
 
#include <prismspf/core/dependency_extents.h>
#include <prismspf/core/group_solution_handler.h>
#include <prismspf/core/initial_conditions.h>
#include <prismspf/core/pde_operator_base.h>
#include <prismspf/core/solve_group.h>
#include <prismspf/core/system_wide.h>
#include <prismspf/core/type_enums.h>
#include <prismspf/core/types.h>
 
#include <prismspf/solvers/solve_context.h>
 
#include <prismspf/user_inputs/user_input_parameters.h>
 
#include <prismspf/config.h>
 
#include <memory>
 
PRISMS_PF_BEGIN_NAMESPACE
 
template <unsigned int dim, unsigned int degree, typename number>
class SolverBase
{
public:
  SolverBase(SolveGroup                               _solve_group,
             const SolveContext<dim, degree, number> &_solve_context)
    : solve_group(std::move(_solve_group))
    , solve_context(&_solve_context)
    , solutions(solve_group, solve_context->get_field_attributes())
  {}
 
  virtual ~SolverBase() = default;
 
  SolverBase(const SolverBase &solver_base) = delete;
 
  SolverBase &
  operator=(const SolverBase &solver_base) = delete;
 
  SolverBase(SolverBase &&solver_base) noexcept = delete;
 
  SolverBase &
  operator=(SolverBase &&solver_base) noexcept = delete;
 
  virtual void
  init(const std::list<DependencyMap> &all_dependeny_sets)
  {
    DependencyExtents extents(solve_group.field_indices, all_dependeny_sets);
    solutions.init(solve_context->get_dof_manager(),
                   solve_context->get_constraint_manager(),
                   extents.oldest_age);
 
    // Apply constraints.
    solutions.apply_constraints();
  }
 
  virtual void
  reinit()
  {
    // Apply constraints.
    solutions.reinit(solve_context->get_dof_manager(),
                     solve_context->get_constraint_manager());
    solutions.apply_constraints();
  }
 
  virtual void
  solve_level([[maybe_unused]] unsigned int relative_level = 0)
  {}
 
  virtual void
  solve()
  {
    if (solve_context->get_simulation_timer().get_increment() == 0 &&
        solve_group.solve_timing == SolveTiming::Primary)
      {
        // Set the initial condition
        set_initial_condition();
        return;
      }
    this->solve_level(0);
    if (solve_context->get_simulation_timer().get_increment() == 0)
      {
        solutions.apply_initial_condition_for_old_fields();
      }
  }
 
  virtual void
  update()
  {
    solutions.update();
  }
 
  virtual void
  update_ghosts()
  {
    solutions.update_ghosts();
  }
 
  void
  prepare_for_solution_transfer()
  {
    solutions.prepare_for_solution_transfer();
  }
 
  void
  execute_solution_transfer()
  {
    solutions.execute_solution_transfer();
  }
 
  virtual void
  print()
  {}
 
  void
  set_initial_condition()
  {
    for (const auto &global_index : solve_group.field_indices)
      {
        if (solve_context->get_user_inputs()
              .load_ic_parameters.get_read_initial_conditions_from_file())
          {
            const auto &initial_condition_parameters =
              solve_context->get_user_inputs().load_ic_parameters;
            for (const auto &initial_condition_file :
                 initial_condition_parameters.get_initial_condition_files())
              {
                auto name_it =
                  std::find(initial_condition_file.simulation_variable_names.begin(),
                            initial_condition_file.simulation_variable_names.end(),
                            solve_context->get_field_attributes()[global_index].name);
                if (name_it != initial_condition_file.simulation_variable_names.end())
                  {
                    solutions.get_solution_vector(global_index).zero_out_ghost_values();
                    dealii::VectorTools::interpolate(
                      SystemWide<dim, degree>::mapping,
                      solve_context->get_dof_manager().get_field_dof_handler(
                        global_index),
                      ReadInitialCondition<dim, number>(
                        *name_it,
                        solve_context->get_field_attributes()[global_index].field_type,
                        initial_condition_file,
                        solve_context->get_user_inputs().spatial_discretization),
                      solutions.get_solution_vector(global_index));
                  }
              }
          }
        else
          {
            solutions.get_solution_vector(global_index).zero_out_ghost_values();
            dealii::VectorTools::interpolate(
              SystemWide<dim, degree>::mapping,
              solve_context->get_dof_manager().get_field_dof_handler(global_index),
              InitialCondition<dim, degree, number>(
                global_index,
                solve_context->get_field_attributes()[global_index].field_type,
                solve_context->get_pde_operator()),
              solutions.get_solution_vector(global_index));
          }
        solutions.apply_initial_condition_for_old_fields();
      }
  }
 
  [[nodiscard]] const GroupSolutionHandler<dim, number> &
  get_solution_manager() const
  {
    return solutions;
  }
 
  [[nodiscard]] GroupSolutionHandler<dim, number> &
  get_solution_manager()
  {
    return solutions;
  }
 
  [[nodiscard]] const SolveGroup &
  get_solve_group() const
  {
    return solve_group;
  }
 
protected:
  SolveGroup solve_group;
 
  const SolveContext<dim, degree, number> *solve_context;
 
  GroupSolutionHandler<dim, number> solutions;
 
  std::vector<SolverBase<dim, degree, number> *> aux_solvers;
};
 
PRISMS_PF_END_NAMESPACE