Skip to content

Commit

Permalink
fix: heat + ghost elimination (#112)
Browse files Browse the repository at this point in the history 
- Fix heat_1D demo
- Possibility to eliminate the ghosts from the linear sytem
  • Loading branch information
@pmatalon
pmatalon authored Jul 13, 2023
1 parent b3de2cb commit 5e43087
Show file tree
Hide file tree
Showing 12 changed files with 425 additions and 150 deletions.
74 changes: 37 additions & 37 deletions include/samurai/interface.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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<coord_index_t, dim>;
using cell_t = typename samurai::Cell<coord_index_t, dim>;

Stencil<2, dim> interface_stencil = in_out_stencil<dim>(direction);
auto interface_it = make_stencil_iterator(mesh, interface_stencil);
Expand Down Expand Up @@ -56,12 +56,12 @@ namespace samurai
int direction_index_int = find(comput_stencil, direction);
std::size_t direction_index = static_cast<std::size_t>(direction_index_int);

Stencil<comput_stencil_size, Mesh::dim> 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<std::size_t>(minus_direction_index_int);
auto minus_comput_stencil_it = make_stencil_iterator(mesh, minus_comput_stencil);
Stencil<comput_stencil_size, dim> 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<std::size_t>(minus_direction_index_int);
auto minus_comput_stencil_it = make_stencil_iterator(mesh, minus_comput_stencil);

for_each_level(
mesh,
Expand All @@ -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<Cell, 2> interface_cells;
std::array<cell_t, 2> 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();

Expand All @@ -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<Cell, 2> interface_cells;
std::array<cell_t, 2> interface_cells;
interface_cells[0] = minus_comput_stencil_it.cells()[minus_direction_index];
interface_cells[1] = coarse_it.cells()[0];

Expand All @@ -141,15 +141,15 @@ namespace samurai
Vector direction,
const Stencil<comput_stencil_size, Mesh::dim>& 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<coord_index_t, dim>;
using cell_t = typename samurai::Cell<coord_index_t, dim>;

Stencil<2, dim> interface_stencil = in_out_stencil<dim>(direction);
auto interface_it = make_stencil_iterator(mesh, interface_stencil);
Expand All @@ -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<mesh_interval_t>(
intersect,
Expand All @@ -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();
}
Expand All @@ -190,12 +190,12 @@ namespace samurai
int direction_index_int = find(comput_stencil, direction);
std::size_t direction_index = static_cast<std::size_t>(direction_index_int);

Stencil<comput_stencil_size, Mesh::dim> 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<std::size_t>(minus_direction_index_int);
auto minus_comput_stencil_it = make_stencil_iterator(mesh, minus_comput_stencil);
Stencil<comput_stencil_size, dim> 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<std::size_t>(minus_direction_index_int);
auto minus_comput_stencil_it = make_stencil_iterator(mesh, minus_comput_stencil);

for_each_level(
mesh,
Expand All @@ -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<mesh_interval_t>(
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<Cell, 2> interface_cells;
std::array<cell_t, 2> 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)
Expand All @@ -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<mesh_interval_t>(
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<Cell, 2> interface_cells;
std::array<cell_t, 2> 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)
Expand Down
2 changes: 1 addition & 1 deletion include/samurai/numeric/gauss_legendre.hpp
View file
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
#pragma once
#include "../cell.hpp"
#include "../utils.hpp"
#include "../static_algorithm.hpp"

namespace samurai
{
Expand Down
10 changes: 5 additions & 5 deletions include/samurai/petsc/block_assembly.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -42,11 +42,11 @@ namespace samurai
});
}

template <class OperatorType>
static Assembly<OperatorType> to_assembly(OperatorType& op)
{
return Assembly<OperatorType>(op);
}
// template <class OperatorType>
// static Assembly<OperatorType> to_assembly(OperatorType& op)
// {
// return Assembly<OperatorType>(op);
// }

auto& block_operator()
{
Expand Down
Loading

0 comments on commit 5e43087

Please sign in to comment.