diff --git a/include/samurai/interface.hpp b/include/samurai/interface.hpp index 4b48c9e7f..3abe562f7 100644 --- a/include/samurai/interface.hpp +++ b/include/samurai/interface.hpp @@ -21,7 +21,7 @@ namespace samurai using mesh_id_t = typename Mesh::mesh_id_t; using mesh_interval_t = typename Mesh::mesh_interval_t; using coord_index_t = typename Mesh::config::interval_t::coord_index_t; - using Cell = typename samurai::Cell; + using cell_t = typename samurai::Cell; Stencil<2, dim> interface_stencil = in_out_stencil(direction); auto interface_it = make_stencil_iterator(mesh, interface_stencil); @@ -56,12 +56,12 @@ namespace samurai int direction_index_int = find(comput_stencil, direction); std::size_t direction_index = static_cast(direction_index_int); - Stencil reversed = xt::eval(xt::flip(comput_stencil, 0)); - auto minus_comput_stencil = xt::eval(-reversed); - auto minus_direction = xt::eval(-direction); - int minus_direction_index_int = find(minus_comput_stencil, minus_direction); - std::size_t minus_direction_index = static_cast(minus_direction_index_int); - auto minus_comput_stencil_it = make_stencil_iterator(mesh, minus_comput_stencil); + Stencil reversed = xt::eval(xt::flip(comput_stencil, 0)); + auto minus_comput_stencil = xt::eval(-reversed); + auto minus_direction = xt::eval(-direction); + int minus_direction_index_int = find(minus_comput_stencil, minus_direction); + std::size_t minus_direction_index = static_cast(minus_direction_index_int); + auto minus_comput_stencil_it = make_stencil_iterator(mesh, minus_comput_stencil); for_each_level( mesh, @@ -81,17 +81,17 @@ namespace samurai fine_intersect, [&](auto fine_mesh_interval) { - mesh_interval_t coarse_mesh_interval(level, fine_mesh_interval.i / 2, fine_mesh_interval.index / 2); + mesh_interval_t coarse_mesh_interval(level, fine_mesh_interval.i >> 1, fine_mesh_interval.index >> 1); comput_stencil_it.init(fine_mesh_interval); coarse_it.init(coarse_mesh_interval); for (std::size_t ii = 0; ii < fine_mesh_interval.i.size(); ++ii) { - std::array interface_cells; + std::array interface_cells; interface_cells[0] = coarse_it.cells()[0]; - interface_cells[1] = comput_stencil_it.cells()[direction_index]; + f(interface_cells, comput_stencil_it.cells()); comput_stencil_it.move_next(); @@ -111,14 +111,14 @@ namespace samurai fine_intersect, [&](auto fine_mesh_interval) { - mesh_interval_t coarse_mesh_interval(level, fine_mesh_interval.i / 2, fine_mesh_interval.index / 2); + mesh_interval_t coarse_mesh_interval(level, fine_mesh_interval.i >> 1, fine_mesh_interval.index >> 1); minus_comput_stencil_it.init(fine_mesh_interval); coarse_it.init(coarse_mesh_interval); for (std::size_t ii = 0; ii < fine_mesh_interval.i.size(); ++ii) { - std::array interface_cells; + std::array interface_cells; interface_cells[0] = minus_comput_stencil_it.cells()[minus_direction_index]; interface_cells[1] = coarse_it.cells()[0]; @@ -141,15 +141,15 @@ namespace samurai Vector direction, const Stencil& comput_stencil, GetFluxCoeffsFunc get_flux_coeffs, - GetCell1CoeffsFunc get_cell1_coeffs, - GetCell2CoeffsFunc get_cell2_coeffs, + GetCell1CoeffsFunc get_left_cell_coeffs, + GetCell2CoeffsFunc get_right_cell_coeffs, Func&& f) { static constexpr std::size_t dim = Mesh::dim; using mesh_id_t = typename Mesh::mesh_id_t; using mesh_interval_t = typename Mesh::mesh_interval_t; using coord_index_t = typename Mesh::config::interval_t::coord_index_t; - using Cell = typename samurai::Cell; + using cell_t = typename samurai::Cell; Stencil<2, dim> interface_stencil = in_out_stencil(direction); auto interface_it = make_stencil_iterator(mesh, interface_stencil); @@ -163,10 +163,10 @@ namespace samurai auto shifted_cells = translate(cells, -direction); auto intersect = intersection(cells, shifted_cells); - auto h = cell_length(level); - auto flux_coeffs = get_flux_coeffs(h); - auto cell1_coeffs = get_cell1_coeffs(flux_coeffs, h, h); - auto cell2_coeffs = get_cell2_coeffs(flux_coeffs, h, h); + auto h = cell_length(level); + auto flux_coeffs = get_flux_coeffs(h); + auto left_cell_coeffs = get_left_cell_coeffs(flux_coeffs, h, h); + auto right_cell_coeffs = get_right_cell_coeffs(flux_coeffs, h, h); for_each_meshinterval( intersect, @@ -176,7 +176,7 @@ namespace samurai comput_stencil_it.init(mesh_interval); for (std::size_t ii = 0; ii < mesh_interval.i.size(); ++ii) { - f(interface_it.cells(), comput_stencil_it.cells(), cell1_coeffs, cell2_coeffs); + f(interface_it.cells(), comput_stencil_it.cells(), left_cell_coeffs, right_cell_coeffs); interface_it.move_next(); comput_stencil_it.move_next(); } @@ -190,12 +190,12 @@ namespace samurai int direction_index_int = find(comput_stencil, direction); std::size_t direction_index = static_cast(direction_index_int); - Stencil reversed = xt::eval(xt::flip(comput_stencil, 0)); - auto minus_comput_stencil = xt::eval(-reversed); - auto minus_direction = xt::eval(-direction); - int minus_direction_index_int = find(minus_comput_stencil, minus_direction); - std::size_t minus_direction_index = static_cast(minus_direction_index_int); - auto minus_comput_stencil_it = make_stencil_iterator(mesh, minus_comput_stencil); + Stencil reversed = xt::eval(xt::flip(comput_stencil, 0)); + auto minus_comput_stencil = xt::eval(-reversed); + auto minus_direction = xt::eval(-direction); + int minus_direction_index_int = find(minus_comput_stencil, minus_direction); + std::size_t minus_direction_index = static_cast(minus_direction_index_int); + auto minus_comput_stencil_it = make_stencil_iterator(mesh, minus_comput_stencil); for_each_level( mesh, @@ -215,25 +215,25 @@ namespace samurai auto shifted_fine_cells = translate(fine_cells, -direction); auto fine_intersect = intersection(coarse_cells, shifted_fine_cells).on(level + 1); - auto cell1_coeffs = get_cell1_coeffs(flux_coeffs, h_lp1, h_l); - auto cell2_coeffs = get_cell2_coeffs(flux_coeffs, h_lp1, h_lp1); + auto left_cell_coeffs = get_left_cell_coeffs(flux_coeffs, h_lp1, h_l); + auto right_cell_coeffs = get_right_cell_coeffs(flux_coeffs, h_lp1, h_lp1); for_each_meshinterval( fine_intersect, [&](auto fine_mesh_interval) { - mesh_interval_t coarse_mesh_interval(level, fine_mesh_interval.i / 2, fine_mesh_interval.index / 2); + mesh_interval_t coarse_mesh_interval(level, fine_mesh_interval.i >> 1, fine_mesh_interval.index >> 1); comput_stencil_it.init(fine_mesh_interval); coarse_it.init(coarse_mesh_interval); for (std::size_t ii = 0; ii < fine_mesh_interval.i.size(); ++ii) { - std::array interface_cells; + std::array interface_cells; interface_cells[0] = coarse_it.cells()[0]; - interface_cells[1] = comput_stencil_it.cells()[direction_index]; - f(interface_cells, comput_stencil_it.cells(), cell1_coeffs, cell2_coeffs); + + f(interface_cells, comput_stencil_it.cells(), left_cell_coeffs, right_cell_coeffs); comput_stencil_it.move_next(); if (ii % 2 == 1) @@ -248,25 +248,25 @@ namespace samurai auto shifted_fine_cells = translate(fine_cells, direction); auto fine_intersect = intersection(coarse_cells, shifted_fine_cells).on(level + 1); - auto cell1_coeffs = get_cell1_coeffs(flux_coeffs, h_lp1, h_lp1); - auto cell2_coeffs = get_cell2_coeffs(flux_coeffs, h_lp1, h_l); + auto left_cell_coeffs = get_left_cell_coeffs(flux_coeffs, h_lp1, h_lp1); + auto right_cell_coeffs = get_right_cell_coeffs(flux_coeffs, h_lp1, h_l); for_each_meshinterval( fine_intersect, [&](auto fine_mesh_interval) { - mesh_interval_t coarse_mesh_interval(level, fine_mesh_interval.i / 2, fine_mesh_interval.index / 2); + mesh_interval_t coarse_mesh_interval(level, fine_mesh_interval.i >> 1, fine_mesh_interval.index >> 1); minus_comput_stencil_it.init(fine_mesh_interval); coarse_it.init(coarse_mesh_interval); for (std::size_t ii = 0; ii < fine_mesh_interval.i.size(); ++ii) { - std::array interface_cells; + std::array interface_cells; interface_cells[0] = minus_comput_stencil_it.cells()[minus_direction_index]; interface_cells[1] = coarse_it.cells()[0]; - f(interface_cells, minus_comput_stencil_it.cells(), cell1_coeffs, cell2_coeffs); + f(interface_cells, minus_comput_stencil_it.cells(), left_cell_coeffs, right_cell_coeffs); minus_comput_stencil_it.move_next(); if (ii % 2 == 1) diff --git a/include/samurai/numeric/gauss_legendre.hpp b/include/samurai/numeric/gauss_legendre.hpp index 302778824..00ad32c93 100644 --- a/include/samurai/numeric/gauss_legendre.hpp +++ b/include/samurai/numeric/gauss_legendre.hpp @@ -1,6 +1,6 @@ #pragma once #include "../cell.hpp" -#include "../utils.hpp" +#include "../static_algorithm.hpp" namespace samurai { diff --git a/include/samurai/petsc/block_assembly.hpp b/include/samurai/petsc/block_assembly.hpp index e92f82169..92d42646d 100644 --- a/include/samurai/petsc/block_assembly.hpp +++ b/include/samurai/petsc/block_assembly.hpp @@ -42,11 +42,11 @@ namespace samurai }); } - template - static Assembly to_assembly(OperatorType& op) - { - return Assembly(op); - } + // template + // static Assembly to_assembly(OperatorType& op) + // { + // return Assembly(op); + // } auto& block_operator() { diff --git a/include/samurai/petsc/fv/FV_scheme_assembly.hpp b/include/samurai/petsc/fv/FV_scheme_assembly.hpp index a89b24ecb..abb9b6977 100644 --- a/include/samurai/petsc/fv/FV_scheme_assembly.hpp +++ b/include/samurai/petsc/fv/FV_scheme_assembly.hpp @@ -33,6 +33,7 @@ namespace samurai using mesh_id_t = typename mesh_t::mesh_id_t; using interval_t = typename mesh_t::interval_t; using field_value_type = typename field_t::value_type; // double + using coord_index_t = typename interval_t::coord_index_t; static constexpr std::size_t dim = field_t::dim; static constexpr std::size_t field_size = field_t::size; static constexpr std::size_t output_field_size = cfg_t::output_field_size; @@ -41,12 +42,15 @@ namespace samurai static constexpr std::size_t bdry_stencil_size = bdry_cfg_t::stencil_size; static constexpr std::size_t nb_bdry_ghosts = bdry_cfg_t::nb_ghosts; static constexpr DirichletEnforcement dirichlet_enfcmt = bdry_cfg_t::dirichlet_enfcmt; + using cell_t = samurai::Cell; using dirichlet_t = Dirichlet; using neumann_t = Neumann; using directional_bdry_config_t = DirectionalBoundaryConfig; + static constexpr bool ghost_elimination_enabled = false; + protected: const Scheme* m_scheme; @@ -55,6 +59,11 @@ namespace samurai InsertMode m_current_insert_mode = INSERT_VALUES; std::vector m_is_row_empty; + // Ghost recursion + using cell_coeff_pair_t = std::pair; + using CellLinearCombination = std::vector; + std::map m_ghost_recursion; + public: explicit FVSchemeAssembly(const Scheme& scheme) @@ -62,7 +71,7 @@ namespace samurai , m_unknown(&scheme.unknown()) { this->set_name(scheme.name()); - reset(); + // reset(); } auto& scheme() const @@ -76,6 +85,72 @@ namespace samurai // std::cout << "reset " << this->name() << ", rows = " << matrix_rows() << std::endl; m_is_row_empty.resize(static_cast(matrix_rows())); std::fill(m_is_row_empty.begin(), m_is_row_empty.end(), true); + + if constexpr (ghost_elimination_enabled) + { + m_ghost_recursion = ghost_recursion(); + } + } + + auto ghost_recursion() + { + static constexpr PetscInt number_of_children = (1 << dim); + + std::map recursion; + + for_each_projection_ghost_and_children_cells( + mesh(), + [&](auto, std::size_t ghost, const std::array& children) + { + CellLinearCombination linear_comb; + for (auto child : children) + { + linear_comb.emplace_back(child, 1. / number_of_children); + } + recursion.insert({ghost, linear_comb}); + }); + + for_each_prediction_ghost(mesh(), + [&](auto& ghost) + { + recursion.insert({ghost.index, prediction_linear_combination(ghost)}); + }); + + bool changes = true; + while (changes) + { + changes = false; + for (auto& pair : recursion) + { + // auto ghost = pair.first; + auto& linear_comb = pair.second; + + CellLinearCombination new_linear_comb; + + // for (auto& c : linear_comb) + for (auto& [cell, coeff] : linear_comb) + { + auto it = recursion.find(cell); + if (it == recursion.end()) // it's a cell + { + new_linear_comb.emplace_back(cell, coeff); + } + else // it's a ghost + { + auto& ghost_lin_comb = it->second; + for (auto& [cell2, coeff2] : ghost_lin_comb) + { + new_linear_comb.emplace_back(cell2, coeff * coeff2); + } + changes = true; + } + } + + linear_comb = new_linear_comb; + } + } + + return recursion; } void set_unknown(field_t& unknown) @@ -580,7 +655,15 @@ namespace samurai { for (unsigned int field_i = 0; field_i < output_field_size; ++field_i) { - nnz[static_cast(row_index(ghost, field_i))] = proj_stencil_size; + if constexpr (ghost_elimination_enabled) + { + nnz[static_cast(row_index(ghost, field_i))] = static_cast( + m_ghost_recursion.at(ghost.index).size()); + } + else + { + nnz[static_cast(row_index(ghost, field_i))] = proj_stencil_size; + } } }); } @@ -595,11 +678,19 @@ namespace samurai static constexpr std::size_t pred_stencil_size = 1 + ce_pow(2 * prediction_order + 1, dim); for_each_prediction_ghost(mesh(), - [&](auto& ghost) + [&](const auto& ghost) { for (unsigned int field_i = 0; field_i < output_field_size; ++field_i) { - nnz[static_cast(row_index(ghost, field_i))] = pred_stencil_size; + if constexpr (ghost_elimination_enabled) + { + nnz[static_cast(row_index(ghost, field_i))] = static_cast( + m_ghost_recursion.at(ghost.index).size()); + } + else + { + nnz[static_cast(row_index(ghost, field_i))] = pred_stencil_size; + } } }); } @@ -629,51 +720,67 @@ namespace samurai public: - void assemble_projection(Mat& A) override + void assemble_projection([[maybe_unused]] Mat& A) override { - static constexpr PetscInt number_of_children = (1 << dim); - - // double scalar = 1; - // if constexpr (is_FluxBasedScheme_v) - // { - // scalar = scheme().scalar(); - // } + if constexpr (!ghost_elimination_enabled) + { + static constexpr PetscInt number_of_children = (1 << dim); - for_each_projection_ghost_and_children_cells( - mesh(), - [&](auto level, PetscInt ghost, const std::array& children) - { - double h = cell_length(level); - double scaling = 1. / (h * h); - for (unsigned int field_i = 0; field_i < output_field_size; ++field_i) + for_each_projection_ghost_and_children_cells( + mesh(), + [&](auto level, PetscInt ghost, const std::array& children) { - PetscInt ghost_index = row_index(ghost, field_i); - MatSetValue(A, ghost_index, ghost_index, scaling, INSERT_VALUES); - for (unsigned int i = 0; i < number_of_children; ++i) + double h = cell_length(level); + double scaling = 1. / (h * h); + for (unsigned int field_i = 0; field_i < output_field_size; ++field_i) { - MatSetValue(A, ghost_index, col_index(children[i], field_i), -scaling / number_of_children, INSERT_VALUES); + PetscInt ghost_index = row_index(ghost, field_i); + MatSetValue(A, ghost_index, ghost_index, scaling, INSERT_VALUES); + for (unsigned int i = 0; i < number_of_children; ++i) + { + MatSetValue(A, ghost_index, col_index(children[i], field_i), -scaling / number_of_children, INSERT_VALUES); + } + set_is_row_not_empty(ghost_index); } - set_is_row_not_empty(ghost_index); - } - }); + }); + } } - void assemble_prediction(Mat& A) override + void assemble_prediction([[maybe_unused]] Mat& A) override { - static_assert(dim >= 1 && dim <= 3, - "assemble_prediction() is not implemented for " - "this dimension."); + if constexpr (!ghost_elimination_enabled) + { + static_assert(dim >= 1 && dim <= 3, "assemble_prediction() is not implemented for this dimension."); + if constexpr (dim == 1) + { + assemble_prediction_1D(A); + } + else if constexpr (dim == 2) + { + assemble_prediction_2D(A); + } + else if constexpr (dim == 3) + { + assemble_prediction_3D(A); + } + } + } + + template + auto prediction_linear_combination(const Cell& ghost) + { + static_assert(dim >= 1 && dim <= 3, "prediction_linear_combination() is not implemented for this dimension."); if constexpr (dim == 1) { - assemble_prediction_1D(A); + return prediction_linear_combination_1D(ghost); } else if constexpr (dim == 2) { - assemble_prediction_2D(A); + return prediction_linear_combination_2D(ghost); } else if constexpr (dim == 3) { - assemble_prediction_3D(A); + return prediction_linear_combination_3D(ghost); } } @@ -701,7 +808,7 @@ namespace samurai MatSetValue(A, ghost_index, ghost_index, scaling, INSERT_VALUES); auto ii = ghost.indices(0); - auto ig = ii / 2; + auto ig = ii >> 1; double isign = (ii & 1) ? -1 : 1; auto interpx = samurai::interp_coeffs<2 * prediction_order + 1>(isign); @@ -726,6 +833,33 @@ namespace samurai }); } + template + auto prediction_linear_combination_1D(const Cell& ghost) + { + using index_t = int; + + std::vector linear_comb; + + auto ii = ghost.indices(0); + auto ig = ii >> 1; + double isign = (ii & 1) ? -1 : 1; + + auto interpx = samurai::interp_coeffs<2 * prediction_order + 1>(isign); + + auto parent_index = this->mesh().get_index(ghost.level - 1, ig); + linear_comb.emplace_back(parent_index, 1.); + + for (std::size_t ci = 0; ci < interpx.size(); ++ci) + { + if (ci != prediction_order) + { + auto coarse_cell_index = this->mesh().get_index(ghost.level - 1, ig + static_cast(ci - prediction_order)); + linear_comb.emplace_back(coarse_cell_index, interpx[ci]); + } + } + return linear_comb; + } + void assemble_prediction_2D(Mat& A) { using index_t = int; @@ -741,9 +875,9 @@ namespace samurai MatSetValue(A, ghost_index, ghost_index, scaling, INSERT_VALUES); auto ii = ghost.indices(0); - auto ig = ii / 2; + auto ig = ii >> 1; auto j = ghost.indices(1); - auto jg = j / 2; + auto jg = j >> 1; double isign = (ii & 1) ? -1 : 1; double jsign = (j & 1) ? -1 : 1; @@ -775,6 +909,42 @@ namespace samurai }); } + template + auto prediction_linear_combination_2D(const Cell& ghost) + { + using index_t = int; + + std::vector linear_comb; + + auto ii = ghost.indices(0); + auto ig = ii >> 1; + auto j = ghost.indices(1); + auto jg = j >> 1; + double isign = (ii & 1) ? -1 : 1; + double jsign = (j & 1) ? -1 : 1; + + auto interpx = samurai::interp_coeffs<2 * prediction_order + 1>(isign); + auto interpy = samurai::interp_coeffs<2 * prediction_order + 1>(jsign); + + auto parent_index = this->mesh().get_index(ghost.level - 1, ig, jg); + linear_comb.emplace_back(parent_index, 1.); + + for (std::size_t ci = 0; ci < interpx.size(); ++ci) + { + for (std::size_t cj = 0; cj < interpy.size(); ++cj) + { + if (ci != prediction_order || cj != prediction_order) + { + auto coarse_cell_index = this->mesh().get_index(ghost.level - 1, + ig + static_cast(ci - prediction_order), + jg + static_cast(cj - prediction_order)); + linear_comb.emplace_back(coarse_cell_index, interpx[ci] * interpy[cj]); + } + } + } + return linear_comb; + } + void assemble_prediction_3D(Mat& A) { using index_t = int; @@ -790,11 +960,11 @@ namespace samurai MatSetValue(A, ghost_index, ghost_index, scaling, INSERT_VALUES); auto ii = ghost.indices(0); - auto ig = ii / 2; + auto ig = ii >> 1; auto j = ghost.indices(1); - auto jg = j / 2; + auto jg = j >> 1; auto k = ghost.indices(2); - auto kg = k / 2; + auto kg = k >> 1; double isign = (ii & 1) ? -1 : 1; double jsign = (j & 1) ? -1 : 1; double ksign = (k & 1) ? -1 : 1; @@ -832,6 +1002,50 @@ namespace samurai }); } + template + auto prediction_linear_combination_3D(const Cell& ghost) + { + using index_t = int; + + std::vector linear_comb; + + auto ii = ghost.indices(0); + auto ig = ii >> 1; + auto j = ghost.indices(1); + auto jg = j >> 1; + auto k = ghost.indices(2); + auto kg = k >> 1; + double isign = (ii & 1) ? -1 : 1; + double jsign = (j & 1) ? -1 : 1; + double ksign = (k & 1) ? -1 : 1; + + auto interpx = samurai::interp_coeffs<2 * prediction_order + 1>(isign); + auto interpy = samurai::interp_coeffs<2 * prediction_order + 1>(jsign); + auto interpz = samurai::interp_coeffs<2 * prediction_order + 1>(ksign); + + auto parent_index = this->mesh().get_index(ghost.level - 1, ig, jg, kg); + linear_comb.emplace_back(parent_index, 1.); + + for (std::size_t ci = 0; ci < interpx.size(); ++ci) + { + for (std::size_t cj = 0; cj < interpy.size(); ++cj) + { + for (std::size_t ck = 0; ck < interpz.size(); ++ck) + { + if (ci != prediction_order || cj != prediction_order || ck != prediction_order) + { + auto coarse_cell_index = this->mesh().get_index(ghost.level - 1, + ig + static_cast(ci - prediction_order), + jg + static_cast(cj - prediction_order), + kg + static_cast(ck - prediction_order)); + linear_comb.emplace_back(coarse_cell_index, interpx[ci] * interpy[cj] * interpz[ck]); + } + } + } + } + return linear_comb; + } + public: bool matrix_is_symmetric() const override diff --git a/include/samurai/petsc/fv/flux_based_scheme_assembly.hpp b/include/samurai/petsc/fv/flux_based_scheme_assembly.hpp index 0d2af726b..08ca210c4 100644 --- a/include/samurai/petsc/fv/flux_based_scheme_assembly.hpp +++ b/include/samurai/petsc/fv/flux_based_scheme_assembly.hpp @@ -21,6 +21,7 @@ namespace samurai public: + using base_class::ghost_elimination_enabled; using base_class::mesh; using base_class::scheme; using base_class::set_current_insert_mode; @@ -60,14 +61,39 @@ namespace samurai mesh(), scheme_coeffs_dir.flux.direction, scheme_coeffs_dir.flux.stencil, - [&](auto& interface_cells, auto&) + [&](auto& interface_cells, auto& comput_cells) { for (unsigned int field_i = 0; field_i < output_field_size; ++field_i) { for (unsigned int field_j = 0; field_j < field_size; ++field_j) { - nnz[static_cast(this->row_index(interface_cells[0], field_i))] += stencil_size * field_size; - nnz[static_cast(this->row_index(interface_cells[1], field_i))] += stencil_size * field_size; + if constexpr (ghost_elimination_enabled) + { + for (std::size_t c = 0; c < stencil_size; ++c) + { + auto it_ghost = this->m_ghost_recursion.find(comput_cells[c].index); + if (it_ghost == this->m_ghost_recursion.end()) + { + nnz[static_cast(this->row_index(interface_cells[0], field_i))] += field_size; + nnz[static_cast(this->row_index(interface_cells[1], field_i))] += field_size; + } + else + { + auto& linear_comb = it_ghost->second; + nnz[static_cast(this->row_index(interface_cells[0], field_i))] += linear_comb.size() + * field_size; + nnz[static_cast(this->row_index(interface_cells[1], field_i))] += linear_comb.size() + * field_size; + } + } + } + else + { + nnz[static_cast(this->row_index(interface_cells[0], field_i))] += stencil_size + * field_size; + nnz[static_cast(this->row_index(interface_cells[1], field_i))] += stencil_size + * field_size; + } } } }); @@ -132,35 +158,67 @@ namespace samurai scheme_coeffs_dir.flux.direction, scheme_coeffs_dir.flux.stencil, scheme_coeffs_dir.flux.get_flux_coeffs, - scheme_coeffs_dir.get_cell1_coeffs, - scheme_coeffs_dir.get_cell2_coeffs, - [&](auto& interface_cells, auto& comput_cells, auto& cell1_coeffs, auto& cell2_coeffs) + scheme_coeffs_dir.get_left_cell_coeffs, + scheme_coeffs_dir.get_right_cell_coeffs, + [&](auto& interface_cells, auto& comput_cells, auto& left_cell_coeffs, auto& right_cell_coeffs) { for (unsigned int field_i = 0; field_i < output_field_size; ++field_i) { - auto interface_cell1_row = this->row_index(interface_cells[0], field_i); - auto interface_cell2_row = this->row_index(interface_cells[1], field_i); + auto left_cell_row = this->row_index(interface_cells[0], field_i); + auto right_cell_row = this->row_index(interface_cells[1], field_i); for (unsigned int field_j = 0; field_j < field_size; ++field_j) { for (std::size_t c = 0; c < stencil_size; ++c) { - auto comput_cell_col = col_index(comput_cells[c], field_j); - double cell1_coeff = scheme().cell_coeff(cell1_coeffs, c, field_i, field_j); - double cell2_coeff = scheme().cell_coeff(cell2_coeffs, c, field_i, field_j); - // if (cell1_coeff != 0) - // { - MatSetValue(A, interface_cell1_row, comput_cell_col, cell1_coeff, ADD_VALUES); - // } - // if (cell2_coeff != 0) - // { - MatSetValue(A, interface_cell2_row, comput_cell_col, cell2_coeff, ADD_VALUES); - // } - // MatSetValue(A, interface_cell1_row, interface_cell1_row, 0, ADD_VALUES); - // MatSetValue(A, interface_cell2_row, interface_cell2_row, 0, ADD_VALUES); + double left_cell_coeff = scheme().cell_coeff(left_cell_coeffs, c, field_i, field_j); + double right_cell_coeff = scheme().cell_coeff(right_cell_coeffs, c, field_i, field_j); + + if constexpr (ghost_elimination_enabled) + { + auto it_ghost = this->m_ghost_recursion.find(comput_cells[c].index); + if (it_ghost == this->m_ghost_recursion.end()) + { + auto comput_cell_col = col_index(comput_cells[c], field_j); + // if (left_cell_coeff != 0) + // { + MatSetValue(A, left_cell_row, comput_cell_col, left_cell_coeff, ADD_VALUES); + // } + // if (right_cell_coeff != 0) + // { + MatSetValue(A, right_cell_row, comput_cell_col, right_cell_coeff, ADD_VALUES); + // } + // MatSetValue(A, left_cell_row, left_cell_row, 0, ADD_VALUES); + // MatSetValue(A, right_cell_row, right_cell_row, 0, ADD_VALUES); + } + else + { + auto& linear_comb = it_ghost->second; + for (auto& [cell, coeff] : linear_comb) + { + auto comput_cell_col = col_index(static_cast(cell), field_j); + MatSetValue(A, left_cell_row, comput_cell_col, left_cell_coeff * coeff, ADD_VALUES); + MatSetValue(A, right_cell_row, comput_cell_col, right_cell_coeff * coeff, ADD_VALUES); + } + } + } + else + { + auto comput_cell_col = col_index(comput_cells[c], field_j); + // if (left_cell_coeff != 0) + // { + MatSetValue(A, left_cell_row, comput_cell_col, left_cell_coeff, ADD_VALUES); + // } + // if (right_cell_coeff != 0) + // { + MatSetValue(A, right_cell_row, comput_cell_col, right_cell_coeff, ADD_VALUES); + // } + // MatSetValue(A, left_cell_row, left_cell_row, 0, ADD_VALUES); + // MatSetValue(A, right_cell_row, right_cell_row, 0, ADD_VALUES); + } } } - set_is_row_not_empty(interface_cell1_row); - set_is_row_not_empty(interface_cell2_row); + set_is_row_not_empty(left_cell_row); + set_is_row_not_empty(right_cell_row); } }); @@ -168,7 +226,7 @@ namespace samurai scheme_coeffs_dir.flux.direction, scheme_coeffs_dir.flux.stencil, scheme_coeffs_dir.flux.get_flux_coeffs, - scheme_coeffs_dir.get_cell1_coeffs, + scheme_coeffs_dir.get_left_cell_coeffs, [&](auto& interface_cells, auto& comput_cells, auto& coeffs) { for (unsigned int field_i = 0; field_i < output_field_size; ++field_i) @@ -197,7 +255,7 @@ namespace samurai opposite_direction, opposite_stencil, scheme_coeffs_dir.flux.get_flux_coeffs, - scheme_coeffs_dir.get_cell2_coeffs, + scheme_coeffs_dir.get_right_cell_coeffs, [&](auto& interface_cells, auto& comput_cells, auto& coeffs) { for (unsigned int field_i = 0; field_i < output_field_size; ++field_i) diff --git a/include/samurai/schemes/fv/explicit_flux_based_scheme.hpp b/include/samurai/schemes/fv/explicit_flux_based_scheme.hpp index 085157130..f0689a60c 100644 --- a/include/samurai/schemes/fv/explicit_flux_based_scheme.hpp +++ b/include/samurai/schemes/fv/explicit_flux_based_scheme.hpp @@ -1,4 +1,5 @@ #pragma once +// #include "../../petsc/fv/flux_based_scheme_assembly.hpp" #include "../explicit_scheme.hpp" #include "flux_based_scheme.hpp" @@ -42,11 +43,11 @@ namespace samurai * Implementation by matrix-vector multiplication */ // Mat A; - // auto assembly = make_assembly(scheme()); + // auto assembly = petsc::make_assembly(scheme()); // assembly.create_matrix(A); // assembly.assemble_matrix(A); - // Vec vec_f = create_petsc_vector_from(f); - // Vec vec_res = create_petsc_vector_from(result); + // Vec vec_f = petsc::create_petsc_vector_from(f); + // Vec vec_res = petsc::create_petsc_vector_from(result); // MatMult(A, vec_f, vec_res); // return result; @@ -59,9 +60,9 @@ namespace samurai scheme_coeffs_dir.flux.direction, scheme_coeffs_dir.flux.stencil, scheme_coeffs_dir.flux.get_flux_coeffs, - scheme_coeffs_dir.get_cell1_coeffs, - scheme_coeffs_dir.get_cell2_coeffs, - [&](auto& interface_cells, auto& comput_cells, auto& cell1_coeffs, auto& cell2_coeffs) + scheme_coeffs_dir.get_left_cell_coeffs, + scheme_coeffs_dir.get_right_cell_coeffs, + [&](auto& interface_cells, auto& comput_cells, auto& left_cell_coeffs, auto& right_cell_coeffs) { for (std::size_t field_i = 0; field_i < output_field_size; ++field_i) { @@ -69,11 +70,11 @@ namespace samurai { for (std::size_t c = 0; c < stencil_size; ++c) { - double cell1_coeff = scheme().cell_coeff(cell1_coeffs, c, field_i, field_j); - double cell2_coeff = scheme().cell_coeff(cell2_coeffs, c, field_i, field_j); - field_value(result, interface_cells[0], field_i) += cell1_coeff + double left_cell_coeff = scheme().cell_coeff(left_cell_coeffs, c, field_i, field_j); + double right_cell_coeff = scheme().cell_coeff(right_cell_coeffs, c, field_i, field_j); + field_value(result, interface_cells[0], field_i) += left_cell_coeff * field_value(f, comput_cells[c], field_j); - field_value(result, interface_cells[1], field_i) += cell2_coeff + field_value(result, interface_cells[1], field_i) += right_cell_coeff * field_value(f, comput_cells[c], field_j); } } @@ -85,7 +86,7 @@ namespace samurai scheme_coeffs_dir.flux.direction, scheme_coeffs_dir.flux.stencil, scheme_coeffs_dir.flux.get_flux_coeffs, - scheme_coeffs_dir.get_cell1_coeffs, + scheme_coeffs_dir.get_left_cell_coeffs, [&](auto& interface_cells, auto& comput_cells, auto& coeffs) { for (std::size_t field_i = 0; field_i < output_field_size; ++field_i) @@ -109,7 +110,7 @@ namespace samurai opposite_direction, opposite_stencil, scheme_coeffs_dir.flux.get_flux_coeffs, - scheme_coeffs_dir.get_cell2_coeffs, + scheme_coeffs_dir.get_right_cell_coeffs, [&](auto& interface_cells, auto& comput_cells, auto& coeffs) { for (std::size_t field_i = 0; field_i < output_field_size; ++field_i) diff --git a/include/samurai/schemes/fv/flux_based_scheme.hpp b/include/samurai/schemes/fv/flux_based_scheme.hpp index 83a993a34..5b7d9e0a1 100644 --- a/include/samurai/schemes/fv/flux_based_scheme.hpp +++ b/include/samurai/schemes/fv/flux_based_scheme.hpp @@ -36,7 +36,7 @@ namespace samurai }; template - auto normal_grad_order2(Vector& direction) + auto normal_grad_order1(Vector& direction) { static constexpr std::size_t dim = Field::dim; static constexpr std::size_t field_size = Field::size; @@ -65,7 +65,7 @@ namespace samurai } template - auto normal_grad_order2() + auto normal_grad_order1() { static constexpr std::size_t dim = Field::dim; using flux_computation_t = NormalFluxComputation; @@ -74,7 +74,7 @@ namespace samurai std::array normal_fluxes; for (std::size_t d = 0; d < dim; ++d) { - normal_fluxes[d] = normal_grad_order2(xt::view(directions, d)); + normal_fluxes[d] = normal_grad_order1(xt::view(directions, d)); } return normal_fluxes; } @@ -90,10 +90,11 @@ namespace samurai using coeff_matrix_t = typename detail::LocalMatrix::Type; using cell_coeffs_t = std::array; using flux_coeffs_t = typename flux_computation_t::flux_coeffs_t; // std::array; + using cell_coeffs_func_t = std::function; flux_computation_t flux; - std::function get_cell1_coeffs; - std::function get_cell2_coeffs; + cell_coeffs_func_t get_left_cell_coeffs; + cell_coeffs_func_t get_right_cell_coeffs; }; /** diff --git a/include/samurai/schemes/fv/operators/diffusion_FV.hpp b/include/samurai/schemes/fv/operators/diffusion_FV.hpp index e5f9c831d..52b45535d 100644 --- a/include/samurai/schemes/fv/operators/diffusion_FV.hpp +++ b/include/samurai/schemes/fv/operators/diffusion_FV.hpp @@ -52,7 +52,7 @@ namespace samurai * Conclusion: the contribution of the face is just the flux received as a parameter, multiplied by |F|/|T|. */ template - static auto get_laplacian_coeffs_cell1(std::array& flux_coeffs, double h_face, double h_cell) + static auto get_laplacian_coeffs_left_cell(std::array& flux_coeffs, double h_face, double h_cell) { double face_measure = pow(h_face, dim - 1); double cell_measure = pow(h_cell, dim); @@ -94,18 +94,19 @@ namespace samurai /** * |-------|-------| - * | cell1 | cell2 | + * | left | right | + * | cell | cell | * |-------|-------| * -------> * normal flux */ - // How the flux is computed in this direction: here, Grad.n = (u2-u1)/h - coeffs.flux = normal_grad_order2(direction); - // Coefficients of the scheme for cell1, in function of the flux - coeffs.get_cell1_coeffs = [](std::array& flux_coeffs, double h_face, double h_cell) + // How the flux is computed in this direction: here, Grad.n = (uR-uL)/h + coeffs.flux = normal_grad_order1(direction); + // Coefficients of the scheme for the left cell, in function of the flux + coeffs.get_left_cell_coeffs = [](std::array& flux_coeffs, double h_face, double h_cell) { - auto cell_coeffs = get_laplacian_coeffs_cell1(flux_coeffs, h_face, h_cell); + auto cell_coeffs = get_laplacian_coeffs_left_cell(flux_coeffs, h_face, h_cell); // We multiply by -1 because we implement the operator -Lap for (auto& coeff : cell_coeffs) { @@ -113,8 +114,8 @@ namespace samurai } return cell_coeffs; }; - // Coefficients of the scheme for cell2, in function of the flux - coeffs.get_cell2_coeffs = get_laplacian_coeffs_cell1; + // Coefficients of the scheme for the right, in function of the flux + coeffs.get_right_cell_coeffs = get_laplacian_coeffs_left_cell; }); return coeffs_by_fluxes; } diff --git a/include/samurai/schemes/fv/operators/divergence_FV.hpp b/include/samurai/schemes/fv/operators/divergence_FV.hpp index 9c214582b..9948eab52 100644 --- a/include/samurai/schemes/fv/operators/divergence_FV.hpp +++ b/include/samurai/schemes/fv/operators/divergence_FV.hpp @@ -106,9 +106,9 @@ namespace samurai auto& coeffs = coeffs_by_fluxes[d]; DirectionVector direction = xt::view(directions, d); - coeffs.flux = normal_grad_order2(direction); - coeffs.get_cell1_coeffs = average; - coeffs.get_cell2_coeffs = minus_average; + coeffs.flux = normal_grad_order1(direction); + coeffs.get_left_cell_coeffs = average; + coeffs.get_right_cell_coeffs = minus_average; }); return coeffs_by_fluxes; } diff --git a/include/samurai/schemes/fv/operators/gradient_FV.hpp b/include/samurai/schemes/fv/operators/gradient_FV.hpp index 9f99d7f59..0992544cc 100644 --- a/include/samurai/schemes/fv/operators/gradient_FV.hpp +++ b/include/samurai/schemes/fv/operators/gradient_FV.hpp @@ -89,9 +89,9 @@ namespace samurai auto& coeffs = coeffs_by_fluxes[d]; DirectionVector direction = xt::view(directions, d); - coeffs.flux = normal_grad_order2(direction); - coeffs.get_cell1_coeffs = half_flux_in_direction; - coeffs.get_cell2_coeffs = half_flux_in_direction; + coeffs.flux = normal_grad_order1(direction); + coeffs.get_left_cell_coeffs = half_flux_in_direction; + coeffs.get_right_cell_coeffs = half_flux_in_direction; }); return coeffs_by_fluxes; } diff --git a/include/samurai/schemes/fv/operators/identity_FV.hpp b/include/samurai/schemes/fv/operators/identity_FV.hpp index 8922860bb..6b7afcdef 100644 --- a/include/samurai/schemes/fv/operators/identity_FV.hpp +++ b/include/samurai/schemes/fv/operators/identity_FV.hpp @@ -67,22 +67,22 @@ namespace samurai if (d == 0) { flux.get_flux_coeffs = flux_coefficients<0>; - flux.get_cell1_coeffs = [&](auto&, double, double) + flux.get_left_cell_coeffs = [&](auto&, double, double) { CellCoeffs coeffs; coeffs[0] = eye(); coeffs[1] = zeros(); return coeffs; }; - flux.get_cell2_coeffs = get_zero_coeffs; + flux.get_right_cell_coeffs = get_zero_coeffs; } if constexpr (dim >= 2) { if (d == 1) { flux.get_flux_coeffs = flux_coefficients<1>; - flux.get_cell1_coeffs = get_zero_coeffs; - flux.get_cell2_coeffs = get_zero_coeffs; + flux.get_left_cell_coeffs = get_zero_coeffs; + flux.get_right_cell_coeffs = get_zero_coeffs; } } if constexpr (dim >= 3) @@ -90,8 +90,8 @@ namespace samurai if (d == 2) { flux.get_flux_coeffs = flux_coefficients<2>; - flux.get_cell1_coeffs = get_zero_coeffs; - flux.get_cell2_coeffs = get_zero_coeffs; + flux.get_left_cell_coeffs = get_zero_coeffs; + flux.get_right_cell_coeffs = get_zero_coeffs; } } } diff --git a/include/samurai/schemes/fv/scheme_operators.hpp b/include/samurai/schemes/fv/scheme_operators.hpp index c5dfd8b8b..89203092d 100644 --- a/include/samurai/schemes/fv/scheme_operators.hpp +++ b/include/samurai/schemes/fv/scheme_operators.hpp @@ -61,12 +61,12 @@ namespace samurai static_for<0, dim>::apply( [&](auto integral_constant_d) { - static constexpr int d = decltype(integral_constant_d)::value; - scalar_x_fluxes[d].get_cell1_coeffs = [&](auto& flux_coeffs, double h_face, double h_cell) + static constexpr int d = decltype(integral_constant_d)::value; + scalar_x_fluxes[d].get_left_cell_coeffs = [&](auto& flux_coeffs, double h_face, double h_cell) { return this->scalar_x_get_cell1_coeffs(flux_coeffs, h_face, h_cell); }; - scalar_x_fluxes[d].get_cell2_coeffs = [&](auto& flux_coeffs, double h_face, double h_cell) + scalar_x_fluxes[d].get_right_cell_coeffs = [&](auto& flux_coeffs, double h_face, double h_cell) { return this->scalar_x_get_cell2_coeffs(flux_coeffs, h_face, h_cell); }; @@ -77,7 +77,7 @@ namespace samurai template auto scalar_x_get_cell1_coeffs(flux_coeffs_t& flux_coeffs, double h_face, double h_cell) const { - auto coeffs = m_scheme.coefficients()[d].get_cell1_coeffs(flux_coeffs, h_face, h_cell); + auto coeffs = m_scheme.coefficients()[d].get_left_cell_coeffs(flux_coeffs, h_face, h_cell); for (auto& coeff : coeffs) { coeff *= m_scalar; @@ -88,7 +88,7 @@ namespace samurai template auto scalar_x_get_cell2_coeffs(flux_coeffs_t& flux_coeffs, double h_face, double h_cell) const { - auto coeffs = m_scheme.coefficients()[d].get_cell2_coeffs(flux_coeffs, h_face, h_cell); + auto coeffs = m_scheme.coefficients()[d].get_right_cell_coeffs(flux_coeffs, h_face, h_cell); for (auto& coeff : coeffs) { coeff *= m_scalar;