16 #include "DataStructures/DataBox/PrefixHelpers.hpp"
18 #include "DataStructures/DataBox/Tag.hpp"
19 #include "DataStructures/DataVector.hpp"
25 #include "Elliptic/DiscontinuousGalerkin/DgElementArray.hpp"
26 #include "Helpers/ParallelAlgorithms/LinearSolver/LinearSolverAlgorithmTestHelpers.hpp"
27 #include "IO/Observer/Actions/RegisterWithObservers.hpp"
28 #include "IO/Observer/ObserverComponent.hpp"
29 #include "IO/Observer/Tags.hpp"
30 #include "NumericalAlgorithms/Convergence/HasConverged.hpp"
31 #include "Parallel/Actions/SetupDataBox.hpp"
32 #include "Parallel/Actions/TerminatePhase.hpp"
35 #include "Parallel/InitializationFunctions.hpp"
36 #include "Parallel/Invoke.hpp"
38 #include "Parallel/ParallelComponentHelpers.hpp"
39 #include "Parallel/PhaseDependentActionList.hpp"
40 #include "Parallel/Reduction.hpp"
41 #include "ParallelAlgorithms/DiscontinuousGalerkin/InitializeDomain.hpp"
42 #include "ParallelAlgorithms/Initialization/Actions/RemoveOptionsAndTerminatePhase.hpp"
43 #include "ParallelAlgorithms/Initialization/MutateAssign.hpp"
44 #include "ParallelAlgorithms/LinearSolver/Actions/MakeIdentityIfSkipped.hpp"
53 namespace helpers = LinearSolverAlgorithmTestHelpers;
58 namespace OptionTags {
65 static constexpr
Options::String help =
"The linear operator A to invert.";
71 static constexpr
Options::String help =
"The source b in the equation Ax=b.";
75 static constexpr
Options::String help =
"The solution x in the equation Ax=b";
83 using option_tags = tmpl::list<OptionTags::LinearOperator>;
85 static constexpr
bool pass_metavariables =
false;
89 linear_operator) noexcept {
90 return linear_operator;
97 using option_tags = tmpl::list<OptionTags::Source>;
99 static constexpr
bool pass_metavariables =
false;
108 using option_tags = tmpl::list<OptionTags::ExpectedResult>;
110 static constexpr
bool pass_metavariables =
false;
113 return expected_result;
120 static std::string name() noexcept {
return "ScalarField"; }
125 using operator_applied_to_fields_tag =
132 size_t get_index(
const ElementId<1>& element_id) noexcept {
133 return element_id.segment_ids()[0].index();
144 template <
typename OperandTag>
147 template <
typename OperandTag>
149 using const_global_cache_tags = tmpl::list<LinearOperator>;
151 template <
typename DbTagsList,
typename... InboxTags,
typename Metavariables,
152 typename ActionList,
typename ParallelComponent>
154 db::DataBox<DbTagsList>& box,
162 const ParallelComponent*
const ) noexcept {
163 const size_t element_index = get_index(element_id);
164 const auto& operator_matrices = get<LinearOperator>(box);
165 const auto number_of_elements = operator_matrices.size();
166 const auto& linear_operator =
gsl::at(operator_matrices, element_index);
167 const auto number_of_grid_points = linear_operator.columns();
168 const auto& operand = get<OperandTag>(box);
170 typename OperandTag::type operator_applied_to_operand{
171 number_of_grid_points * number_of_elements};
172 dgemv_(
'N', linear_operator.rows(), linear_operator.columns(), 1,
173 linear_operator.data(), linear_operator.spacing(), operand.data(), 1,
174 0, operator_applied_to_operand.data(), 1);
176 Parallel::contribute_to_reduction<CollectOperatorAction<OperandTag>>(
177 Parallel::ReductionData<
179 operator_applied_to_operand},
180 Parallel::get_parallel_component<ParallelComponent>(
cache)[element_id],
181 Parallel::get_parallel_component<ParallelComponent>(
cache));
185 return {std::move(box),
true};
189 template <
typename OperandTag>
191 using local_operator_applied_to_operand_tag =
194 template <
typename ParallelComponent,
typename DbTagsList,
195 typename Metavariables,
typename ScalarFieldOperandTag,
197 DbTagsList, local_operator_applied_to_operand_tag>> =
nullptr>
198 static void apply(db::DataBox<DbTagsList>& box,
201 const Variables<tmpl::list<ScalarFieldOperandTag>>&
202 Ap_global_data) noexcept {
203 const size_t element_index = get_index(element_id);
206 const auto number_of_grid_points = get<LinearOperator>(box)[0].columns();
207 const auto& Ap_global = get<ScalarFieldOperandTag>(Ap_global_data).get();
209 std::copy(Ap_global.begin() +
210 static_cast<int>(element_index * number_of_grid_points),
212 static_cast<int>((element_index + 1) * number_of_grid_points),
214 db::mutate<local_operator_applied_to_operand_tag>(
216 [&Ap_local, &number_of_grid_points](
auto Ap) noexcept {
217 *Ap =
typename local_operator_applied_to_operand_tag::type{
218 number_of_grid_points};
223 Parallel::get_parallel_component<ParallelComponent>(
cache)[element_id]
224 .perform_algorithm(
true);
229 template <
typename OptionsGroup>
231 using const_global_cache_tags =
232 tmpl::list<ExpectedResult, helpers::ExpectedConvergenceReason>;
234 template <
typename DbTagsList,
typename... InboxTags,
typename Metavariables,
235 typename ActionList,
typename ParallelComponent>
237 db::DataBox<DbTagsList>& box,
245 const ParallelComponent*
const ) noexcept {
246 const size_t element_index = get_index(element_id);
247 const auto& has_converged =
248 get<Convergence::Tags::HasConverged<OptionsGroup>>(box);
251 get<helpers::ExpectedConvergenceReason>(box));
252 const auto& expected_result =
253 gsl::at(get<ExpectedResult>(box), element_index);
254 const auto& result = get<ScalarFieldTag>(box).get();
255 for (
size_t i = 0; i < expected_result.size(); i++) {
258 return {std::move(box),
true};
263 using const_global_cache_tags = tmpl::list<Source>;
265 using simple_tags = tmpl::list<fields_tag, sources_tag>;
266 using compute_tags = tmpl::list<>;
267 template <
typename DbTagsList,
typename... InboxTags,
typename Metavariables,
268 typename ActionList,
typename ParallelComponent>
269 static auto apply(db::DataBox<DbTagsList>& box,
273 const ParallelComponent*
const ) noexcept {
274 const size_t element_index = get_index(element_id);
275 const auto& source =
gsl::at(get<Source>(box), element_index);
276 const size_t num_points = source.size();
277 ::Initialization::mutate_assign<simple_tags>(
279 typename sources_tag::type{source});
281 return std::make_tuple(std::move(box));
287 template <
typename Preconditioner>
288 struct run_preconditioner_impl {
290 tmpl::list<ComputeOperatorAction<typename Preconditioner::fields_tag>,
292 typename Preconditioner::operand_tag>>,
297 struct run_preconditioner_impl<void> {
298 using type = tmpl::list<>;
301 template <
typename Preconditioner>
302 using run_preconditioner =
303 typename run_preconditioner_impl<Preconditioner>::type;
307 template <
typename Metavariables,
308 typename LinearSolverType =
typename Metavariables::linear_solver,
309 typename PreconditionerType =
typename Metavariables::preconditioner>
310 using initialization_actions =
311 tmpl::list<Actions::SetupDataBox, dg::Actions::InitializeDomain<1>,
312 InitializeElement,
typename LinearSolverType::initialize_element,
313 ComputeOperatorAction<fields_tag>,
314 helpers::detail::init_preconditioner<PreconditionerType>,
317 template <
typename Metavariables,
318 typename LinearSolverType =
typename Metavariables::linear_solver,
319 typename PreconditionerType =
typename Metavariables::preconditioner>
320 using register_actions =
321 tmpl::list<
typename LinearSolverType::register_element,
322 helpers::detail::register_preconditioner<PreconditionerType>,
325 template <
typename Metavariables,
326 typename LinearSolverType =
typename Metavariables::linear_solver,
327 typename PreconditionerType =
typename Metavariables::preconditioner>
328 using solve_actions = tmpl::list<
329 typename LinearSolverType::template solve<tmpl::list<
330 detail::run_preconditioner<PreconditionerType>,
331 ComputeOperatorAction<typename LinearSolverType::operand_tag>>>,
334 template <
typename Metavariables,
335 typename LinearSolverType =
typename Metavariables::linear_solver,
336 typename PreconditionerType =
typename Metavariables::preconditioner>
338 tmpl::list<TestResult<typename LinearSolverType::options_group>>;
340 template <
typename Metavariables>
345 Metavariables::Phase::Initialization,
346 initialization_actions<Metavariables>>,
348 Metavariables::Phase::RegisterWithObserver,
349 register_actions<Metavariables>>,
351 Metavariables::Phase::PerformLinearSolve,
352 solve_actions<Metavariables>>,
354 Metavariables::Phase::TestResult,
355 test_actions<Metavariables>>>>;
357 template <
typename Metavariables>
358 using component_list =
360 typename Metavariables::linear_solver>,
361 helpers::detail::get_component_list<
362 typename Metavariables::preconditioner>>,
363 ElementArray<Metavariables>,