Skip to content

Commit e97e1f4

Browse files
ddemidovddemidov
committed
[wip] [skip ci] Enable amgcl/opencl in trilinos application
1 parent dc926f4 commit e97e1f4

4 files changed

Lines changed: 251 additions & 83 deletions

File tree

‎applications/trilinos_application/CMakeLists.txt‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -50,6 +50,7 @@ set( KRATOS_TRILINOS_APPLICATION_SOURCES
5050
${CMAKE_CURRENT_SOURCE_DIR}/custom_factories/trilinos_linear_solver_factory.cpp
5151
${CMAKE_SOURCE_DIR}/applications/FluidDynamicsApplication/fluid_dynamics_application_variables.cpp #TODO: this should REALLY NOT BE HERE
5252
${CMAKE_CURRENT_SOURCE_DIR}/custom_utilities/mpi_normal_calculation_utilities.cpp;
53+
${CMAKE_CURRENT_SOURCE_DIR}/amgcl_mpi_solver_impl.cpp;
5354
)
5455

5556
## Kratos tests sources. Enabled by default

‎applications/trilinos_application/amgcl_mpi_solver_impl.cpp‎

Lines changed: 244 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,244 @@
1+
#ifndef KRATOS_AMGCL_MPI_SOLVE_FUNCTIONS_H
2+
#define KRATOS_AMGCL_MPI_SOLVE_FUNCTIONS_H
3+
4+
#include <boost/range/iterator_range.hpp>
5+
#include <boost/property_tree/ptree.hpp>
6+
7+
#include <amgcl/adapter/crs_tuple.hpp>
8+
#include <amgcl/adapter/epetra.hpp>
9+
#include <amgcl/adapter/ublas.hpp>
10+
#include <amgcl/adapter/zero_copy.hpp>
11+
#include <amgcl/adapter/block_matrix.hpp>
12+
#include <amgcl/backend/builtin.hpp>
13+
#include <amgcl/value_type/static_matrix.hpp>
14+
#include <amgcl/solver/runtime.hpp>
15+
16+
#include <amgcl/mpi/util.hpp>
17+
#include <amgcl/mpi/make_solver.hpp>
18+
#include <amgcl/mpi/amg.hpp>
19+
#include <amgcl/mpi/coarsening/runtime.hpp>
20+
#include <amgcl/mpi/relaxation/runtime.hpp>
21+
#include <amgcl/mpi/direct_solver/runtime.hpp>
22+
#include <amgcl/mpi/partition/runtime.hpp>
23+
24+
#ifdef KRATOS_HAVE_OPENCL
25+
# include <amgcl/backend/vexcl.hpp>
26+
# include <amgcl/backend/vexcl_static_matrix.hpp>
27+
#endif
28+
29+
#include "Epetra_FECrsMatrix.h"
30+
#include "Epetra_FEVector.h"
31+
#include "trilinos_space.h"
32+
33+
namespace Kratos
34+
{
35+
36+
#ifdef KRATOS_HAVE_OPENCL
37+
vex::Context& vexcl_context();
38+
39+
template <int TBlockSize>
40+
void register_vexcl_static_matrix_type();
41+
#endif
42+
43+
// Spacialization of AMGCLScalarSolve for distribued systems.
44+
template <class TSparseSpaceType>
45+
typename std::enable_if<TSparseSpaceType::IsDistributed(), void>::type
46+
AMGCLScalarSolve(
47+
typename TSparseSpaceType::MatrixType& rA,
48+
typename TSparseSpaceType::VectorType& rX,
49+
typename TSparseSpaceType::VectorType& rB,
50+
typename TSparseSpaceType::IndexType& rIterationNumber,
51+
double& rResidual,
52+
const boost::property_tree::ptree &amgclParams,
53+
int verbosity_level,
54+
bool use_opencl
55+
)
56+
{
57+
#ifdef KRATOS_HAVE_OPENCL
58+
if (use_opencl && vexcl_context()) {
59+
auto &ctx = vexcl_context();
60+
61+
typedef amgcl::backend::vexcl<double> Backend;
62+
63+
typedef
64+
amgcl::mpi::make_solver<
65+
amgcl::mpi::amg<
66+
Backend,
67+
amgcl::runtime::mpi::coarsening::wrapper<Backend>,
68+
amgcl::runtime::mpi::relaxation::wrapper<Backend>,
69+
amgcl::runtime::mpi::direct::solver<double>,
70+
amgcl::runtime::mpi::partition::wrapper<Backend>
71+
>,
72+
amgcl::runtime::solver::wrapper
73+
>
74+
Solver;
75+
76+
Backend::params bprm;
77+
bprm.q = ctx;
78+
79+
Solver solve(MPI_COMM_WORLD, amgcl::adapter::map(rA), amgclParams, bprm);
80+
81+
std::size_t n = rA.NumMyRows();
82+
83+
vex::vector<double> b(ctx, n, rB.Values());
84+
vex::vector<double> x(ctx, n, rX.Values());
85+
86+
std::tie(rIterationNumber, rResidual) = solve(b, x);
87+
88+
vex::copy(x.begin(), x.end(), rX.Values());
89+
} else
90+
#endif
91+
{
92+
typedef amgcl::backend::builtin<double> Backend;
93+
94+
typedef
95+
amgcl::mpi::make_solver<
96+
amgcl::mpi::amg<
97+
Backend,
98+
amgcl::runtime::mpi::coarsening::wrapper<Backend>,
99+
amgcl::runtime::mpi::relaxation::wrapper<Backend>,
100+
amgcl::runtime::mpi::direct::solver<double>,
101+
amgcl::runtime::mpi::partition::wrapper<Backend>
102+
>,
103+
amgcl::runtime::solver::wrapper
104+
>
105+
Solver;
106+
107+
Solver solve(MPI_COMM_WORLD, amgcl::adapter::map(rA), amgclParams);
108+
109+
std::size_t n = rA.NumMyRows();
110+
111+
auto b_range = boost::make_iterator_range(rB.Values(), rB.Values() + n);
112+
auto x_range = boost::make_iterator_range(rX.Values(), rX.Values() + n);
113+
114+
std::tie(rIterationNumber, rResidual) = solve(b_range, x_range);
115+
}
116+
}
117+
118+
// Spacialization of AMGCLBlockSolve for distribued systems.
119+
template <int TBlockSize, class TSparseSpaceType>
120+
typename std::enable_if<TSparseSpaceType::IsDistributed(), void>::type
121+
AMGCLBlockSolve(
122+
typename TSparseSpaceType::MatrixType & rA,
123+
typename TSparseSpaceType::VectorType& rX,
124+
typename TSparseSpaceType::VectorType& rB,
125+
typename TSparseSpaceType::IndexType& rIterationNumber,
126+
double& rResidual,
127+
boost::property_tree::ptree amgclParams,
128+
int verbosity_level,
129+
bool use_opencl
130+
)
131+
{
132+
amgclParams.put("precond.coarsening.aggr.block_size",1);
133+
134+
typedef amgcl::static_matrix<double, TBlockSize, TBlockSize> val_type;
135+
typedef amgcl::static_matrix<double, TBlockSize, 1> rhs_type;
136+
137+
std::size_t n = rA.RowMap().NumMyElements();
138+
std::size_t nb = n / TBlockSize;
139+
140+
#ifdef KRATOS_HAVE_OPENCL
141+
if (use_opencl && vexcl_context()) {
142+
auto &ctx = vexcl_context();
143+
register_vexcl_static_matrix_type<TBlockSize>();
144+
145+
typedef amgcl::backend::vexcl<val_type> Backend;
146+
147+
typedef
148+
amgcl::mpi::make_solver<
149+
amgcl::mpi::amg<
150+
Backend,
151+
amgcl::runtime::mpi::coarsening::wrapper<Backend>,
152+
amgcl::runtime::mpi::relaxation::wrapper<Backend>,
153+
amgcl::runtime::mpi::direct::solver<val_type>,
154+
amgcl::runtime::mpi::partition::wrapper<Backend>
155+
>,
156+
amgcl::runtime::solver::wrapper
157+
>
158+
Solver;
159+
160+
typename Backend::params bprm;
161+
bprm.q = ctx;
162+
163+
Solver solve(
164+
MPI_COMM_WORLD,
165+
amgcl::adapter::block_matrix<val_type>(amgcl::adapter::map(rA)),
166+
amgclParams, bprm
167+
);
168+
169+
auto b_begin = reinterpret_cast<const rhs_type*>(rB.Values());
170+
auto x_begin = reinterpret_cast<rhs_type*>(rX.Values());
171+
172+
vex::vector<rhs_type> x(ctx, nb, x_begin);
173+
vex::vector<rhs_type> b(ctx, nb, b_begin);
174+
175+
std::tie(rIterationNumber, rResidual) = solve(b, x);
176+
177+
vex::copy(x.begin(), x.end(), x_begin);
178+
} else
179+
#endif
180+
{
181+
typedef amgcl::backend::builtin<val_type> Backend;
182+
183+
typedef
184+
amgcl::mpi::make_solver<
185+
amgcl::mpi::amg<
186+
Backend,
187+
amgcl::runtime::mpi::coarsening::wrapper<Backend>,
188+
amgcl::runtime::mpi::relaxation::wrapper<Backend>,
189+
amgcl::runtime::mpi::direct::solver<val_type>,
190+
amgcl::runtime::mpi::partition::wrapper<Backend>
191+
>,
192+
amgcl::runtime::solver::wrapper
193+
>
194+
Solver;
195+
196+
Solver solve(
197+
MPI_COMM_WORLD,
198+
amgcl::adapter::block_matrix<val_type>(amgcl::adapter::map(rA)),
199+
amgclParams
200+
);
201+
202+
auto b_begin = reinterpret_cast<const rhs_type*>(rB.Values());
203+
auto x_begin = reinterpret_cast<rhs_type*>(rX.Values());
204+
205+
auto b_range = boost::make_iterator_range(b_begin, b_begin + n / TBlockSize);
206+
auto x_range = boost::make_iterator_range(x_begin, x_begin + n / TBlockSize);
207+
208+
std::tie(rIterationNumber, rResidual) = solve(b_range, x_range);
209+
}
210+
}
211+
212+
// Exlplicit instantiations:
213+
template void AMGCLScalarSolve< TrilinosSpace<Epetra_FECrsMatrix, Epetra_FEVector> >(
214+
TrilinosSpace<Epetra_FECrsMatrix, Epetra_FEVector>::MatrixType& rA,
215+
TrilinosSpace<Epetra_FECrsMatrix, Epetra_FEVector>::VectorType& rX,
216+
TrilinosSpace<Epetra_FECrsMatrix, Epetra_FEVector>::VectorType& rB,
217+
TrilinosSpace<Epetra_FECrsMatrix, Epetra_FEVector>::IndexType& rIterationNumber,
218+
double& rResidual,
219+
const boost::property_tree::ptree &amgclParams,
220+
int verbosity_level,
221+
bool use_opencl
222+
);
223+
224+
#define INSTANTIATE_BLOCK_SOLVER(B) \
225+
template void AMGCLBlockSolve<B, TrilinosSpace<Epetra_FECrsMatrix, Epetra_FEVector> >( \
226+
TrilinosSpace<Epetra_FECrsMatrix, Epetra_FEVector>::MatrixType& rA, \
227+
TrilinosSpace<Epetra_FECrsMatrix, Epetra_FEVector>::VectorType& rX, \
228+
TrilinosSpace<Epetra_FECrsMatrix, Epetra_FEVector>::VectorType& rB, \
229+
TrilinosSpace<Epetra_FECrsMatrix, Epetra_FEVector>::IndexType& rIterationNumber, \
230+
double& rResidual, \
231+
boost::property_tree::ptree amgclParams, \
232+
int verbosity_level, \
233+
bool use_opencl \
234+
)
235+
236+
INSTANTIATE_BLOCK_SOLVER(2);
237+
INSTANTIATE_BLOCK_SOLVER(3);
238+
INSTANTIATE_BLOCK_SOLVER(4);
239+
240+
#undef INSTANTIATE_BLOCK_SOLVER
241+
242+
} // namespace Kratos
243+
244+
#endif

‎applications/trilinos_application/external_includes/amgcl_mpi_solve_functions.h‎

Lines changed: 5 additions & 83 deletions
Original file line numberDiff line numberDiff line change
@@ -1,24 +1,7 @@
11
#ifndef KRATOS_AMGCL_MPI_SOLVE_FUNCTIONS_H
22
#define KRATOS_AMGCL_MPI_SOLVE_FUNCTIONS_H
33

4-
#include <boost/range/iterator_range.hpp>
5-
6-
#include <amgcl/adapter/crs_tuple.hpp>
7-
#include <amgcl/adapter/epetra.hpp>
8-
#include <amgcl/adapter/ublas.hpp>
9-
#include <amgcl/adapter/zero_copy.hpp>
10-
#include <amgcl/adapter/block_matrix.hpp>
11-
#include <amgcl/backend/builtin.hpp>
12-
#include <amgcl/value_type/static_matrix.hpp>
13-
#include <amgcl/solver/runtime.hpp>
14-
15-
#include <amgcl/mpi/util.hpp>
16-
#include <amgcl/mpi/make_solver.hpp>
17-
#include <amgcl/mpi/amg.hpp>
18-
#include <amgcl/mpi/coarsening/runtime.hpp>
19-
#include <amgcl/mpi/relaxation/runtime.hpp>
20-
#include <amgcl/mpi/direct_solver/runtime.hpp>
21-
#include <amgcl/mpi/partition/runtime.hpp>
4+
#include <boost/property_tree/ptree.hpp>
225

236
namespace Kratos
247
{
@@ -34,33 +17,8 @@ AMGCLScalarSolve(
3417
double& rResidual,
3518
const boost::property_tree::ptree &amgclParams,
3619
int verbosity_level,
37-
bool /*use_opencl*/
38-
)
39-
{
40-
typedef amgcl::backend::builtin<double> Backend;
41-
42-
typedef
43-
amgcl::mpi::make_solver<
44-
amgcl::mpi::amg<
45-
Backend,
46-
amgcl::runtime::mpi::coarsening::wrapper<Backend>,
47-
amgcl::runtime::mpi::relaxation::wrapper<Backend>,
48-
amgcl::runtime::mpi::direct::solver<double>,
49-
amgcl::runtime::mpi::partition::wrapper<Backend>
50-
>,
51-
amgcl::runtime::solver::wrapper
52-
>
53-
Solver;
54-
55-
Solver solve(MPI_COMM_WORLD, amgcl::adapter::map(rA), amgclParams);
56-
57-
std::size_t n = rA.NumMyRows();
58-
59-
auto b_range = boost::make_iterator_range(rB.Values(), rB.Values() + n);
60-
auto x_range = boost::make_iterator_range(rX.Values(), rX.Values() + n);
61-
62-
std::tie(rIterationNumber, rResidual) = solve(b_range, x_range);
63-
}
20+
bool use_opencl
21+
);
6422

6523
// Spacialization of AMGCLBlockSolve for distribued systems.
6624
template <int TBlockSize, class TSparseSpaceType>
@@ -73,44 +31,8 @@ AMGCLBlockSolve(
7331
double& rResidual,
7432
boost::property_tree::ptree amgclParams,
7533
int verbosity_level,
76-
bool /*use_opencl*/
77-
)
78-
{
79-
amgclParams.put("precond.coarsening.aggr.block_size",1);
80-
81-
typedef amgcl::static_matrix<double, TBlockSize, TBlockSize> val_type;
82-
typedef amgcl::static_matrix<double, TBlockSize, 1> rhs_type;
83-
typedef amgcl::backend::builtin<val_type> Backend;
84-
85-
std::size_t n = rA.RowMap().NumMyElements();
86-
87-
typedef
88-
amgcl::mpi::make_solver<
89-
amgcl::mpi::amg<
90-
Backend,
91-
amgcl::runtime::mpi::coarsening::wrapper<Backend>,
92-
amgcl::runtime::mpi::relaxation::wrapper<Backend>,
93-
amgcl::runtime::mpi::direct::solver<val_type>,
94-
amgcl::runtime::mpi::partition::wrapper<Backend>
95-
>,
96-
amgcl::runtime::solver::wrapper
97-
>
98-
Solver;
99-
100-
Solver solve(
101-
MPI_COMM_WORLD,
102-
amgcl::adapter::block_matrix<val_type>(amgcl::adapter::map(rA)),
103-
amgclParams
104-
);
105-
106-
auto b_begin = reinterpret_cast<const rhs_type*>(rB.Values());
107-
auto x_begin = reinterpret_cast<rhs_type*>(rX.Values());
108-
109-
auto b_range = boost::make_iterator_range(b_begin, b_begin + n / TBlockSize);
110-
auto x_range = boost::make_iterator_range(x_begin, x_begin + n / TBlockSize);
111-
112-
std::tie(rIterationNumber, rResidual) = solve(b_range, x_range);
113-
}
34+
bool use_opencl
35+
);
11436

11537
} // namespace Kratos
11638

‎external_libraries/amgcl/adapter/epetra.hpp‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -38,6 +38,7 @@ THE SOFTWARE.
3838
#include <Epetra_CrsMatrix.h>
3939
#include <Epetra_IntVector.h>
4040
#include <Epetra_Import.h>
41+
#include <Epetra_Comm.h>
4142

4243
#include <amgcl/backend/interface.hpp>
4344

0 commit comments

Comments
 (0)