Skip to content

Commit

Permalink
fix: Fix some minor issues in the flux-based schemes (#217)
Browse files Browse the repository at this point in the history 
## Description

- Fix the `for_each_boundary_interface()` function with minimal
parameters (it wasn't compiling).
- Fix the PETSc insert mode of the projection/prediction coefficients
(implicit context).
- The scheme is now copied into the `Assembly` object (instead of
storing a reference).
- Possibility to disable the computation of boundary fluxes.

## How has this been tested?
Very quickly.

## Code of Conduct
By submitting this PR, you agree to follow our [Code of
Conduct](https://github.com/hpc-maths/samurai/blob/master/docs/CODE_OF_CONDUCT.md)
- [x] I agree to follow this project's Code of Conduct
  • Loading branch information
@pmatalon
pmatalon authored Jun 18, 2024
1 parent 51ad39c commit c15b06d
Show file tree
Hide file tree
Showing 14 changed files with 136 additions and 67 deletions.
2 changes: 2 additions & 0 deletions conda/mpi-environment.yml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -13,4 +13,6 @@ dependencies:
- pytest
- h5py
- libboost-mpi
- libboost-devel
- libboost-headers
- hdf5=*=mpi*
2 changes: 1 addition & 1 deletion demos/FiniteVolume/manual_block_matrix_assembly.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -237,6 +237,7 @@ int main(int argc, char* argv[])
// Diffusion operators for the electrolyte and the solid
double D_e = 1; // diffusion coefficient
auto diff_e = samurai::make_diffusion_order2<decltype(u_e)>(D_e);
diff_e.include_boundary_fluxes(false);
double D_s = 2;
auto diff_s = samurai::make_diffusion_order2<decltype(u_s)>(D_s);
auto id = samurai::make_identity<decltype(u_e)>();
Expand Down Expand Up @@ -265,7 +266,6 @@ int main(int argc, char* argv[])
// Disable the assembly of the BC for the diffusion operators
assembly.get<0, 0>().include_bc(false);
assembly.get<2, 2>().include_bc(false);
assembly.get<0, 0>().include_boundary_fluxes(false);
assembly.set_diag_value_for_useless_ghosts(9);

// Set the unknowns of the system (even if you don't want the solve it).
Expand Down
2 changes: 1 addition & 1 deletion include/samurai/interface.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -502,7 +502,7 @@ namespace samurai
DirectionVector<Mesh::dim> direction;
direction.fill(0);
direction[d] = 1;
for_each_boundary_interface<run_type, get_type>(mesh, direction, std::forward<Func>(f));
for_each_boundary_interface__both_directions<run_type, get_type>(mesh, direction, std::forward<Func>(f));
}
}

Expand Down
8 changes: 4 additions & 4 deletions include/samurai/petsc/block_assembly.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -19,13 +19,13 @@ namespace samurai

private:

const block_operator_t* m_block_operator;
block_operator_t m_block_operator;
std::tuple<Assembly<Operators>...> m_assembly_ops;

public:

explicit BlockAssembly(const block_operator_t& block_op)
: m_block_operator(&block_op)
: m_block_operator(block_op)
, m_assembly_ops(transform(block_op.operators(),
[](const auto& op)
{
Expand All @@ -50,12 +50,12 @@ namespace samurai

auto& block_operator()
{
return *m_block_operator;
return m_block_operator;
}

auto& block_operator() const
{
return *m_block_operator;
return m_block_operator;
}

template <class Func>
Expand Down
43 changes: 29 additions & 14 deletions include/samurai/petsc/fv/FV_scheme_assembly.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -55,7 +55,7 @@ namespace samurai

protected:

const Scheme* m_scheme;
Scheme m_scheme;
field_t* m_unknown = nullptr;
std::size_t m_n_cells = 0;
std::vector<bool> m_is_row_empty;
Expand All @@ -69,14 +69,14 @@ namespace samurai
public:

explicit FVSchemeAssembly(const Scheme& scheme)
: m_scheme(&scheme)
: m_scheme(scheme)
{
this->set_name(scheme.name());
}

auto& scheme() const
{
return *m_scheme;
return m_scheme;
}

void reset() override
Expand Down Expand Up @@ -623,7 +623,11 @@ namespace samurai
INSERT_VALUES);
if (error)
{
std::cerr << scheme().name() << ": failure to insert diagonal coefficient at ("
<< m_row_shift + static_cast<PetscInt>(i) << ", " << m_col_shift + static_cast<PetscInt>(i)
<< "), i.e. (" << i << ", " << i << ") in the block." << std::endl;
assert(false);
exit(EXIT_FAILURE);
}
// m_is_row_empty[i] = false;
}
Expand Down Expand Up @@ -756,10 +760,21 @@ namespace samurai
for (unsigned int field_i = 0; field_i < output_field_size; ++field_i)
{
PetscInt ghost_index = row_index(ghost, field_i);
MatSetValue(A, ghost_index, ghost_index, scaling, INSERT_VALUES);
MatSetValue(A, ghost_index, ghost_index, scaling, current_insert_mode());
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);
auto error = MatSetValue(A,
ghost_index,
col_index(children[i], field_i),
-scaling / number_of_children,
current_insert_mode());
if (error)
{
std::cerr << scheme().name() << ": failure to insert projection coefficient at (" << ghost_index
<< ", " << m_col_shift + col_index(children[i], field_i) << ")." << std::endl;
assert(false);
exit(EXIT_FAILURE);
}
}
set_is_row_not_empty(ghost_index);
}
Expand Down Expand Up @@ -824,7 +839,7 @@ namespace samurai
for (unsigned int field_i = 0; field_i < field_size; ++field_i)
{
PetscInt ghost_index = this->row_index(ghost, field_i);
MatSetValue(A, ghost_index, ghost_index, scaling, INSERT_VALUES);
MatSetValue(A, ghost_index, ghost_index, scaling, current_insert_mode());

auto ii = ghost.indices(0);
auto ig = ii >> 1;
Expand All @@ -833,7 +848,7 @@ namespace samurai
auto interpx = samurai::interp_coeffs<2 * prediction_order + 1>(isign);

auto parent_index = this->col_index(static_cast<PetscInt>(this->mesh().get_index(ghost.level - 1, ig)), field_i);
MatSetValue(A, ghost_index, parent_index, -scaling, INSERT_VALUES);
MatSetValue(A, ghost_index, parent_index, -scaling, current_insert_mode());

for (std::size_t ci = 0; ci < interpx.size(); ++ci)
{
Expand All @@ -844,7 +859,7 @@ namespace samurai
static_cast<PetscInt>(
this->mesh().get_index(ghost.level - 1, ig + static_cast<coord_index_t>(ci - prediction_order))),
field_i);
MatSetValue(A, ghost_index, coarse_cell_index, scaling * value, INSERT_VALUES);
MatSetValue(A, ghost_index, coarse_cell_index, scaling * value, current_insert_mode());
}
}
set_is_row_not_empty(ghost_index);
Expand Down Expand Up @@ -889,7 +904,7 @@ namespace samurai
for (unsigned int field_i = 0; field_i < field_size; ++field_i)
{
PetscInt ghost_index = this->row_index(ghost, field_i);
MatSetValue(A, ghost_index, ghost_index, scaling, INSERT_VALUES);
MatSetValue(A, ghost_index, ghost_index, scaling, current_insert_mode());

auto ii = ghost.indices(0);
auto ig = ii >> 1;
Expand All @@ -903,7 +918,7 @@ namespace samurai

auto parent_index = this->col_index(static_cast<PetscInt>(this->mesh().get_index(ghost.level - 1, ig, jg)),
field_i);
MatSetValue(A, ghost_index, parent_index, -scaling, INSERT_VALUES);
MatSetValue(A, ghost_index, parent_index, -scaling, current_insert_mode());

for (std::size_t ci = 0; ci < interpx.size(); ++ci)
{
Expand All @@ -917,7 +932,7 @@ namespace samurai
ig + static_cast<coord_index_t>(ci - prediction_order),
jg + static_cast<coord_index_t>(cj - prediction_order))),
field_i);
MatSetValue(A, ghost_index, coarse_cell_index, scaling * value, INSERT_VALUES);
MatSetValue(A, ghost_index, coarse_cell_index, scaling * value, current_insert_mode());
}
}
}
Expand Down Expand Up @@ -971,7 +986,7 @@ namespace samurai
for (unsigned int field_i = 0; field_i < field_size; ++field_i)
{
PetscInt ghost_index = this->row_index(ghost, field_i);
MatSetValue(A, ghost_index, ghost_index, scaling, INSERT_VALUES);
MatSetValue(A, ghost_index, ghost_index, scaling, current_insert_mode());

auto ii = ghost.indices(0);
auto ig = ii >> 1;
Expand All @@ -989,7 +1004,7 @@ namespace samurai

auto parent_index = this->col_index(static_cast<PetscInt>(this->mesh().get_index(ghost.level - 1, ig, jg, kg)),
field_i);
MatSetValue(A, ghost_index, parent_index, -scaling, INSERT_VALUES);
MatSetValue(A, ghost_index, parent_index, -scaling, current_insert_mode());

for (std::size_t ci = 0; ci < interpx.size(); ++ci)
{
Expand All @@ -1007,7 +1022,7 @@ namespace samurai
jg + static_cast<coord_index_t>(cj - prediction_order),
kg + static_cast<coord_index_t>(ck - prediction_order))),
field_i);
MatSetValue(A, ghost_index, coarse_cell_index, scaling * value, INSERT_VALUES);
MatSetValue(A, ghost_index, coarse_cell_index, scaling * value, current_insert_mode());
}
}
}
Expand Down
1 change: 1 addition & 0 deletions include/samurai/petsc/fv/flux_based_scheme_assembly.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -47,6 +47,7 @@ namespace samurai
: base_class(s)
{
set_current_insert_mode(ADD_VALUES);
m_include_boundary_fluxes = s.include_boundary_fluxes();
}

void include_boundary_fluxes(bool include)
Expand Down
1 change: 1 addition & 0 deletions include/samurai/petsc/zero_block_assembly.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -22,6 +22,7 @@ namespace samurai
exit(EXIT_FAILURE);
}
this->fit_block_dimensions(true);
this->set_name("0");
}

void sparsity_pattern_scheme(std::vector<PetscInt>&) const override
Expand Down
6 changes: 6 additions & 0 deletions include/samurai/schemes/block_operator.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -29,6 +29,12 @@ namespace samurai
return m_operators;
}

template <std::size_t row, std::size_t col>
auto& get()
{
return std::get<row * cols + col>(m_operators);
}

template <class Func>
void for_each_operator(Func&& f)
{
Expand Down
37 changes: 20 additions & 17 deletions include/samurai/schemes/fv/flux_based/explicit_flux_based_scheme__lin_het.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -74,30 +74,33 @@ namespace samurai
});

// Boundary interfaces
scheme().for_each_boundary_interface_and_coeffs(
input_field,
[&](const auto& cell, const auto& comput_cells, auto& coeffs)
{
for (std::size_t field_i = 0; field_i < output_field_size; ++field_i)
if (scheme().include_boundary_fluxes())
{
scheme().for_each_boundary_interface_and_coeffs(
input_field,
[&](const auto& cell, const auto& comput_cells, auto& coeffs)
{
for (std::size_t field_j = 0; field_j < field_size; ++field_j)
for (std::size_t field_i = 0; field_i < output_field_size; ++field_i)
{
for (std::size_t c = 0; c < stencil_size; ++c)
for (std::size_t field_j = 0; field_j < field_size; ++field_j)
{
#ifdef SAMURAI_CHECK_NAN
if (std::isnan(field_value(input_field, comput_cells[c], field_j)))
for (std::size_t c = 0; c < stencil_size; ++c)
{
std::cerr << "NaN detected when computing the flux on the boundary interfaces: " << comput_cells[c]
<< std::endl;
assert(false);
}
#ifdef SAMURAI_CHECK_NAN
if (std::isnan(field_value(input_field, comput_cells[c], field_j)))
{
std::cerr << "NaN detected when computing the flux on the boundary interfaces: " << comput_cells[c]
<< std::endl;
assert(false);
}
#endif
double coeff = this->scheme().cell_coeff(coeffs, c, field_i, field_j);
field_value(output_field, cell, field_i) += coeff * field_value(input_field, comput_cells[c], field_j);
double coeff = this->scheme().cell_coeff(coeffs, c, field_i, field_j);
field_value(output_field, cell, field_i) += coeff * field_value(input_field, comput_cells[c], field_j);
}
}
}
}
});
});
}
}
};

Expand Down
49 changes: 27 additions & 22 deletions include/samurai/schemes/fv/flux_based/explicit_flux_based_scheme__lin_hom.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -265,43 +265,48 @@ namespace samurai
});

// Boundary interfaces
scheme().template for_each_boundary_interface_and_coeffs<Run::Parallel, Get::Intervals>(
input_field,
[&](auto& cell, auto& stencil, auto& coeffs)
{
for (std::size_t field_i = 0; field_i < output_field_size; ++field_i)
if (scheme().include_boundary_fluxes())
{
scheme().template for_each_boundary_interface_and_coeffs<Run::Parallel, Get::Intervals>(
input_field,
[&](auto& cell, auto& stencil, auto& coeffs)
{
for (std::size_t field_j = 0; field_j < field_size; ++field_j)
for (std::size_t field_i = 0; field_i < output_field_size; ++field_i)
{
for (std::size_t c = 0; c < stencil_size; ++c)
for (std::size_t field_j = 0; field_j < field_size; ++field_j)
{
#ifdef SAMURAI_CHECK_NAN
if (std::isnan(field_value(input_field, stencil.cells()[c], field_j)))
for (std::size_t c = 0; c < stencil_size; ++c)
{
std::cerr << "NaN detected when computing the flux on the boundary interfaces: " << stencil.cells()[c]
<< std::endl;
assert(false);
}
#ifdef SAMURAI_CHECK_NAN
if (std::isnan(field_value(input_field, stencil.cells()[c], field_j)))
{
std::cerr
<< "NaN detected when computing the flux on the boundary interfaces: " << stencil.cells()[c]
<< std::endl;
assert(false);
}
#endif
auto coeff = this->scheme().cell_coeff(coeffs, c, field_i, field_j);
// field_value(output_field, cell, field_i) += coeff * field_value(input_field, stencil[c], field_j);
auto coeff = this->scheme().cell_coeff(coeffs, c, field_i, field_j);
// field_value(output_field, cell, field_i) += coeff * field_value(input_field, stencil[c], field_j);

auto cell_index_init = cell.index;
auto comput_index_init = stencil.cells()[c].index;
auto cell_index_init = cell.index;
auto comput_index_init = stencil.cells()[c].index;

using index_t = decltype(cell_index_init);
using index_t = decltype(cell_index_init);

// clang-format off

// clang-format off
#pragma omp simd
for (index_t ii = 0; ii < static_cast<index_t>(stencil.interval().size()); ++ii)
{
field_value(output_field, cell_index_init + ii, field_i) += coeff * field_value(input_field, comput_index_init + ii, field_j);
}
// clang-format on
// clang-format on
}
}
}
}
});
});
}
}
};

Expand Down
19 changes: 11 additions & 8 deletions include/samurai/schemes/fv/flux_based/explicit_flux_based_scheme__nonlin.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -47,15 +47,18 @@ namespace samurai
});

// Boundary interfaces
scheme().template for_each_boundary_interface<Run::Parallel>( // We need the 'template' keyword...
input_field,
[&](const auto& cell, auto& contrib)
{
for (std::size_t field_i = 0; field_i < output_field_size; ++field_i)
if (scheme().include_boundary_fluxes())
{
scheme().template for_each_boundary_interface<Run::Parallel>( // We need the 'template' keyword...
input_field,
[&](const auto& cell, auto& contrib)
{
field_value(output_field, cell, field_i) += this->scheme().flux_value_cmpnent(contrib, field_i);
}
});
for (std::size_t field_i = 0; field_i < output_field_size; ++field_i)
{
field_value(output_field, cell, field_i) += this->scheme().flux_value_cmpnent(contrib, field_i);
}
});
}
}
};
} // end namespace samurai
Loading

0 comments on commit c15b06d

Please sign in to comment.