FEAT 3
Finite Element Analysis Toolbox
Loading...
Searching...
No Matches
frosch.cpp
1#include <kernel/base_header.hpp>
2#include <string>
3
4#if defined(FEAT_HAVE_TRILINOS)
5
6// preprocessor macro sanity checks
7#ifdef FROSCH_SEQUENTIAL
8# error FROSCH_SEQUENTIAL does not make sense
9#endif
10#ifdef FEAT_HAVE_MPI
11# ifndef FROSCH_HAVE_MPI
12# define FROSCH_HAVE_MPI 1
13# endif
14#else
15# error FROSch without MPI does not make sense.
16#endif // FEAT_HAVE_MPI
17
18#define FROSCH_TIMER_DETAILS 2
19
20#include <ShyLU_DDFROSch_config.h>
21
22#include <mpi.h>
23
24#include <Teuchos_GlobalMPISession.hpp>
25#include <Teuchos_CommandLineProcessor.hpp>
26#include <Teuchos_XMLParameterListCoreHelpers.hpp>
27#include <Teuchos_StackedTimer.hpp>
28
29// Thyra includes
30#include <Thyra_LinearOpWithSolveBase.hpp>
31#include <Thyra_LinearOpWithSolveFactoryHelpers.hpp>
32#include <Thyra_SolveSupportTypes.hpp>
33#include <Thyra_TpetraLinearOp.hpp>
34#include <Thyra_TpetraMultiVector.hpp>
35#include <Thyra_TpetraThyraWrappers.hpp>
36#include <Thyra_TpetraVector.hpp>
37#include <Thyra_VectorBase.hpp>
38#include <Thyra_VectorStdOps.hpp>
39#ifdef HAVE_SHYLU_DDFROSCH_EPETRA
40#include <Thyra_EpetraLinearOp.hpp>
41#endif
42#include <Thyra_VectorSpaceBase_decl.hpp>
43#include <Thyra_VectorSpaceBase_def.hpp>
44
45// Stratimikos includes
46#include <Stratimikos_DefaultLinearSolverBuilder.hpp>
47#include <Stratimikos_FROSch_def.hpp>
48
49#include <Tpetra_Core.hpp>
50#include <Tpetra_MultiVector.hpp>
51#include <Tpetra_CrsMatrix.hpp>
52
53// FROSCH thyra includes
54#include <Thyra_FROSchLinearOp_def.hpp>
55#include <Thyra_FROSchFactory_def.hpp>
56#include <FROSch_Tools_def.hpp>
57
58// undefine assert to define it in feat, which was previously defined in trilinos
59#undef ASSERT
60
61// includes, FEAT
62#include <kernel/lafem/dense_vector.hpp>
63#include <kernel/util/simple_arg_parser.hpp>
64#include <kernel/util/dist_file_io.hpp>
65#include <kernel/util/hash.hpp>
66
67// includes, system
68#include <vector>
69
70using DataType = double;
71using IndexType = FEAT::Index;
72
73namespace FEAT
74{
75 namespace Solver
76 {
77 namespace Trilinos
78 {
104 class FROSchParameterList
105 {
106 public:
107 /* Use MUMPS or UMFPACK or ILU */
108 enum DirectSolver {
109 KLU, /* Trilinos internal solver: works always */
110 MUMPS, /* Trilinos needs to be compiled with Amesos2 and Mumps */
111 UMFPACK, /* Trilinos needs to be compiled with Amesos2 and Umfpack */
112 ILU /* Trilinos needs to be compiled with Ifpack2 */
113 };
114 /* Normally, use TWOLEVELBLOCK */
115 enum Preconditioner {
116 ONELEVEL, /* one-level overlapping schwarz; just for testing */
117 TWOLEVEL, /* works only for pressure Poisson and two-level neither for multi-level pressure Poisson nor for the saddle point problem */
118 TWOLEVELBLOCK /* works for pressure Poisson (two-level and multi-level) and the saddle point problem (two-level and multi-level) */
119 };
120
121 enum ProblemType {
122 PRESSUREPOISSON,
123 SADDLEPOINT
124 };
125
126 enum PartitionType {
127 PHG, /* Parallel hyper-graph, Zoltan */
128 PARMETIS, /* Trilinos needs to be compiled with Parmetis */
129 BLOCK, /* Should work always (not good) */
130 SCOTCH, /* Scotch */
131 PTSCOTCH /* Scotch */
132 };
133
134 enum PartitionApproach {
135 PARTITION,
136 REPARTITION /* probably the best, but I'm not sure */
137 };
138
139 enum CombineOverlap {
140 FULL, /* produces a symmetric preconditioner */
141 RESTRICTED, /* __ATTENTION:__ produces a non-symmetric preconditioner. As default: Restricted with (F)GMRES should work better than FULL with PCG. */
142 AVERAGING /* __ATTENTION:__ produces a non-symmetric preconditioner. Should work better than FULL with PCG. */
143 };
144
145 enum CombineLevel {
146 ADDITIVE,
147 MULTIPLICATIVE
148 };
149
150 /* interface partition of unity
151 *
152 * Normally:
153 * velocity:
154 * In 2D use GDSW.
155 * In 3D use GDSWSTAR.
156 * pressure:
157 * Use GDSW, but for discont. pressure it does not matter which is choosen
158 *
159 */
160 enum IPOU {
161 GDSW,
162 GDSWSTAR, /* Only for three dimensional problems */
163 RGDSW
164 };
165
191 explicit FROSchParameterList(const Dist::Comm& comm_, const int dim_, const FROSchParameterList::ProblemType problem_, const int nlevels_ = 1, const std::string& name_ = "FEAT-FROSch Parameterlist");
192
213 explicit FROSchParameterList(const Dist::Comm& comm_, const int dim_, const std::string& xml_file_, const std::string& name_ = "Parameter List");
214 ~FROSchParameterList();
215
216 void read_from_xml_str(const std::string &xml_str_);
217 void read_from_xml_file(const std::string &xml_file_);
218
219 /* Print parameter list to rank zero */
220 void print() const;
221
222 /* return a pointer to the core */
223 void* get_core() const;
224
225 void create_core();
226 void delete_core();
227
228 int get_dim() const;
229 int get_dpe_velo() const;
230 int get_dpe_pres() const;
231 int get_nlevels() const;
232 bool get_use_timer() const;
233 bool get_print() const;
234
235 void set_dim(const int dim);
236 void set_nlevels(const int nlevels);
237 void set_problem_type(const ProblemType problem_);
238 void set_use_timer(const bool use_timer_);
239 void set_print(const bool print_);
240
241 void set_precond(const Preconditioner precond_);
242 void set_coarse_solver(const DirectSolver solver_);
243
244 void set_gsteps(const int gsteps_);
245 void set_gsteps(const std::vector<int>& gsteps_);
246
247 void set_overlaps(const int overlap_);
248 void set_overlaps(const std::vector<int>& overlaps_);
249
250 void set_combine_overlap(const CombineOverlap combine_mode_);
251 void set_combine_overlap(const std::vector<CombineOverlap>& combine_mode_);
252
253 void set_combine_lvl(const CombineLevel combine_mode_);
254 void set_combine_lvl(const std::vector<CombineLevel>& combine_lvl_);
255
256 void set_solvers(const DirectSolver solver_ao_, const DirectSolver solver_ext_);
257 void set_solvers(const std::vector<DirectSolver>& solvers_ao_, const std::vector<DirectSolver>& solvers_ext_);
258 void set_solvers_ao(const DirectSolver solver_);
259 void set_solvers_ao(const std::vector<DirectSolver>& solvers_);
260 void set_solvers_ext(const DirectSolver solver_);
261 void set_solvers_ext(const std::vector<DirectSolver>& solvers_);
262
263 void set_paritioner(const std::vector<int>& subregions_, const PartitionType parti_type_, const PartitionApproach parti_approach_, const double parti_imbl_);
264 void set_paritioner(const std::vector<int>& subregions_, const std::vector<PartitionType>& parti_types_, const std::vector<PartitionApproach>& parti_approach_, const std::vector<double>& parti_imbl_);
265
266 void set_subregions(const std::vector<int>& subregions_);
267
268 void set_parti_types(const PartitionType parti_type_);
269 void set_parti_types(const std::vector<PartitionType>& parti_types_);
270
271 void set_parti_approach(const PartitionApproach parti_approach_);
272 void set_parti_approach(const std::vector<PartitionApproach>& parti_approach_);
273
274 void set_parti_imbl(const double parti_imbl_);
275 void set_parti_imbl(const std::vector<double>& parti_imbl_);
276
277 void set_excludes(const bool velo_pres_, const bool pres_velo_);
278 void set_excludes(const std::vector<bool>& velo_pres_, const std::vector<bool>& pres_velo_);
279 void set_excludes_velo_pres(const bool velo_pres_);
280 void set_excludes_velo_pres(const std::vector<bool>& velo_pres_);
281 void set_excludes_pres_velo(const bool pres_velo_);
282 void set_excludes_pres_velo(const std::vector<bool>& pres_velo_);
283
284 void set_ipous(const IPOU ipou_velo_, const IPOU ipou_pres_);
285 void set_ipous(const std::vector<IPOU>& ipous_velo_, const std::vector<IPOU>& ipous_pres_);
286 void set_ipous_velo(const IPOU ipou_);
287 void set_ipous_velo(const std::vector<IPOU>& ipous_);
288 void set_ipous_pres(const IPOU ipou_);
289 void set_ipous_pres(const std::vector<IPOU>& ipous_);
290
291 void set_cspace(const bool cspace_velo_, const bool cspace_pres_);
292 void set_cspace(const std::vector<bool>& cspace_velo_, const std::vector<bool>& cspace_pres_);
293 void set_cspace_velo(const bool cspace_);
294 void set_cspace_velo(const std::vector<bool>& cspace_);
295 void set_cspace_pres(const bool cspace_);
296 void set_cspace_pres(const std::vector<bool>& cspace_);
297
298 void set_reuse_sf(const bool rsf_ao_, const bool rsf_cm_);
299 void set_reuse_sf(const std::vector<bool> &rsf_ao_, const std::vector<bool> &rsf_cm_);
300 void set_reuse_sf_ao(const bool rsf_ao_);
301 void set_reuse_sf_ao(const std::vector<bool> &rsf_ao_);
302 void set_reuse_sf_cm(const bool rsf_cm_);
303 void set_reuse_sf_cm(const std::vector<bool> &rsf_cm_);
304
305 void set_reuse_coarse(const bool reuse_cm_, const bool reuse_cb_);
306 void set_reuse_coarse(const std::vector<bool> &reuse_cm_, const std::vector<bool> &reuse_cb_);
307 void set_reuse_coarse_matrix(const bool reuse_cm_);
308 void set_reuse_coarse_matrix(const std::vector<bool> &reuse_cm_);
309 void set_reuse_coarse_basis(const bool reuse_cb_);
310 void set_reuse_coarse_basis(const std::vector<bool> &reuse_cb_);
311
312 void set_phi_dropping_threshold(const double phi_dt_);
313 void set_phi_dropping_threshold(const std::vector<double> &phi_dt_);
314
315 void set_verbose(const bool verbose_);
316
317 bool parse_args(SimpleArgParser& args);
318
319 protected:
322 void* _core;
323
324 /* for parsing of arguments */
325 const Dist::Comm& _comm;
326
327 /* spatial dimensios */
328 int _dimension;
329 /* dofs per entity (velocity) */
330 int _dpe_velo;
331 /* dofs per pressure (velocity) */
332 int _dpe_pres;
333
334 /* Number of levels (Domain Decomposition levels) before the coarse problem is solved
335 *
336 * For example:
337 * The two-level preconditioner has _nlevels = 1
338 * The three-level preconditioner has _nlevels = 2
339 * */
340 int _nlevels;
341
342 ProblemType _problem;
343 Preconditioner _preconditioner;
344 DirectSolver _solver_coarse;
345
346 std::vector<int> _gsteps;
347
348 std::vector<int> _overlaps;
349 std::vector<DirectSolver> _solvers_ao;
350 std::vector<DirectSolver> _solvers_ext;
351
352 /* __ATTENTION:__ the last entry has to be one */
353 std::vector<int> _subregions;
354 std::vector<PartitionType> _partition_types;
355 std::vector<PartitionApproach> _partition_approach;
356 std::vector<double> _partition_imbl_tol;
357
358 /* exclude blocks for the coarse space:
359 *
360 * _exclude_velo_pres: exclude the velocity-pressure coupling
361 * _exclude_pres_velo: exclude the pressure-velocity coupling
362 */
363 std::vector<bool> _exclude_pres_velo;
364 std::vector<bool> _exclude_velo_pres;
365 /*
366 * use block for coarse space. Normally, set this both to true.
367 */
368 std::vector<bool> _cspace_velo;
369 std::vector<bool> _cspace_pres;
370 /* interface partition of unity for blocks */
371 std::vector<IPOU> _ipous_velo;
372 std::vector<IPOU> _ipous_pres;
373
374 /* combine mode of overlap */
375 std::vector<CombineOverlap> _combine_ov;
376
377 /* combine mode of lvl */
378 std::vector<CombineLevel> _combine_lvl;
379
380 /* reuse symbolic factorization: algebraic overlapping Operator */
381 std::vector<bool> _ao_rsf;
382 /* reuse symbolic factorization: coarse matrix */
383 std::vector<bool> _cm_rsf;
384
385 /* reuse the complete coarse matrix */
386 std::vector<bool> _cm_r;
387 /* reuse the complete coarse basis */
388 std::vector<bool> _cb_r;
389 /* Phi: Dropping Threshold ==> default 1.e-8.
390 * The higher the dropping threshold the more Krylov iterations
391 * but the faster the construction of the coarse basis */
392 std::vector<double> _phi_dt;
393
394 std::string _name;
395
396 /* print the parameter list after construction */
397 bool _print_list;
398 /* use Trilinos timers for FROSch solver */
399 bool _use_timer;
400
401 /* verbosity flag for the FROSch/Trilinos class */
402 bool _verbose;
403
404 void init_defaults();
405
406 }; // FROSchParameterList
407
408 inline void parameterlist_set_level_id(const int levelid, Teuchos::ParameterList &params)
409 {
410 params.set("Level ID", (unsigned int)levelid);
411 }
412
413 inline void parameterlist_set_default_block_id(const int defblockid, Teuchos::ParameterList &params)
414 {
415 params.set("Default BlockId Subdomain Graph", defblockid);
416 }
417
418 inline void parameterlist_set_dim(const int dimension, Teuchos::ParameterList &params)
419 {
420 params.set("Dimension", dimension);
421 }
422
423 inline void parameterlist_set_overlap(const int overlap, Teuchos::ParameterList &params)
424 {
425 params.set("Overlap", overlap);
426 }
427
428 inline void parameterlist_set_precond(const FROSchParameterList::Preconditioner precond, const bool reuse, Teuchos::ParameterList &params)
429 {
430 if(precond == FROSchParameterList::TWOLEVEL)
431 {
432 params.set("FROSch Preconditioner Type", "TwoLevelPreconditioner");
433 } else if(precond == FROSchParameterList::TWOLEVELBLOCK)
434 {
435 params.set("FROSch Preconditioner Type", "TwoLevelBlockPreconditioner");
436 } else if(precond == FROSchParameterList::ONELEVEL)
437 {
438 params.set("FROSch Preconditioner Type", "OneLevelPreconditioner");
439 } else
440 {
441 XASSERTM(false, "Unkown FROSchParameterList::Preconditioner type");
442 }
443 params.set("Recycling", reuse);
444 }
445
446 inline void parameterlist_set_combine_lvl(const FROSchParameterList::CombineLevel combine_lvl, Teuchos::ParameterList &params)
447 {
448 if(combine_lvl == FROSchParameterList::ADDITIVE)
449 params.set("Level Combination", "Additive");
450 else if(combine_lvl == FROSchParameterList::MULTIPLICATIVE)
451 params.set("Level Combination", "Multiplicative");
452 else
453 XASSERTM(false, "Unkown FROSchParameterList::CombineLevel type");
454 }
455
456 inline void parameterlist_set_verbose(const bool verbose, Teuchos::ParameterList &params)
457 {
458 params.set("Verbose", verbose);
459 }
460
461 inline void parameterlist_header(const int dim, const int overlap, const int levelid, const int defblockid, const FROSchParameterList::Preconditioner precond, const std::string& nullspace, const FROSchParameterList::CombineLevel combine_lvl, const bool reuse, const bool verbose, Teuchos::ParameterList &params)
462 {
463 parameterlist_set_dim(dim, params);
464 parameterlist_set_overlap(overlap, params);
465 parameterlist_set_default_block_id(defblockid, params);
466 parameterlist_set_level_id(levelid, params);
467 parameterlist_set_precond(precond, reuse, params);
468 parameterlist_set_combine_lvl(combine_lvl, params);
469 parameterlist_set_verbose(verbose, params);
470 params.set("Null Space Type", nullspace);
471 params.set("OverlappingOperator Type", "AlgebraicOverlappingOperator");
472 params.set("CoarseOperator Type", "IPOUHarmonicCoarseOperator");
473 }
474
475 inline void parameterlist_set_solver(const FROSchParameterList::DirectSolver solver, Teuchos::ParameterList &params)
476 {
477 if(solver == FROSchParameterList::ILU)
478 {
479 params.set("SolverType", "Ifpack2");
480 params.set("Solver", "RILUK");
481 params.set("Ifpack2", Teuchos::ParameterList("Ifpack2"));
482 params.sublist("Ifpack2").set("fact: iluk level-of-fill", 0);
483 }else
484 {
485 params.set("SolverType", "Amesos2");
486 if(solver == FROSchParameterList::KLU)
487 params.set("Solver", "Klu");
488 else if(solver == FROSchParameterList::MUMPS)
489 params.set("Solver", "Mumps");
490 else if(solver == FROSchParameterList::UMFPACK)
491 params.set("Solver", "Umfpack");
492 else
493 XASSERTM(false, "Unkown FROSchParameterList::DirectSolver type");
494 }
495 }
496
497 inline void parameterlist_set_aoo(const FROSchParameterList::DirectSolver solver, Teuchos::ParameterList &params)
498 {
499 params.set("AlgebraicOverlappingOperator", Teuchos::ParameterList("AlgebraicOverlappingOperator"));
500 if(solver == FROSchParameterList::MUMPS)
501 params.sublist("AlgebraicOverlappingOperator").set("Reuse: Symbolic Factorization", false);
502 params.sublist("AlgebraicOverlappingOperator").set("Solver", Teuchos::ParameterList("Solver"));
503 parameterlist_set_solver(solver, params.sublist("AlgebraicOverlappingOperator").sublist("Solver"));
504 }
505
506 inline void parameterlist_set_extsolv(const FROSchParameterList::DirectSolver solver, Teuchos::ParameterList &params)
507 {
508 params.set("ExtensionSolver", Teuchos::ParameterList("ExtensionSolver"));
509 parameterlist_set_solver(solver, params.sublist("ExtensionSolver"));
510 }
511
512 inline void parameterlist_set_coarsesolv(const FROSchParameterList::DirectSolver solver, Teuchos::ParameterList &params)
513 {
514 params.set("CoarseSolver", Teuchos::ParameterList("CoarseSolver"));
515 parameterlist_set_solver(solver, params.sublist("CoarseSolver"));
516 }
517
518 inline void parameterlist_combine_overlap(const FROSchParameterList::CombineOverlap mode, Teuchos::ParameterList &params)
519 {
520 if(mode == FROSchParameterList::FULL)
521 params.set("Combine Values in Overlap", "Full");
522 else if(mode == FROSchParameterList::RESTRICTED)
523 params.set("Combine Values in Overlap", "Restricted");
524 else if(mode == FROSchParameterList::AVERAGING)
525 params.set("Combine Values in Overlap", "Averaging");
526 else
527 XASSERTM(false, "Unkown FROSchParameterList::DirectSolver type");
528 }
529
530 inline void parameterlist_set_zoltan2(const FROSchParameterList::PartitionType parti_type, const FROSchParameterList::PartitionApproach parti_app, const double parti_imbl, Teuchos::ParameterList &params)
531 {
532 params.set("Zoltan2 Parameter", Teuchos::ParameterList("Zoltan2 Parameter"));
533 if(parti_type == FROSchParameterList::PARMETIS)
534 params.sublist("Zoltan2 Parameter").set("algorithm", "parmetis");
535 else if(parti_type == FROSchParameterList::SCOTCH)
536 params.sublist("Zoltan2 Parameter").set("algorithm", "scotch");
537 else if(parti_type == FROSchParameterList::PTSCOTCH)
538 params.sublist("Zoltan2 Parameter").set("algorithm", "ptscotch");
539 else if(parti_type == FROSchParameterList::BLOCK)
540 params.sublist("Zoltan2 Parameter").set("algorithm", "block");
541 else if(parti_type == FROSchParameterList::PHG)
542 params.sublist("Zoltan2 Parameter").set("algorithm", "phg");
543 else
544 XASSERTM(false, "Unkown FROSchParameterList::PartitionType type");
545
546 /* __TODO:__ Funktioniert das auch für "block" oder nur für "parmetis" */
547 params.sublist("Zoltan2 Parameter").set("imbalance_tolerance", parti_imbl);
548 if(parti_app == FROSchParameterList::PARTITION)
549 params.sublist("Zoltan2 Parameter").set("partitioning_approach", "partition");
550 else if(parti_app == FROSchParameterList::REPARTITION)
551 params.sublist("Zoltan2 Parameter").set("partitioning_approach", "repartition");
552 else
553 XASSERTM(false, "Unkown FROSchParameterList::PartitionApproach type");
554 }
555
556 inline void parameterlist_set_distribution(const int subreg, const FROSchParameterList::PartitionType parti_type, const FROSchParameterList::PartitionApproach parti_app, const double parti_imbl, const int gsteps, Teuchos::ParameterList &params)
557 {
558 params.set("Distribution", Teuchos::ParameterList("Distribution"));
559 params.sublist("Distribution").set("NumProcs", subreg);
560 if(subreg == (int)1)
561 {
562 params.sublist("Distribution").set("Type", "linear");
563 }else
564 {
565 params.sublist("Distribution").set("Type", "ZoltanDual");
566 parameterlist_set_zoltan2(parti_type, parti_app, parti_imbl, params.sublist("Distribution"));
567 }
568 params.sublist("Distribution").set("Factor", (double)1.0);
569 params.sublist("Distribution").set("GatheringSteps", gsteps);
570 params.sublist("Distribution").set("Gathering Communication", Teuchos::ParameterList("Gathering Communication"));
571 params.sublist("Distribution").sublist("Gathering Communication").set("Send type", "Send");
572 }
573
574 inline void parameterlist_set_ipou(const FROSchParameterList::IPOU ipou, const bool verbose, Teuchos::ParameterList &params)
575 {
576 std::string ipou_str;
577 if(ipou == FROSchParameterList::GDSW)
578 ipou_str = "GDSW";
579 else if(ipou == FROSchParameterList::GDSWSTAR)
580 ipou_str = "GDSWStar";
581 else if(ipou == FROSchParameterList::RGDSW)
582 ipou_str = "RGDSW";
583 else
584 XASSERTM(false, "Unkown FROSchParameterList::IPOU type");
585 params.set("InterfacePartitionOfUnity", Teuchos::ParameterList("InterfacePartitionOfUnity"));
586 params.sublist("InterfacePartitionOfUnity").set("Verbose", verbose);
587 params.sublist("InterfacePartitionOfUnity").set("Type", ipou_str);
588 params.sublist("InterfacePartitionOfUnity").set(ipou_str, Teuchos::ParameterList(ipou_str));
589 params.sublist("InterfacePartitionOfUnity").sublist(ipou_str).set("Verbose", verbose);
590 params.sublist("InterfacePartitionOfUnity").sublist(ipou_str).set("Type", "Full");
591 params.sublist("InterfacePartitionOfUnity").sublist(ipou_str).set("Distance Function", "Constant");
592 params.sublist("InterfacePartitionOfUnity").sublist(ipou_str).set("Custom", Teuchos::ParameterList("Custom"));
593 params.sublist("InterfacePartitionOfUnity").sublist(ipou_str).sublist("Custom").set("Roots", true);
594
595 /* This suppresses the output for the constant partition of unity */
596 params.set("PartitionOfUnity", Teuchos::ParameterList("PartitionOfUnity"));
597 params.sublist("PartitionOfUnity").set("Constant", Teuchos::ParameterList("Constant"));
598 params.sublist("PartitionOfUnity").sublist("Constant").set("Verbose", verbose);
599 }
600
601 inline void parameterlist_set_block(const int blockid, const bool use_cspace, const int exclude, const FROSchParameterList::IPOU ipou, const bool verbose, Teuchos::ParameterList &params)
602 {
603 std::string blockid_str = std::to_string(blockid);
604 std::string excl_str = std::to_string(exclude);
605
606 params.set(blockid_str, Teuchos::ParameterList(blockid_str));
607 params.sublist(blockid_str).set("Use For Coarse Space", use_cspace);
608 if(verbose)
609 params.sublist(blockid_str).set("Verbosity", "All");
610 else
611 params.sublist(blockid_str).set("Verbosity", "None");
612 if(exclude > 0)
613 params.sublist(blockid_str).set("Exclude", excl_str);
614 params.sublist(blockid_str).set("Rotations", false);
615 parameterlist_set_ipou(ipou, verbose, params.sublist(blockid_str));
616 }
617
618 inline void parameterlist_set_reuse_coarse(const bool rsf_cm, const bool reuse_cm, const bool reuse_cb, Teuchos::ParameterList &params)
619 {
620 params.set("Reuse: Coarse Matrix Symbolic Factorization", rsf_cm);
621 params.set("Reuse: Coarse Matrix", reuse_cm);
622 params.set("Reuse: Coarse Basis", reuse_cb);
623 }
624
625 FROSchParameterList::FROSchParameterList(const Dist::Comm& comm_, const int dim_, const FROSchParameterList::ProblemType problem_, const int nlevels_, const std::string& name_) :
626 _core(nullptr),
627 _comm(comm_),
628 _dimension(dim_),
629 _dpe_velo(dim_),
630 _dpe_pres(dim_+1),
631 _nlevels(nlevels_),
632 _problem(problem_),
633 _name(name_),
634 _print_list(false),
635 _use_timer(false)
636 {
637 XASSERTM(dim_ == 2 || dim_ == 3, "Only implemented for dim == 2 or dim == 3");
638 init_defaults();
639 }
640
641 FROSchParameterList::FROSchParameterList(const Dist::Comm& comm_, const int dim_, const std::string& xml_file_, const std::string& name_) :
642 _core(nullptr),
643 _comm(comm_),
644 _dimension(dim_),
645 _dpe_velo(dim_),
646 _dpe_pres(dim_+1),
647 _nlevels(1),
648 _name(name_),
649 _print_list(false),
650 _use_timer(false)
651 {
652 read_from_xml_file(xml_file_);
653 }
654
655 FROSchParameterList::~FROSchParameterList()
656 {
657 delete_core();
658 }
659
660 void FROSchParameterList::read_from_xml_str(const std::string& xml_str_)
661 {
662 delete_core();
663
664 Teuchos::ParameterList *param = new Teuchos::ParameterList(_name);
665 Teuchos::Ptr<Teuchos::ParameterList> ptr_param(param);
666 Teuchos::updateParametersFromXmlString(xml_str_, ptr_param);
667
668 _core = (void*)ptr_param.get();
669 }
670
671 void FROSchParameterList::read_from_xml_file(const std::string& xml_file_)
672 {
673 std::stringstream xml_str;
674 DistFileIO::read_common(xml_str, xml_file_, _comm);
675 read_from_xml_str(xml_str.str());
676 }
677
678 void FROSchParameterList::init_defaults()
679 {
680 set_overlaps(1);
681 set_gsteps(1);
682
683 /* Es fehlen noch die Partitionsparameter für nlevel > 1 */
684 if(_nlevels == 1)
685 set_subregions(std::vector<int>(1, 1));
686
687 set_parti_types(FROSchParameterList::PARMETIS);
688 set_parti_approach(FROSchParameterList::REPARTITION);
689 set_parti_imbl(1.1);
690 if(_dimension == 2)
691 set_ipous(GDSW, GDSW);
692 else
693 set_ipous(GDSWSTAR, GDSW);
694
695 set_solvers(FROSchParameterList::KLU, FROSchParameterList::KLU);
696 _solver_coarse = FROSchParameterList::KLU;
697 _preconditioner = FROSchParameterList::TWOLEVELBLOCK;
698 set_cspace(true, true);
699 set_excludes(false, false);
700 /* __ATTENTION:__ Is produces a non-symmetric preconditioner */
701 set_combine_overlap(FROSchParameterList::RESTRICTED);
702 set_combine_lvl(FROSchParameterList::ADDITIVE);
703
704 set_reuse_sf(false, false);
705 set_reuse_coarse(false, false);
706 set_phi_dropping_threshold(1.e-8);
707
708 set_verbose(false);
709 }
710
711 void FROSchParameterList::print() const
712 {
713 if(_comm.rank() == 0)
714 {
715 XASSERTM(_core != nullptr, "_core == nullptr");
716 (reinterpret_cast<Teuchos::ParameterList*>(_core))->print(std::cout, 4);
717 }
718 }
719
720 void* FROSchParameterList::get_core() const
721 {
722 return _core;
723 }
724
725 void FROSchParameterList::create_core()
726 {
727 if(_core != nullptr)
728 {
729 /* core is already created */
730 return;
731 }
732
733 XASSERTM(_nlevels == int(_subregions.size()), "number of subregions is not consistent with number of levels.");
734 XASSERTM(_subregions.back() == 1, "last entry of subregions is not 1.");
735
736 std::string precond_str;
737 if(_preconditioner == FROSchParameterList::TWOLEVEL)
738 precond_str = "TwoLevelPreconditioner";
739 else if(_preconditioner == FROSchParameterList::TWOLEVELBLOCK)
740 precond_str = "TwoLevelBlockPreconditioner";
741 else if(_preconditioner == FROSchParameterList::ONELEVEL)
742 precond_str = "OneLevelPreconditioner";
743 else
744 XASSERTM(false, "Unkown FROSchParameterList::Preconditioner type");
745
746 Teuchos::ParameterList *params = new Teuchos::ParameterList(_name);
747 params->set("Preconditioner Type", "FROSch");
748 params->set("Preconditioner Types", Teuchos::ParameterList("Preconditioner Types"));
749 params->sublist("Preconditioner Types").set("FROSch", Teuchos::ParameterList("FROSch"));
750
751 int overlap, gsteps;
752 DirectSolver solver_ao, solver_ext;
753 CombineOverlap mode;
754 IPOU ipou_velo, ipou_pres;
755 bool cspace_velo, cspace_pres, excl_velo_pres, excl_pres_velo;
756 bool rsf_ao, rsf_cm, reuse_cm, reuse_cb;
757 double phi_dt;
758 PartitionType partitype;
759 PartitionApproach partiapp;
760 double partiimbl;
761 CombineLevel combine_lvl;
762 bool reuse;
763 bool verbose;
764
765 Teuchos::ParameterList *tmp_params = &(params->sublist("Preconditioner Types").sublist("FROSch"));
766 for(int i = 0; i < _nlevels; ++i)
767 {
768 overlap = _overlaps.size() == 1 ? _overlaps.at(0) : _overlaps.at(i);
769 solver_ao = _solvers_ao.size() == 1 ? _solvers_ao.at(0) : _solvers_ao.at(i);
770 solver_ext = _solvers_ext.size() == 1 ? _solvers_ext.at(0) : _solvers_ext.at(i);
771 mode = _combine_ov.size() == 1 ? _combine_ov.at(0) : _combine_ov.at(i);
772 combine_lvl = _combine_ov.size() == 1 ? _combine_lvl.at(0) : _combine_lvl.at(i);
773 ipou_velo = _ipous_velo.size() == 1 ? _ipous_velo.at(0) : _ipous_velo.at(i);
774 ipou_pres = _ipous_pres.size() == 1 ? _ipous_pres.at(0) : _ipous_pres.at(i);
775 cspace_velo = _cspace_velo.size() == 1 ? _cspace_velo.at(0) : _cspace_velo.at(i);
776 cspace_pres = _cspace_pres.size() == 1 ? _cspace_pres.at(0) : _cspace_pres.at(i);
777 excl_velo_pres = _exclude_velo_pres.size() == 1 ? _exclude_velo_pres.at(0) : _exclude_velo_pres.at(i);
778 excl_pres_velo = _exclude_pres_velo.size() == 1 ? _exclude_pres_velo.at(0) : _exclude_pres_velo.at(i);
779 partitype = _partition_types.size() == 1 ? _partition_types.at(0) : _partition_types.at(i);
780 partiapp = _partition_approach.size() == 1 ? _partition_approach.at(0) : _partition_approach.at(i);
781 partiimbl = _partition_imbl_tol.size() == 1 ? _partition_imbl_tol.at(0) : _partition_imbl_tol.at(i);
782 gsteps = _gsteps.size() == 1 ? _gsteps.at(0) : _gsteps.at(i);
783 rsf_ao = _ao_rsf.size() == 1 ? _ao_rsf.at(0) : _ao_rsf.at(i);
784 rsf_cm = _cm_rsf.size() == 1 ? _cm_rsf.at(0) : _cm_rsf.at(i);
785 reuse_cm = _cm_r.size() == 1 ? _cm_r.at(0) : _cm_r.at(i);
786 reuse_cb = _cb_r.size() == 1 ? _cb_r.at(0) : _cb_r.at(i);
787 phi_dt = _phi_dt.size() == 1 ? _phi_dt.at(0) : _phi_dt.at(i);
788
789 reuse = rsf_ao || rsf_cm || reuse_cb || reuse_cm;
790 verbose = _verbose;
791
792 parameterlist_header(_dimension, overlap, i+1, 1, _preconditioner, (i == 0 ? "Input" : "Laplace"), combine_lvl, reuse, verbose, *tmp_params);
793 parameterlist_set_aoo(solver_ao, *tmp_params);
794 parameterlist_combine_overlap(mode, tmp_params->sublist("AlgebraicOverlappingOperator"));
795 tmp_params->sublist("AlgebraicOverlappingOperator").set("Reuse: Symbolic Factorization", rsf_ao);
796 parameterlist_set_verbose(verbose, tmp_params->sublist("AlgebraicOverlappingOperator"));
797
798 tmp_params->set("IPOUHarmonicCoarseOperator", Teuchos::ParameterList("IPOUHarmonicCoarseOperator"));
799 tmp_params->sublist("IPOUHarmonicCoarseOperator").set("Blocks", Teuchos::ParameterList("Blocks"));
800 parameterlist_set_reuse_coarse(rsf_cm, reuse_cm, reuse_cb, tmp_params->sublist("IPOUHarmonicCoarseOperator"));
801 tmp_params->sublist("IPOUHarmonicCoarseOperator").set("Phi: Dropping Threshold", phi_dt);
802 parameterlist_set_verbose(verbose, tmp_params->sublist("IPOUHarmonicCoarseOperator"));
803
804 if(_problem == FROSchParameterList::PRESSUREPOISSON)
805 {
806 /* for pressure poisson: we always use the pressure in the coarse space and do not exclude the velocity (not existent) */
807 /* pressure block */
808 parameterlist_set_block(1, true, -1, ipou_pres, _verbose, tmp_params->sublist("IPOUHarmonicCoarseOperator").sublist("Blocks"));
809 }else if(_problem == FROSchParameterList::SADDLEPOINT)
810 {
811 /* velocity block */
812 parameterlist_set_block(1, cspace_velo, (excl_velo_pres ? 2 : -1), ipou_velo, _verbose, tmp_params->sublist("IPOUHarmonicCoarseOperator").sublist("Blocks"));
813 /* pressure block */
814 parameterlist_set_block(2, cspace_pres, (excl_pres_velo ? 1 : -1), ipou_pres, _verbose, tmp_params->sublist("IPOUHarmonicCoarseOperator").sublist("Blocks"));
815
816 }else
817 XASSERTM(false, "Unkown FROSchParameterList::ProblemType type");
818
819 parameterlist_set_extsolv(solver_ext, tmp_params->sublist("IPOUHarmonicCoarseOperator"));
820 parameterlist_set_distribution(_subregions.at(i), partitype, partiapp, partiimbl, gsteps, tmp_params->sublist("IPOUHarmonicCoarseOperator"));
821
822 /* next level */
823 tmp_params->sublist("IPOUHarmonicCoarseOperator").set("CoarseSolver", Teuchos::ParameterList("CoarseSolver"));
824
825 if(i+1 < _nlevels)
826 {
827 /* there is another level */
828 tmp_params->sublist("IPOUHarmonicCoarseOperator").sublist("CoarseSolver").set("SolverType", "FROSchPreconditioner");
829 tmp_params->sublist("IPOUHarmonicCoarseOperator").sublist("CoarseSolver").set("Solver", precond_str);
830 tmp_params->sublist("IPOUHarmonicCoarseOperator").sublist("CoarseSolver").set("FROSchPreconditioner", Teuchos::ParameterList("FROSchPreconditioner"));
831 tmp_params->sublist("IPOUHarmonicCoarseOperator").sublist("CoarseSolver").sublist("FROSchPreconditioner").set(precond_str, Teuchos::ParameterList(precond_str));
832
833 tmp_params = &(tmp_params->sublist("IPOUHarmonicCoarseOperator").sublist("CoarseSolver").sublist("FROSchPreconditioner").sublist(precond_str));
834 }else
835 {
836 /* add coarse solver */
837 parameterlist_set_solver(_solver_coarse, tmp_params->sublist("IPOUHarmonicCoarseOperator").sublist("CoarseSolver"));
838 }
839 }
840
841 _core = (void*) params;
842
843 if(_print_list)
844 {
845 print();
846 }
847
848 }
849
850 void FROSchParameterList::delete_core()
851 {
852 if(_core != nullptr)
853 {
854 // /* DEBUG */
855 // if(_print_list)
856 // print();
857 delete reinterpret_cast<Teuchos::ParameterList*>(_core);
858 }
859 }
860
861 int FROSchParameterList::get_dim() const
862 {
863 return _dimension;
864 }
865
866 int FROSchParameterList::get_dpe_velo() const
867 {
868 return _dpe_velo;
869 }
870
871 int FROSchParameterList::get_dpe_pres() const
872 {
873 return _dpe_pres;
874 }
875
876 int FROSchParameterList::get_nlevels() const
877 {
878 return _nlevels;
879 }
880
881 bool FROSchParameterList::get_use_timer() const
882 {
883 return _use_timer;
884 }
885
886 bool FROSchParameterList::get_print() const
887 {
888 return _print_list;
889 }
890
891 void FROSchParameterList::set_print(const bool print_)
892 {
893 _print_list = print_;
894 }
895
896 void FROSchParameterList::set_use_timer(const bool use_timer_)
897 {
898 _use_timer = use_timer_;
899 }
900
901 void FROSchParameterList::set_nlevels(const int nlevels_)
902 {
903 _nlevels = nlevels_;
904 }
905
906 void FROSchParameterList::set_dim(const int dim_)
907 {
908 XASSERTM(dim_ == 2 || dim_ == 3, "Only implemented for dim_ == 2 or dim_ == 3");
909 _dimension = dim_;
910 _dpe_velo = dim_;
911 _dpe_pres = dim_ + 1;
912 }
913
914 void FROSchParameterList::set_problem_type(const ProblemType problem_)
915 {
916 _problem = problem_;
917 }
918
919 void FROSchParameterList::set_precond(const Preconditioner precond_)
920 {
921 _preconditioner = precond_;
922 }
923
924 void FROSchParameterList::set_coarse_solver(const DirectSolver solver_)
925 {
926 _solver_coarse = solver_;
927 }
928
929 void FROSchParameterList::set_gsteps(const int gsteps_)
930 {
931 _gsteps.assign(_nlevels, gsteps_);
932 }
933
934 void FROSchParameterList::set_gsteps(const std::vector<int>& gsteps_)
935 {
936 XASSERT(int(gsteps_.size()) == _nlevels);
937 _gsteps.assign(gsteps_.begin(), gsteps_.end());
938 }
939
940 void FROSchParameterList::set_overlaps(const int overlap_)
941 {
942 _overlaps.assign(_nlevels, overlap_);
943 }
944
945 void FROSchParameterList::set_overlaps(const std::vector<int>& overlaps_)
946 {
947 XASSERT(int(overlaps_.size()) == _nlevels);
948 _overlaps.assign(overlaps_.begin(), overlaps_.end());
949 }
950
951 void FROSchParameterList::set_combine_overlap(const CombineOverlap combine_mode_)
952 {
953 _combine_ov.assign(_nlevels, combine_mode_);
954 }
955
956 void FROSchParameterList::set_combine_overlap(const std::vector<CombineOverlap>& combine_mode_)
957 {
958 XASSERT(int(combine_mode_.size()) == _nlevels);
959 _combine_ov.assign(combine_mode_.begin(), combine_mode_.end());
960 }
961
962 void FROSchParameterList::set_combine_lvl(const CombineLevel combine_lvl_)
963 {
964 _combine_lvl.assign(_nlevels, combine_lvl_);
965 }
966
967 void FROSchParameterList::set_combine_lvl(const std::vector<CombineLevel>& combine_lvl_)
968 {
969 XASSERT(int(_combine_lvl.size()) == _nlevels);
970 _combine_lvl.assign(combine_lvl_.begin(), combine_lvl_.end());
971 }
972
973 void FROSchParameterList::set_solvers(const DirectSolver solver_ao_, const DirectSolver solver_ext_)
974 {
975 set_solvers_ao(solver_ao_);
976 set_solvers_ext(solver_ext_);
977 }
978
979 void FROSchParameterList::set_solvers(const std::vector<DirectSolver>& solver_ao_, const std::vector<DirectSolver>& solver_ext_)
980 {
981 set_solvers_ao(solver_ao_);
982 set_solvers_ext(solver_ext_);
983 }
984
985 void FROSchParameterList::set_solvers_ao(const DirectSolver solver_)
986 {
987 _solvers_ao.assign(_nlevels, solver_);
988 }
989
990 void FROSchParameterList::set_solvers_ao(const std::vector<DirectSolver>& solvers_)
991 {
992 XASSERT(int(solvers_.size()) == _nlevels);
993 _solvers_ao.assign(solvers_.begin(), solvers_.end());
994 }
995
996 void FROSchParameterList::set_solvers_ext(const DirectSolver solver_)
997 {
998 _solvers_ext.assign(_nlevels, solver_);
999 }
1000
1001 void FROSchParameterList::set_solvers_ext(const std::vector<DirectSolver>& solvers_)
1002 {
1003 XASSERT(int(solvers_.size()) == _nlevels);
1004 _solvers_ext.assign(solvers_.begin(), solvers_.end());
1005 }
1006
1007 void FROSchParameterList::set_paritioner(const std::vector<int>& subregions_, const PartitionType parti_type_, const PartitionApproach parti_approach_, const double parti_imbl_)
1008 {
1009 set_subregions(subregions_);
1010 set_parti_types(parti_type_);
1011 set_parti_approach(parti_approach_);
1012 set_parti_imbl(parti_imbl_);
1013 }
1014
1015 void FROSchParameterList::set_paritioner(const std::vector<int>& subregions_, const std::vector<PartitionType>& parti_types_, const std::vector<PartitionApproach>& parti_approach_, const std::vector<double>& parti_imbl_)
1016 {
1017 set_subregions(subregions_);
1018 set_parti_types(parti_types_);
1019 set_parti_approach(parti_approach_);
1020 set_parti_imbl(parti_imbl_);
1021 }
1022
1023 void FROSchParameterList::set_subregions(const std::vector<int>& subregions_)
1024 {
1025 XASSERT(int(subregions_.size()) == _nlevels);
1026 XASSERTM(int(subregions_.back()) == int(1), "Last entry has to be one.");
1027 _subregions.assign(subregions_.begin(), subregions_.end());
1028 }
1029
1030 void FROSchParameterList::set_parti_types(const PartitionType parti_type_)
1031 {
1032 _partition_types.assign(_nlevels, parti_type_);
1033 }
1034
1035 void FROSchParameterList::set_parti_types(const std::vector<PartitionType>& parti_types_)
1036 {
1037 XASSERT(int(parti_types_.size()) == _nlevels);
1038 _partition_types.assign(parti_types_.begin(), parti_types_.end());
1039 }
1040
1041 void FROSchParameterList::set_parti_approach(const PartitionApproach parti_approach_)
1042 {
1043 _partition_approach.assign(_nlevels, parti_approach_);
1044 }
1045
1046 void FROSchParameterList::set_parti_approach(const std::vector<PartitionApproach>& parti_approach_)
1047 {
1048 XASSERT(int(parti_approach_.size()) == _nlevels);
1049 _partition_approach.assign(parti_approach_.begin(), parti_approach_.end());
1050 }
1051
1052 void FROSchParameterList::set_parti_imbl(const double parti_imbl_)
1053 {
1054 _partition_imbl_tol.assign(_nlevels, parti_imbl_);
1055 }
1056
1057 void FROSchParameterList::set_parti_imbl(const std::vector<double>& parti_imbl_)
1058 {
1059 XASSERT(int(parti_imbl_.size()) == _nlevels);
1060 _partition_imbl_tol.assign(parti_imbl_.begin(), parti_imbl_.end());
1061 }
1062
1063 void FROSchParameterList::set_excludes(const bool velo_pres_, const bool pres_velo_)
1064 {
1065 _exclude_velo_pres.assign(_nlevels, velo_pres_);
1066 _exclude_pres_velo.assign(_nlevels, pres_velo_);
1067 }
1068
1069 void FROSchParameterList::set_excludes(const std::vector<bool>& velo_pres_, const std::vector<bool>& pres_velo_)
1070 {
1071 XASSERT(int(velo_pres_.size()) == _nlevels);
1072 XASSERT(int(pres_velo_.size()) == _nlevels);
1073 _exclude_velo_pres.assign(velo_pres_.begin(), velo_pres_.end());
1074 _exclude_pres_velo.assign(pres_velo_.begin(), pres_velo_.end());
1075 }
1076
1077 void FROSchParameterList::set_excludes_velo_pres(const bool velo_pres_)
1078 {
1079 _exclude_velo_pres.assign(_nlevels, velo_pres_);
1080 }
1081
1082 void FROSchParameterList::set_excludes_velo_pres(const std::vector<bool>& velo_pres_)
1083 {
1084 XASSERT(int(velo_pres_.size()) == _nlevels);
1085 _exclude_velo_pres.assign(velo_pres_.begin(), velo_pres_.end());
1086 }
1087
1088 void FROSchParameterList::set_excludes_pres_velo(const bool pres_velo_)
1089 {
1090 _exclude_pres_velo.assign(_nlevels, pres_velo_);
1091 }
1092
1093 void FROSchParameterList::set_excludes_pres_velo(const std::vector<bool>& pres_velo_)
1094 {
1095 XASSERT(int(pres_velo_.size()) == _nlevels);
1096 _exclude_pres_velo.assign(pres_velo_.begin(), pres_velo_.end());
1097 }
1098
1099 void FROSchParameterList::set_ipous(const IPOU ipou_velo_, const IPOU ipou_pres_)
1100 {
1101 set_ipous_velo(ipou_velo_);
1102 set_ipous_pres(ipou_pres_);
1103 }
1104
1105 void FROSchParameterList::set_ipous(const std::vector<IPOU>& ipous_velo_, const std::vector<IPOU>& ipous_pres_)
1106 {
1107 set_ipous_velo(ipous_velo_);
1108 set_ipous_pres(ipous_pres_);
1109 }
1110
1111 void FROSchParameterList::set_ipous_velo(const IPOU ipou_)
1112 {
1113 _ipous_velo.assign(_nlevels, ipou_);
1114 }
1115
1116 void FROSchParameterList::set_ipous_velo(const std::vector<IPOU>& ipous_)
1117 {
1118 XASSERT(int(ipous_.size()) == _nlevels);
1119 _ipous_velo.assign(ipous_.begin(), ipous_.end());
1120 }
1121
1122 void FROSchParameterList::set_ipous_pres(const IPOU ipou_)
1123 {
1124 _ipous_pres.assign(_nlevels, ipou_);
1125 }
1126
1127 void FROSchParameterList::set_ipous_pres(const std::vector<IPOU>& ipous_)
1128 {
1129 XASSERT(int(ipous_.size()) == _nlevels);
1130 _ipous_pres.assign(ipous_.begin(), ipous_.end());
1131 }
1132
1133 void FROSchParameterList::set_cspace(const bool cspace_velo_, const bool cspace_pres_)
1134 {
1135 set_cspace_velo(cspace_velo_);
1136 set_cspace_pres(cspace_pres_);
1137 }
1138
1139 void FROSchParameterList::set_cspace(const std::vector<bool>& cspace_velo_, const std::vector<bool>& cspace_pres_)
1140 {
1141 set_cspace_velo(cspace_velo_);
1142 set_cspace_pres(cspace_pres_);
1143 }
1144
1145 void FROSchParameterList::set_cspace_velo(const bool cspace_)
1146 {
1147 _cspace_velo.assign(_nlevels, cspace_);
1148 }
1149
1150 void FROSchParameterList::set_cspace_velo(const std::vector<bool>& cspace_)
1151 {
1152 XASSERT(int(cspace_.size()) == _nlevels);
1153 _cspace_velo.assign(cspace_.begin(), cspace_.end());
1154 }
1155
1156 void FROSchParameterList::set_cspace_pres(const bool cspace_)
1157 {
1158 _cspace_pres.assign(_nlevels, cspace_);
1159 }
1160
1161 void FROSchParameterList::set_cspace_pres(const std::vector<bool>& cspace_)
1162 {
1163 XASSERT(int(cspace_.size()) == _nlevels);
1164 _cspace_pres.assign(cspace_.begin(), cspace_.end());
1165 }
1166
1167 void FROSchParameterList::set_reuse_sf(const bool rsf_ao_, const bool rsf_cm_)
1168 {
1169 set_reuse_sf_ao(rsf_ao_);
1170 set_reuse_sf_cm(rsf_cm_);
1171 }
1172
1173 void FROSchParameterList::set_reuse_sf(const std::vector<bool> &rsf_ao_, const std::vector<bool> &rsf_cm_)
1174 {
1175 set_reuse_sf_ao(rsf_ao_);
1176 set_reuse_sf_cm(rsf_cm_);
1177 }
1178
1179 void FROSchParameterList::set_reuse_sf_ao(const bool rsf_ao_)
1180 {
1181 _ao_rsf.assign(_nlevels, rsf_ao_);
1182 }
1183
1184 void FROSchParameterList::set_reuse_sf_ao(const std::vector<bool> &rsf_ao_)
1185 {
1186 XASSERT(int(rsf_ao_.size()) == _nlevels);
1187 _ao_rsf.assign(rsf_ao_.begin(), rsf_ao_.end());
1188 }
1189
1190 void FROSchParameterList::set_reuse_sf_cm(const bool rsf_cm_)
1191 {
1192 _cm_rsf.assign(_nlevels, rsf_cm_);
1193 }
1194
1195 void FROSchParameterList::set_reuse_sf_cm(const std::vector<bool> &rsf_cm_)
1196 {
1197 XASSERT(int(rsf_cm_.size()) == _nlevels);
1198 _cm_rsf.assign(rsf_cm_.begin(), rsf_cm_.end());
1199 }
1200
1201 void FROSchParameterList::set_reuse_coarse(const bool reuse_cm_, const bool reuse_cb_)
1202 {
1203 set_reuse_coarse_matrix(reuse_cm_);
1204 set_reuse_coarse_basis(reuse_cb_);
1205 }
1206
1207 void FROSchParameterList::set_reuse_coarse(const std::vector<bool> &reuse_cm_, const std::vector<bool> &reuse_cb_)
1208 {
1209 set_reuse_coarse_matrix(reuse_cm_);
1210 set_reuse_coarse_basis(reuse_cb_);
1211 }
1212
1213 void FROSchParameterList::set_reuse_coarse_matrix(const bool reuse_cm_)
1214 {
1215 _cm_r.assign(_nlevels, reuse_cm_);
1216 }
1217
1218 void FROSchParameterList::set_reuse_coarse_matrix(const std::vector<bool> &reuse_cm_)
1219 {
1220 XASSERT(int(reuse_cm_.size()) == _nlevels);
1221 _cm_r.assign(reuse_cm_.begin(), reuse_cm_.end());
1222 }
1223
1224 void FROSchParameterList::set_reuse_coarse_basis(const bool rsf_cm_)
1225 {
1226 _cb_r.assign(_nlevels, rsf_cm_);
1227 }
1228
1229 void FROSchParameterList::set_reuse_coarse_basis(const std::vector<bool> &reuse_cb_)
1230 {
1231 XASSERT(int(reuse_cb_.size()) == _nlevels);
1232 _cb_r.assign(reuse_cb_.begin(), reuse_cb_.end());
1233 }
1234
1235 void FROSchParameterList::set_phi_dropping_threshold(const double phi_dt_)
1236 {
1237 _phi_dt.assign(_nlevels, phi_dt_);
1238 }
1239
1240 void FROSchParameterList::set_phi_dropping_threshold(const std::vector<double> &phi_dt_)
1241 {
1242 XASSERT(int(phi_dt_.size()) == _nlevels);
1243 _phi_dt.assign(phi_dt_.begin(), phi_dt_.end());
1244 }
1245
1246 void FROSchParameterList::set_verbose(const bool verbose_)
1247 {
1248 _verbose = verbose_;
1249 }
1250
1251 inline bool parse_fpl_levels(const Dist::Comm &comm, SimpleArgParser &args, FROSchParameterList &params)
1252 {
1253 auto it = args.query("frosch-nlevels");
1254 if(it != nullptr)
1255 {
1256 if(it->second.size() == 1)
1257 {
1258 int nlevels = std::stoi(it->second.front());
1259 if(nlevels < 1)
1260 {
1261 comm.print("ERROR: given number of levels is smaller then 1");
1262 }
1263 params.set_nlevels(nlevels);
1264 }else
1265 {
1266 comm.print("ERROR: multiple levels given");
1267 return false;
1268 }
1269 }
1270 return true;
1271 }
1272
1273 inline bool parse_fpl_dim(const Dist::Comm& comm, SimpleArgParser& args, FROSchParameterList &params)
1274 {
1275 auto it = args.query("frosch-dim");
1276 if(it != nullptr)
1277 {
1278 if(it->second.size() == 1)
1279 {
1280 int dim = std::stoi(it->second.front());
1281 if(dim != 2 && dim != 3)
1282 {
1283 comm.print("ERROR: given dimension has to be 2 or 3");
1284 }
1285 params.set_dim(dim);
1286 }else
1287 {
1288 comm.print("ERROR: multiple dimensions given");
1289 return false;
1290 }
1291 }
1292
1293 return true;
1294 }
1295
1296 inline bool parse_fpl_overlaps(const Dist::Comm& comm, SimpleArgParser& args, FROSchParameterList &params)
1297 {
1298 auto it = args.query("frosch-overlap");
1299 if(it != nullptr)
1300 {
1301 if(it->second.size() == 1u)
1302 {
1303 auto ov = std::stoi(it->second.front());
1304 if(ov < 0 )
1305 {
1306 comm.print("ERROR: overlap is smaller than 0");
1307 return false;
1308 }
1309 params.set_overlaps((int)ov);
1310 }else if(int(it->second.size()) == params.get_nlevels())
1311 {
1312 std::vector<int> overlaps;
1313 overlaps.reserve(params.get_nlevels());
1314 for(const auto& t : it->second)
1315 {
1316 auto ov = std::stoi(t);
1317 if(ov < 0)
1318 {
1319 comm.print("ERROR: overlap is smaller than 0");
1320 return false;
1321 }
1322 overlaps.push_back(ov);
1323 }
1324 params.set_overlaps(overlaps);
1325 }else
1326 {
1327 comm.print("ERROR: not matching number of overlaps: nlevels: " + std::to_string(params.get_nlevels()) + " number of given overlaps: " + std::to_string(it->second.size()) + "'");
1328 return false;
1329 }
1330 }
1331 return true;
1332 }
1333
1334 inline bool parse_fpl_subreg(const Dist::Comm& comm, SimpleArgParser& args, FROSchParameterList &params)
1335 {
1336 auto it = args.query("frosch-subreg");
1337 if(it != nullptr)
1338 {
1339 if(int(it->second.size()) == params.get_nlevels())
1340 {
1341 std::vector<int> subreg;
1342 subreg.reserve(params.get_nlevels());
1343 for(const auto& t : it->second)
1344 {
1345 auto sr = std::stoi(t);
1346 if(sr < 0)
1347 {
1348 /* __TODO:__ Eigentlich müsste man hier testen auf "< 1", aber der letzte Eintrag muss 1 sein */
1349 comm.print("ERROR: subreg is smaller than 0");
1350 return false;
1351 }
1352 subreg.push_back(sr);
1353 }
1354 if(subreg.back() != 1)
1355 {
1356 comm.print("ERROR: last entry of subreg is not 1: " + std::to_string(subreg.back()));
1357 }
1358 params.set_subregions(subreg);
1359 }else
1360 {
1361 comm.print("ERROR: not matching number of subregions: nlevels: " + std::to_string(params.get_nlevels()) + " number of given subregions: " + std::to_string(it->second.size()) + "'");
1362 return false;
1363 }
1364 }
1365 return true;
1366 }
1367
1368 inline bool parse_fpl_gsteps(const Dist::Comm& comm, SimpleArgParser& args, FROSchParameterList &params)
1369 {
1370 auto it = args.query("frosch-gstep");
1371 if(it != nullptr)
1372 {
1373 if(it->second.size() == 1u)
1374 {
1375 auto gsteps = std::stoi(it->second.front());
1376 if(gsteps < 1 )
1377 {
1378 comm.print("ERROR: gstep is smaller than 1");
1379 return false;
1380 }
1381 params.set_gsteps((int)gsteps);
1382 }else if(int(it->second.size()) == params.get_nlevels())
1383 {
1384 std::vector<int> gsteps;
1385 gsteps.reserve(params.get_nlevels());
1386 for(const auto& t : it->second)
1387 {
1388 auto gs = std::stoi(t);
1389 if(gs < 1)
1390 {
1391 comm.print("ERROR: gstep is smaller than 0");
1392 return false;
1393 }
1394 gsteps.push_back(gs);
1395 }
1396 params.set_gsteps(gsteps);
1397 }else
1398 {
1399 comm.print("ERROR: not matching number of gsteps: nlevels: " + std::to_string(params.get_nlevels()) + " number of given gsteps: " + std::to_string(it->second.size()) + "'");
1400 return false;
1401 }
1402 }
1403 return true;
1404 }
1405
1406 inline bool parse_fpl_parti_type(const Dist::Comm& comm, SimpleArgParser& args, FROSchParameterList &params)
1407 {
1408 auto it = args.query("frosch-parti-type");
1409 if(it != nullptr)
1410 {
1411 if(it->second.size() == 1u)
1412 {
1413 FROSchParameterList::PartitionType parti_type;
1414 std::string type = it->second.front();
1415 std::transform(type.begin(), type.end(), type.begin(), [](unsigned char c){return std::tolower(c);});
1416
1417 if(!type.compare("parmetis"))
1418 {
1419 parti_type = Solver::Trilinos::FROSchParameterList::PARMETIS;
1420 }else if(!type.compare("scotch"))
1421 {
1422 parti_type = Solver::Trilinos::FROSchParameterList::SCOTCH;
1423 }else if(!type.compare("ptscotch"))
1424 {
1425 parti_type = Solver::Trilinos::FROSchParameterList::PTSCOTCH;
1426 }else if(!type.compare("block"))
1427 {
1428 parti_type = Solver::Trilinos::FROSchParameterList::BLOCK;
1429 }else if(!type.compare("phg"))
1430 {
1431 parti_type = Solver::Trilinos::FROSchParameterList::PHG;
1432 }else
1433 {
1434 comm.print("ERROR: unknown partitioner type '" + it->second.front() + "'");
1435 return false;
1436 }
1437 params.set_parti_types(parti_type);
1438 }else if(int(it->second.size()) == params.get_nlevels())
1439 {
1440 std::vector<FROSchParameterList::PartitionType> parti_types;
1441 parti_types.reserve(params.get_nlevels());
1442 for(const auto& t : it->second)
1443 {
1444 std::string type = t;
1445 std::transform(type.begin(), type.end(), type.begin(), [](unsigned char c){return std::tolower(c);});
1446
1447 if(!type.compare("parmetis"))
1448 {
1449 parti_types.push_back(Solver::Trilinos::FROSchParameterList::PARMETIS);
1450 }else if(!type.compare("scotch"))
1451 {
1452 parti_types.push_back(Solver::Trilinos::FROSchParameterList::SCOTCH);
1453 }else if(!type.compare("ptscotch"))
1454 {
1455 parti_types.push_back(Solver::Trilinos::FROSchParameterList::PTSCOTCH);
1456 }else if(!type.compare("block"))
1457 {
1458 parti_types.push_back(Solver::Trilinos::FROSchParameterList::BLOCK);
1459 }else if(!type.compare("phg"))
1460 {
1461 parti_types.push_back(Solver::Trilinos::FROSchParameterList::PHG);
1462 }else
1463 {
1464 comm.print("ERROR: unknown partitioner type '" + t + "'");
1465 return false;
1466 }
1467 }
1468 params.set_parti_types(parti_types);
1469 }else
1470 {
1471 comm.print("ERROR: not matching number of parti-type: nlevels: " + std::to_string(params.get_nlevels()) + " number of given parti-type: " + std::to_string(it->second.size()) + "'");
1472 return false;
1473 }
1474 }
1475 return true;
1476 }
1477
1478 inline bool parse_fpl_parti_approach(const Dist::Comm& comm, SimpleArgParser& args, FROSchParameterList &params)
1479 {
1480 auto it = args.query("frosch-parti-approach");
1481 if(it != nullptr)
1482 {
1483 if(it->second.size() == 1u)
1484 {
1485 FROSchParameterList::PartitionApproach parti_app;
1486 std::string type = it->second.front();
1487 std::transform(type.begin(), type.end(), type.begin(), [](unsigned char c){return std::tolower(c);});
1488
1489 if(!type.compare("partition"))
1490 {
1491 parti_app = FROSchParameterList::PARTITION;
1492 }else if(!type.compare("repartition"))
1493 {
1494 parti_app = FROSchParameterList::REPARTITION;
1495 }else
1496 {
1497 comm.print("ERROR: unknown partitioner approach '" + it->second.front() + "'");
1498 return false;
1499 }
1500 params.set_parti_approach(parti_app);
1501 }else if(int(it->second.size()) == params.get_nlevels())
1502 {
1503 std::vector<FROSchParameterList::PartitionApproach> parti_apps;
1504 parti_apps.reserve(params.get_nlevels());
1505 for(const auto& t : it->second)
1506 {
1507 std::string type = t;
1508 std::transform(type.begin(), type.end(), type.begin(), [](unsigned char c){return std::tolower(c);});
1509
1510 if(!type.compare("partition"))
1511 {
1512 parti_apps.push_back(FROSchParameterList::PARTITION);
1513 }else if(!type.compare("repartitionblock"))
1514 {
1515 parti_apps.push_back(FROSchParameterList::REPARTITION);
1516 }else
1517 {
1518 comm.print("ERROR: unknown partitioner approach '" + t + "'");
1519 return false;
1520 }
1521 }
1522 params.set_parti_approach(parti_apps);
1523 }else
1524 {
1525 comm.print("ERROR: not matching number of parti-approach: nlevels: " + std::to_string(params.get_nlevels()) + " number of given parti-approach: " + std::to_string(it->second.size()) + "'");
1526 return false;
1527 }
1528 }
1529 return true;
1530 }
1531
1532 inline bool parse_fpl_parti_imbl(const Dist::Comm& comm, SimpleArgParser& args, FROSchParameterList &params)
1533 {
1534 auto it = args.query("frosch-parti-imbl");
1535 if(it != nullptr)
1536 {
1537 if(it->second.size() == 1u)
1538 {
1539 double parti_imbl = std::stod(it->second.front());
1540 if(parti_imbl < 1.)
1541 {
1542 /* __TODO:__ Stimmt das? */
1543 comm.print("ERROR: unknown partition imblance smaller than 1.: '" + it->second.front() + "'");
1544 return false;
1545 }
1546 params.set_parti_imbl(parti_imbl);
1547 }else if(int(it->second.size()) == params.get_nlevels())
1548 {
1549 std::vector<double> parti_imbl;
1550 parti_imbl.reserve(params.get_nlevels());
1551 for(const auto& t : it->second)
1552 {
1553 double imbl = std::stod(t);
1554 if(imbl < 1.)
1555 {
1556 /* __TODO:__ Stimmt das? */
1557 comm.print("ERROR: unknown partition imblance smaller than 1.: '" + t + "'");
1558 return false;
1559 }
1560 parti_imbl.push_back(imbl);
1561 }
1562 params.set_parti_imbl(parti_imbl);
1563 }else
1564 {
1565 comm.print("ERROR: not matching number of parti-imbl: nlevels: " + std::to_string(params.get_nlevels()) + " number of given parti-imbl: " + std::to_string(it->second.size()) + "'");
1566 return false;
1567 }
1568 }
1569 return true;
1570 }
1571
1572 inline bool parse_fpl_ipou_pres(const Dist::Comm& comm, SimpleArgParser& args, FROSchParameterList &params)
1573 {
1574 auto it = args.query("frosch-ipou-pres");
1575 if(it != nullptr)
1576 {
1577 if(it->second.size() == 1u)
1578 {
1579 FROSchParameterList::IPOU ipou;
1580 std::string type = it->second.front();
1581 std::transform(type.begin(), type.end(), type.begin(), [](unsigned char c){return std::tolower(c);});
1582
1583 if(!type.compare("gdsw"))
1584 {
1585 ipou = FROSchParameterList::GDSW;
1586 }else if(!type.compare("gdswstar"))
1587 {
1588 ipou = FROSchParameterList::GDSWSTAR;
1589 }else if(!type.compare("rgdsw"))
1590 {
1591 ipou = FROSchParameterList::RGDSW;
1592 }else
1593 {
1594 comm.print("ERROR: unknown ipou pressure type '" + it->second.front() + "'");
1595 return false;
1596 }
1597 params.set_ipous_pres(ipou);
1598 }else if(int(it->second.size()) == params.get_nlevels())
1599 {
1600 std::vector<FROSchParameterList::IPOU> ipous;
1601 ipous.reserve(params.get_nlevels());
1602 for(const auto& t : it->second)
1603 {
1604 std::string type = t;
1605 std::transform(type.begin(), type.end(), type.begin(), [](unsigned char c){return std::tolower(c);});
1606
1607 if(!type.compare("gdsw"))
1608 {
1609 ipous.push_back(FROSchParameterList::GDSW);
1610 }else if(!type.compare("gdswstar"))
1611 {
1612 ipous.push_back(FROSchParameterList::GDSWSTAR);
1613 }else if(!type.compare("rgdsw"))
1614 {
1615 ipous.push_back(FROSchParameterList::RGDSW);
1616 }else
1617 {
1618 comm.print("ERROR: unknown ipou pressure type '" + t + "'");
1619 return false;
1620 }
1621 }
1622 params.set_ipous_pres(ipous);
1623 }else
1624 {
1625 comm.print("ERROR: not matching number of ipou-pres: nlevels: " + std::to_string(params.get_nlevels()) + " number of given ipou-pres: " + std::to_string(it->second.size()) + "'");
1626 return false;
1627 }
1628 }
1629 return true;
1630 }
1631
1632 inline bool parse_fpl_ipou_velo(const Dist::Comm& comm, SimpleArgParser& args, FROSchParameterList &params)
1633 {
1634 auto it = args.query("frosch-ipou-velo");
1635 if(it != nullptr)
1636 {
1637 if(it->second.size() == 1u)
1638 {
1639 FROSchParameterList::IPOU ipou;
1640 std::string type = it->second.front();
1641 std::transform(type.begin(), type.end(), type.begin(), [](unsigned char c){return std::tolower(c);});
1642
1643 if(!type.compare("gdsw"))
1644 {
1645 ipou = FROSchParameterList::GDSW;
1646 }else if(!type.compare("gdswstar"))
1647 {
1648 ipou = FROSchParameterList::GDSWSTAR;
1649 }else if(!type.compare("rgdsw"))
1650 {
1651 ipou = FROSchParameterList::RGDSW;
1652 }else
1653 {
1654 comm.print("ERROR: unknown ipou velocity type '" + it->second.front() + "'");
1655 return false;
1656 }
1657 params.set_ipous_velo(ipou);
1658 }else if(int(it->second.size()) == params.get_nlevels())
1659 {
1660 std::vector<FROSchParameterList::IPOU> ipous;
1661 ipous.reserve(params.get_nlevels());
1662 for(const auto& t : it->second)
1663 {
1664 std::string type = t;
1665 std::transform(type.begin(), type.end(), type.begin(), [](unsigned char c){return std::tolower(c);});
1666
1667 if(!type.compare("gdsw"))
1668 {
1669 ipous.push_back(FROSchParameterList::GDSW);
1670 }else if(!type.compare("gdswstar"))
1671 {
1672 ipous.push_back(FROSchParameterList::GDSWSTAR);
1673 }else if(!type.compare("rgdsw"))
1674 {
1675 ipous.push_back(FROSchParameterList::RGDSW);
1676 }else
1677 {
1678 comm.print("ERROR: unknown ipou velo type '" + t + "'");
1679 return false;
1680 }
1681 }
1682 params.set_ipous_pres(ipous);
1683 }else
1684 {
1685 comm.print("ERROR: not matching number of ipou-velo: nlevels: " + std::to_string(params.get_nlevels()) + " number of given ipou-velo: " + std::to_string(it->second.size()) + "'");
1686 return false;
1687 }
1688 }
1689 return true;
1690 }
1691
1692 inline bool parse_fpl_solver_ao(const Dist::Comm& comm, SimpleArgParser& args, FROSchParameterList &params)
1693 {
1694 auto it = args.query("frosch-solver-ao");
1695 if(it != nullptr)
1696 {
1697 if(it->second.size() == 1u)
1698 {
1699 FROSchParameterList::DirectSolver solver;
1700 std::string type = it->second.front();
1701 std::transform(type.begin(), type.end(), type.begin(), [](unsigned char c){return std::tolower(c);});
1702
1703 if(!type.compare("klu"))
1704 {
1705 solver = FROSchParameterList::KLU;
1706 }else if(!type.compare("mumps"))
1707 {
1708 solver = FROSchParameterList::MUMPS;
1709 }else if(!type.compare("umfpack"))
1710 {
1711 solver = FROSchParameterList::UMFPACK;
1712 }else
1713 {
1714 comm.print("ERROR: unknown solver-ao type '" + it->second.front() + "'");
1715 return false;
1716 }
1717 params.set_solvers_ao(solver);
1718 }else if(int(it->second.size()) == params.get_nlevels())
1719 {
1720 std::vector<FROSchParameterList::DirectSolver> solvers;
1721 solvers.reserve(params.get_nlevels());
1722 for(const auto& t : it->second)
1723 {
1724 std::string type = t;
1725 std::transform(type.begin(), type.end(), type.begin(), [](unsigned char c){return std::tolower(c);});
1726
1727 if(!type.compare("klu"))
1728 {
1729 solvers.push_back(FROSchParameterList::KLU);
1730 }else if(!type.compare("mumps"))
1731 {
1732 solvers.push_back(FROSchParameterList::MUMPS);
1733 }else if(!type.compare("umfpack"))
1734 {
1735 solvers.push_back(FROSchParameterList::UMFPACK);
1736 }else
1737 {
1738 comm.print("ERROR: unknown solver-ao type '" + t + "'");
1739 return false;
1740 }
1741 }
1742 params.set_solvers_ao(solvers);
1743 }else
1744 {
1745 comm.print("ERROR: not matching number of solver-ao: nlevels: " + std::to_string(params.get_nlevels()) + " number of given solver-ao: " + std::to_string(it->second.size()) + "'");
1746 return false;
1747 }
1748 }
1749 return true;
1750 }
1751
1752 inline bool parse_fpl_solver_ext(const Dist::Comm& comm, SimpleArgParser& args, FROSchParameterList &params)
1753 {
1754 auto it = args.query("frosch-solver-ext");
1755 if(it != nullptr)
1756 {
1757 if(it->second.size() == 1u)
1758 {
1759 FROSchParameterList::DirectSolver solver;
1760 std::string type = it->second.front();
1761 std::transform(type.begin(), type.end(), type.begin(), [](unsigned char c){return std::tolower(c);});
1762
1763 if(!type.compare("klu"))
1764 {
1765 solver = FROSchParameterList::KLU;
1766 }else if(!type.compare("mumps"))
1767 {
1768 solver = FROSchParameterList::MUMPS;
1769 }else if(!type.compare("umfpack"))
1770 {
1771 solver = FROSchParameterList::UMFPACK;
1772 }else
1773 {
1774 comm.print("ERROR: unknown solver-ext type '" + it->second.front() + "'");
1775 return false;
1776 }
1777 params.set_solvers_ext(solver);
1778 }else if(int(it->second.size()) == params.get_nlevels())
1779 {
1780 std::vector<FROSchParameterList::DirectSolver> solvers;
1781 solvers.reserve(params.get_nlevels());
1782 for(const auto& t : it->second)
1783 {
1784 std::string type = t;
1785 std::transform(type.begin(), type.end(), type.begin(), [](unsigned char c){return std::tolower(c);});
1786
1787 if(!type.compare("klu"))
1788 {
1789 solvers.push_back(FROSchParameterList::KLU);
1790 }else if(!type.compare("mumps"))
1791 {
1792 solvers.push_back(FROSchParameterList::MUMPS);
1793 }else if(!type.compare("umfpack"))
1794 {
1795 solvers.push_back(FROSchParameterList::UMFPACK);
1796 }else
1797 {
1798 comm.print("ERROR: unknown solver-ext type '" + t + "'");
1799 return false;
1800 }
1801 }
1802 params.set_solvers_ext(solvers);
1803 }else
1804 {
1805 comm.print("ERROR: not matching number of solver-ext: nlevels: " + std::to_string(params.get_nlevels()) + " number of given solver-ext: " + std::to_string(it->second.size()) + "'");
1806 return false;
1807 }
1808 }
1809 return true;
1810 }
1811
1812 inline bool parse_fpl_solver_coarse(const Dist::Comm& comm, SimpleArgParser& args, FROSchParameterList &params)
1813 {
1814 auto it = args.query("frosch-solver-coarse");
1815 if(it != nullptr)
1816 {
1817 if(it->second.size() == 1u)
1818 {
1819 FROSchParameterList::DirectSolver solver;
1820 std::string type = it->second.front();
1821 std::transform(type.begin(), type.end(), type.begin(), [](unsigned char c){return std::tolower(c);});
1822
1823 if(!type.compare("klu"))
1824 {
1825 solver = FROSchParameterList::KLU;
1826 }else if(!type.compare("mumps"))
1827 {
1828 solver = FROSchParameterList::MUMPS;
1829 }else if(!type.compare("umfpack"))
1830 {
1831 solver = FROSchParameterList::UMFPACK;
1832 }else
1833 {
1834 comm.print("ERROR: unknown solver-coarse type '" + it->second.front() + "'");
1835 return false;
1836 }
1837 params.set_coarse_solver(solver);
1838 }else
1839 {
1840 comm.print("ERROR: multiple coarse solvers are given");
1841 return false;
1842 }
1843 }
1844 return true;
1845 }
1846
1847 inline bool parse_fpl_precond(const Dist::Comm& comm, SimpleArgParser& args, FROSchParameterList &params)
1848 {
1849 auto it = args.query("frosch-precond-type");
1850 if(it != nullptr)
1851 {
1852 if(it->second.size() == 1u)
1853 {
1854 FROSchParameterList::Preconditioner precond;
1855 std::string type = it->second.front();
1856 std::transform(type.begin(), type.end(), type.begin(), [](unsigned char c){return std::tolower(c);});
1857
1858 if(!type.compare("twolevel"))
1859 {
1860 precond = FROSchParameterList::TWOLEVEL;
1861 }else if(!type.compare("twolevelblock"))
1862 {
1863 precond = FROSchParameterList::TWOLEVELBLOCK;
1864 }else if(!type.compare("onelevel"))
1865 {
1866 precond = FROSchParameterList::ONELEVEL;
1867 }else
1868 {
1869 comm.print("ERROR: unknown precond-type '" + it->second.front() + "'");
1870 return false;
1871 }
1872 params.set_precond(precond);
1873 }else
1874 {
1875 comm.print("ERROR: multiple preconditioner types are given");
1876 return false;
1877 }
1878 }
1879 return true;
1880 }
1881
1882 inline bool parse_fpl_cspace_pres(const Dist::Comm& comm, SimpleArgParser& args, FROSchParameterList &params)
1883 {
1884 auto it = args.query("frosch-cspace-pres");
1885 if(it != nullptr)
1886 {
1887 if(it->second.size() == 1u)
1888 {
1889 bool cspace;
1890 std::string type = it->second.front();
1891 std::transform(type.begin(), type.end(), type.begin(), [](unsigned char c){return std::tolower(c);});
1892
1893 if(!type.compare("true"))
1894 {
1895 cspace = true;
1896 }else if(!type.compare("false"))
1897 {
1898 cspace = false;
1899 }else
1900 {
1901 comm.print("ERROR: unknown cspace-pres type '" + it->second.front() + "'");
1902 return false;
1903 }
1904 params.set_cspace_pres(cspace);
1905 }else if(int(it->second.size()) == params.get_nlevels())
1906 {
1907 std::vector<bool> cspaces;
1908 cspaces.reserve(params.get_nlevels());
1909 for(const auto& t : it->second)
1910 {
1911 std::string type = t;
1912 std::transform(type.begin(), type.end(), type.begin(), [](unsigned char c){return std::tolower(c);});
1913
1914 if(!type.compare("true"))
1915 {
1916 cspaces.push_back(true);
1917 }else if(!type.compare("false"))
1918 {
1919 cspaces.push_back(false);
1920 }else
1921 {
1922 comm.print("ERROR: unknown cspace-pres type '" + t + "'");
1923 return false;
1924 }
1925 }
1926 params.set_cspace_pres(cspaces);
1927 }else
1928 {
1929 comm.print("ERROR: not matching number of cspace-pres: nlevels: " + std::to_string(params.get_nlevels()) + " number of given cspace-pres: " + std::to_string(it->second.size()) + "'");
1930 return false;
1931 }
1932 }
1933 return true;
1934 }
1935
1936 inline bool parse_fpl_cspace_velo(const Dist::Comm& comm, SimpleArgParser& args, FROSchParameterList &params)
1937 {
1938 auto it = args.query("frosch-cspace-velo");
1939 if(it != nullptr)
1940 {
1941 if(it->second.size() == 1u)
1942 {
1943 bool cspace;
1944 std::string type = it->second.front();
1945 std::transform(type.begin(), type.end(), type.begin(), [](unsigned char c){return std::tolower(c);});
1946
1947 if(!type.compare("true"))
1948 {
1949 cspace = true;
1950 }else if(!type.compare("false"))
1951 {
1952 cspace = false;
1953 }else
1954 {
1955 comm.print("ERROR: unknown cspace-velo type '" + it->second.front() + "'");
1956 return false;
1957 }
1958 params.set_cspace_velo(cspace);
1959 }else if(int(it->second.size()) == params.get_nlevels())
1960 {
1961 std::vector<bool> cspaces;
1962 cspaces.reserve(params.get_nlevels());
1963 for(const auto& t : it->second)
1964 {
1965 std::string type = t;
1966 std::transform(type.begin(), type.end(), type.begin(), [](unsigned char c){return std::tolower(c);});
1967
1968 if(!type.compare("true"))
1969 {
1970 cspaces.push_back(true);
1971 }else if(!type.compare("false"))
1972 {
1973 cspaces.push_back(false);
1974 }else
1975 {
1976 comm.print("ERROR: unknown cspace-velo type '" + t + "'");
1977 return false;
1978 }
1979 }
1980 params.set_cspace_velo(cspaces);
1981 }else
1982 {
1983 comm.print("ERROR: not matching number of cspace-velo: nlevels: " + std::to_string(params.get_nlevels()) + " number of given cspace-velo: " + std::to_string(it->second.size()) + "'");
1984 return false;
1985 }
1986 }
1987 return true;
1988 }
1989
1990 inline bool parse_fpl_exclude_velo_pres(const Dist::Comm& comm, SimpleArgParser& args, FROSchParameterList &params)
1991 {
1992 auto it = args.query("frosch-exclude-velo-pres");
1993 if(it != nullptr)
1994 {
1995 if(it->second.size() == 1u)
1996 {
1997 bool excl;
1998 std::string type = it->second.front();
1999 std::transform(type.begin(), type.end(), type.begin(), [](unsigned char c){return std::tolower(c);});
2000
2001 if(!type.compare("true"))
2002 {
2003 excl = true;
2004 }else if(!type.compare("false"))
2005 {
2006 excl = false;
2007 }else
2008 {
2009 comm.print("ERROR: unknown exclude-velo-pres type '" + it->second.front() + "'");
2010 return false;
2011 }
2012 params.set_excludes_velo_pres(excl);
2013 }else if(int(it->second.size()) == params.get_nlevels())
2014 {
2015 std::vector<bool> excl;
2016 excl.reserve(params.get_nlevels());
2017 for(const auto& t : it->second)
2018 {
2019 std::string type = t;
2020 std::transform(type.begin(), type.end(), type.begin(), [](unsigned char c){return std::tolower(c);});
2021
2022 if(!type.compare("true"))
2023 {
2024 excl.push_back(true);
2025 }else if(!type.compare("false"))
2026 {
2027 excl.push_back(false);
2028 }else
2029 {
2030 comm.print("ERROR: unknown exclude-velo-pres type '" + t + "'");
2031 return false;
2032 }
2033 }
2034 params.set_excludes_velo_pres(excl);
2035 }else
2036 {
2037 comm.print("ERROR: not matching number of exclude-velo-pres: nlevels: " + std::to_string(params.get_nlevels()) + " number of given exclude-velo-pres: " + std::to_string(it->second.size()) + "'");
2038 return false;
2039 }
2040 }
2041 return true;
2042 }
2043
2044 inline bool parse_fpl_exclude_pres_velo(const Dist::Comm& comm, SimpleArgParser& args, FROSchParameterList &params)
2045 {
2046 auto it = args.query("frosch-exclude-pres-velo");
2047 if(it != nullptr)
2048 {
2049 if(it->second.size() == 1u)
2050 {
2051 bool excl;
2052 std::string type = it->second.front();
2053 std::transform(type.begin(), type.end(), type.begin(), [](unsigned char c){return std::tolower(c);});
2054
2055 if(!type.compare("true"))
2056 {
2057 excl = true;
2058 }else if(!type.compare("false"))
2059 {
2060 excl = false;
2061 }else
2062 {
2063 comm.print("ERROR: unknown exclude-velo-pres type '" + it->second.front() + "'");
2064 return false;
2065 }
2066 params.set_excludes_pres_velo(excl);
2067 }else if(int(it->second.size()) == params.get_nlevels())
2068 {
2069 std::vector<bool> excl;
2070 excl.reserve(params.get_nlevels());
2071 for(const auto& t : it->second)
2072 {
2073 std::string type = t;
2074 std::transform(type.begin(), type.end(), type.begin(), [](unsigned char c){return std::tolower(c);});
2075
2076 if(!type.compare("true"))
2077 {
2078 excl.push_back(true);
2079 }else if(!type.compare("false"))
2080 {
2081 excl.push_back(false);
2082 }else
2083 {
2084 comm.print("ERROR: unknown exclude-pres-velo type '" + t + "'");
2085 return false;
2086 }
2087 }
2088 params.set_excludes_pres_velo(excl);
2089 }else
2090 {
2091 comm.print("ERROR: not matching number of exclude-pres-velo: nlevels: " + std::to_string(params.get_nlevels()) + " number of given exclude-pres-velo: " + std::to_string(it->second.size()) + "'");
2092 return false;
2093 }
2094 }
2095 return true;
2096 }
2097
2098 inline bool parse_fpl_print(const Dist::Comm& comm, SimpleArgParser& args, FROSchParameterList &params)
2099 {
2100 auto it = args.query("frosch-print-list");
2101 if(it != nullptr)
2102 {
2103 if(it->second.size() == 1u)
2104 {
2105 bool print;
2106 std::string type = it->second.front();
2107 std::transform(type.begin(), type.end(), type.begin(), [](unsigned char c){return std::tolower(c);});
2108 if(!type.compare("true"))
2109 {
2110 print = true;
2111 }else if(!type.compare("false"))
2112 {
2113 print = false;
2114 }else
2115 {
2116 comm.print("ERROR: unknown print-list '" + it->second.front() + "'");
2117 return false;
2118 }
2119 params.set_print(print);
2120 }else
2121 {
2122 comm.print("ERROR: more than one print-list arg given: number of args: " + std::to_string(it->second.size()) + "'");
2123 return false;
2124 }
2125 }
2126 return true;
2127 }
2128
2129 inline bool parse_fpl_verbose(const Dist::Comm& comm, SimpleArgParser& args, FROSchParameterList &params)
2130 {
2131 auto it = args.query("frosch-verbose");
2132 if(it != nullptr)
2133 {
2134 if(it->second.size() == 1)
2135 {
2136 bool verbose;
2137 std::string type = it->second.front();
2138 std::transform(type.begin(), type.end(), type.begin(), [](unsigned char c){return std::tolower(c);});
2139 if(!type.compare("true"))
2140 {
2141 verbose = true;
2142 }else if(!type.compare("false"))
2143 {
2144 verbose = false;
2145 }else
2146 {
2147 comm.print("ERROR: unknown verbose '" + it->second.front() + "'");
2148 return false;
2149 }
2150 params.set_verbose(verbose);
2151 }else
2152 {
2153 comm.print("ERROR: more than one verbose arg given: number of args: " + std::to_string(it->second.size()) + "'");
2154 return false;
2155 }
2156 }
2157 return true;
2158 }
2159
2160 inline bool parse_fpl_timer(const Dist::Comm& comm, SimpleArgParser& args, FROSchParameterList &params)
2161 {
2162 auto it = args.query("frosch-use-timer");
2163 if(it != nullptr)
2164 {
2165 if(it->second.size() == 1u)
2166 {
2167 bool timer;
2168 std::string type = it->second.front();
2169 std::transform(type.begin(), type.end(), type.begin(), [](unsigned char c){return std::tolower(c);});
2170
2171 if(!type.compare("true"))
2172 {
2173 timer = true;
2174 }else if(!type.compare("false"))
2175 {
2176 timer = false;
2177 }else
2178 {
2179 comm.print("ERROR: unknown use-timer '" + it->second.front() + "'");
2180 return false;
2181 }
2182 params.set_use_timer(timer);
2183 }else
2184 {
2185 comm.print("ERROR: more than one use-timer arg given: number of args: " + std::to_string(it->second.size()) + "'");
2186 return false;
2187 }
2188 }
2189 return true;
2190 }
2191
2192 inline bool parse_fpl_combine(const Dist::Comm& comm, SimpleArgParser& args, FROSchParameterList &params)
2193 {
2194 auto it = args.query("frosch-combine-overlap");
2195 if(it != nullptr)
2196 {
2197 if(it->second.size() == 1u)
2198 {
2199 FROSchParameterList::CombineOverlap mode;
2200 std::string type = it->second.front();
2201 std::transform(type.begin(), type.end(), type.begin(), [](unsigned char c){return std::tolower(c);});
2202
2203 if(!type.compare("full"))
2204 {
2205 mode = FROSchParameterList::FULL;
2206 }else if(!type.compare("restricted"))
2207 {
2208 mode = FROSchParameterList::RESTRICTED;
2209 }else if(!type.compare("averaging"))
2210 {
2211 mode = FROSchParameterList::AVERAGING;
2212 }else
2213 {
2214 comm.print("ERROR: unknown combine-overlap type '" + it->second.front() + "'");
2215 return false;
2216 }
2217 params.set_combine_overlap(mode);
2218 }else if(int(it->second.size()) == params.get_nlevels())
2219 {
2220 std::vector<FROSchParameterList::CombineOverlap> modes;
2221 modes.reserve(params.get_nlevels());
2222 for(const auto& t : it->second)
2223 {
2224 std::string type = t;
2225 std::transform(type.begin(), type.end(), type.begin(), [](unsigned char c){return std::tolower(c);});
2226
2227 if(!type.compare("full"))
2228 {
2229 modes.push_back(FROSchParameterList::FULL);
2230 }else if(!type.compare("restricted"))
2231 {
2232 modes.push_back(FROSchParameterList::RESTRICTED);
2233 }else if(!type.compare("umfpack"))
2234 {
2235 modes.push_back(FROSchParameterList::AVERAGING);
2236 }else
2237 {
2238 comm.print("ERROR: unknown combine-overlap type '" + t + "'");
2239 return false;
2240 }
2241 }
2242 params.set_combine_overlap(modes);
2243 }else
2244 {
2245 comm.print("ERROR: not matching number of combine-overlap: nlevels: " + std::to_string(params.get_nlevels()) + " number of given combine-overlap: " + std::to_string(it->second.size()) + "'");
2246 return false;
2247 }
2248 }
2249 return true;
2250 }
2251
2252 inline bool parse_fpl_combine_lvl(const Dist::Comm& comm, SimpleArgParser& args, FROSchParameterList &params)
2253 {
2254 auto it = args.query("frosch-combine-lvl");
2255 if(it != nullptr)
2256 {
2257 if(it->second.size() == 1u)
2258 {
2259 FROSchParameterList::CombineLevel mode;
2260 std::string type = it->second.front();
2261 std::transform(type.begin(), type.end(), type.begin(), [](unsigned char c){return std::tolower(c);});
2262
2263 if(!type.compare("additive"))
2264 {
2265 mode = FROSchParameterList::ADDITIVE;
2266 }else if(!type.compare("multiplicative"))
2267 {
2268 mode = FROSchParameterList::MULTIPLICATIVE;
2269 }else
2270 {
2271 comm.print("ERROR: unknown combine-lvl type '" + it->second.front() + "'");
2272 return false;
2273 }
2274 params.set_combine_lvl(mode);
2275 }else if(int(it->second.size()) == params.get_nlevels())
2276 {
2277 std::vector<FROSchParameterList::CombineLevel> modes;
2278 modes.reserve(params.get_nlevels());
2279 for(const auto& t : it->second)
2280 {
2281 std::string type = t;
2282 std::transform(type.begin(), type.end(), type.begin(), [](unsigned char c){return std::tolower(c);});
2283
2284 if(!type.compare("additive"))
2285 {
2286 modes.push_back(FROSchParameterList::ADDITIVE);
2287 }else if(!type.compare("multiplicative"))
2288 {
2289 modes.push_back(FROSchParameterList::MULTIPLICATIVE);
2290 }else
2291 {
2292 comm.print("ERROR: unknown combine-lvl type '" + t + "'");
2293 return false;
2294 }
2295 }
2296 params.set_combine_lvl(modes);
2297 }else
2298 {
2299 comm.print("ERROR: not matching number of combine-lvl: nlevels: " + std::to_string(params.get_nlevels()) + " number of given combine-lvl: " + std::to_string(it->second.size()) + "'");
2300 return false;
2301 }
2302 }
2303 return true;
2304 }
2305
2306 inline bool parse_fpl_reuse_sf_ao(const Dist::Comm &comm, SimpleArgParser &args, FROSchParameterList &params)
2307 {
2308 auto it = args.query("frosch-reuse-sf-ao");
2309 if(it != nullptr)
2310 {
2311 if(it->second.size() == 1u)
2312 {
2313 bool rsf_ao;
2314 std::string type = it->second.front();
2315 std::transform(type.begin(), type.end(), type.begin(), [](unsigned char c){return std::tolower(c);});
2316
2317 if(!type.compare("true"))
2318 {
2319 rsf_ao = true;
2320 }else if(!type.compare("false"))
2321 {
2322 rsf_ao = false;
2323 }else
2324 {
2325 comm.print("ERROR: unknown reuse-sf-ao '" + it->second.front() + "'");
2326 return false;
2327 }
2328 params.set_reuse_sf_ao(rsf_ao);
2329 }else if(int(it->second.size()) == params.get_nlevels())
2330 {
2331 std::vector<bool> rsf_ao;
2332 rsf_ao.reserve(params.get_nlevels());
2333 for(const auto& t : it->second)
2334 {
2335 std::string type = t;
2336 std::transform(type.begin(), type.end(), type.begin(), [](unsigned char c){return std::tolower(c);});
2337
2338 if(!type.compare("true"))
2339 {
2340 rsf_ao.push_back(true);
2341 }else if(!type.compare("false"))
2342 {
2343 rsf_ao.push_back(false);
2344 }else
2345 {
2346 comm.print("ERROR: unknown reuse-sf-ao type '" + t + "'");
2347 return false;
2348 }
2349 }
2350 params.set_reuse_sf_ao(rsf_ao);
2351 }else
2352 {
2353 comm.print("ERROR: not matching number of reuse-sf-ao: nlevels: " + std::to_string(params.get_nlevels()) + " number of given reuse-sf-ao: " + std::to_string(it->second.size()) + "'");
2354 return false;
2355 }
2356 }
2357 return true;
2358 }
2359
2360 inline bool parse_fpl_reuse_sf_cm(const Dist::Comm &comm, SimpleArgParser &args, FROSchParameterList &params)
2361 {
2362 auto it = args.query("frosch-reuse-sf-cm");
2363 if(it != nullptr)
2364 {
2365 if(it->second.size() == 1u)
2366 {
2367 bool rsf_cm;
2368 std::string type = it->second.front();
2369 std::transform(type.begin(), type.end(), type.begin(), [](unsigned char c){return std::tolower(c);});
2370
2371 if(!type.compare("true"))
2372 {
2373 rsf_cm = true;
2374 }else if(!type.compare("false"))
2375 {
2376 rsf_cm = false;
2377 }else
2378 {
2379 comm.print("ERROR: unknown reuse-sf-cm '" + it->second.front() + "'");
2380 return false;
2381 }
2382 params.set_reuse_sf_cm(rsf_cm);
2383 }else if(int(it->second.size()) == params.get_nlevels())
2384 {
2385 std::vector<bool> rsf_cm;
2386 rsf_cm.reserve(params.get_nlevels());
2387 for(const auto& t : it->second)
2388 {
2389 std::string type = t;
2390 std::transform(type.begin(), type.end(), type.begin(), [](unsigned char c){return std::tolower(c);});
2391
2392 if(!type.compare("true"))
2393 {
2394 rsf_cm.push_back(true);
2395 }else if(!type.compare("false"))
2396 {
2397 rsf_cm.push_back(false);
2398 }else
2399 {
2400 comm.print("ERROR: unknown reuse-sf-cm type '" + t + "'");
2401 return false;
2402 }
2403 }
2404 params.set_reuse_sf_cm(rsf_cm);
2405 }else
2406 {
2407 comm.print("ERROR: not matching number of reuse-sf-cm: nlevels: " + std::to_string(params.get_nlevels()) + " number of given reuse-sf-cm: " + std::to_string(it->second.size()) + "'");
2408 return false;
2409 }
2410 }
2411 return true;
2412 }
2413
2414 inline bool parse_fpl_reuse_coarse_matrix(const Dist::Comm &comm, SimpleArgParser &args, FROSchParameterList &params)
2415 {
2416 auto it = args.query("frosch-reuse-coarse-matrix");
2417 if(it != nullptr)
2418 {
2419 if(it->second.size() == 1u)
2420 {
2421 bool reuse_cm;
2422 std::string type = it->second.front();
2423 std::transform(type.begin(), type.end(), type.begin(), [](unsigned char c){return std::tolower(c);});
2424
2425 if(!type.compare("true"))
2426 {
2427 reuse_cm = true;
2428 }else if(!type.compare("false"))
2429 {
2430 reuse_cm = false;
2431 }else
2432 {
2433 comm.print("ERROR: unknown reuse-coarse-matrix '" + it->second.front() + "'");
2434 return false;
2435 }
2436 params.set_reuse_coarse_matrix(reuse_cm);
2437 }else if(int(it->second.size()) == params.get_nlevels())
2438 {
2439 std::vector<bool> reuse_cm;
2440 reuse_cm.reserve(params.get_nlevels());
2441 for(const auto& t : it->second)
2442 {
2443 std::string type = t;
2444 std::transform(type.begin(), type.end(), type.begin(), [](unsigned char c){return std::tolower(c);});
2445
2446 if(!type.compare("true"))
2447 {
2448 reuse_cm.push_back(true);
2449 }else if(!type.compare("false"))
2450 {
2451 reuse_cm.push_back(false);
2452 }else
2453 {
2454 comm.print("ERROR: unknown reuse-coarse-matrix type '" + t + "'");
2455 return false;
2456 }
2457 }
2458 params.set_reuse_coarse_matrix(reuse_cm);
2459 }else
2460 {
2461 comm.print("ERROR: not matching number of reuse-coarse-matrix: nlevels: " + std::to_string(params.get_nlevels()) + " number of given reuse-coarse-matrix: " + std::to_string(it->second.size()) + "'");
2462 return false;
2463 }
2464 }
2465 return true;
2466 }
2467
2468 inline bool parse_fpl_reuse_coarse_basis(const Dist::Comm &comm, SimpleArgParser &args, FROSchParameterList &params)
2469 {
2470 auto it = args.query("frosch-reuse-coarse-basis");
2471 if(it != nullptr)
2472 {
2473 if(it->second.size() == 1u)
2474 {
2475 bool reuse_cb;
2476 std::string type = it->second.front();
2477 std::transform(type.begin(), type.end(), type.begin(), [](unsigned char c){return std::tolower(c);});
2478
2479 if(!type.compare("true"))
2480 {
2481 reuse_cb = true;
2482 }else if(!type.compare("false"))
2483 {
2484 reuse_cb = false;
2485 }else
2486 {
2487 comm.print("ERROR: unknown reuse-coarse-basis '" + it->second.front() + "'");
2488 return false;
2489 }
2490 params.set_reuse_coarse_basis(reuse_cb);
2491 }else if(int(it->second.size()) == params.get_nlevels())
2492 {
2493 std::vector<bool> reuse_cb;
2494 reuse_cb.reserve(params.get_nlevels());
2495 for(const auto& t : it->second)
2496 {
2497 std::string type = t;
2498 std::transform(type.begin(), type.end(), type.begin(), [](unsigned char c){return std::tolower(c);});
2499
2500 if(!type.compare("true"))
2501 {
2502 reuse_cb.push_back(true);
2503 }else if(!type.compare("false"))
2504 {
2505 reuse_cb.push_back(false);
2506 }else
2507 {
2508 comm.print("ERROR: unknown reuse-coarse-basis type '" + t + "'");
2509 return false;
2510 }
2511 }
2512 params.set_reuse_coarse_basis(reuse_cb);
2513 }else
2514 {
2515 comm.print("ERROR: not matching number of reuse-coarse-basis: nlevels: " + std::to_string(params.get_nlevels()) + " number of given reuse-coarse-basis: " + std::to_string(it->second.size()) + "'");
2516 return false;
2517 }
2518 }
2519 return true;
2520 }
2521
2522 inline bool parse_fpl_phi_dt(const Dist::Comm &comm, SimpleArgParser &args, FROSchParameterList &params)
2523 {
2524 auto it = args.query("frosch-phi-dropping-threshold");
2525 if(it != nullptr)
2526 {
2527 if(it->second.size() == 1u)
2528 {
2529 double phi_dt = std::stod(it->second.front());
2530 if (phi_dt < 0.)
2531 {
2532 /* __TODO:__ Stimmt das? */
2533 comm.print("ERROR: unknown dropping threshold smaller than 0.: '" + it->second.front() + "'");
2534 return false;
2535 }
2536 params.set_phi_dropping_threshold(phi_dt);
2537 }else if(int(it->second.size()) == params.get_nlevels())
2538 {
2539 std::vector<double> phi_dt;
2540 phi_dt.reserve(params.get_nlevels());
2541 for(const auto& t : it->second)
2542 {
2543 double dt = std::stod(t);
2544 if (dt < 0.)
2545 {
2546 /* __TODO:__ Stimmt das? */
2547 comm.print("ERROR: unknown partition imblance smaller than 0.: '" + t + "'");
2548 return false;
2549 }
2550 phi_dt.push_back(dt);
2551 }
2552 params.set_phi_dropping_threshold(phi_dt);
2553 }else
2554 {
2555 comm.print("ERROR: not matching number of phi-dropping-threshold: nlevels: " + std::to_string(params.get_nlevels()) + " number of given phi-dropping-threshold: " + std::to_string(it->second.size()) + "'");
2556 return false;
2557 }
2558 }
2559 return true;
2560 }
2561
2562 bool FROSchParameterList::parse_args(SimpleArgParser& args)
2563 {
2564 {
2565 auto it = args.query("frosch-plist");
2566 if(it != nullptr)
2567 {
2568 if(it->second.size() == 1u)
2569 {
2570 read_from_xml_file(it->second.front());
2571 /* given parameter list wins */
2572 return true;
2573 }else
2574 {
2575 _comm.print("ERROR: multiple XML files are given");
2576 return false;
2577 }
2578 }
2579 }
2580
2581 if(!parse_fpl_print(_comm, args, *this))
2582 {
2583 return false;
2584 }
2585 if(!parse_fpl_verbose(_comm, args, *this))
2586 {
2587 return false;
2588 }
2589 if(!parse_fpl_timer(_comm, args, *this))
2590 {
2591 return false;
2592 }
2593 if(!parse_fpl_levels(_comm, args, *this))
2594 {
2595 return false;
2596 }
2597 if(!parse_fpl_dim(_comm, args, *this))
2598 {
2599 return false;
2600 }
2601 if(!parse_fpl_overlaps(_comm, args, *this))
2602 {
2603 return false;
2604 }
2605 if(!parse_fpl_subreg(_comm, args, *this))
2606 {
2607 return false;
2608 }
2609 if(!parse_fpl_gsteps(_comm, args, *this))
2610 {
2611 return false;
2612 }
2613 if(!parse_fpl_parti_type(_comm, args, *this))
2614 {
2615 return false;
2616 }
2617 if(!parse_fpl_parti_approach(_comm, args, *this))
2618 {
2619 return false;
2620 }
2621 if(!parse_fpl_parti_imbl(_comm, args, *this))
2622 {
2623 return false;
2624 }
2625 if(!parse_fpl_ipou_pres(_comm, args, *this))
2626 {
2627 return false;
2628 }
2629 if(!parse_fpl_ipou_velo(_comm, args, *this))
2630 {
2631 return false;
2632 }
2633 if(!parse_fpl_solver_ao(_comm, args, *this))
2634 {
2635 return false;
2636 }
2637 if(!parse_fpl_solver_ext(_comm, args, *this))
2638 {
2639 return false;
2640 }
2641 if(!parse_fpl_solver_coarse(_comm, args, *this))
2642 {
2643 return false;
2644 }
2645 if(!parse_fpl_precond(_comm, args, *this))
2646 {
2647 return false;
2648 }
2649 if(!parse_fpl_cspace_pres(_comm, args, *this))
2650 {
2651 return false;
2652 }
2653 if(!parse_fpl_cspace_velo(_comm, args, *this))
2654 {
2655 return false;
2656 }
2657 if(!parse_fpl_exclude_velo_pres(_comm, args, *this))
2658 {
2659 return false;
2660 }
2661 if(!parse_fpl_exclude_pres_velo(_comm, args, *this))
2662 {
2663 return false;
2664 }
2665 if(!parse_fpl_combine(_comm, args, *this))
2666 {
2667 return false;
2668 }
2669 if(!parse_fpl_combine_lvl(_comm, args, *this))
2670 {
2671 return false;
2672 }
2673 if (!parse_fpl_reuse_sf_ao(_comm, args, *this))
2674 {
2675 return false;
2676 }
2677 if (!parse_fpl_reuse_sf_cm(_comm, args, *this))
2678 {
2679 return false;
2680 }
2681 if (!parse_fpl_reuse_coarse_matrix(_comm, args, *this))
2682 {
2683 return false;
2684 }
2685 if (!parse_fpl_reuse_coarse_basis(_comm, args, *this))
2686 {
2687 return false;
2688 }
2689 if (!parse_fpl_phi_dt(_comm, args, *this))
2690 {
2691 return false;
2692 }
2693
2694 return true;
2695 }
2696
2700 class TpetraCoreBase
2701 {
2702 public:
2703 // some typedefs
2704 using UN = unsigned;
2705 using SC = double;
2706 using LO = int;
2707 using GO = FROSch::DefaultGlobalOrdinal;
2708 using NO = Tpetra::KokkosClassic::DefaultNode::DefaultNodeType;
2709
2710 using AGS = Teuchos::Array<GO>::size_type;
2711 using ARCPS = Teuchos::ArrayRCP<std::size_t>::size_type;
2712 using TGS = Tpetra::global_size_t;
2713 using VGS = std::vector<GO>::size_type;
2714
2715 //----------------------------------------
2716 Teuchos::RCP<const Teuchos::Comm<int>> _comm;
2717
2718 const LO _num_owned_dofs, _num_owned_nonzeros;
2719 const GO _my_dof_offset, _num_global_dofs;
2720 const int _dpe_pres;
2722 // note: TGS = Tpetra::global_size_t should usually have the same size as std::int64_t
2723 Teuchos::ArrayRCP<std::int64_t> _row_ptr;
2725 Teuchos::ArrayRCP<TGS> _num_nze;
2727 // note: GO = FROSch::DefaultGlobalOrdinal should usually have the same size as std::int32_t
2728 Teuchos::ArrayRCP<std::int32_t> _col_idx;
2730 Teuchos::ArrayRCP<SC> _mat_val;
2731
2733 Teuchos::RCP<const Tpetra::Map<LO, GO, NO>> _map;
2734
2736 Teuchos::RCP<Tpetra::MultiVector<SC, LO,GO,NO>> _vec_def, _vec_cor;
2737 Teuchos::RCP<Tpetra::CrsMatrix<SC,LO,GO,NO>> _matrix;
2738
2740 Teuchos::RCP<Teuchos::ParameterList> _params;
2741
2742 bool _use_timer;
2743 Teuchos::RCP<Teuchos::StackedTimer> _stacked_timer;
2744
2745 explicit TpetraCoreBase(const void* comm, Index num_global_dofs, Index my_dof_offset,
2746 Index num_owned_dofs, Index num_nonzeros, const Trilinos::FROSchParameterList & params) :
2747 _comm(Teuchos::rcp(new Teuchos::MpiComm<int>(*reinterpret_cast<const MPI_Comm*>(comm)))),
2748 _num_owned_dofs(static_cast<LO>(num_owned_dofs)),
2749 _num_owned_nonzeros(static_cast<LO>(num_nonzeros)),
2750 _my_dof_offset(static_cast<GO>(my_dof_offset)),
2751 _num_global_dofs(static_cast<GO>(num_global_dofs)),
2752 _dpe_pres(params.get_dpe_pres()),
2753 _row_ptr(static_cast<ARCPS>(_num_owned_dofs+1u), static_cast<std::int64_t>(0)),
2754 _num_nze(static_cast<ARCPS>(_num_owned_dofs), static_cast<TGS>(0)),
2755 _col_idx(static_cast<ARCPS>(_num_owned_nonzeros), static_cast<std::int32_t>(0)),
2756 _mat_val(static_cast<ARCPS>(_num_owned_nonzeros), static_cast<SC>(0.)),
2757 _map(nullptr),
2758 _vec_def(nullptr),
2759 _vec_cor(nullptr),
2760 _matrix(nullptr),
2761 _params(Teuchos::null),
2762 _use_timer(params.get_use_timer()),
2763 _stacked_timer(Teuchos::rcp(new Teuchos::StackedTimer("TpetraCoreBase")))
2764 {
2765 void* pcore = params.get_core();
2766 if(pcore != nullptr)
2767 {
2768 _params = Teuchos::rcp(reinterpret_cast<Teuchos::ParameterList*>(pcore), false);
2769 }
2770
2771 if(_use_timer)
2772 {
2773 _comm->barrier();
2774 Teuchos::TimeMonitor::setStackedTimer(_stacked_timer);
2775 }
2776 }
2777
2778 virtual ~TpetraCoreBase()
2779 {
2780 if(_use_timer)
2781 {
2782 _comm->barrier();
2783 _stacked_timer->stop("TpetraCoreBase");
2784 Teuchos::StackedTimer::OutputOptions options;
2785 options.output_fraction = options.output_histogram = options.output_minmax = true;
2786 Teuchos::RCP<Teuchos::FancyOStream> out = Teuchos::VerboseObjectBase::getDefaultOStream();
2787 _stacked_timer->report(*out, _comm, options);
2788
2789 }
2790 }
2791
2792 virtual void init_symbolic() = 0;
2793
2794 virtual void init_numeric() = 0;
2795
2796 virtual std::int64_t* get_row_ptr()
2797 {
2798 return this->_row_ptr.get();
2799 }
2800
2801 virtual std::int32_t* get_col_idx()
2802 {
2803 return this->_col_idx.get();
2804 }
2805
2806 virtual double* get_mat_val()
2807 {
2808 return this->_mat_val.get();
2809 }
2810
2811 virtual double* get_vec_def()
2812 {
2813 return this->_vec_def->getDataNonConst(0).get();
2814 }
2815
2816 virtual double* get_vec_cor()
2817 {
2818 return this->_vec_cor->getDataNonConst(0).get();
2819 }
2820
2821 virtual void format_vec_cor()
2822 {
2823 return this->_vec_cor->putScalar(0.0);
2824 }
2825 }; // class TpetraCore
2826
2827 void set_parameter_list(void* core, const FROSchParameterList& params)
2828 {
2829 void* pcore = params.get_core();
2830 if(pcore == nullptr)
2831 {
2832 reinterpret_cast<TpetraCoreBase*>(core)->_params = Teuchos::null;
2833 }
2834 else
2835 {
2836 reinterpret_cast<TpetraCoreBase*>(core)->_params
2837 = Teuchos::rcp(reinterpret_cast<Teuchos::ParameterList*>(pcore), false);
2838 }
2839 }
2840
2841 void destroy_core(void* core)
2842 {
2843 delete reinterpret_cast<TpetraCoreBase*>(core);
2844 }
2845
2846 void init_symbolic(void* core)
2847 {
2848 reinterpret_cast<TpetraCoreBase*>(core)->init_symbolic();
2849 }
2850
2851 void init_numeric(void* core)
2852 {
2853 reinterpret_cast<TpetraCoreBase*>(core)->init_numeric();
2854 }
2855
2856 std::int64_t* get_row_ptr(void* core)
2857 {
2858 return reinterpret_cast<TpetraCoreBase*>(core)->get_row_ptr();
2859 }
2860
2861 std::int32_t* get_col_idx(void* core)
2862 {
2863 return reinterpret_cast<TpetraCoreBase*>(core)->get_col_idx();
2864 }
2865
2866 double* get_mat_val(void* core)
2867 {
2868 return reinterpret_cast<TpetraCoreBase*>(core)->get_mat_val();
2869 }
2870
2871 double* get_vec_def(void* core)
2872 {
2873 return reinterpret_cast<TpetraCoreBase*>(core)->get_vec_def();
2874 }
2875
2876 double* get_vec_cor(void* core)
2877 {
2878 return reinterpret_cast<TpetraCoreBase*>(core)->get_vec_cor();
2879 }
2880
2881 void format_vec_cor(void* core)
2882 {
2883 return reinterpret_cast<TpetraCoreBase*>(core)->format_vec_cor();
2884 }
2885
2889 class TpetraCoreScalar : public TpetraCoreBase
2890 {
2891 public:
2892 // some typedefs
2893 typedef TpetraCoreBase BaseClass;
2894 using UN = unsigned;
2895 using SC = double;
2896 using LO = int;
2897 using GO = FROSch::DefaultGlobalOrdinal;
2898 using NO = Tpetra::KokkosClassic::DefaultNode::DefaultNodeType;
2899
2900 using AGS = Teuchos::Array<GO>::size_type;
2901 using ARCPS = Teuchos::ArrayRCP<std::size_t>::size_type;
2902 using TGS = Tpetra::global_size_t;
2903 using VGS = std::vector<GO>::size_type;
2904
2905 explicit TpetraCoreScalar(const void* comm, Index num_global_dofs, Index my_dof_offset,
2906 Index num_owned_dofs, Index num_nonzeros, const Trilinos::FROSchParameterList & params) :
2907 TpetraCoreBase(comm, num_global_dofs, my_dof_offset, num_owned_dofs, num_nonzeros, params)
2908 {
2909 }
2910
2911 virtual void init_symbolic() override
2912 {
2913 if(_use_timer)
2914 {
2915 _stacked_timer->start("TpetraCoreScalar::init_core");
2916 _stacked_timer->start("TpetraCoreScalar::init_core::build_map");
2917 }
2918
2920 _num_nze = Teuchos::ArrayRCP<TGS>(static_cast<ARCPS>(_num_owned_dofs), static_cast<TGS>(0));
2921
2923 for(LO i(0); i < _num_owned_dofs; ++i)
2924 _num_nze[i] = static_cast<TGS>(_row_ptr[i+1] - _row_ptr[i]);
2925
2927
2933 std::vector<GO> element_list_vec(static_cast<VGS>(_num_owned_dofs));
2934 for(IndexType i = 0; i < static_cast<IndexType>(_num_owned_dofs); ++i)
2935 element_list_vec.at(i) = static_cast<GO>(_my_dof_offset + i);
2936
2937 GO num_total_dofs = 0;
2938 GO my_dof_offset_tmp = static_cast<GO>(_num_owned_dofs);
2939 Teuchos::reduceAll(*_comm, Teuchos::REDUCE_SUM, 1, &my_dof_offset_tmp, &num_total_dofs);
2940
2941 const GO _index_base = 0;
2942 const Teuchos::ArrayView<const GO> element_list = Teuchos::arrayViewFromVector(element_list_vec);
2943 _map = Teuchos::rcp(new Tpetra::Map(static_cast<TGS>(num_total_dofs), element_list, _index_base, _comm));
2944
2945 if(_use_timer)
2946 {
2947 _stacked_timer->stop("TpetraCoreScalar::init_core::build_map");
2948 _stacked_timer->start("TpetraCoreScalar::init_core::build_system_matrix");
2949 }
2950
2952 auto crsgraph = Teuchos::rcp(new Tpetra::CrsGraph<LO, GO, NO >(_map, _num_nze()));
2953 const TGS max_nnz_per_row = *std::max_element(_num_nze.begin(), _num_nze.end());
2954 std::vector<GO> col_buffer(max_nnz_per_row); // Reusable buffer
2955
2957
2961 AGS entries_start = 0;
2962 for(GO i = 0; i < static_cast<GO>(_num_owned_dofs); ++i)
2963 {
2964 const auto num_entries = static_cast<AGS>(_num_nze[i]);
2965
2966 for(Teuchos::Array<GO>::size_type j = 0; j < num_entries; ++j)
2967 col_buffer[j] = static_cast<GO>(this->_col_idx[entries_start + j]);
2968
2969 crsgraph->insertGlobalIndices(static_cast<GO>(_my_dof_offset + i), int(num_entries), col_buffer.data());
2970 entries_start += num_entries;
2971 }
2972 crsgraph->fillComplete();
2973 _matrix = Teuchos::rcp(new Tpetra::CrsMatrix<SC, LO, GO, NO>(crsgraph));
2974 if(_use_timer)
2975 {
2976 _stacked_timer->stop("TpetraCoreScalar::init_core::build_system_matrix");
2977 }
2978
2979 // create vectors
2980 const LO num_vectors = 1;
2981 _vec_def = Teuchos::rcp(new Tpetra::MultiVector<SC, LO, GO, NO>(_map, num_vectors));
2982 _vec_cor = Teuchos::rcp(new Tpetra::MultiVector<SC, LO, GO, NO>(_map, num_vectors));
2983
2984 // set entries for TwoLevelBlockPreconditionier
2985 if (! sublist(sublist(_params, "Preconditioner Types"), "FROSch")->get("FROSch Preconditioner Type","TwoLevelPreconditioner").compare("TwoLevelBlockPreconditioner"))
2986 {
2987 const int NumberOfBlocks = 1;
2988 Teuchos::ArrayRCP<Teuchos::RCP<Xpetra::Map<LO, GO, NO>>> XRepeatedMaps(NumberOfBlocks);
2989 auto xcrsgraph = Teuchos::rcp(new const Xpetra::TpetraCrsGraph<LO, GO, NO>(Teuchos::rcp_const_cast<Tpetra::CrsGraph<LO, GO, NO >>(_matrix->getCrsGraph ())));
2990 XRepeatedMaps[0] = FROSch::BuildRepeatedMapNonConst<LO, GO, NO>(xcrsgraph);
2991 sublist(sublist(_params, "Preconditioner Types"), "FROSch")->set("Repeated Map Vector", XRepeatedMaps);
2992
2993 Teuchos::ArrayRCP<FROSch::DofOrdering> dofOrderings(NumberOfBlocks);
2994 dofOrderings[0] = FROSch::NodeWise;
2995 sublist(sublist(_params, "Preconditioner Types"), "FROSch")->set("DofOrdering Vector", dofOrderings);
2996
2997 Teuchos::ArrayRCP<unsigned> dofsPerNodeVector(NumberOfBlocks);
2998 dofsPerNodeVector[0] = _dpe_pres;
2999 sublist(sublist(_params, "Preconditioner Types"), "FROSch")->set("DofsPerNode Vector", dofsPerNodeVector);
3000
3001 /* NULL SPACE */
3002 Teuchos::ArrayRCP<Teuchos::RCP<Xpetra::MultiVector<SC,LO,GO,NO>>> nullSpaceBasis(NumberOfBlocks);
3003 nullSpaceBasis[0] = Xpetra::MultiVectorFactory<SC,LO,GO,NO>::Build(XRepeatedMaps[0].getConst(), 1);
3004 unsigned i = 0;
3005 for (unsigned j = 0; j < XRepeatedMaps[i]->getLocalNumElements(); j += dofsPerNodeVector[i])
3006 {
3007 nullSpaceBasis[0]->getDataNonConst(i)[XRepeatedMaps[i]->getLocalElement(XRepeatedMaps[i]->getGlobalElement(j))] = Teuchos::ScalarTraits<SC>::one();
3008 }
3009 sublist(sublist(_params, "Preconditioner Types"), "FROSch")->set("Null Space Vector", nullSpaceBasis);
3010 }
3011 else if(! sublist(sublist(_params, "Preconditioner Types"), "FROSch")->get("FROSch Preconditioner Type","TwoLevelPreconditioner").compare("TwoLevelPreconditioner"))
3012 {
3013 /* NULL SPACE */
3014 auto xcrsgraph = Teuchos::rcp(new const Xpetra::TpetraCrsGraph<LO, GO, NO>(Teuchos::rcp_const_cast<Tpetra::CrsGraph<LO, GO, NO >>(_matrix->getCrsGraph ())));
3015 Teuchos::RCP<Xpetra::Map<LO, GO, NO>> XRepeatedMap = FROSch::BuildRepeatedMapNonConst<LO, GO, NO>(xcrsgraph);
3016 sublist(sublist(_params, "Preconditioner Types"), "FROSch")->set("Repeated Map", XRepeatedMap);
3017
3018 Teuchos::RCP<Xpetra::MultiVector<SC,LO,GO,NO>> nullSpaceBasis = Xpetra::MultiVectorFactory<SC,LO,GO,NO>::Build(XRepeatedMap.getConst(), 1);
3019 for (unsigned j = 0; j < XRepeatedMap->getLocalNumElements(); j += _dpe_pres)
3020 {
3021 nullSpaceBasis->getDataNonConst(0)[XRepeatedMap->getLocalElement(XRepeatedMap->getGlobalElement(j))] = Teuchos::ScalarTraits<SC>::one();
3022 }
3023 sublist(sublist(_params, "Preconditioner Types"), "FROSch")->set("Null Space", nullSpaceBasis);
3024 }
3025
3026 if(_use_timer)
3027 {
3028 _stacked_timer->stop("TpetraCoreScalar::init_core");
3029 }
3030 }
3031
3032 virtual void init_numeric() override
3033 {
3034 if(this->_use_timer)
3035 {
3036 this->_stacked_timer->start("TpetraCoreScalar::SetValues");
3037 }
3038
3039 Teuchos::ArrayView<const LO> Indices;
3040
3041 _matrix->resumeFill();
3042 AGS entries_start = 0;
3043
3044 for(LO i(0); i < _num_owned_dofs; ++i)
3045 {
3046 const auto num_entries = static_cast<Teuchos::Array<SC>::size_type>(_num_nze[i]);
3047 // copy vals to Teuchos::Array
3048 Teuchos::Array<SC> vals_ar(num_entries);
3049 for(Teuchos::Array<SC>::size_type j = 0; j < num_entries; ++j)
3050 vals_ar[j] = _mat_val[static_cast<VGS>(entries_start + j)];
3051
3052 Teuchos::Array<GO> cols(num_entries);
3053 for(Teuchos::Array<GO>::size_type j = 0; j < num_entries; ++j)
3054 cols[j] = static_cast<GO>(_col_idx[entries_start + j]);
3055
3056 _matrix->replaceGlobalValues(_my_dof_offset + static_cast<GO>(i),
3057 cols.view(0, static_cast<AGS>(num_entries)),
3058 vals_ar.view(0, static_cast<Teuchos::Array<SC>::size_type>(num_entries)));
3059
3060 entries_start += num_entries;
3061 }
3062
3063 _matrix->fillComplete();
3064
3065 if(_use_timer)
3066 {
3067 _stacked_timer->stop("TpetraCoreScalar::SetValues");
3068 }
3069 }
3070 }; // class TpetraCoreScalar
3071
3072 void* create_core_scalar(const void* comm, Index num_global_dofs, Index my_dof_offset,
3073 Index num_owned_dofs, Index num_nonzeros, const FROSchParameterList& params)
3074 {
3075 return new TpetraCoreScalar(comm, num_global_dofs, my_dof_offset, num_owned_dofs, num_nonzeros, params);
3076 }
3077
3083 class TpetraCoreStokes : public TpetraCoreBase
3084 {
3085 public:
3086 // some typedefs
3087 typedef TpetraCoreBase BaseClass;
3088 using UN = unsigned;
3089 using SC = double;
3090 using LO = int;
3091 using GO = FROSch::DefaultGlobalOrdinal;
3092 using NO = Tpetra::KokkosClassic::DefaultNode::DefaultNodeType;
3093
3094 using AGS = Teuchos::Array<GO>::size_type;
3095 using ARCPS = Teuchos::ArrayRCP<std::size_t>::size_type;
3096 using TGS = Tpetra::global_size_t;
3097 using VGS = std::vector<GO>::size_type;
3098
3099 //----------------------------------------
3100 const int _dpe_velo;
3101 // number of velocity/pressure dofs owned by this process
3102 const IndexType _num_owned_velo_dofs;
3103 const IndexType _num_owned_pres_dofs; // <- given to constructor
3104 // index of first velocity/pressure dof owned by this process
3105 const IndexType _first_owned_velo_dof;
3106 const IndexType _first_owned_pres_dof;
3107 // total number of velocity/pressure dofs
3108 const IndexType _num_global_velo_dofs;
3109 const IndexType _num_global_pres_dofs;
3111 Teuchos::ArrayRCP<std::int32_t> _col_idx_mapped;
3112
3113 Teuchos::ArrayRCP<GO> _all_global_dof_offset;
3114 Teuchos::ArrayRCP<IndexType> _all_num_owned_velo_dofs, _all_num_owned_pres_dofs;
3115 Teuchos::ArrayRCP<IndexType> _all_first_owned_velo_dof, _all_first_owned_pres_dof;
3116
3117 explicit TpetraCoreStokes(const void* comm,
3118 Index num_global_dofs, Index my_dof_offset,
3119 Index num_owned_dofs, Index num_nonzeros,
3120 IndexType num_owned_velo_dofs, IndexType num_owned_pres_dofs,
3121 IndexType first_owned_velo_dof, IndexType first_owned_pres_dof,
3122 IndexType num_global_velo_dofs, IndexType num_global_pres_dofs,
3123 const Trilinos::FROSchParameterList & params) :
3124 BaseClass(comm, num_global_dofs, my_dof_offset, num_owned_dofs, num_nonzeros, params),
3125 _dpe_velo(params.get_dpe_velo()),
3126 _num_owned_velo_dofs(num_owned_velo_dofs),
3127 _num_owned_pres_dofs(num_owned_pres_dofs),
3128 _first_owned_velo_dof(first_owned_velo_dof),
3129 _first_owned_pres_dof(first_owned_pres_dof),
3130 _num_global_velo_dofs(num_global_velo_dofs),
3131 _num_global_pres_dofs(num_global_pres_dofs),
3132 _col_idx_mapped(static_cast<ARCPS>(_num_owned_nonzeros), 0),
3133 _all_global_dof_offset(static_cast<ARCPS>(this->_comm->getSize())),
3134 _all_num_owned_velo_dofs(static_cast<ARCPS>(this->_comm->getSize())),
3135 _all_num_owned_pres_dofs(static_cast<ARCPS>(this->_comm->getSize())),
3136 _all_first_owned_velo_dof(static_cast<ARCPS>(this->_comm->getSize())),
3137 _all_first_owned_pres_dof(static_cast<ARCPS>(this->_comm->getSize()))
3138 {
3139 }
3140
3141 virtual void init_symbolic() override
3142 {
3143 if(_use_timer)
3144 {
3145 _comm->barrier();
3146 Teuchos::TimeMonitor::setStackedTimer(this->_stacked_timer);
3147 _stacked_timer->start("TpetraCoreStokes::init_symbolic");
3148 }
3149
3150 // ceate a deep copy of the original column indices array, as it will be mapped later on
3151 this->_col_idx_mapped.assign(this->_col_idx.begin(), this->_col_idx.end());
3152
3153 if(_use_timer)
3154 {
3155 _stacked_timer->start("TpetraCoreStokes::init_symbolic::gather");
3156 }
3158 Teuchos::gatherAll(*(_comm), static_cast<int>(1), &_my_dof_offset, static_cast<int>(1*_comm->getSize()), _all_global_dof_offset.getRawPtr());
3159 Teuchos::gatherAll(*(_comm), static_cast<int>(1), &_num_owned_velo_dofs, static_cast<int>(1*_comm->getSize()), _all_num_owned_velo_dofs.getRawPtr());
3160 Teuchos::gatherAll(*(_comm), static_cast<int>(1), &_num_owned_pres_dofs, static_cast<int>(1*_comm->getSize()), _all_num_owned_pres_dofs.getRawPtr());
3161 Teuchos::gatherAll(*(_comm), static_cast<int>(1), &_first_owned_velo_dof, static_cast<int>(1*_comm->getSize()), _all_first_owned_velo_dof.getRawPtr());
3162 Teuchos::gatherAll(*(_comm), static_cast<int>(1), &_first_owned_pres_dof, static_cast<int>(1*_comm->getSize()), _all_first_owned_pres_dof.getRawPtr());
3163 if(_use_timer)
3164 {
3165 _stacked_timer->stop("TpetraCoreStokes::init_symbolic::gather");
3166 _stacked_timer->start("TpetraCoreStokes::init_symbolic::build_map");
3167 }
3168
3170 _num_nze = Teuchos::ArrayRCP<TGS>(static_cast<ARCPS>(_num_owned_dofs), static_cast<TGS>(0));
3171
3173 for(LO i(0); i < _num_owned_dofs; ++i)
3174 _num_nze[i] = static_cast<TGS>(this->_row_ptr[i+1] - this->_row_ptr[i]);
3175
3189 const GO _index_base = 0;
3190 std::vector<GO> element_list_vec(static_cast<VGS>(_num_owned_dofs));
3191 for(IndexType i = 0; i < static_cast<IndexType>(_num_owned_velo_dofs); ++i)
3192 element_list_vec.at(i) = static_cast<GO>(_first_owned_velo_dof + i);
3193 for(IndexType i = 0; i < static_cast<IndexType>(_num_owned_pres_dofs); ++i)
3194 element_list_vec.at(_num_owned_velo_dofs + i) = static_cast<GO>(_num_global_velo_dofs + _first_owned_pres_dof + i);
3195 const Teuchos::ArrayView<const GO> element_list = Teuchos::arrayViewFromVector(element_list_vec);
3196
3197 _map = Teuchos::rcp(new Tpetra::Map(static_cast<TGS>(_num_global_velo_dofs + _num_global_pres_dofs), element_list, _index_base, _comm));
3198
3199 if(_use_timer)
3200 {
3201 _stacked_timer->stop("TpetraCoreStokes::init_symbolic::build_map");
3202 _stacked_timer->start("TpetraCoreStokes::init_symbolic::build_repeated_maps");
3203 }
3204 const int NumberOfBlocks = 2;
3205
3206 // Pre-compute rank boundaries for O(1) access
3207 std::vector<GO> rank_start(_comm->getSize());
3208 std::vector<GO> rank_velo_end(_comm->getSize());
3209 std::vector<GO> rank_total_end(_comm->getSize());
3210 for (int rank = 0; rank < _comm->getSize(); ++rank)
3211 {
3212 rank_start[rank] = _all_global_dof_offset[rank];
3213 rank_velo_end[rank] = _all_global_dof_offset[rank] + GO(_all_num_owned_velo_dofs[rank]);
3214 rank_total_end[rank] = _all_global_dof_offset[rank] + GO(_all_num_owned_velo_dofs[rank]) + GO(_all_num_owned_pres_dofs[rank]);
3215 }
3216
3218 {
3219 if(_use_timer)
3220 {
3221 _stacked_timer->start("TpetraCoreStokes::init_symbolic::build_repeated_maps::velocity_graph");
3222 }
3223
3224 std::vector<GO> element_list_vec_velo(static_cast<VGS>(_num_owned_velo_dofs));
3225 for(IndexType i = 0; i < static_cast<IndexType>(_num_owned_velo_dofs); ++i)
3226 element_list_vec_velo.at(i) = static_cast<GO>(_first_owned_velo_dof + i);
3227 const Teuchos::ArrayView<const GO> element_list_velo = Teuchos::arrayViewFromVector(element_list_vec_velo);
3228 Teuchos::RCP<const Tpetra::Map<LO, GO, NO>> _map_velo = Teuchos::rcp(new Tpetra::Map(static_cast<TGS>(_num_global_velo_dofs), element_list_velo, _index_base, _comm));
3229
3230 auto crsgraph = Teuchos::rcp(new Tpetra::CrsGraph<LO, GO, NO >(_map_velo, _num_nze.view(0, static_cast<AGS>(_num_owned_velo_dofs))));
3231 const TGS max_nnz_per_row = *std::max_element(_num_nze.begin(), _num_nze.end());
3232 std::vector<GO> col_buffer(max_nnz_per_row); // Reusable buffer
3233
3235 for(IndexType row(0); row < _num_owned_velo_dofs; ++row)
3236 {
3237 Teuchos::Array<GO> cols(static_cast<GO>(_num_nze[static_cast<ARCPS>(row)]));
3238 GO counter = 0;
3239
3240 // Start scanning from rank 0 for each row
3241 int current_rank = 0;
3242
3243 for(auto colidx = this->_row_ptr[row]; colidx < this->_row_ptr[row+1]; ++colidx)
3244 {
3245 GO col = static_cast<GO>(this->_col_idx_mapped[colidx]);
3246
3247 // Advance rank until we find one that might contain this column
3248 while(current_rank < int(_comm->getSize()) && col >= rank_total_end[current_rank])
3249 {
3250 current_rank++;
3251 }
3252
3253 // Check if current rank contains this column
3254 if(current_rank < int(_comm->getSize()) && col >= rank_start[current_rank])
3255 {
3256 if(col < rank_velo_end[current_rank])
3257 {
3258 // Velocity DOF
3259 GO mapped_dof = GO(_all_first_owned_velo_dof[current_rank]) + GO(col - rank_start[current_rank]);
3260 cols[counter++] = mapped_dof;
3261 }
3262 }
3263 // Note: current_rank is NOT reset - it continues from where it left off
3264 }
3265
3266 crsgraph->insertGlobalIndices(static_cast<GO>(_first_owned_velo_dof + row),
3267 cols.view(0, static_cast<AGS>(counter)));
3268 }
3269 crsgraph->fillComplete();
3270
3271 if(_use_timer)
3272 {
3273 _stacked_timer->stop("TpetraCoreStokes::init_symbolic::build_repeated_maps::velocity_graph");
3274 }
3275
3276 Teuchos::ArrayRCP<Teuchos::RCP<Tpetra::Map<LO, GO, NO>>> RepeatedMaps(NumberOfBlocks);
3277 Teuchos::ArrayRCP<Teuchos::RCP<Xpetra::Map<LO, GO, NO>>> XRepeatedMaps(NumberOfBlocks);
3278
3279 const GO _index_base_tmp = 0;
3280 /* pressure map */
3281 std::vector<GO> element_list_vec_tmp(static_cast<VGS>(_num_owned_pres_dofs));
3282 for(IndexType i = 0; i < static_cast<IndexType>(_num_owned_pres_dofs); ++i)
3283 element_list_vec_tmp.at(i) = static_cast<GO>(_first_owned_pres_dof + i);
3284 const Teuchos::ArrayView<const GO> element_list_tmp = Teuchos::arrayViewFromVector(element_list_vec_tmp);
3285 RepeatedMaps[1] = Teuchos::rcp(new Tpetra::Map(static_cast<TGS>(_num_owned_pres_dofs), element_list_tmp, _index_base_tmp, _comm));
3286
3287 auto xcrsgraph = Teuchos::rcp(new Xpetra::TpetraCrsGraph<LO, GO, NO>(crsgraph));
3288 XRepeatedMaps[0] = FROSch::BuildRepeatedMapNonConst<LO, GO, NO>(xcrsgraph);
3289 // __ATTENTION:__ This works only for a discontinuous pressure
3290 XRepeatedMaps[1] = Teuchos::rcp(new Xpetra::TpetraMap<LO, GO, NO>(RepeatedMaps[1]));
3291
3292 sublist(sublist(_params, "Preconditioner Types"), "FROSch")->set("Repeated Map Vector", XRepeatedMaps);
3293
3295 Teuchos::ArrayRCP<FROSch::DofOrdering> dofOrderings(NumberOfBlocks);
3296 dofOrderings[0] = FROSch::NodeWise;
3297 dofOrderings[1] = FROSch::NodeWise;
3298 sublist(sublist(_params, "Preconditioner Types"), "FROSch")->set("DofOrdering Vector", dofOrderings);
3299
3300 Teuchos::ArrayRCP<unsigned> dofsPerNodeVector(NumberOfBlocks);
3301 dofsPerNodeVector[0] = _dpe_velo;
3302 dofsPerNodeVector[1] = _dpe_pres;
3303 sublist(sublist(_params, "Preconditioner Types"), "FROSch")->set("DofsPerNode Vector", dofsPerNodeVector);
3304
3305 /* Null Space */
3306 Teuchos::ArrayRCP<Teuchos::RCP<Xpetra::MultiVector<SC,LO,GO,NO>>> nullSpaceBasis(NumberOfBlocks);
3307 Teuchos::ArrayRCP<Teuchos::RCP<const Xpetra::Map<LO, GO, NO>>> XRepeatedMapsTmp(NumberOfBlocks);
3308 XRepeatedMapsTmp[0] = XRepeatedMaps[0].getConst();
3309 XRepeatedMapsTmp[1] = XRepeatedMaps[1].getConst();
3310 Teuchos::RCP<const Xpetra::Map<LO, GO, NO>> repeatedMap = FROSch::MergeMaps(XRepeatedMapsTmp);
3311
3312 nullSpaceBasis[0] = Xpetra::MultiVectorFactory<SC,LO,GO,NO>::Build(repeatedMap, dofsPerNodeVector[0]);
3313 nullSpaceBasis[1] = Xpetra::MultiVectorFactory<SC,LO,GO,NO>::Build(repeatedMap, 1);
3314 for(unsigned i = 0; i < dofsPerNodeVector[0]; i++)
3315 for (unsigned j = i; j < XRepeatedMaps[0]->getLocalNumElements(); j += dofsPerNodeVector[0]) {
3316 nullSpaceBasis[0]->getDataNonConst(i)[XRepeatedMaps[0]->getLocalElement(XRepeatedMaps[0]->getGlobalElement(j))] = Teuchos::ScalarTraits<SC>::one();
3317 }
3318 LO offset = LO(XRepeatedMaps[0]->getLocalNumElements());
3319 for (unsigned j = 0; j < XRepeatedMaps[1]->getLocalNumElements(); j += dofsPerNodeVector[1]) {
3320 nullSpaceBasis[1]->getDataNonConst(0)[offset + XRepeatedMaps[1]->getLocalElement(XRepeatedMaps[1]->getGlobalElement(j))] = Teuchos::ScalarTraits<SC>::one();
3321 }
3322
3323 sublist(sublist(_params, "Preconditioner Types"), "FROSch")->set("Null Space Vector", nullSpaceBasis);
3324 }
3325 if(_use_timer)
3326 {
3327 _stacked_timer->stop("TpetraCoreStokes::init_symbolic::build_repeated_maps");
3328 _stacked_timer->start("TpetraCoreStokes::init_symbolic::build_crsgraph_system");
3329 }
3330
3332 GO counter_col_idx = 0;
3333
3335 auto crsgraph = Teuchos::rcp(new Tpetra::CrsGraph<LO, GO, NO >(_map, _num_nze.view(0, static_cast<AGS>(_num_owned_dofs))));
3336
3338 for(IndexType row(0); row < _num_owned_velo_dofs; ++row)
3339 {
3340 Teuchos::Array<GO> cols(static_cast<GO>(_num_nze[static_cast<ARCPS>(row)]));
3341 GO counter = 0;
3342
3343 // Start scanning from rank 0 for each row
3344 int current_rank = 0;
3345
3346 for(auto colidx = this->_row_ptr[row]; colidx < this->_row_ptr[row+1]; ++colidx)
3347 {
3348 GO col = static_cast<GO>(this->_col_idx_mapped[colidx]);
3349
3350 // Advance rank until we find one that might contain this column
3351 while(current_rank < int(_comm->getSize()) && col >= rank_total_end[current_rank])
3352 {
3353 current_rank++;
3354 }
3355
3356 // Check if current rank contains this column
3357 if(current_rank < int(_comm->getSize()) && col >= rank_start[current_rank])
3358 {
3359 if(col < rank_velo_end[current_rank])
3360 {
3361 // Velocity DOF
3362 GO mapped_dof = GO(_all_first_owned_velo_dof[current_rank]) + GO(col - rank_start[current_rank]);
3363 cols[counter++] = mapped_dof;
3364 this->_col_idx_mapped[counter_col_idx++] = mapped_dof;
3365 }
3366 else if (col < rank_total_end[current_rank])
3367 {
3368 // Pressure DOF
3369 GO pres_offset = col - rank_velo_end[current_rank];
3370 GO mapped_dof = GO(_num_global_velo_dofs) + GO(_all_first_owned_pres_dof[current_rank]) + GO(pres_offset);
3371 cols[counter++] = mapped_dof;
3372 this->_col_idx_mapped[counter_col_idx++] = mapped_dof;
3373 }
3374 }
3375 // Note: current_rank is NOT reset - it continues from where it left off
3376 }
3377
3378 crsgraph->insertGlobalIndices(static_cast<GO>(_first_owned_velo_dof + row), cols.view(0, static_cast<AGS>(counter)));
3379 }
3381 for(IndexType row(_num_owned_velo_dofs); row < (_num_owned_velo_dofs + _num_owned_pres_dofs); ++row)
3382 {
3383 Teuchos::Array<GO> cols(static_cast<GO>(_num_nze[static_cast<ARCPS>(row)]));
3384 GO counter = 0;
3385
3386 // Start scanning from rank 0 for each row
3387 int current_rank = 0;
3388
3389 for(auto colidx = this->_row_ptr[row]; colidx < this->_row_ptr[row+1]; ++colidx)
3390 {
3391 GO col = this->_col_idx_mapped[colidx];
3392
3393 // Advance rank until we find one that might contain this column
3394 while(current_rank < int(_comm->getSize()) && col >= rank_total_end[current_rank])
3395 {
3396 current_rank++;
3397 }
3398
3399 // Check if current rank contains this column
3400 if(current_rank < int(_comm->getSize()) && col >= rank_start[current_rank])
3401 {
3402
3403 if(col < rank_velo_end[current_rank])
3404 {
3405 // Velocity DOF
3406 GO mapped_dof = GO(_all_first_owned_velo_dof[current_rank]) + GO(col - rank_start[current_rank]);
3407 cols[counter++] = mapped_dof;
3408 this->_col_idx_mapped[counter_col_idx++] = mapped_dof;
3409 }
3410 else if(col < rank_total_end[current_rank])
3411 {
3412 // Pressure DOF
3413 GO pres_offset = col - rank_velo_end[current_rank];
3414 GO mapped_dof = GO(_num_global_velo_dofs) + GO(_all_first_owned_pres_dof[current_rank]) + GO(pres_offset);
3415 cols[counter++] = mapped_dof;
3416 this->_col_idx_mapped[counter_col_idx++] = mapped_dof;
3417 }
3418 }
3419 }
3420 crsgraph->insertGlobalIndices(static_cast<GO>(_num_global_velo_dofs + _first_owned_pres_dof + (row - _num_owned_velo_dofs)),
3421 cols.view(0, static_cast<AGS>(counter)));
3422 }
3423 crsgraph->fillComplete();
3424
3425 if(_use_timer)
3426 {
3427 _stacked_timer->stop("TpetraCoreStokes::init_symbolic::build_crsgraph_system");
3428 _stacked_timer->start("TpetraCoreStokes::init_symbolic::build_system_matrix");
3429 }
3430
3432 _matrix = Teuchos::rcp(new Tpetra::CrsMatrix<SC, LO, GO, NO>(crsgraph, _params));
3433
3434 // create vectors
3435 const LO num_vectors = 1;
3436 _vec_def = Teuchos::rcp(new Tpetra::MultiVector<SC, LO, GO, NO>(_map, num_vectors));
3437 _vec_cor = Teuchos::rcp(new Tpetra::MultiVector<SC, LO, GO, NO>(_map, num_vectors));
3438
3439 if(_use_timer)
3440 {
3441 _stacked_timer->stop("TpetraCoreStokes::init_symbolic::build_system_matrix");
3442 _stacked_timer->stop("TpetraCoreStokes::init_symbolic");
3443 }
3444 }
3445
3446 virtual void init_numeric() override
3447 {
3448 if(_use_timer)
3449 _stacked_timer->start("TpetraCoreStokes::init_numeric");
3450
3451 Teuchos::ArrayView<const LO> Indices;
3452
3453 _matrix->resumeFill();
3454
3455 AGS entries_start = 0;
3457 for(IndexType row(0); row < _num_owned_velo_dofs; ++row)
3458 {
3459 const auto num_entries = static_cast<Teuchos::Array<SC>::size_type>(_num_nze[static_cast<AGS>(row)]);
3460
3461 // copy vals to Teuchos::Array
3462 Teuchos::Array<SC> vals_ar(num_entries);
3463 for(Teuchos::Array<SC>::size_type j = 0; j < num_entries; ++j)
3464 vals_ar[j] = static_cast<SC>(this->_mat_val[static_cast<VGS>(entries_start + j)]);
3465
3466 Teuchos::Array<GO> cols(num_entries);
3467 for(Teuchos::Array<GO>::size_type j = 0; j < num_entries; ++j)
3468 cols[j] = static_cast<GO>(this->_col_idx_mapped[entries_start + j]);
3469
3470 _matrix->replaceGlobalValues(static_cast<GO>(_first_owned_velo_dof + row),
3471 cols.view(0, static_cast<AGS>(num_entries)),
3472 vals_ar.view(0, static_cast<Teuchos::Array<SC>::size_type>(num_entries)));
3473
3474 entries_start += num_entries;
3475 }
3476
3478 for(IndexType row(0); row < _num_owned_pres_dofs; ++row)
3479 {
3480 const auto num_entries = static_cast<Teuchos::Array<SC>::size_type>(_num_nze[static_cast<AGS>(_num_owned_velo_dofs+row)]);
3481
3482 // copy vals to Teuchos::Array
3483 Teuchos::Array<SC> vals_ar(num_entries);
3484 for(Teuchos::Array<SC>::size_type j = 0; j < num_entries; ++j)
3485 vals_ar[j] = static_cast<SC>(this->_mat_val[static_cast<VGS>(entries_start + j)]);
3486
3487 Teuchos::Array<GO> cols(num_entries);
3488 for(Teuchos::Array<GO>::size_type j = 0; j < num_entries; ++j)
3489 cols[j] = static_cast<GO>(this->_col_idx_mapped[entries_start + j]);
3490
3491 _matrix->replaceGlobalValues(static_cast<GO>(_num_global_velo_dofs + _first_owned_pres_dof + row),
3492 cols.view(0, static_cast<AGS>(num_entries)),
3493 vals_ar.view(0, static_cast<Teuchos::Array<SC>::size_type>(num_entries)));
3494
3495 entries_start += num_entries;
3496 }
3497
3498 _matrix->fillComplete();
3499
3500 if(this->_use_timer)
3501 {
3502 this->_stacked_timer->stop("TpetraCoreStokes::init_numeric");
3503 }
3504 }
3505 };
3506
3507 void* create_core_stokes(const void* comm, Index num_global_dofs, Index my_dof_offset, Index num_owned_dofs, Index num_nonzeros,
3508 IndexType num_owned_velo_dofs, IndexType num_owned_pres_dofs, IndexType first_owned_velo_dof,
3509 IndexType first_owned_pres_dof, IndexType num_global_velo_dofs, IndexType num_global_pres_dofs, const Trilinos::FROSchParameterList & params)
3510 {
3511 return new TpetraCoreStokes(comm, num_global_dofs, my_dof_offset, num_owned_dofs, num_nonzeros, num_owned_velo_dofs, num_owned_pres_dofs, first_owned_velo_dof, first_owned_pres_dof, num_global_velo_dofs, num_global_pres_dofs, params);
3512 }
3513
3514 /* *********************************************************************************************************** */
3515 /* *********************************************************************************************************** */
3516 /* *********************************************************************************************************** */
3517
3518 typedef struct FROSchPrecond
3519 {
3520 using SC = TpetraCoreBase::SC;
3521 Teuchos::RCP<Thyra::PreconditionerFactoryBase<SC>> pfbFactory;
3522 Teuchos::RCP<const Thyra::LinearOpBase<SC>> K_thyra;
3523 Teuchos::RCP<Thyra::PreconditionerBase<SC>> ThyraPrec;
3524 } FROSchPrecond;
3525
3526 //----------------------------------------
3527 void* create_frosch(void* cr)
3528 {
3529 using SC = TpetraCoreBase::SC;
3530 using LO = TpetraCoreBase::LO;
3531 using GO = TpetraCoreBase::GO;
3532 using NO = TpetraCoreBase::NO;
3533
3534 TpetraCoreBase *core = reinterpret_cast<TpetraCoreBase*>(cr);
3535
3536 Stratimikos::DefaultLinearSolverBuilder linearSolverBuilder;
3537 Stratimikos::enableFROSch<SC, LO, GO, NO>(linearSolverBuilder);
3538 // set parameter list
3539 linearSolverBuilder.setParameterList(core->_params);
3540
3541 Teuchos::RCP<Thyra::PreconditionerFactoryBase<SC>> pfbFactory = linearSolverBuilder.createPreconditioningStrategy("");
3542
3543 if(core->_use_timer)
3544 {
3545 Teuchos::RCP<Teuchos::FancyOStream> out = Teuchos::VerboseObjectBase::getDefaultOStream();
3546 pfbFactory->setOStream(out);
3547 pfbFactory->setVerbLevel(Teuchos::VERB_HIGH);
3548 }
3549
3550 // system matrix and right-hand side to Thyra interface
3551 Teuchos::RCP<const Tpetra::RowMatrix<SC, LO, GO, NO>> tpRowMat = Teuchos::rcp_dynamic_cast<const Tpetra::RowMatrix<SC, LO, GO, NO>>(core->_matrix, true);
3552 Teuchos::RCP<const Tpetra::Operator<SC, LO, GO, NO>> tpOperator = Teuchos::rcp_dynamic_cast<const Tpetra::Operator<SC, LO, GO, NO>> (tpRowMat, true);
3553 // thyraOp = Thyra::createConstLinearOp(tpOperator);
3554
3555 // // Xpetra::CrsMatrixWrap<SC, LO, GO, NO> K(core->_matrix);
3556 // // Teuchos::RCP<const Thyra::LinearOpBase<SC>> K_thyra = FROSch::ThyraUtils<SC,LO,GO,NO>::toThyra(K.getCrsMatrix());
3557 // Teuchos::RCP<const Thyra::LinearOpBase<SC>> K_thyra = Thyra::createConstLinearOp(tpOperator);
3558 // Teuchos::RCP<Thyra::PreconditionerBase<SC>> ThyraPrec = Thyra::prec(*pfbFactory, K_thyra);
3559 // Teuchos::RCP<Thyra::LinearOpBase<SC>> LinearPrecOp = ThyraPrec->getNonconstUnspecifiedPrecOp();
3560
3561 // /**
3562 // * It is a bit weird to use a raw pointer to a smart pointer.
3563 // * It would be nicer to use a pointer to Thyra::LinearOpWithSolveBase<SC>, but
3564 // * such things are not intended by the Thyra interface
3565 // * or at least I'm to stupid to achieve this with a little effort.
3566 // **/
3567 // Teuchos::RCP<Thyra::LinearOpBase<SC>> *precond = new Teuchos::RCP<Thyra::LinearOpBase<SC>>();
3568 // *precond = LinearPrecOp;
3569
3570 FROSchPrecond *precond = new FROSchPrecond;
3571 precond->pfbFactory = pfbFactory;
3572 precond->K_thyra = Thyra::createConstLinearOp(tpOperator);
3573 precond->ThyraPrec = Thyra::prec(*pfbFactory, precond->K_thyra);
3574
3575 return precond;
3576 }
3577
3578 void reinit_frosch(void *pre)
3579 {
3580 using SC = TpetraCoreBase::SC;
3581 FROSchPrecond *precond = reinterpret_cast<FROSchPrecond *>(pre);
3582 Thyra::initializePrec<SC>(*(precond->pfbFactory), precond->K_thyra, precond->ThyraPrec.ptr());
3583 }
3584
3585 void destroy_frosch(void *pre)
3586 {
3587 // using SC = TpetraCoreBase::SC;
3588 // Teuchos::RCP<Thyra::LinearOpBase<SC>> *precond = reinterpret_cast<Teuchos::RCP<Thyra::LinearOpBase<SC>> *>(pre);
3589 // delete precond;
3590 FROSchPrecond *precond = reinterpret_cast<FROSchPrecond*>(pre);
3591 delete precond;
3592 }
3593
3594 void solve_frosch(void* cr, void* pre)
3595 {
3596 using SC = TpetraCoreBase::SC;
3597 using LO = TpetraCoreBase::LO;
3598 using GO = TpetraCoreBase::GO;
3599 using NO = TpetraCoreBase::NO;
3600
3601 TpetraCoreBase *core = reinterpret_cast<TpetraCoreBase *>(cr);
3602 FROSchPrecond *fprecond = reinterpret_cast<FROSchPrecond *>(pre);
3603 // Teuchos::RCP<Thyra::LinearOpBase<SC>> *precond = reinterpret_cast<Teuchos::RCP<Thyra::LinearOpBase<SC>> *>(pre);
3604
3605 auto thyTpMap = Thyra::tpetraVectorSpace<SC, LO, GO, NO>(core->_vec_cor->getMap());
3606 auto thyDomMap = Thyra::tpetraVectorSpace<SC, LO, GO, NO>(Tpetra::createLocalMapWithNode<LO, GO, NO>(core->_vec_cor->getNumVectors(), core->_vec_cor->getMap()->getComm()) );
3607 auto thyra_X = Teuchos::rcp(new Thyra::TpetraMultiVector<SC, LO, GO, NO>());
3608 thyra_X->initialize(thyTpMap, thyDomMap, core->_vec_cor);
3609
3610 thyTpMap = Thyra::tpetraVectorSpace<SC, LO, GO, NO>(core->_vec_def->getMap());
3611 thyDomMap = Thyra::tpetraVectorSpace<SC, LO, GO, NO>(Tpetra::createLocalMapWithNode<LO, GO, NO>(core->_vec_def->getNumVectors(), core->_vec_def->getMap()->getComm()) );
3612 auto thyra_B = Teuchos::rcp(new Thyra::TpetraMultiVector<SC, LO, GO, NO>());
3613 thyra_B->initialize(thyTpMap, thyDomMap, core->_vec_def);
3614
3615 //Teuchos::RCP<Thyra::MultiVectorBase<SC>> thyra_X = Teuchos::rcp_const_cast<Thyra::MultiVectorBase<SC> >(FROSch::ThyraUtils<SC,LO,GO,NO>::toThyraMultiVector(core->_vec_cor));
3616 //Teuchos::RCP<const Thyra::MultiVectorBase<SC>> thyra_B = FROSch::ThyraUtils<SC,LO,GO,NO>::toThyraMultiVector((core->_vec_def));
3617
3618 // Thyra::apply<SC>(*(*precond), Thyra::NOTRANS, *(thyra_B.getConst()), thyra_X.ptr());
3619 Teuchos::RCP<Thyra::LinearOpBase<SC>> LinearPrecOp = fprecond->ThyraPrec->getNonconstUnspecifiedPrecOp();
3620 Thyra::apply<SC>(*LinearPrecOp, Thyra::NOTRANS, *(thyra_B.getConst()), thyra_X.ptr());
3621 }
3622 } // namespace Trilinos
3623 } // namespace Solver
3624} // namespace FEAT
3625
3626#else
3627// insert dummy function to suppress linker warnings
3628void dummy_frosch_function() {}
3629#endif // FEAT_HAVE_TRILINOS
XASSERT
#define XASSERT(expr)
Assertion macro definition.
Definition: assertion.hpp:262
XASSERTM
#define XASSERTM(expr, msg)
Assertion macro definition with custom message.
Definition: assertion.hpp:263
base_header.hpp
FEAT Kernel base header.
FEAT
FEAT namespace.
Definition: adjactor.hpp:12
FEAT::Index
std::uint64_t Index
Index data type.
Definition: base_header.hpp:122