1#ifndef DUNE_COPASI_MODEL_MAKE_DIFFUSION_REACTION_STEP_OPERATOR_HH
2#define DUNE_COPASI_MODEL_MAKE_DIFFUSION_REACTION_STEP_OPERATOR_HH
7#if HAVE_METIS && DUNE_COPASI_CONCURRENT_ASSEMBLY
8#include <dune/pdelab/common/partition/metis.hh>
11#include <dune/pdelab/common/partition/simple.hh>
12#include <dune/pdelab/common/execution.hh>
13#include <dune/pdelab/operator/operator.hh>
15#include <spdlog/spdlog.h>
21template<
class LocalBasisTraits,
24 class ResidualQuantity,
26 PDELab::Concept::Basis Basis,
27 Dune::Concept::GridView CellDataGridView =
typename Basis::EntitySet,
28 class CellDataType =
double>
29[[nodiscard]]
inline static std::unique_ptr<PDELab::Operator<Coefficients, Coefficients>>
30make_step_operator(
const ParameterTree& config,
33 std::shared_ptr<
const FunctorFactory<Basis::EntitySet::dimension>> functor_factory,
34 std::shared_ptr<
const CellData<CellDataGridView, CellDataType>> grid_cell_data =
nullptr)
37 std::unique_ptr<PDELab::Operator<Coefficients, Coefficients>> one_step;
38 const auto& assembly_cfg = config.sub(
"assembly");
40 auto make_one_step_op = [&, functor_factory, grid_cell_data]<
class ExecutionPolicy, PDELab::Concept::Basis OperatorBasis>(
41 ExecutionPolicy execution_policy, const OperatorBasis& operator_basis) {
42 spdlog::info(
"Creating mass/stiffness local operator");
43 const auto& time_step_cfg = config.sub(
"time_step_operator");
45 LocalOperator<OperatorBasis, LocalBasisTraits, CellDataGridView, CellDataType, ExecutionPolicy>;
46 LocalOperator
const stiff_lop(operator_basis,
52 LocalOperator
const mass_lop(operator_basis,
58 return Dune::Copasi::make_step_operator<Coefficients, Residual, ResidualQuantity, TimeQuantity>(
59 time_step_cfg, operator_basis, mass_lop, stiff_lop);
62 auto entity_set = basis.entitySet();
63 std::size_t part_patches =
64 assembly_cfg.get(
"partition.patches",
static_cast<std::size_t
>(entity_set.size(0) / 40));
66#if DUNE_COPASI_CONCURRENT_ASSEMBLY
67 auto make_concurrent_one_step_op =
68 [&]<
class ExecutionPolicy>(const ExecutionPolicy execution_policy) {
69 auto entity_set = basis.entitySet();
70 auto default_part_type =
76 auto part_type = assembly_cfg.get(
"partition.type", default_part_type);
77 auto part_coloring = assembly_cfg.get(
"partition.coloring",
"None");
78 spdlog::info(
" Concurrent Partitioning:");
79 spdlog::info(
" Type: {}", part_type);
80 spdlog::info(
" Patches: {}", part_patches);
81 spdlog::info(
" Halo: {}", halo);
82 spdlog::info(
" Coloring: {}", part_coloring);
83 if (part_type ==
"METIS") {
85 if (part_coloring ==
"None") {
86 return make_one_step_op(
89 PDELab::EntitySetPartitioner::Metis{ entity_set, part_patches, halo },
90 TypeTree::treePath()));
91 }
else if (part_coloring ==
"DSatur") {
92 return make_one_step_op(
94 basis.subSpace(PDELab::EntitySetPartitioner::MetisColored{ entity_set, part_patches, halo },
95 TypeTree::treePath()));
97 throw format_exception(IOError{},
"Not known coloring algorithm '{}' known", part_coloring);
100 throw format_exception(IOError{},
"'METIS' partitioner is not available on this executable");
102 }
else if (part_type ==
"Simple") {
103 if (part_coloring ==
"None") {
104 return make_one_step_op(
107 PDELab::EntitySetPartitioner::Simple{ entity_set, part_patches, halo },
108 TypeTree::treePath()));
109 }
else if (part_coloring ==
"DSatur") {
110 return make_one_step_op(
112 basis.subSpace(PDELab::EntitySetPartitioner::SimpleColored{ entity_set, part_patches, halo },
113 TypeTree::treePath()));
115 throw format_exception(IOError{},
"Not known coloring algorithm '{}' known", part_coloring);
118 throw format_exception(IOError{},
"Not known partition algorithm '{}' known", part_type);
122 auto exec_policy = assembly_cfg.get(
"type",
"concurrent");
124 auto exec_policy = assembly_cfg.get(
"type",
"sequential");
127 spdlog::info(
"Creating time-step operator with '{}' execution policy", exec_policy);
128 if (exec_policy !=
"sequential" and part_patches < std::thread::hardware_concurrency()) {
129 exec_policy =
"sequential";
130 spdlog::warn(
"Too few patches '{}', execution policy of assembler is switched to 'sequential'", part_patches);
133 if (exec_policy ==
"sequential") {
134 one_step = make_one_step_op(PDELab::Execution::seq, basis);
135#if DUNE_COPASI_CONCURRENT_ASSEMBLY
136 }
else if (exec_policy ==
"concurrent") {
137 one_step = make_concurrent_one_step_op(PDELab::Execution::par);
140 throw format_exception(IOError{},
"Not known execution '{}' policy known", exec_policy);
Definition: factory.hh:28
auto format_exception(Exception &&e, fmt::format_string< Args... > format, Args &&... args)
Definition: exceptions.hh:23