combined_system.hpp

Go to the documentation of this file.

// Copyright (c) Lawrence Livermore National Security, LLC and
// other Smith Project Developers. See the top-level LICENSE file for
// details.
//
// SPDX-License-Identifier: (BSD-3-Clause)
 
#pragma once
 
#include <vector>
#include <memory>
#include <tuple>
#include "smith/differentiable_numerics/system_base.hpp"
#include "smith/differentiable_numerics/field_store.hpp"
 
namespace smith {
 
struct CombinedSystem : public SystemBase {
  std::vector<std::shared_ptr<SystemBase>> subsystems;  
 
  using SystemBase::SystemBase;
 
  template <typename SubSystemType, typename MaterialType>
  void setMaterialOn(std::shared_ptr<SubSystemType> sub, const MaterialType& material, const std::string& domain_name)
  {
    sub->setMaterial(material, domain_name);
  }
};
 
template <typename... SubSystems>
auto combineSystems(std::shared_ptr<SubSystems>... subs)
{
  static_assert(sizeof...(subs) >= 2, "combineSystems requires at least two sub-systems");
 
  auto first_sub = std::get<0>(std::forward_as_tuple(subs...));
  auto field_store = first_sub->field_store;
 
  auto combined = std::make_shared<CombinedSystem>();
  combined->field_store = field_store;
 
  int max_stagger_iters = 1;
  bool exact_staggered_steps = false;
 
  std::vector<std::shared_ptr<SystemBase>> post_solve_systems;
  std::vector<std::vector<size_t>> subsystem_global_block_indices;
 
  (
      [&](auto& sub) {
        std::vector<size_t> global_block_indices;
        combined->subsystems.push_back(sub);
        for (auto& wf : sub->weak_forms) {
          global_block_indices.push_back(combined->weak_forms.size());
          combined->weak_forms.push_back(wf);
        }
        subsystem_global_block_indices.push_back(global_block_indices);
        if (sub->solver) {
          max_stagger_iters = std::max(max_stagger_iters, sub->solver->maxStaggeredIterations());
          exact_staggered_steps = exact_staggered_steps || sub->solver->exactStaggeredSteps();
        }
        for (auto& cz_sys : sub->cycle_zero_systems) {
          combined->cycle_zero_systems.push_back(cz_sys);
        }
        post_solve_systems.insert(post_solve_systems.end(), sub->post_solve_systems.begin(),
                                  sub->post_solve_systems.end());
      }(subs),
      ...);
 
  combined->solver = std::make_shared<SystemSolver>(max_stagger_iters, exact_staggered_steps);
  for (size_t i = 0; i < combined->subsystems.size(); ++i) {
    const auto& sub = combined->subsystems[i];
    SLIC_ERROR_IF(!sub->solver, "Combined subsystem must have a solver");
    combined->solver->appendStagesWithBlockMapping(*sub->solver, subsystem_global_block_indices[i]);
  }
 
  combined->post_solve_systems = post_solve_systems;
 
  return combined;
}
 
template <typename... SubSystems>
std::shared_ptr<SystemBase> combineSystems(std::shared_ptr<SystemSolver> solver, std::shared_ptr<SubSystems>... subs)
{
  static_assert(sizeof...(subs) >= 2, "combineSystems requires at least two sub-systems");
 
  auto field_store = std::get<0>(std::forward_as_tuple(subs...))->field_store;
 
  std::vector<std::shared_ptr<WeakForm>> wfs;
  std::vector<std::shared_ptr<SystemBase>> cycle_zero_systems;
  std::vector<std::shared_ptr<SystemBase>> post_solve_systems;
 
  (
      [&](auto& sub) {
        for (auto& wf : sub->weak_forms) {
          wfs.push_back(wf);
        }
        for (auto& cz_sys : sub->cycle_zero_systems) {
          cycle_zero_systems.push_back(cz_sys);
        }
        post_solve_systems.insert(post_solve_systems.end(), sub->post_solve_systems.begin(),
                                  sub->post_solve_systems.end());
      }(subs),
      ...);
 
  auto combined = std::make_shared<SystemBase>(field_store, solver, wfs);
  combined->cycle_zero_systems = cycle_zero_systems;
  combined->post_solve_systems = post_solve_systems;
 
  return combined;
}
 
}  // namespace smith