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_block.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>
 
PRISMS_PF_BEGIN_NAMESPACE
 
template <unsigned int dim, unsigned int degree, typename number>
class SolverBase
{
public:
  SolverBase(SolveBlock                               _solve_block,
             const SolveContext<dim, degree, number> &_solve_context)
    : solve_block(std::move(_solve_block))
    , solve_context(&_solve_context)
    , solutions(solve_block,
                solve_context->get_field_attributes(),
                solve_context->get_matrix_free_manager())
  {}

  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<SolveBlock> &all_solve_blocks)
  {
    NewDependencyExtents extents(solve_block.field_indices, all_solve_blocks);
    solutions.init(extents);
 
    // Apply constraints.
    solutions.apply_constraints();
 
    if (solutions.num_levels() > 0)
      {
        solutions.reinit_mg_transfer(solve_context->get_dof_manager());
      }
  }

  virtual void
  reinit()
  {
    // Apply constraints.
    solutions.reinit();
    solutions.apply_constraints();
    if (solutions.num_levels() > 0)
      {
        solutions.reinit_mg_transfer(solve_context->get_dof_manager());
        solutions.mg_transfer_down(
          solve_context->get_dof_manager(),
          solve_context->get_user_inputs().spatial_discretization.global_refinement,
          true);
      }
  }

  virtual void
  solve_impl()
  {}

  virtual void
  solve()
  {
    if (solve_context->get_simulation_timer().get_increment() == 0 &&
        solve_block.solve_timing == SolveTiming::Primary)
      {
        // Set the initial condition
        set_initial_condition();
      }
    else
      {
        this->solve_impl();
      }
    if (solve_context->get_simulation_timer().get_increment() == 0)
      {
        solutions.apply_initial_condition_for_old_fields();
      }
    if (solutions.num_levels() > 0)
      {
        solutions.mg_transfer_down(
          solve_context->get_dof_manager(),
          solve_context->get_user_inputs().spatial_discretization.global_refinement,
          false);
      }
  }

  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_block.field_indices)
      {
        bool initialized_from_file = false;
 
        // First, try to find this variable in IC files
        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())
              {
                // Found in file - read from file
                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));
 
                initialized_from_file = true;
                break; // Stop searching once found
              }
          }
        // If not found in files, use programmatic IC
        if (!initialized_from_file)
          {
            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 SolveBlock &
  get_solve_block() const
  {
    return solve_block;
  }
 
protected:
  SolveBlock solve_block;

  const SolveContext<dim, degree, number> *solve_context;

  GroupSolutionHandler<dim, number> solutions;
 
  std::vector<SolverBase<dim, degree, number> *> aux_solvers;
};
 
PRISMS_PF_END_NAMESPACE