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

Jump to bottom

feat: new flux-based schemes #117

Merged
merged 11 commits into from
Jul 26, 2023
Merged

feat: new flux-based schemes #117

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
11 commits
Select commit Hold shift + click to select a range
40a7405
New fluxes
pmatalon Jul 17, 2023
8ce444b
for_each_boundary_interface() now also handles the opposite direction
pmatalon Jul 17, 2023
70b04c0
for_each_interface() refactoring
pmatalon Jul 17, 2023
9433785
Fix make_bc when the bc is a function
pmatalon Jul 18, 2023
74f5f5d
for_each_boundary_interface() refactoring
pmatalon Jul 18, 2023
b7d1d0c
Comments
pmatalon Jul 18, 2023
9c9be5a
Grad and Div built from 'average_quantity' flux
pmatalon Jul 20, 2023
19f0171
Add some useful methods in block_assembly to apply the matrix to a ve…
pmatalon Jul 20, 2023
373a71b
Rename operator variables in Stokes
pmatalon Jul 20, 2023
cb21673
Some renaming
pmatalon Jul 20, 2023
2a00760
Add "this->" to fix bug in gcc
pmatalon Jul 20, 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
148 changes: 89 additions & 59 deletions demos/FiniteVolume/stokes_2d.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -70,32 +70,6 @@ void configure_LU_solver(Solver& solver)
{
#if defined(PETSC_HAVE_SUPERLU)
PCFactorSetMatSolverType(pc, MATSOLVERSUPERLU); // (equiv. '-pc_factor_mat_solver_type superlu')
// Configure the following options:
// -mat_superlu_rowperm LargeDiag -mat_superlu_colperm NATURAL -mat_superlu_diagpivotthresh 0.5 -mat_superlu_iterrefine SINGLE
PetscBool set = PETSC_FALSE;
PetscOptionsHasName(NULL, NULL, "-mat_superlu_rowperm", &set);
if (!set)
{
PetscOptionsSetValue(NULL, "-mat_superlu_rowperm", "LargeDiag");
}

PetscOptionsHasName(NULL, NULL, "-mat_superlu_colperm", &set);
if (!set)
{
PetscOptionsSetValue(NULL, "-mat_superlu_colperm", "NATURAL");
}

PetscOptionsHasName(NULL, NULL, "-mat_superlu_diagpivotthresh", &set);
if (!set)
{
PetscOptionsSetValue(NULL, "-mat_superlu_diagpivotthresh", "0.5");
}

PetscOptionsHasName(NULL, NULL, "-mat_superlu_iterrefine", &set);
if (!set)
{
PetscOptionsSetValue(NULL, "-mat_superlu_iterrefine", "SINGLE");
}
#endif
}
else if (use_mumps)
Expand All @@ -104,8 +78,9 @@ void configure_LU_solver(Solver& solver)
PCFactorSetMatSolverType(pc, MATSOLVERMUMPS); // (equiv. '-pc_factor_mat_solver_type mumps')
#endif
}
KSPSetFromOptions(ksp); // KSP and PC overwritten by user value if needed
// If SuperLU is not installed, you can use: -ksp_type gmres -pc_type ilu
// KSP and PC overwritten by user value if needed
KSPSetFromOptions(ksp);
// If neither SuperLU nor MUMPS is installed, you can try: -ksp_type gmres -pc_type ilu
}

//
Expand Down Expand Up @@ -179,7 +154,7 @@ int main(int argc, char* argv[])
using Mesh = samurai::MRMesh<Config>;
using mesh_id_t = typename Mesh::mesh_id_t;
static constexpr bool is_soa = false;
static constexpr bool monolithic = false;
static constexpr bool monolithic = true;

//----------------//
// Parameters //
Expand Down Expand Up @@ -278,13 +253,13 @@ int main(int argc, char* argv[])
// Stokes operator
// | Diff Grad |
// | -Div 0 |
auto diff_v = samurai::make_diffusion_FV(velocity);
auto grad_p = samurai::make_gradient_FV(pressure);
auto div_v = samurai::make_divergence_FV(velocity);
auto zero_p = samurai::make_zero_operator_FV(pressure);
auto diff = samurai::make_diffusion_FV(velocity);
auto grad = samurai::make_gradient_FV(pressure);
auto div = samurai::make_divergence_FV(velocity);
auto zero_op = samurai::make_zero_operator_FV(pressure);

auto stokes = samurai::make_block_operator<2, 2>(diff_v, grad_p,
-div_v, zero_p);
auto stokes = samurai::make_block_operator<2, 2>(diff, grad,
-div, zero_op);

// clang-format on

Expand Down Expand Up @@ -433,17 +408,17 @@ int main(int argc, char* argv[])
// Stokes operator
// | Diff Grad |
// | -Div 0 |
auto diff_v = diff_coeff * samurai::make_diffusion_FV(velocity_np1);
auto grad_p = samurai::make_gradient_FV(pressure_np1);
auto div_v = samurai::make_divergence_FV(velocity_np1);
auto zero_p = samurai::make_zero_operator_FV(pressure_np1);
auto id_v = samurai::make_identity_FV(velocity_np1);
auto diff = diff_coeff * samurai::make_diffusion_FV(velocity_np1);
auto grad = samurai::make_gradient_FV(pressure_np1);
auto div = samurai::make_divergence_FV(velocity_np1);
auto zero_op = samurai::make_zero_operator_FV(pressure_np1);
auto id = samurai::make_identity_FV(velocity_np1);

// Stokes with backward Euler
// | I + dt*Diff dt*Grad |
// | -Div 0 |
auto stokes = samurai::make_block_operator<2, 2>(id_v + dt * diff_v, dt * grad_p,
-div_v, zero_p);
auto stokes = samurai::make_block_operator<2, 2>(id + dt * diff, dt * grad,
-div, zero_op);
// clang-format on

// Linear solver
Expand Down Expand Up @@ -521,7 +496,7 @@ int main(int argc, char* argv[])
{
if (dt_has_changed)
{
stokes = samurai::make_block_operator<2, 2>(id_v + dt * diff_v, dt * grad_p, -div_v, zero_p);
stokes = samurai::make_block_operator<2, 2>(id + dt * diff, dt * grad, -div, zero_op);
}
stokes_solver = samurai::petsc::make_solver<monolithic>(stokes);
configure_solver(stokes_solver);
Expand Down Expand Up @@ -558,25 +533,28 @@ int main(int argc, char* argv[])
std::cout.precision(2);
std::cout << ", L2-error: " << std::scientific << error;

// Divergence
auto div_velocity = div(velocity);

// Save the result
std::string suffix = (nfiles != 1) ? fmt::format("_ite_{}", nsave++) : "";
samurai::save(path, fmt::format("{}{}", filename, suffix), velocity.mesh(), velocity);
samurai::save(path, fmt::format("{}{}", filename, suffix), velocity.mesh(), velocity, div_velocity);

if (min_level != max_level)
{
// Reconstruction on the finest level
samurai::update_ghost_mr(velocity);
auto velocity_recons = samurai::reconstruction(velocity);
// Error
double error_recons = L2_error(velocity_recons,
[&](auto& coord)
{
return exact_velocity(t_n, coord);
});
double error_recons = samurai::L2_error(velocity_recons,
[&](auto& coord)
{
return exact_velocity(t_n, coord);
});
std::cout.precision(2);
std::cout << ", L2-error (recons): " << std::scientific << error_recons;
// Save
samurai::save(path, fmt::format("{}_recons{}", filename, suffix), velocity_recons.mesh(), velocity_recons);
samurai::save(path, fmt::format("{}_recons{}", filename, suffix), velocity_recons.mesh(), velocity_recons, div_velocity);
}
std::cout << std::endl;
}
Expand Down Expand Up @@ -629,17 +607,17 @@ int main(int argc, char* argv[])
// Stokes operator
// | Diff Grad |
// | -Div 0 |
auto diff_v = diff_coeff * samurai::make_diffusion_FV(velocity_np1);
auto grad_p = samurai::make_gradient_FV(pressure_np1);
auto div_v = samurai::make_divergence_FV(velocity_np1);
auto zero_p = samurai::make_zero_operator_FV(pressure_np1);
auto id_v = samurai::make_identity_FV(velocity_np1);
auto diff = diff_coeff * samurai::make_diffusion_FV(velocity_np1);
auto grad = samurai::make_gradient_FV(pressure_np1);
auto div = samurai::make_divergence_FV(velocity_np1);
auto zero_op = samurai::make_zero_operator_FV(pressure_np1);
auto id = samurai::make_identity_FV(velocity_np1);

// Stokes with backward Euler
// | I + dt*Diff dt*Grad |
// | -Div 0 |
auto stokes = samurai::make_block_operator<2, 2>(id_v + dt * diff_v, dt * grad_p,
-div_v, zero_p);
auto stokes = samurai::make_block_operator<2, 2>(id + dt * diff, dt * grad,
-div, zero_op);
// clang-format on

auto MRadaptation = samurai::make_MRAdapt(velocity);
Expand Down Expand Up @@ -696,7 +674,7 @@ int main(int argc, char* argv[])
{
if (dt_has_changed)
{
stokes = samurai::make_block_operator<2, 2>(id_v + dt * diff_v, dt * grad_p, -div_v, zero_p);
stokes = samurai::make_block_operator<2, 2>(id + dt * diff, dt * grad, -div, zero_op);
}
stokes_solver = samurai::petsc::make_solver<monolithic>(stokes);
configure_solver(stokes_solver);
Expand All @@ -718,11 +696,63 @@ int main(int argc, char* argv[])
{
samurai::update_ghost_mr(velocity);
auto velocity_recons = samurai::reconstruction(velocity);
auto div_velocity = div(velocity);

std::string suffix = (nfiles != 1) ? fmt::format("_ite_{}", nsave++) : "";
samurai::save(path, fmt::format("{}{}", filename, suffix), velocity.mesh(), velocity);
samurai::save(path, fmt::format("{}{}", filename, suffix), velocity.mesh(), velocity, div_velocity);
samurai::save(path, fmt::format("{}_recons{}", filename, suffix), velocity_recons.mesh(), velocity_recons);
}

// srand(time(NULL));
// auto x_velocity = samurai::make_field<dim, is_soa>("x_velocity", mesh);
// auto x_pressure = samurai::make_field<1, is_soa>("x_pressure", mesh);
// samurai::for_each_cell(mesh[mesh_id_t::reference],
// [&](auto cell)
// {
// x_velocity[cell] = rand() % 10 + 1;
// x_pressure[cell] = rand() % 10 + 1;
// });
// auto x = stokes.tie_unknowns(x_velocity, x_pressure);

// auto monolithicAssembly = samurai::petsc::make_assembly<true>(stokes);
// Mat monolithicA;
// monolithicAssembly.create_matrix(monolithicA);
// monolithicAssembly.assemble_matrix(monolithicA);
// Vec mono_x = monolithicAssembly.create_applicable_vector(x); // copy
// auto result_velocity_mono = samurai::make_field<dim, is_soa>("result_velocity", mesh);
// auto result_pressure_mono = samurai::make_field<1, is_soa>("result_pressure", mesh);
// auto result_mono = stokes.tie_rhs(result_velocity_mono, result_pressure_mono);
// Vec mono_result = monolithicAssembly.create_rhs_vector(result_mono); // copy
// MatMult(monolithicA, mono_x, mono_result);
// monolithicAssembly.update_result_fields(mono_result, result_mono); // copy

// auto nestedAssembly = samurai::petsc::make_assembly<false>(stokes);
// Mat nestedA;
// nestedAssembly.create_matrix(nestedA);
// nestedAssembly.assemble_matrix(nestedA);
// Vec nest_x = nestedAssembly.create_applicable_vector(x);
// auto result_velocity_nest = samurai::make_field<dim, is_soa>("result_velocity", mesh);
// auto result_pressure_nest = samurai::make_field<1, is_soa>("result_pressure", mesh);
// auto result_nest = stokes.tie_rhs(result_velocity_nest, result_pressure_nest);
// Vec nest_result = nestedAssembly.create_rhs_vector(result_nest);
// MatMult(nestedA, nest_x, nest_result);

// std::cout << std::setprecision(15);
// std::cout << std::fixed;
// samurai::for_each_cell(
// mesh[mesh_id_t::reference],
// [&](auto cell)
// {
// std::cout << round(result_velocity_mono[cell][0] * 1.e8) / 1.e8 << std::endl;
// if (round(result_velocity_mono[cell][0] * 1.e5) / 1.e5 != round(result_velocity_nest[cell][0] * 1.e5) / 1.e5)
// {
// std::cout << result_velocity_mono[cell][0] << " =? " << result_velocity_nest[cell][0] << std::endl;
// }
// if (result_pressure_mono[cell] != result_pressure_nest[cell])
// {
// std::cout << result_pressure_mono[cell] << " =? " << result_pressure_nest[cell] << std::endl;
// }
// });
}
}
else
Expand Down
4 changes: 1 addition & 3 deletions include/samurai/bc.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -710,9 +710,7 @@ namespace samurai
}

template <template <class> class bc_type, class Field>
auto make_bc(Field& field,
const std::function<detail::return_type_t<typename Field::value_type, Field::size>(
const xt::xtensor_fixed<typename Field::mesh_t::interval_t::value_t, xt::xshape<Field::dim>>&)>& func)
auto make_bc(Field& field, typename FunctionBc<Field::dim, typename Field::value_type, Field::size>::function_t func)
{
using value_t = typename Field::value_type;
constexpr std::size_t dim = Field::dim;
Expand Down
Loading