Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

feat: remove field object from operators #124

Merged
merged 2 commits into from
Sep 12, 2023
Merged

feat: remove field object from operators #124

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
362617f
Remove field object from operators
pmatalon Sep 11, 2023
d6e949c
Fix compil issue
pmatalon Sep 11, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion demos/FiniteVolume/burgers.cpp
View file
Open in desktop
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
Open in desktop
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
Open in desktop
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
Open in desktop
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
Open in desktop
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
Open in desktop
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
Open in desktop
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
Open in desktop
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