16 #include <deal.II/lac/trilinos_index_access.h> 17 #include <deal.II/lac/trilinos_precondition.h> 19 #ifdef DEAL_II_WITH_TRILINOS 20 # if DEAL_II_TRILINOS_VERSION_GTE(11, 14, 0) 22 # include <deal.II/lac/sparse_matrix.h> 23 # include <deal.II/lac/trilinos_index_access.h> 24 # include <deal.II/lac/trilinos_sparse_matrix.h> 25 # include <deal.II/lac/vector.h> 27 # include <Epetra_MultiVector.h> 29 # include <MueLu_EpetraOperator.hpp> 30 # include <MueLu_MLParameterListInterpreter.hpp> 31 # include <Teuchos_ParameterList.hpp> 32 # include <Teuchos_RCP.hpp> 33 # include <ml_MultiLevelPreconditioner.h> 34 # include <ml_include.h> 36 DEAL_II_NAMESPACE_OPEN
42 const unsigned int n_cycles,
44 const double aggregation_threshold,
45 const std::vector<std::vector<bool>> &constant_modes,
46 const unsigned int smoother_sweeps,
47 const unsigned int smoother_overlap,
48 const bool output_details,
49 const char * smoother_type,
50 const char * coarse_type)
54 , aggregation_threshold(aggregation_threshold)
55 , constant_modes(constant_modes)
56 , smoother_sweeps(smoother_sweeps)
57 , smoother_overlap(smoother_overlap)
58 , output_details(output_details)
59 , smoother_type(smoother_type)
60 , coarse_type(coarse_type)
67 # ifdef DEAL_II_WITH_64BIT_INDICES 70 "PreconditionAMGMueLu does not support 64bit-indices!"));
98 Teuchos::ParameterList parameter_list;
100 if (additional_data.
elliptic ==
true)
101 ML_Epetra::SetDefaults(
"SA", parameter_list);
104 ML_Epetra::SetDefaults(
"NSSA", parameter_list);
105 parameter_list.set(
"aggregation: block scaling",
true);
110 parameter_list.set(
"aggregation: type",
"Uncoupled");
112 parameter_list.set(
"smoother: type", additional_data.
smoother_type);
113 parameter_list.set(
"coarse: type", additional_data.
coarse_type);
115 parameter_list.set(
"smoother: sweeps",
117 parameter_list.set(
"cycle applications",
118 static_cast<int>(additional_data.
n_cycles));
119 if (additional_data.
w_cycle ==
true)
120 parameter_list.set(
"prec type",
"MGW");
122 parameter_list.set(
"prec type",
"MGV");
124 parameter_list.set(
"smoother: Chebyshev alpha", 10.);
125 parameter_list.set(
"smoother: ifpack overlap",
127 parameter_list.set(
"aggregation: threshold",
129 parameter_list.set(
"coarse: max size", 2000);
132 parameter_list.set(
"ML output", 10);
134 parameter_list.set(
"ML output", 0);
136 const Epetra_Map &domain_map = matrix.OperatorDomainMap();
138 const size_type constant_modes_dimension =
140 Epetra_MultiVector distributed_constant_modes(
141 domain_map, constant_modes_dimension > 0 ? constant_modes_dimension : 1);
142 std::vector<double> dummy(constant_modes_dimension);
144 if (constant_modes_dimension > 0)
147 const bool constant_modes_are_global =
150 constant_modes_are_global ? n_rows :
152 const size_type my_size = domain_map.NumMyElements();
153 if (constant_modes_are_global ==
false)
154 Assert(n_relevant_rows == my_size,
157 distributed_constant_modes)),
160 distributed_constant_modes)));
162 (void)n_relevant_rows;
167 for (
size_type d = 0; d < constant_modes_dimension; ++d)
168 for (
size_type row = 0; row < my_size; ++row)
170 TrilinosWrappers::types::int_type global_row_id =
171 constant_modes_are_global ?
174 distributed_constant_modes[d][row] =
178 parameter_list.set(
"null space: type",
"pre-computed");
179 parameter_list.set(
"null space: dimension",
180 distributed_constant_modes.NumVectors());
182 parameter_list.set(
"null space: vectors",
183 distributed_constant_modes.Values());
187 parameter_list.set(
"null space: vectors", dummy.data());
197 Teuchos::ParameterList &muelu_parameters)
206 Teuchos::ParameterList &muelu_parameters)
213 using node = KokkosClassic::DefaultNode::DefaultNodeType;
218 Teuchos::RCP<Epetra_CrsMatrix> rcp_matrix =
219 Teuchos::rcpFromRef(*(const_cast<Epetra_CrsMatrix *>(&matrix)));
220 Teuchos::RCP<Xpetra::CrsMatrix<double, int, int, node>> muelu_crs_matrix =
221 Teuchos::rcp(
new Xpetra::EpetraCrsMatrix(rcp_matrix));
222 Teuchos::RCP<Xpetra::Matrix<double, int, int, node>> muelu_matrix =
224 new Xpetra::CrsMatrixWrap<double, int, int, node>(muelu_crs_matrix));
227 Teuchos::RCP<MueLu::HierarchyManager<double, int, int, node>>
229 hierarchy_factory = Teuchos::rcp(
230 new MueLu::MLParameterListInterpreter<double, int, int, node>(
232 Teuchos::RCP<MueLu::Hierarchy<double, int, int, node>> hierarchy =
233 hierarchy_factory->CreateHierarchy();
234 hierarchy->GetLevel(0)->Set(
"A", muelu_matrix);
235 hierarchy_factory->SetupHierarchy(*hierarchy);
239 preconditioner = std::make_shared<MueLu::EpetraOperator>(hierarchy);
244 template <
typename number>
247 const ::SparseMatrix<number> &deal_ii_sparse_matrix,
249 const double drop_tolerance,
250 const ::SparsityPattern * use_this_sparsity)
253 const size_type n_rows = deal_ii_sparse_matrix.m();
260 static_cast<TrilinosWrappers::types::int_type
>(n_rows), 0,
communicator);
267 deal_ii_sparse_matrix,
289 unsigned int memory =
sizeof(*this);
305 const ::SparsityPattern *);
310 const ::SparsityPattern *);
314 DEAL_II_NAMESPACE_CLOSE
316 # endif // DEAL_II_TRILINOS_VERSION_GTE(11,14,0) 317 #endif // DEAL_II_WITH_TRILINOS unsigned int smoother_sweeps
void initialize(const SparseMatrix &matrix, const AdditionalData &additional_data=AdditionalData())
::types::global_dof_index size_type
const Epetra_CrsMatrix & trilinos_matrix() const
TrilinosWrappers::types::int_type global_length(const Epetra_MultiVector &vector)
size_type memory_consumption() const
const char * smoother_type
std::shared_ptr< Epetra_Map > vector_distributor
#define AssertThrow(cond, exc)
static::ExceptionBase & ExcMessage(std::string arg1)
TrilinosWrappers::types::int_type n_global_rows(const Epetra_CrsGraph &graph)
double aggregation_threshold
#define Assert(cond, exc)
static::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
std::vector< std::vector< bool > > constant_modes
AdditionalData(const bool elliptic=true, const unsigned int n_cycles=1, const bool w_cyle=false, const double aggregation_threshold=1e-4, const std::vector< std::vector< bool >> &constant_modes=std::vector< std::vector< bool >>(0), const unsigned int smoother_sweeps=2, const unsigned int smoother_overlap=0, const bool output_details=false, const char *smoother_type="Chebyshev", const char *coarse_type="Amesos-KLU")
~PreconditionAMGMueLu() override
std::shared_ptr< SparseMatrix > trilinos_matrix
TrilinosWrappers::types::int_type global_index(const Epetra_BlockMap &map, const ::types::global_dof_index i)
std::shared_ptr< Epetra_Operator > preconditioner
unsigned int smoother_overlap
Epetra_MpiComm communicator