Skip to content

Commit

Permalink
feat: remove field object from operators (#124)
Browse files Browse the repository at this point in the history 
Field object removed from the operators. Now only the field type is
used.
  • Loading branch information
@pmatalon
pmatalon authored Sep 12, 2023
1 parent 25431d0 commit 8c11e00
Show file tree
Hide file tree
Showing 21 changed files with 254 additions and 157 deletions.
2 changes: 1 addition & 1 deletion demos/FiniteVolume/burgers.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -208,7 +208,7 @@ int main_dim(int argc, char* argv[])
}

double cst = dim == 1 ? 0.5 : 1; // if dim == 1, we want f(u) = (1/2)*u^2
auto conv = cst * samurai::make_convection(u);
auto conv = cst * samurai::make_convection<decltype(u)>();

/**
* The following is another implementation of the convection operator (here in 1D):
Expand Down
4 changes: 2 additions & 2 deletions demos/FiniteVolume/heat.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -169,8 +169,8 @@ int main(int argc, char* argv[])
samurai::make_bc<samurai::Neumann>(u, 0.);
samurai::make_bc<samurai::Neumann>(unp1, 0.);

auto diff = diff_coeff * samurai::make_diffusion(u); // diffusion = -Laplacian
auto id = samurai::make_identity(u);
auto diff = diff_coeff * samurai::make_diffusion<decltype(u)>(); // diffusion = -Laplacian
auto id = samurai::make_identity<decltype(u)>();

//--------------------//
// Time iteration //
Expand Down
61 changes: 37 additions & 24 deletions demos/FiniteVolume/stokes_2d.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -233,6 +233,9 @@ int main(int argc, char* argv[])
auto velocity = samurai::make_field<dim, is_soa>("velocity", mesh);
auto pressure = samurai::make_field<1, is_soa>("pressure", mesh);

using VelocityField = decltype(velocity);
using PressureField = decltype(pressure);

// Boundary conditions
samurai::make_bc<samurai::Dirichlet>(velocity,
[](const auto&, const auto& coord)
Expand All @@ -258,10 +261,10 @@ int main(int argc, char* argv[])
// Stokes operator
// | Diff Grad |
// | -Div 0 |
auto diff = samurai::make_diffusion(velocity);
auto grad = samurai::make_gradient(pressure);
auto div = samurai::make_divergence(velocity);
auto zero_op = samurai::make_zero_operator(pressure);
auto diff = samurai::make_diffusion<VelocityField>();
auto grad = samurai::make_gradient<PressureField>();
auto div = samurai::make_divergence<VelocityField>();
auto zero_op = samurai::make_zero_operator<PressureField>();

auto stokes = samurai::make_block_operator<2, 2>(diff, grad,
-div, zero_op);
Expand All @@ -285,11 +288,11 @@ int main(int argc, char* argv[])
// Linear solver
std::cout << "Solving Stokes system..." << std::endl;
auto stokes_solver = samurai::petsc::make_solver<monolithic>(stokes);

stokes_solver.set_unknowns(velocity, pressure);
configure_solver(stokes_solver);

auto unknowns = stokes.tie_unknowns(velocity, pressure);
auto rhs = stokes.tie_rhs(f, zero);
stokes_solver.solve(unknowns, rhs);
stokes_solver.solve(f, zero);
std::cout << stokes_solver.iterations() << " iterations" << std::endl << std::endl;

// Error
Expand Down Expand Up @@ -397,6 +400,10 @@ int main(int argc, char* argv[])
auto velocity = samurai::make_field<dim, is_soa>("velocity", mesh);
auto velocity_np1 = samurai::make_field<dim, is_soa>("velocity_np1", mesh);
auto pressure_np1 = samurai::make_field<1, is_soa>("pressure_np1", mesh);

using VelocityField = decltype(velocity);
using PressureField = decltype(pressure_np1);

// Right-hand side
auto rhs = samurai::make_field<dim, is_soa>("rhs", mesh);
auto zero = samurai::make_field<1, is_soa>("zero", mesh);
Expand All @@ -418,11 +425,11 @@ int main(int argc, char* argv[])
// Stokes operator
// | Diff Grad |
// | -Div 0 |
auto diff = diff_coeff * samurai::make_diffusion(velocity_np1);
auto grad = samurai::make_gradient(pressure_np1);
auto div = samurai::make_divergence(velocity_np1);
auto zero_op = samurai::make_zero_operator(pressure_np1);
auto id = samurai::make_identity(velocity_np1);
auto diff = diff_coeff * samurai::make_diffusion<VelocityField>();
auto grad = samurai::make_gradient<PressureField>();
auto div = samurai::make_divergence<VelocityField>();
auto zero_op = samurai::make_zero_operator<PressureField>();
auto id = samurai::make_identity<VelocityField>();

// Stokes with backward Euler
// | I + dt*Diff dt*Grad |
Expand All @@ -433,6 +440,8 @@ int main(int argc, char* argv[])

// Linear solver
auto stokes_solver = samurai::petsc::make_solver<monolithic>(stokes);

stokes_solver.set_unknowns(velocity_np1, pressure_np1);
configure_solver(stokes_solver);

// Initial condition
Expand Down Expand Up @@ -512,6 +521,7 @@ int main(int argc, char* argv[])
stokes = samurai::make_block_operator<2, 2>(id + dt * diff, dt * grad, -div, zero_op);
}
stokes_solver = samurai::petsc::make_solver<monolithic>(stokes);
stokes_solver.set_unknowns(velocity_np1, pressure_np1);
configure_solver(stokes_solver);
}

Expand All @@ -527,9 +537,7 @@ int main(int argc, char* argv[])
rhs.fill(0);
rhs = velocity + dt * f;
zero.fill(0);
auto unknowns = stokes.tie_unknowns(velocity_np1, pressure_np1);
auto tied_rhs = stokes.tie_rhs(rhs, zero);
stokes_solver.solve(unknowns, tied_rhs);
stokes_solver.solve(rhs, zero);

// Prepare next step
std::swap(velocity.array(), velocity_np1.array());
Expand Down Expand Up @@ -597,6 +605,10 @@ int main(int argc, char* argv[])
auto velocity = samurai::make_field<dim, is_soa>("velocity", mesh);
auto velocity_np1 = samurai::make_field<dim, is_soa>("velocity_np1", mesh);
auto pressure_np1 = samurai::make_field<1, is_soa>("pressure_np1", mesh);

using VelocityField = decltype(velocity);
using PressureField = decltype(pressure_np1);

// Right-hand side
auto rhs = samurai::make_field<dim, is_soa>("rhs", mesh);
auto zero = samurai::make_field<1, is_soa>("zero", mesh);
Expand Down Expand Up @@ -627,12 +639,12 @@ int main(int argc, char* argv[])
// Stokes operator
// | Diff Grad |
// | -Div 0 |
auto diff = diff_coeff * samurai::make_diffusion(velocity_np1);
auto grad = samurai::make_gradient(pressure_np1);
auto conv = samurai::make_convection(velocity);
auto div = samurai::make_divergence(velocity_np1);
auto zero_op = samurai::make_zero_operator(pressure_np1);
auto id = samurai::make_identity(velocity_np1);
auto diff = diff_coeff * samurai::make_diffusion<VelocityField>();
auto grad = samurai::make_gradient<PressureField>();
auto conv = samurai::make_convection<VelocityField>();
auto div = samurai::make_divergence<VelocityField>();
auto zero_op = samurai::make_zero_operator<PressureField>();
auto id = samurai::make_identity<VelocityField>();

// Stokes with backward Euler
// | I + dt*Diff dt*Grad |
Expand All @@ -645,6 +657,8 @@ int main(int argc, char* argv[])

// Linear solver
auto stokes_solver = samurai::petsc::make_solver<monolithic>(stokes);

stokes_solver.set_unknowns(velocity_np1, pressure_np1);
configure_solver(stokes_solver);

// Time iteration
Expand Down Expand Up @@ -702,6 +716,7 @@ int main(int argc, char* argv[])
stokes = samurai::make_block_operator<2, 2>(id + dt * diff, dt * grad, -div, zero_op);
}
stokes_solver = samurai::petsc::make_solver<monolithic>(stokes);
stokes_solver.set_unknowns(velocity_np1, pressure_np1);
configure_solver(stokes_solver);
}

Expand All @@ -710,9 +725,7 @@ int main(int argc, char* argv[])
auto conv_v = conv(velocity);
rhs = velocity - dt * conv_v;
zero.fill(0);
auto unknowns = stokes.tie_unknowns(velocity_np1, pressure_np1);
auto block_rhs = stokes.tie_rhs(rhs, zero);
stokes_solver.solve(unknowns, block_rhs);
stokes_solver.solve(rhs, zero);

// Prepare next step
std::swap(velocity.array(), velocity_np1.array());
Expand Down
9 changes: 4 additions & 5 deletions demos/highorder/main.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -27,8 +27,7 @@ class HighOrderDiffusion : public samurai::CellBasedScheme<HighOrderDiffusion<Fi
using field_t = Field;
using directional_bdry_config_t = typename base_class::directional_bdry_config_t;

HighOrderDiffusion(Field& unknown)
: base_class(unknown)
HighOrderDiffusion()
{
}

Expand Down Expand Up @@ -116,9 +115,9 @@ class HighOrderDiffusion : public samurai::CellBasedScheme<HighOrderDiffusion<Fi
};

template <class Field>
auto make_high_order_diffusion(Field& f)
auto make_high_order_diffusion()
{
return HighOrderDiffusion<Field>(f);
return HighOrderDiffusion<Field>();
}

int main(int argc, char* argv[])
Expand Down Expand Up @@ -276,7 +275,7 @@ int main(int argc, char* argv[])
});
u.fill(0);

auto diff = make_high_order_diffusion(u);
HighOrderDiffusion<decltype(u)> diff;

auto solver = samurai::petsc::make_solver(diff);
KSP ksp = solver.Ksp();
Expand Down
5 changes: 3 additions & 2 deletions demos/multigrid/main.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -236,8 +236,9 @@ int main(int argc, char* argv[])
// Solve linear system //
//---------------------//

auto diff = samurai::make_diffusion<samurai::DirichletEnforcement::Equation>(solution);
auto diff = samurai::make_diffusion<samurai::DirichletEnforcement::Equation, decltype(solution)>();
auto solver = samurai::petsc::make_solver(diff);
solver.set_unknown(solution);

Timer setup_timer, solve_timer, total_timer;

Expand All @@ -250,7 +251,7 @@ int main(int argc, char* argv[])

std::cout << "Solving..." << std::endl;
solve_timer.Start();
solver.solve(solution, source);
solver.solve(source);
solve_timer.Stop();

total_timer.Stop();
Expand Down
19 changes: 18 additions & 1 deletion include/samurai/petsc/block_assembly.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -86,6 +86,13 @@ namespace samurai
});
}

template <class... Fields>
void set_unknowns(Fields&... unknowns)
{
auto unknown_tuple = block_operator().tie_unknowns(unknowns...);
set_unknown(unknown_tuple);
}

template <class... Fields>
void set_unknown(std::tuple<Fields...>& unknowns)
{
Expand Down Expand Up @@ -116,6 +123,17 @@ namespace samurai
});
}

bool undefined_unknown() const
{
bool undefined = false;
for_each_assembly_op(
[&](auto& op, auto, auto)
{
undefined = undefined || !op.unknown_ptr();
});
return undefined;
}

std::array<std::string, cols> field_names() const
{
std::array<std::string, cols> names;
Expand Down Expand Up @@ -333,7 +351,6 @@ namespace samurai
{
op.set_is_block(true);
});
reset();
}

void reset() override
Expand Down
20 changes: 17 additions & 3 deletions include/samurai/petsc/fv/FV_scheme_assembly.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -71,10 +71,8 @@ namespace samurai

explicit FVSchemeAssembly(const Scheme& scheme)
: m_scheme(&scheme)
, m_unknown(&scheme.unknown())
{
this->set_name(scheme.name());
// reset();
}

auto& scheme() const
Expand All @@ -84,6 +82,12 @@ namespace samurai

void reset() override
{
if (!m_unknown)
{
std::cerr << "Undefined unknown for operator " << scheme().name() << ". Assembly initialization failed!" << std::endl;
assert(false);
exit(EXIT_FAILURE);
}
m_n_cells = mesh().nb_cells();
// std::cout << "reset " << this->name() << ", rows = " << matrix_rows() << std::endl;
m_is_row_empty.resize(static_cast<std::size_t>(matrix_rows()));
Expand Down Expand Up @@ -163,10 +167,20 @@ namespace samurai

auto& unknown() const
{
assert(m_unknown);
assert(m_unknown && "undefined unknown");
return *m_unknown;
}

auto unknown_ptr() const
{
return m_unknown;
}

bool undefined_unknown() const
{
return !m_unknown;
}

auto& mesh() const
{
return unknown().mesh();
Expand Down
17 changes: 13 additions & 4 deletions include/samurai/petsc/fv/scheme_operators_assembly.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -32,7 +32,6 @@ namespace samurai

explicit Assembly(const scheme_t& sum_scheme)
: m_sum_scheme(&sum_scheme)
, m_unknown(&sum_scheme.unknown())
, m_flux_assembly(sum_scheme.flux_scheme())
, m_cell_assembly(sum_scheme.cell_scheme())
{
Expand All @@ -41,13 +40,23 @@ namespace samurai

auto& unknown() const
{
assert(m_unknown);
return *m_unknown;
return m_flux_assembly.unknown();
}

auto unknown_ptr() const
{
return m_flux_assembly.unknown_ptr();
}

bool undefined_unknown() const
{
return !m_flux_assembly.unknown_ptr();
}

void set_unknown(field_t& unknown)
{
m_unknown = &unknown;
m_flux_assembly.set_unknown(unknown);
m_cell_assembly.set_unknown(unknown);
}

auto& scheme() const
Expand Down
Loading

0 comments on commit 8c11e00

Please sign in to comment.