52 std::shared_ptr<StructuredBlockForest> blocks, BlockDataID field_id,
53 BlockDataID tmp_field_id,
unsigned int n_ghost_layers,
54 unsigned int shear_direction,
unsigned int shear_plane_normal,
56 std::function<
double()>
get_shift = []() {
return 0.0; })
57 : m_blocks(std::move(blocks)), m_field_id(field_id),
58 m_tmp_field_id(tmp_field_id), m_n_ghost_layers(uint_c(n_ghost_layers)),
59 m_shear_direction(uint_c(shear_direction)),
60 m_shear_plane_normal(uint_c(shear_plane_normal)),
63 if (m_n_ghost_layers != 1u) {
64 throw std::domain_error(
"The Lees-Edwards sweep is implemented "
65 "for a ghost layer of thickness 1");
67 if (m_shear_plane_normal == 0
u) {
68 m_slab_min = stencil::W;
69 m_slab_max = stencil::E;
70 }
else if (m_shear_plane_normal == 1u) {
71 m_slab_min = stencil::S;
72 m_slab_max = stencil::N;
73 }
else if (m_shear_plane_normal == 2u) {
74 m_slab_min = stencil::B;
75 m_slab_max = stencil::T;
80 return numeric_cast<FloatType>(m_get_pos_offset());
83 FloatType
get_shift()
const {
return numeric_cast<FloatType>(m_get_shift()); }
86 kernel(
block, m_slab_min);
87 kernel(
block, m_slab_max);
91 void kernel(IBlock *
block, stencil::Direction slab_dir) {
93 assert(m_blocks->getNumberOfCells(*
block, m_shear_plane_normal) >= 2u);
94 auto const dir = m_shear_direction;
95 auto const dim = cell_idx_c(m_blocks->getNumberOfCells(*
block, dir));
96 auto const length = numeric_cast<FloatType>(dim);
101 auto field =
block->template getData<FieldType>(m_field_id);
102 auto tmp_field =
block->template getData<FieldType>(m_tmp_field_id);
105 field->getGhostRegion(slab_dir, ci, cell_idx_t{1},
true);
111 auto const prefactor =
112 ((slab_dir == m_slab_max) ? FloatType{-1} : FloatType{1});
114 for (
auto const &&cell : ci) {
117 source1[dir] = cell_idx_c(std::floor(
118 static_cast<FloatType
>(source1[dir]) + offset)) %
120 source1[dir] = cell_idx_c(
static_cast<FloatType
>(source1[dir]) + length);
121 source1[dir] = cell_idx_c(source1[dir] % dim);
124 cell_idx_c(std::ceil(
static_cast<FloatType
>(source2[dir]) + offset)) %
126 source2[dir] = cell_idx_c(
static_cast<FloatType
>(source2[dir]) + length);
127 source2[dir] = cell_idx_c(source2[dir] % dim);
129 for (uint_t
f = 0;
f < FieldType::F_SIZE; ++
f) {
130 tmp_field->get(cell,
f) = field->get(source1,
f) * (1 -
weight) +
131 field->get(source2,
f) *
weight;
133 tmp_field->get(cell, m_shear_direction) -= prefactor * shift;
137 for (
auto const &&cell : ci) {
138 for (uint_t
f = 0;
f < FieldType::F_SIZE; ++
f) {
139 field->get(cell,
f) = tmp_field->get(cell,
f);
145 std::shared_ptr<StructuredBlockForest> m_blocks;
146 BlockDataID m_field_id;
147 BlockDataID m_tmp_field_id;
148 uint_t m_n_ghost_layers;
149 uint_t m_shear_direction;
150 uint_t m_shear_plane_normal;
151 std::function<double()> m_get_pos_offset;
152 std::function<double()> m_get_shift;
153 stencil::Direction m_slab_min;
154 stencil::Direction m_slab_max;
InterpolateAndShiftAtBoundary(std::shared_ptr< StructuredBlockForest > blocks, BlockDataID field_id, BlockDataID tmp_field_id, unsigned int n_ghost_layers, unsigned int shear_direction, unsigned int shear_plane_normal, std::function< double()> get_pos_offset, std::function< double()> get_shift=[]() { return 0.0;})