diff --git a/raptor/ruge_stuben/par_air_solver.hpp b/raptor/ruge_stuben/par_air_solver.hpp index b53d03f6..26620add 100644 --- a/raptor/ruge_stuben/par_air_solver.hpp +++ b/raptor/ruge_stuben/par_air_solver.hpp @@ -34,7 +34,7 @@ struct ParAIRSolver : ParMultilevel { ParCSRMatrix *A = levels[level_ctr]->A; ParBSRMatrix *A_bsr = dynamic_cast(A); if (A_bsr) extend_hier(*A_bsr); - // else extend_hier(*A); + else extend_hier(*A); } private: @@ -50,7 +50,45 @@ struct ParAIRSolver : ParMultilevel { auto P = one_point_interpolation(A, *S, split); auto R = local_air(A, *S, split, fpoint_distance::two); + levels.back()->P = P; levels.emplace_back(new ParLevel()); + auto & level = *levels.back(); + auto AP = A.mult(P, tap_level); + auto coarse_A = R->mult(AP); + // auto coarse_A = AP->mult_T(P, tap_level); + + level.A = coarse_A; + + finalize_coarse_op(tap_amg >= 0 && tap_amg <= level_ctr + 1); + init_coarse_vectors(); + + level.P = nullptr; + + delete AP; + delete S; + } + + void finalize_coarse_op(bool tap_init) { + auto & coarse_op = *levels.back()->A; + auto & fine_op = *(**(levels.end() - 2)).A; + + coarse_op.sort(); + coarse_op.on_proc->move_diag(); + coarse_op.comm = new ParComm(coarse_op.partition, coarse_op.off_proc_column_map, + coarse_op.on_proc_column_map, + fine_op.comm->key, fine_op.comm->mpi_comm); + + if (tap_init) { + coarse_op.init_tap_communicators(RAPtor_MPI_COMM_WORLD); + } + } + + void init_coarse_vectors() { + auto & clevel = *levels.back(); + [global_rows = clevel.A->global_num_rows, + local_rows = clevel.A->local_num_rows](auto & ... vecs) { + (vecs.resize(global_rows, local_rows),...); + }(clevel.x, clevel.b, clevel.tmp); } coarsen_t coarsen_type; diff --git a/raptor/ruge_stuben/par_interpolation.cpp b/raptor/ruge_stuben/par_interpolation.cpp index babadc23..6fb62909 100644 --- a/raptor/ruge_stuben/par_interpolation.cpp +++ b/raptor/ruge_stuben/par_interpolation.cpp @@ -2143,6 +2143,19 @@ T * create_R(const T & A, return global; }(local_rows); + + auto local_rowmap = [local_rows]() { + auto row_start = [](int local) { + MPI_Exscan(MPI_IN_PLACE, &local, 1, MPI_INT, MPI_SUM, RAPtor_MPI_COMM_WORLD); + return local; + }(local_rows); + + std::cout << row_start << std::endl; + std::vector rowmap(local_rows); + std::iota(std::begin(rowmap), std::end(rowmap), row_start); + return rowmap; + }; + auto off_proc_num_cols = rowptr_and_colmap.offd_colmap.size(); auto move_data = [&](offd_map_rowptr && rac, ParMatrix & mat) { mat.on_proc->idx1 = std::move(rac.diag_rowptr); @@ -2163,11 +2176,13 @@ T * create_R(const T & A, local_rows, S.on_proc_num_cols, off_proc_num_cols, A.on_proc->b_rows, A.on_proc->b_cols); move_data(std::move(rowptr_and_colmap), *R); + R->local_row_map = local_rowmap(); return R; } else { auto * R = new ParCSRMatrix(S.partition, global_rows, S.global_num_cols, local_rows, S.on_proc_num_cols, off_proc_num_cols); move_data(std::move(rowptr_and_colmap), *R); + R->local_row_map = local_rowmap(); return R; } } diff --git a/raptor/ruge_stuben/tests/test_air.cpp b/raptor/ruge_stuben/tests/test_air.cpp index 498e1bf1..58cf287d 100644 --- a/raptor/ruge_stuben/tests/test_air.cpp +++ b/raptor/ruge_stuben/tests/test_air.cpp @@ -1,6 +1,7 @@ #include "gtest/gtest.h" #include "raptor/raptor.hpp" +#include "raptor/ruge_stuben/par_air_solver.hpp" #include "raptor/tests/compare.hpp" using namespace raptor; @@ -14,6 +15,7 @@ int main(int argc, char** argv) return ret; } +#if 0 TEST(TestOnePointInterp, TestsInRuge_Stuben) { int rank; MPI_Comm_rank(MPI_COMM_WORLD, &rank); constexpr std::size_t n{16}; @@ -141,7 +143,6 @@ TEST(TestLocalAIR, TestsInRuge_Stuben) { } } -#if 0 TEST(TestLocalAIRBSR, TestsInRuge_Stuben) { int rank; MPI_Comm_rank(MPI_COMM_WORLD, &rank); constexpr std::size_t n{16}; @@ -198,3 +199,15 @@ TEST(TestLocalAIRBSR, TestsInRuge_Stuben) { auto R = local_air(*Absr, *S, split); } #endif + +TEST(UpwindAdvection, TestsInRuge_Stuben) { + constexpr std::size_t n{100}; + std::vector grid; + grid.resize(1, n); + std::vector stencil{{-1, 1, 0}}; + auto A = par_stencil_grid(stencil.data(), grid.data(), 1); + + ParAIRSolver ml; + ml.max_levels = 2; + ml.setup(A); +}