diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index c93366d20..710aa38f9 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -4,8 +4,8 @@ on: [pull_request] # This section cancels previous runs for the same PR concurrency: - group: ${{ github.workflow }}-${{ github.event.pull_request.number || github.ref }} - cancel-in-progress: true + group: ${{ github.workflow }}-${{ github.event.pull_request.number || github.ref }} + cancel-in-progress: true jobs: # @@ -211,13 +211,14 @@ jobs: export LD_LIBRARY_PATH="$CONDA_PREFIX/lib:$LD_LIBRARY_PATH" cd tests pytest -v -s --h5diff + # # Check linux build with mamba environment # ######################################################### linux-mpi-mamba: needs: [pre-commit, cppcheck] - runs-on: ubuntu-latest + runs-on: self-hosted steps: - uses: actions/checkout@v4 @@ -233,29 +234,33 @@ jobs: - name: Mamba and samurai env installation uses: mamba-org/setup-micromamba@v2 + env: + MAMBA_ROOT_PREFIX: ${{ runner.temp }}/micromamba-root with: + download-micromamba: true + micromamba-binary-path: ${{ runner.temp }}/bin/micromamba environment-file: conda/mpi-environment.yml environment-name: samurai-env - cache-environment: true - - # - name: Conda update - # shell: bash -l {0} - # run: | - # conda update --all + generate-run-shell: true + cache-environment: false + cache-downloads: false + post-cleanup: all - name: Petsc installation - shell: bash -l {0} + shell: micromamba-shell {0} + env: + MAMBA_ROOT_PREFIX: ${{ runner.temp }}/micromamba-root run: | - conda install -y petsc pkg-config cxx-compiler + micromamba install -y petsc pkg-config cxx-compiler - name: Conda informations - shell: bash -l {0} + shell: micromamba-shell {0} run: | - conda info - conda list + micromamba info + micromamba list - name: Configure - shell: bash -l {0} + shell: micromamba-shell {0} run: | cmake \ . \ @@ -267,13 +272,12 @@ jobs: -DBUILD_TESTS=ON - name: Build - shell: bash -l {0} + shell: micromamba-shell {0} run: | - cmake --build build --target finite-volume-advection-2d --parallel 4 - cmake --build build --target finite-volume-burgers --parallel 4 + cmake --build build --target finite-volume-advection-2d finite-volume-burgers finite-volume-burgers-os-2d-mpi --parallel 4 - name: MPI test finite-volume-advection-2d - shell: bash -l {0} + shell: micromamba-shell {0} run: | cd build set -e # Stop on first failure @@ -289,7 +293,7 @@ jobs: python ../python/compare.py FV_advection_2d_size_1_ite_ FV_advection_2d_size_9_ite_ --start 1 --end 19 - name: MPI test finite-volume-burgers - shell: bash -l {0} + shell: micromamba-shell {0} run: | cd build set -e # Stop on first failure @@ -304,6 +308,18 @@ jobs: python ../python/compare.py burgers_2D_size_1_ite_ burgers_2D_size_4_ite_ --start 1 --end 19 python ../python/compare.py burgers_2D_size_1_ite_ burgers_2D_size_6_ite_ --start 1 --end 19 + - name: MPI test burgers OSMP 2D (weak scaling) + shell: micromamba-shell {0} + run: | + cd build + set -e # Stop on first failure + mpiexec -n 2 ./demos/FiniteVolume/finite-volume-burgers-os-2d-mpi --npx=2 --npy=1 --min-level=0 --max-level=8 --Tf=6 --nfiles 1 + mpiexec -n 2 ./demos/FiniteVolume/finite-volume-burgers-os-2d-mpi --npx=2 --npy=1 --min-level=0 --max-level=8 --Tf=6 --nfiles 1 + mpiexec -n 4 ./demos/FiniteVolume/finite-volume-burgers-os-2d-mpi --npx=2 --npy=2 --min-level=0 --max-level=8 --Tf=6 --nfiles 1 + mpiexec -n 8 ./demos/FiniteVolume/finite-volume-burgers-os-2d-mpi --npx=4 --npy=2 --min-level=0 --max-level=8 --Tf=6 --nfiles 1 + mpiexec -n 8 ./demos/FiniteVolume/finite-volume-burgers-os-2d-mpi --npx=2 --npy=4 --min-level=0 --max-level=8 --Tf=6 --nfiles 1 + mpiexec -n 9 ./demos/FiniteVolume/finite-volume-burgers-os-2d-mpi --npx=3 --npy=3 --min-level=0 --max-level=8 --Tf=6 --nfiles 1 + macos-mamba: needs: [pre-commit, cppcheck] runs-on: macos-14 diff --git a/demos/FiniteVolume/CMakeLists.txt b/demos/FiniteVolume/CMakeLists.txt index 8ff4b2bf9..b86353eac 100644 --- a/demos/FiniteVolume/CMakeLists.txt +++ b/demos/FiniteVolume/CMakeLists.txt @@ -25,6 +25,13 @@ burgers_mra.cpp:finite-volume-burgers-mra burgers_os.cpp:finite-volume-burgers-os ) +# List of MPI demos (only if WITH_MPI is enabled) +if(${WITH_MPI}) + list(APPEND STANDARD_DEMOS + burgers_os_2d_mpi.cpp:finite-volume-burgers-os-2d-mpi + ) +endif() + # Create executables with PETSc if(${WITH_PETSC}) message(STATUS "Building demos with PETSc support") diff --git a/demos/FiniteVolume/burgers_os_2d_mpi.cpp b/demos/FiniteVolume/burgers_os_2d_mpi.cpp new file mode 100644 index 000000000..a157141c2 --- /dev/null +++ b/demos/FiniteVolume/burgers_os_2d_mpi.cpp @@ -0,0 +1,278 @@ +// Copyright 2018-2025 the samurai's authors +// SPDX-License-Identifier: BSD-3-Clause + +#include +#include +#include +#include +#include +#include + +#include +namespace fs = std::filesystem; + +#include + +#include "convection_nonlinear_osmp.hpp" + +template +void save(const fs::path& path, const std::string& filename, const Field& u, const std::string& suffix = "") +{ + auto mesh = u.mesh(); + auto level_ = samurai::make_scalar_field("level", mesh); + + if (!fs::exists(path)) + { + fs::create_directory(path); + } + + samurai::for_each_cell(mesh, + [&](const auto& cell) + { + level_[cell] = cell.level; + }); + + samurai::save(path, fmt::format("{}{}", filename, suffix), mesh, u, level_); + // samurai::save(path, fmt::format("{}_full_{}", filename, suffix), {true, true}, mesh, u, level_); +} + +void check_diff(auto& mesh) +{ + using mesh_id_t = typename std::decay_t::mesh_id_t; + + mpi::communicator world; + + auto my_min_indices = mesh.subdomain().min_indices(); + + bool different = false; + for (const auto& neighbour : mesh.mpi_neighbourhood()) + { + auto neighbour_min_indices = neighbour.mesh.subdomain().min_indices(); + + xt::xtensor_fixed> translation{my_min_indices[0] - neighbour_min_indices[0], + my_min_indices[1] - neighbour_min_indices[1]}; + + samurai::for_each_level(mesh, + [&](auto level) + { + auto set = samurai::difference(mesh[mesh_id_t::cells][level], + samurai::translate(neighbour.mesh[mesh_id_t::cells][level], + translation >> (mesh.subdomain().level() - level))); + different = !set.empty(); + }); + if (different) + { + break; + } + } + + if (mpi::all_reduce(world, different, std::logical_or<>())) + { + for (const auto& neighbour : mesh.mpi_neighbourhood()) + { + auto neighbour_min_indices = neighbour.mesh.subdomain().min_indices(); + + xt::xtensor_fixed> translation{my_min_indices[0] - neighbour_min_indices[0], + my_min_indices[1] - neighbour_min_indices[1]}; + + samurai::for_each_level(mesh, + [&](auto level) + { + auto set = samurai::difference(mesh[mesh_id_t::cells][level], + samurai::translate(neighbour.mesh[mesh_id_t::cells][level], + translation >> (mesh.subdomain().level() - level))); + set( + [&](const auto& i, const auto& index) + { + std::cout << "Difference found !! " << level << " " << i << " " << index << " for domain " + << world.rank() << " with subdomain " << neighbour.rank << "\n"; + different = true; + }); + }); + } + samurai::save("diff_mesh", mesh); + throw std::runtime_error("Difference found between subdomains"); + } +} + +auto get_box(const xt::xtensor_fixed>& min_corner, const xt::xtensor_fixed>& max_corner, int npx) +{ + mpi::communicator world; + + const xt::xtensor_fixed> pcoords{static_cast(world.rank() % npx), static_cast(world.rank() / npx)}; + + auto length = max_corner - min_corner; + + return samurai::Box(min_corner + pcoords * length, max_corner + pcoords * length); +} + +int main(int argc, char* argv[]) +{ + static constexpr std::size_t dim = 2; + using Config = samurai::MRConfig; + + auto& app = samurai::initialize("Finite volume example for the linear convection equation", argc, argv); + + mpi::communicator world; + + std::cout << world.rank() << " / " << world.size() << std::endl; + + std::cout << "------------------------- Burgers 2D with OSMP scheme -------------------------" << std::endl; + + //--------------------// + // Program parameters // + //--------------------// + + // Simulation parameters + xt::xtensor_fixed> min_corner = {-1., -1.}; + xt::xtensor_fixed> max_corner = {1., 1.}; + + // Time integration + double Tf = 0.5; + double dt = 0; + double cfl = 0.95; + double t = 0.; + + // MPI parameters + int npx = 1; + int npy = 1; + + // Multiresolution parameters + std::size_t min_level = 8; + std::size_t max_level = 8; + + // Output parameters + fs::path path = fs::current_path(); + std::string filename = "burgers"; + std::size_t nfiles = 0; + + app.add_option("--min-corner", min_corner, "The min corner of the first box")->capture_default_str()->group("Simulation parameters"); + app.add_option("--max-corner", max_corner, "The max corner of the first box")->capture_default_str()->group("Simulation parameters"); + app.add_option("--Ti", t, "Initial time")->capture_default_str()->group("Simulation parameters"); + app.add_option("--Tf", Tf, "Final time")->capture_default_str()->group("Simulation parameters"); + app.add_option("--dt", dt, "Time step")->capture_default_str()->group("Simulation parameters"); + app.add_option("--cfl", cfl, "The CFL")->capture_default_str()->group("Simulation parameters"); + app.add_option("--npx", npx, "Number of processes in x direction")->capture_default_str()->group("MPI parameters"); + app.add_option("--npy", npy, "Number of processes in y direction")->capture_default_str()->group("MPI parameters"); + app.add_option("--min-level", min_level, "Minimum level of the multiresolution")->capture_default_str()->group("Multiresolution"); + app.add_option("--max-level", max_level, "Maximum level of the multiresolution")->capture_default_str()->group("Multiresolution"); + app.add_option("--path", path, "Output path")->capture_default_str()->group("Output"); + app.add_option("--filename", filename, "File name prefix")->capture_default_str()->group("Output"); + app.add_option("--nfiles", nfiles, "Number of output files")->capture_default_str()->group("Output"); + + app.allow_extras(); + SAMURAI_PARSE(argc, argv); + + if (world.size() != npx * npy) + { + throw std::runtime_error("Number of MPI processes must be equal to npx * npy"); + } + + std::cout << " max_level = " << max_level << " min_level = " << min_level << std::endl; + + double approx_box_tol = 0.05; + double scaling_factor = 1; + + //--------------------// + // Problem definition // + //--------------------// + + std::array periodic; + periodic.fill(true); + auto box = get_box(min_corner, max_corner, npx); + samurai::CellArray<2> cells; + for (std::size_t level = 0; level < cells.max_size; ++level) + { + cells[level].set_origin_point(min_corner); + cells[level].set_scaling_factor(scaling_factor); + } + cells[max_level] = {max_level, box, min_corner, approx_box_tol, scaling_factor}; + samurai::MRMesh mesh{cells, min_level, max_level}; + mesh.periodicity().fill(true); + mesh.find_neighbourhood(); + mesh = {cells, mesh}; + + auto u = samurai::make_scalar_field("u", mesh); + auto unp1 = samurai::make_scalar_field("unp1", mesh); + auto u1 = samurai::make_scalar_field("u1", mesh); + auto u2 = samurai::make_scalar_field("u2", mesh); + + auto middle = xt::eval(0.5 * (box.min_corner() + box.max_corner())); + middle += 0.25; + + // Initial solution + samurai::for_each_cell(mesh, + [&](auto& cell) + { + double coef = 0.5; + double stiffness = 50.; + u[cell] = coef + * std::exp(-stiffness + * ((cell.center(0) - middle(0)) * (cell.center(0) - middle(0)) + + (cell.center(1) - middle(1)) * (cell.center(1) - middle(1)))); + }); + + auto MRadaptation = samurai::make_MRAdapt(u); + auto mra_config = samurai::mra_config().epsilon(1e-3); + MRadaptation(mra_config); + + double dt_save = nfiles == 0 ? dt : Tf / static_cast(nfiles); + std::size_t nsave = 0, nt = 0; + if (nfiles != 1) + { + std::string suffix = (nfiles != 1) ? fmt::format("_level_{}_{}_np_{}_{}_ite_{}", min_level, max_level, npx, npy, nsave) : ""; + save(path, filename, u, suffix); + } + + // Convection operator + xt::xtensor_fixed> velocity = {0.5, 0.5}; + static constexpr std::size_t norder = 2; + auto conv = samurai::make_convection_nonlinear_osmp(dt); + + //--------------------// + // Time iteration // + //--------------------// + + double dx = mesh.cell_length(max_level); + dt = cfl * dx / velocity(0); + + while (t != Tf) + { + // Move to next timestep + t += dt; + if (t > Tf) + { + dt += Tf - t; + t = Tf; + } + std::cout << fmt::format("iteration {}: t = {:.12f}, dt = {}", nt++, t, dt) << std::flush; + + // Mesh adaptation + MRadaptation(mra_config); + check_diff(mesh); + u1.resize(); + u2.resize(); + unp1.resize(); + + u1 = u - 0.5 * dt * conv(0, u); + u2 = u1 - dt * conv(1, u1); + unp1 = u2 - 0.5 * dt * conv(0, u2); + + // u <-- unp1 + std::swap(u.array(), unp1.array()); + + // Save the result + if (nfiles == 0 || t >= static_cast(nsave + 1) * dt_save || t == Tf) + { + std::cout << " (saving results)" << std::flush; + std::string suffix = (nfiles != 1) ? fmt::format("_level_{}_{}_np_{}_{}_ite_{}", min_level, max_level, npx, npy, nsave) : ""; + save(path, filename, u, suffix); + nsave++; + } + + std::cout << std::endl; + } + + samurai::finalize(); + return 0; +} diff --git a/demos/FiniteVolume/convection_nonlinear_osmp.hpp b/demos/FiniteVolume/convection_nonlinear_osmp.hpp new file mode 100644 index 000000000..724837e33 --- /dev/null +++ b/demos/FiniteVolume/convection_nonlinear_osmp.hpp @@ -0,0 +1,327 @@ +#pragma once +#include + +namespace samurai +{ + + template + auto compute_osmp_flux_limiter(xtensor_t& d_alpha, xtensor_t& lambdalpha, xtensor_nu& nu, std::size_t j) + { + // using value_type = typename xtensor_t::value_type; + + static constexpr double zero = 1e-14; + + xt::xtensor_fixed> c_order; + + // value_type flux; + double phi_o; + double phi_lim; + + // Lax-Wendroff + phi_o = d_alpha[j]; + + // 3rd order + if (order >= 2) + { + for (std::size_t l = 0; l < stencil_size - 1; l++) + { + c_order(0, l) = (1. + nu[l]) / 3.; + } + phi_o += -c_order(0, j) * d_alpha[j] + c_order(0, j - 1) * d_alpha[j - 1]; + } + + if (order >= 3) + { + for (std::size_t l = 0; l < stencil_size - 1; l++) + { + c_order(1, l) = c_order(0, l) * (nu[l] - 2) / 4.; + c_order(2, l) = c_order(1, l) * (nu[l] + 2) / 5.; + } + // 4th order + phi_o += c_order(1, j + 1) * d_alpha[j + 1] - 2. * c_order(1, j) * d_alpha[j] + c_order(1, j - 1) * d_alpha[j - 1]; + + // 5th order + phi_o += -(c_order(2, j + 1) * d_alpha[j + 1] - 3. * c_order(2, j) * d_alpha[j] + 3. * c_order(2, j - 1) * d_alpha[j - 1] + - c_order(2, j - 2) * d_alpha[j - 2]); + } + + if (order >= 4) + { + for (std::size_t l = 0; l < stencil_size - 1; l++) + { + c_order(3, l) = c_order(2, l) * (nu[l] - 3) / 6.; + c_order(4, l) = c_order(3, l) * (nu[l] + 3) / 7.; + } + // 6th order + phi_o += c_order(3, j + 2) * d_alpha[j + 2] - 4. * c_order(3, j + 1) * d_alpha[j + 1] + 6. * c_order(3, j) * d_alpha[j] + - 4. * c_order(3, j - 1) * d_alpha[j - 1] + c_order(3, j - 2) * d_alpha[j - 2]; + // 7th order + phi_o += -(c_order(4, j + 2) * d_alpha[j + 2] - 5. * c_order(4, j + 1) * d_alpha[j + 1] + 10. * c_order(4, j) * d_alpha[j] + - 10. * c_order(4, j - 1) * d_alpha[j - 1] + 5. * c_order(4, j - 2) * d_alpha[j - 2] + - c_order(4, j - 3) * d_alpha[j - 3]); + } + + // MP constraint + xt::xtensor_fixed> d2; + xt::xtensor_fixed> djp12; + for (std::size_t i = 1; i < stencil_size - 1; i++) + { + d2[i] = lambdalpha[i] - lambdalpha[i - 1]; + } + + // minmod limiter between i and i+1 + for (std::size_t i = 1; i < stencil_size - 2; i++) + { + double dmm = 0.; + double d4min = 0.; + if (d2[i] * d2[i + 1] < 0.) + { + dmm = 0.; + } + else + { + dmm = std::min(std::abs(d2[i]), std::abs(d2[i + 1])); + if (d2[i] < 0.) + { + dmm = -dmm; + } + } + + if ((4. * d2[i] - d2[i + 1]) * (4. * d2[i + 1] - d2[i]) < 0.) + { + d4min = 0.; + } + else + { + d4min = std::min(std::abs(4. * d2[i] - d2[i + 1]), std::abs(4. * d2[i + 1] - d2[i])); + if ((4. * d2[i] - d2[i + 1]) < 0.) + { + d4min = -d4min; + } + } + + if (d4min * dmm < 0.) + { + djp12[i] = 0.; + } + else + { + djp12[i] = std::min(std::abs(d4min), std::abs(dmm)); + if (d4min < 0) + { + djp12[i] = -djp12[i]; + } + } + } + if (j <= 1) + { + djp12[j - 1] = djp12[j]; + } + + // TVD constraints + double phi_tvdmin = 2. * d_alpha[j - 1] / (nu[j - 1] + zero); + // double phi_tvdmin = 2. * d_alpha[j-1] * (1.-nu[j-1])/ (nu[j-1]*(1.-nu[j])+zero); + double phi_tvdmax = 2. * d_alpha[j] / (1 - nu[j] + zero); + + // MP limiter + double dfo = 0.5 * phi_o; + double dfabs = 0.5 * phi_tvdmax; + double dful = 0.5 * phi_tvdmin; + double dfmd = 0.5 * (dfabs - djp12[j]); + // double dflc = 0.5*(dful + (1.-nu[j-1])*djp12[j-1]/(nu[j]+zero)); + double dflc = 0.5 * (dful + (1. - nu[j - 1]) * djp12[j - 1] / (nu[j - 1] + zero)); + + double dfmin = std::max(std::min(0., std::min(dfabs, dfmd)), std::min(0., std::min(dful, dflc))); + double dfmax = std::min(std::max(0., std::max(dfabs, dfmd)), std::max(0., std::max(dful, dflc))); + + double val = (dfmin - dfo) * (dfmax - dfo); + if (val >= 0.) + { + phi_lim = std::max(0., std::min(phi_tvdmin, std::min(phi_o, phi_tvdmax))); + } + else + { + phi_lim = phi_o; + } + + return phi_lim; + } + + template + auto make_convection_nonlinear_osmp(double& dt) + { + // using field_value_t = typename Field::value_type; + + static constexpr std::size_t dim = Field::dim; + static constexpr std::size_t field_size = Field::n_comp; + static constexpr std::size_t stencil_size = 2 * order; + + using cfg = samurai::FluxConfig; + + samurai::FluxDefinition osmp; + + samurai::static_for<0, dim>::apply( // for each positive Cartesian direction 'd' + [&](auto integral_constant_d) + { + static constexpr int d = decltype(integral_constant_d)::value; + + auto f = [](auto u) -> samurai::FluxValue + { + if constexpr (field_size == 1) + { + return 0.5 * (u * u); + // return u; + } + else + { + return u(d) * u; + } + }; + + // osmp[d].stencil = line_stencil_from(1-static_cast(order)); + + osmp[d].cons_flux_function = [f, &dt](FluxValue& flux, const StencilData& data, const StencilValues& u) + { + static constexpr std::size_t j = order - 1; + + // for (std::size_t l = 0; l < stencil_size; l++) + // { + // std::cout << " Cell : l = " << l << data.cells[l] << std::endl; + // } + // std::cout << " j = " << data.cells[j] << " F = " << f(u[j]) << std::endl; + // std::cout << " U = " << u[j] << std::endl; + + // mean values + xt::xtensor_fixed, xt::xshape> lambda; + for (std::size_t l = 0; l < stencil_size - 1; l++) + { + lambda[l] = 0.5 * (u[l] + u[l + 1]); + // lambda[l] = 1.; + } + + auto dx = data.cell_length; + double sigma = dt / dx; + + // Variation at each interface l+1/2 of the stencil + xt::xtensor_fixed, xt::xshape> dalpha; + for (std::size_t l = 0; l < stencil_size - 1; l++) + { + dalpha[l] = u[l + 1] - u[l]; + } + + // Calculation of flux correction for each component + xt::xtensor_fixed> nu; + xt::xtensor_fixed> unmnudalpha; + xt::xtensor_fixed> lambdadalpha; + + if constexpr (field_size == 1) + { + if (lambda[j] >= 0) + { + for (std::size_t l = 0; l < stencil_size - 1; l++) + { + nu[l] = sigma * std::abs(lambda[l]); + } + // Factorized term of the Lax-Wendroff scheme + for (std::size_t l = 0; l < stencil_size - 1; l++) + { + unmnudalpha[l] = std::abs(lambda[l]) * (1. - nu[l]) * dalpha[l]; + lambdadalpha[l] = lambda[l] * dalpha[l]; + } + flux = f(u[j]); + } + else + { + for (std::size_t l = 0; l < stencil_size - 1; l++) + { + nu[l] = sigma * std::abs(lambda[stencil_size - 2 - l]); + } + // Factorized term of the Lax-Wendroff scheme + for (std::size_t l = 0; l < stencil_size - 1; l++) + { + unmnudalpha[l] = std::abs(lambda[stencil_size - 2 - l]) * (1. - nu[l]) * dalpha[stencil_size - 2 - l]; + lambdadalpha[l] = lambda[stencil_size - 2 - l] * dalpha[stencil_size - 2 - l]; + } + flux = f(u[j + 1]); + } + + // Limiter giving 1st-order Roe scheme + double phi_lim = 0.; + // Accuracy function for high-order approximations up to 7th-order + if (order > 1) + { + // Flux correction + phi_lim = compute_osmp_flux_limiter(unmnudalpha, + lambdadalpha, + nu, + j); + } + + flux += 0.5 * phi_lim; + } + else + { + // For each k-wave + samurai::FluxValue psi; + + for (std::size_t k = 0; k < field_size; k++) + { + if (lambda[j](k) >= 0) + { + for (std::size_t l = 0; l < stencil_size - 1; l++) + { + nu[l] = sigma * std::abs(lambda[l](k)); + } + // Factorized term of the Lax-Wendroff scheme + for (std::size_t l = 0; l < stencil_size - 1; l++) + { + unmnudalpha[l] = std::abs(lambda[l](k)) * (1. - nu[l]) * dalpha[l](k); + lambdadalpha[l] = lambda[l](k) * dalpha[l](k); + } + flux = f(u[j]); + } + else + { + for (std::size_t l = 0; l < stencil_size - 1; l++) + { + nu[l] = (dt / dx) * std::abs(lambda[stencil_size - 2 - l](k)); + } + // Factorized term of the Lax-Wendroff scheme + for (std::size_t l = 0; l < stencil_size - 1; l++) + { + unmnudalpha[l] = std::abs(lambda[stencil_size - 2 - l](k)) * (1. - nu[l]) + * dalpha[stencil_size - 2 - l](k); + lambdadalpha[l] = lambda[stencil_size - 2 - l](k) * dalpha[stencil_size - 2 - l](k); + } + flux = f(u[j + 1]); + } + + // Accuracy function for high-order approximations up to 7th-order + if (order > 1) + { + // Flux correction + double phi_lim = compute_osmp_flux_limiter( + unmnudalpha, + lambdadalpha, + nu, + j); + + psi(k) = phi_lim; + } + } + + // Projection back to the physical space + for (std::size_t k = 0; k < field_size; k++) + { + flux(k) += 0.5 * psi(k); + } + } + }; + }); + + auto scheme = make_flux_based_scheme(osmp); + scheme.set_name("convection"); + return scheme; + } + +} // end namespace samurai diff --git a/include/samurai/algorithm/graduation.hpp b/include/samurai/algorithm/graduation.hpp index e15574ddf..d9efc5c52 100644 --- a/include/samurai/algorithm/graduation.hpp +++ b/include/samurai/algorithm/graduation.hpp @@ -6,8 +6,17 @@ #include #include +#ifdef SAMURAI_WITH_MPI +#include + +#include +#include +namespace mpi = boost::mpi; +#endif + #include "../array_of_interval_and_point.hpp" #include "../cell_flag.hpp" +#include "../concepts.hpp" #include "../mesh.hpp" #include "../stencil.hpp" #include "../subset/node.hpp" @@ -228,36 +237,29 @@ namespace samurai return true; } - template + template void list_interval_to_refine_for_graduation( const size_t grad_width, const CellArray& ca, const LevelCellArray& domain, - [[maybe_unused]] const std::vector>& mpi_neighbourhood, + [[maybe_unused]] const auto& mpi_meshes, const std::array& is_periodic, const std::array& nb_cells_finest_level, std::array, CellArray::max_size>& out) { - const size_t max_level = ca.max_level() + 1; - const size_t min_level = ca.min_level(); - const size_t min_fine_level = (min_level + 2) - 1; // fine_level = max_level, max_level-1, ..., min_level+2. Thus fine_level != - // min_level+2-1 const int max_width = int(grad_width) + 1; for (size_t i = 0; i != max_size; ++i) { out[i].clear(); } -#ifdef SAMURAI_WITH_MPI - mpi::communicator world; - std::vector req; - for (const auto& mpi_neighbor : mpi_neighbourhood) - { - req.push_back(world.isend(mpi_neighbor.rank, mpi_neighbor.rank, ca)); - } -#endif // SAMURAI_WITH_MPI + const auto list_overlapping_intervals = [&](const auto& lhs_ca, const auto& rhs_ca) { + const size_t max_level = lhs_ca.max_level() + 1; + const size_t min_level = lhs_ca.min_level(); + const size_t min_fine_level = (min_level + 2) - 1; // fine_level = max_level, max_level-1, ..., min_level+2. Thus fine_level != + // min_level+2-1 auto apply_refine = [&](const auto& union_func, auto coarse_level, auto& isIntersectionEmpty) { for (int width = 1; isIntersectionEmpty and width != max_width; ++width) @@ -309,14 +311,42 @@ namespace samurai list_overlapping_intervals(ca, ca); #ifdef SAMURAI_WITH_MPI - CellArray neighbor_ca; - for (const auto& mpi_neighbor : mpi_neighbourhood) + for (const auto& mpi_mesh : mpi_meshes) { - world.recv(mpi_neighbor.rank, world.rank(), neighbor_ca); - list_overlapping_intervals(neighbor_ca, ca); + list_overlapping_intervals(mpi_mesh, ca); } +#endif // SAMURAI_WITH_MPI + } + + template + auto update_subdomains_mpi([[maybe_unused]] const mesh_t& mesh, const auto& mpi_neighbourhood) + { + std::vector mpi_meshes(mpi_neighbourhood.size()); +#ifdef SAMURAI_WITH_MPI + mpi::communicator world; + std::vector req; + + boost::mpi::packed_oarchive::buffer_type buffer; + boost::mpi::packed_oarchive oa(world, buffer); + oa << mesh; + + std::transform(mpi_neighbourhood.cbegin(), + mpi_neighbourhood.cend(), + std::back_inserter(req), + [&](const auto& neighbour) + { + return world.isend(neighbour.rank, neighbour.rank, buffer); + }); + + std::size_t index = 0; + for (auto& neighbour : mpi_neighbourhood) + { + world.recv(neighbour.rank, world.rank(), mpi_meshes[index++]); + } + mpi::wait_all(req.begin(), req.end()); #endif // SAMURAI_WITH_MPI + return mpi_meshes; } template @@ -324,6 +354,7 @@ namespace samurai const size_t half_stencil_width, const CellArray& ca, const LevelCellArray& domain, + [[maybe_unused]] const auto& mpi_meshes, const std::array& is_periodic, std::array, CellArray::max_size>& out) { @@ -332,8 +363,14 @@ namespace samurai return; } - const size_t max_level = ca.max_level(); - const size_t min_level = ca.min_level(); + size_t max_level = ca.max_level(); + size_t min_level = ca.min_level(); + + for (const auto& mpi_mesh : mpi_meshes) + { + max_level = std::max(max_level, mpi_mesh.max_level()); + min_level = std::min(min_level, mpi_mesh.min_level()); + } // We want to avoid a flux being computed with ghosts outside of the domain if the cell doesn't touch the boundary, // because we only want to apply the B.C. on the cells that touch the boundary. @@ -368,14 +405,24 @@ namespace samurai // the intersection. When the problem is fixed, remove the two following lines and uncomment the line // below. LevelCellArray translated_boundary(translate(boundaryCells, -i * translation)); - auto refine_subset = intersection(translated_boundary, ca[level - 1]).on(level - 1); - // auto refine_subset = intersection(translate(boundaryCells, -i*translation), ca[level-1]).on(level-1); - - refine_subset( - [&](const auto& x_interval, const auto& yz) - { - out[level - 1].push_back(x_interval, yz); - }); + auto refine_impl = [&](const auto& ca_lm1) + { + auto refine_subset = intersection(translated_boundary, ca_lm1).on(level - 1); + // auto refine_subset = intersection(translate(boundaryCells, -i*translation), ca[level-1]).on(level-1); + + refine_subset( + [&](const auto& x_interval, const auto& yz) + { + out[level - 1].push_back(x_interval, yz); + }); + }; + refine_impl(ca[level - 1]); +#ifdef SAMURAI_WITH_MPI + for (const auto& mpi_mesh : mpi_meshes) + { + refine_impl(mpi_mesh[level - 1]); + } +#endif } } } @@ -393,15 +440,24 @@ namespace samurai auto boundaryCells = difference(ca[level], translate(self(domain).on(level), -translation)); for (size_t i = 1; i != half_stencil_width; ++i) { - auto refine_subset = translate( - intersection(translate(boundaryCells, -i * translation), ca[level + 1]).on(level), - i * translation) - .on(level); - refine_subset( - [&](const auto& x_interval, const auto& yz) - { - out[level].push_back(x_interval, yz); - }); + auto refine_impl = [&](const auto& ca_lp1) + { + auto refine_subset = translate(intersection(translate(boundaryCells, -i * translation), ca_lp1).on(level), + i * translation) + .on(level); + refine_subset( + [&](const auto& x_interval, const auto& yz) + { + out[level].push_back(x_interval, yz); + }); + }; + refine_impl(ca[level + 1]); +#ifdef SAMURAI_WITH_MPI + for (const auto& mpi_mesh : mpi_meshes) + { + refine_impl(mpi_mesh[level + 1]); + } +#endif } } } @@ -419,10 +475,11 @@ namespace samurai const std::array& nb_cells_finest_level, std::array, CellArray::max_size>& out) { - list_interval_to_refine_for_graduation(grad_width, ca, domain, mpi_neighbourhood, is_periodic, nb_cells_finest_level, out); + auto mpi_meshes = update_subdomains_mpi(ca, mpi_neighbourhood); + list_interval_to_refine_for_graduation(grad_width, ca, domain, mpi_meshes, is_periodic, nb_cells_finest_level, out); if (!domain.empty()) { - list_interval_to_refine_for_contiguous_boundary_cells(half_stencil_width, ca, domain, is_periodic, out); + list_interval_to_refine_for_contiguous_boundary_cells(half_stencil_width, ca, domain, mpi_meshes, is_periodic, out); } } @@ -581,7 +638,9 @@ namespace samurai } // end for level // We then create new_ca as ca U ca_add new_ca.clear(); - for (std::size_t level = min_level; level != max_level + 1; ++level) + for (std::size_t level = std::min(ca.min_level(), ca_add_p.min_level()); + level != std::max(ca.max_level(), ca_add_p.max_level()) + 1; + ++level) { auto set = difference(union_(ca[level], ca_add_p[level]), ca_remove_p[level]); set( @@ -705,5 +764,4 @@ namespace samurai } return new_ca; } - } diff --git a/include/samurai/algorithm/update.hpp b/include/samurai/algorithm/update.hpp index 196e8fed5..d9b66ffa6 100644 --- a/include/samurai/algorithm/update.hpp +++ b/include/samurai/algorithm/update.hpp @@ -537,8 +537,8 @@ namespace samurai for (std::size_t level = max_level; level > min_level; --level) { - update_ghost_periodic(level, field, other_fields...); update_ghost_subdomains(level, field, other_fields...); + update_ghost_periodic(level, field, other_fields...); auto set_at_levelm1 = intersection(mesh[mesh_id_t::reference][level], mesh[mesh_id_t::proj_cells][level - 1]).on(level - 1); set_at_levelm1.apply_op(variadic_projection(field, other_fields...)); @@ -548,12 +548,12 @@ namespace samurai if (min_level > 0 && min_level != max_level) { - update_ghost_periodic(min_level - 1, field, other_fields...); update_ghost_subdomains(min_level - 1, field, other_fields...); + update_ghost_periodic(min_level - 1, field, other_fields...); update_outer_ghosts(min_level - 1, field, other_fields...); } - update_ghost_periodic(min_level, field, other_fields...); update_ghost_subdomains(min_level, field, other_fields...); + update_ghost_periodic(min_level, field, other_fields...); for (std::size_t level = min_level + 1; level <= max_level; ++level) { @@ -562,8 +562,8 @@ namespace samurai auto expr = intersection(pred_ghosts, mesh.subdomain(), mesh[mesh_id_t::all_cells][level - 1]).on(level); expr.apply_op(variadic_prediction(field, other_fields...)); - update_ghost_periodic(level, field, other_fields...); update_ghost_subdomains(level, field, other_fields...); + update_ghost_periodic(level, field, other_fields...); } // save(fs::current_path(), "update_ghosts", {true, true}, mesh, field); diff --git a/include/samurai/mesh.hpp b/include/samurai/mesh.hpp index 1fa2192f7..ce8c1ae24 100644 --- a/include/samurai/mesh.hpp +++ b/include/samurai/mesh.hpp @@ -110,6 +110,7 @@ namespace samurai bool is_periodic() const; bool is_periodic(std::size_t d) const; const std::array& periodicity() const; + std::array& periodicity(); // std::vector& neighbouring_ranks(); std::vector& mpi_neighbourhood(); const std::vector& mpi_neighbourhood() const; @@ -148,6 +149,8 @@ namespace samurai const lca_type& corner(const DirectionVector& direction) const; + void find_neighbourhood(); + protected: using derived_type = D; @@ -192,8 +195,6 @@ namespace samurai void update_sub_mesh(); void renumbering(); - void find_neighbourhood(); - void partition_mesh(std::size_t start_level, const Box& global_box); void load_balancing(); void load_transfer(const std::vector& load_fluxes); @@ -272,7 +273,7 @@ namespace samurai update_sub_mesh(); construct_corners(); renumbering(); - update_mesh_neighbour(); + // update_mesh_neighbour(); set_origin_point(origin_point()); set_scaling_factor(scaling_factor()); @@ -310,7 +311,7 @@ namespace samurai update_sub_mesh(); construct_corners(); renumbering(); - update_mesh_neighbour(); + // update_mesh_neighbour(); set_origin_point(domain_builder.origin_point()); set_scaling_factor(scaling_factor_); @@ -335,7 +336,7 @@ namespace samurai partition_mesh(start_level, b); // load_balancing(); #else - this->m_cells[mesh_id_t::cells][start_level] = {start_level, b, approx_box_tol, scaling_factor_}; + m_cells[mesh_id_t::cells][start_level] = {start_level, b, approx_box_tol, scaling_factor_}; #endif construct_subdomain(); @@ -343,7 +344,7 @@ namespace samurai update_sub_mesh(); construct_corners(); renumbering(); - update_mesh_neighbour(); + // update_mesh_neighbour(); set_origin_point(origin_point()); set_scaling_factor(scaling_factor()); @@ -357,7 +358,8 @@ namespace samurai m_periodic.fill(false); assert(min_level <= max_level); - this->m_cells[mesh_id_t::cells] = {cl}; + m_cells[mesh_id_t::cells] = {cl}; + update_meshid_neighbour(mesh_id_t::cells); construct_subdomain(); construct_domain(); @@ -365,7 +367,7 @@ namespace samurai update_sub_mesh(); construct_corners(); renumbering(); - update_mesh_neighbour(); + // update_mesh_neighbour(); set_origin_point(cl.origin_point()); set_scaling_factor(cl.scaling_factor()); @@ -379,7 +381,11 @@ namespace samurai m_periodic.fill(false); assert(min_level <= max_level); - this->m_cells[mesh_id_t::cells] = ca; + m_cells[mesh_id_t::cells] = ca; + update_meshid_neighbour(mesh_id_t::cells); + + set_origin_point(ca.origin_point()); + set_scaling_factor(ca.scaling_factor()); construct_subdomain(); construct_domain(); @@ -388,18 +394,18 @@ namespace samurai construct_corners(); renumbering(); -#ifdef SAMURAI_WITH_MPI - mpi::communicator world; - m_mpi_neighbourhood.clear(); - for (int i = 0; i < world.size(); ++i) - { - if (i != world.rank()) - { - m_mpi_neighbourhood.emplace_back(i); - } - } -#endif - update_mesh_neighbour(); + // #ifdef SAMURAI_WITH_MPI + // mpi::communicator world; + // m_mpi_neighbourhood.clear(); + // for (int i = 0; i < world.size(); ++i) + // { + // if (i != world.rank()) + // { + // m_mpi_neighbourhood.emplace_back(i); + // } + // } + // #endif + // update_mesh_neighbour(); set_origin_point(ca.origin_point()); set_scaling_factor(ca.scaling_factor()); @@ -415,13 +421,14 @@ namespace samurai { m_cells[mesh_id_t::cells] = ca; + update_meshid_neighbour(mesh_id_t::cells); construct_subdomain(); construct_union(); update_sub_mesh(); construct_corners(); renumbering(); - update_mesh_neighbour(); + // update_mesh_neighbour(); set_origin_point(ref_mesh.origin_point()); set_scaling_factor(ref_mesh.scaling_factor()); @@ -437,13 +444,14 @@ namespace samurai { m_cells[mesh_id_t::cells] = {cl, false}; + update_meshid_neighbour(mesh_id_t::cells); construct_subdomain(); construct_union(); update_sub_mesh(); construct_corners(); renumbering(); - update_mesh_neighbour(); + // update_mesh_neighbour(); set_origin_point(ref_mesh.origin_point()); set_scaling_factor(ref_mesh.scaling_factor()); @@ -751,6 +759,12 @@ namespace samurai return m_periodic; } + template + inline auto Mesh_base::periodicity() -> std::array& + { + return m_periodic; + } + template inline auto Mesh_base::mpi_neighbourhood() -> std::vector& { @@ -990,6 +1004,7 @@ namespace samurai find_neighbourhood(); } #endif + update_neighbour_subdomain(); } template diff --git a/include/samurai/mr/mesh.hpp b/include/samurai/mr/mesh.hpp index 70cbae98b..18e32c3a4 100644 --- a/include/samurai/mr/mesh.hpp +++ b/include/samurai/mr/mesh.hpp @@ -8,12 +8,15 @@ #include #include +// #include "../algorithm/graduation.hpp" +#include "../algorithm/graduation.hpp" #include "../box.hpp" #include "../mesh.hpp" #include "../samurai_config.hpp" #include "../stencil.hpp" #include "../subset/apply.hpp" #include "../subset/node.hpp" +#include "../timers.hpp" using namespace xt::placeholders; @@ -25,9 +28,9 @@ namespace samurai cells_and_ghosts = 1, proj_cells = 2, union_cells = 3, - all_cells = 4, + reference = 4, count = 5, - reference = all_cells + all_cells = reference }; template inline void MRMesh::update_sub_mesh_impl() { -#ifdef SAMURAI_WITH_MPI - mpi::communicator world; - // cppcheck-suppress redundantInitialization - auto max_level = mpi::all_reduce(world, this->cells()[mesh_id_t::cells].max_level(), mpi::maximum()); - // cppcheck-suppress redundantInitialization - auto min_level = mpi::all_reduce(world, this->cells()[mesh_id_t::cells].min_level(), mpi::minimum()); - cl_type cell_list; -#else - // cppcheck-suppress redundantInitialization - auto max_level = this->cells()[mesh_id_t::cells].max_level(); - // cppcheck-suppress redundantInitialization - auto min_level = this->cells()[mesh_id_t::cells].min_level(); + times::timers.start("mesh construction"); + cl_type cell_list; -#endif + // Construction of ghost cells // =========================== // // Example with max_stencil_width = 1 // - // level 2 |.|-|-|-|-|.| |-| - // cells - // |.| - // ghost - // cells - // level 1 |...|---|---|...|...|---|---|...| + // level 2 |.|-|-|-|-|-|-|.| |-| cells + // |.| ghost cells // - // level 0 |.......|-------|.......| |.......|-------|.......| + // level 1 |...|---|---|...| |...|---|---|...| // - for_each_interval( - this->cells()[mesh_id_t::cells], - [&](std::size_t level, const auto& interval, const auto& index_yz) - { - lcl_type& lcl = cell_list[level]; - static_nested_loop( - [&](auto stencil) - { - auto index = xt::eval(index_yz + stencil); - lcl[index].add_interval({interval.start - config::max_stencil_width, interval.end + config::max_stencil_width}); - }); - }); + // level 0 |.......|-------|.......| |.......|-------|.......| + // + for_each_interval(this->cells()[mesh_id_t::cells], + [&](std::size_t level, const auto& interval, const auto& index_yz) + { + lcl_type& lcl = cell_list[level]; + static_nested_loop( + [&](auto stencil) + { + auto index = xt::eval(index_yz + stencil); + lcl[index].add_interval({interval.start - config::ghost_width, interval.end + config::ghost_width}); + }); + }); this->cells()[mesh_id_t::cells_and_ghosts] = {cell_list, false}; - // Add cells for the MRA - if (this->max_level() != this->min_level()) + // Manage neighbourhood for MPI + // ============================== + // We do the same with the subdomains of the neighbouring processes + // to be sure that we have all the ghost cells and the cells at the interface + // with the other processes. + for (auto& neighbour : this->mpi_neighbourhood()) { - for (std::size_t level = max_level; level >= ((min_level == 0) ? 1 : min_level); --level) + for_each_level(neighbour.mesh[mesh_id_t::cells], + [&](std::size_t level) + { + lcl_type& lcl = cell_list[level]; + auto set = intersection(nestedExpand(neighbour.mesh[mesh_id_t::cells][level]), + nestedExpand(self(this->subdomain()).on(level))); + set( + [&](const auto& interval, const auto& index) + { + lcl[index].add_interval(interval); + }); + }); + } + + // Manage periodicity + // =================== + // This process is similar to the MPI one, but we need first to compute the directions to move the subdomains + // on the other side of the domain. + std::vector> directions; + if (this->is_periodic()) + { + std::array nb_cells_finest_level; + const auto& min_indices = this->domain().min_indices(); + const auto& max_indices = this->domain().max_indices(); + + for (size_t d = 0; d != max_indices.size(); ++d) { - auto expr = difference(intersection(this->cells()[mesh_id_t::cells_and_ghosts][level], self(this->domain()).on(level)), - this->get_union()[level]); + nb_cells_finest_level[d] = max_indices[d] - min_indices[d]; + } + directions = detail::get_periodic_directions(nb_cells_finest_level, 0, this->periodicity()); + } - expr( - [&](const auto& interval, const auto& index_yz) + auto add_periodic_cells = [&](const auto& subset, auto level) + { + lcl_type& lcl = cell_list[level]; + const int delta_l = int(this->domain().level() - level); + + for (const auto& d : directions) + { + auto set = intersection(nestedExpand(translate(subset, d >> delta_l)), + nestedExpand(self(this->subdomain()).on(level))); + set( + [&](const auto& interval, const auto& index) { - lcl_type& lcl = cell_list[level - 1]; - - static_nested_loop( - [&](auto stencil) - { - auto new_interval = interval >> 1; - lcl[(index_yz >> 1) + stencil].add_interval( - {new_interval.start - config::prediction_order, new_interval.end + config::prediction_order}); - }); + lcl[index].add_interval(interval); }); + } + }; + + // ghost cells are added by translating the mesh on the different periodic directions + // and intersecting with the subdomain + for_each_level(this->cells()[mesh_id_t::cells], + [&](std::size_t level) + { + add_periodic_cells(this->cells()[mesh_id_t::cells][level], level); + }); + + // We do the same with the subdomains of the neighbouring processes + for (auto& neighbour : this->mpi_neighbourhood()) + { + for_each_level(neighbour.mesh[mesh_id_t::cells], + [&](std::size_t level) + { + add_periodic_cells(neighbour.mesh[mesh_id_t::cells][level], level); + }); + } - auto expr_2 = intersection(this->cells()[mesh_id_t::cells][level], this->cells()[mesh_id_t::cells][level]); + // Add cells for the MRA to be able to compute the details at one and two levels below each cells. + // The prediction operator gives the number of ghost cells to add in each direction. + if (this->max_level() != this->min_level()) + { + auto add_prediction_ghosts = [&](auto&& subset_cells, auto&& subset_cells_and_ghosts, auto level) + { + lcl_type& lcl_m1 = cell_list[level - 1]; - expr_2( + subset_cells_and_ghosts( [&](const auto& interval, const auto& index_yz) { - if (level - 1 > 0) - { - lcl_type& lcl = cell_list[level - 2]; - - static_nested_loop( - [&](auto stencil) - { - auto new_interval = interval >> 2; - lcl[(index_yz >> 2) + stencil].add_interval( - {new_interval.start - config::prediction_order, new_interval.end + config::prediction_order}); - }); - } + lcl_m1[index_yz].add_interval(interval); }); - } - this->cells()[mesh_id_t::all_cells] = {cell_list, false}; - - this->update_meshid_neighbour(mesh_id_t::cells_and_ghosts); - this->update_meshid_neighbour(mesh_id_t::reference); - for (auto& neighbour : this->mpi_neighbourhood()) - { - for (std::size_t level = 0; level <= this->max_level(); ++level) + if (level - 1 > 0) { - auto expr = intersection(this->subdomain(), neighbour.mesh[mesh_id_t::reference][level]).on(level); - expr( + lcl_type& lcl_m2 = cell_list[level - 2]; + subset_cells( [&](const auto& interval, const auto& index_yz) { - lcl_type& lcl = cell_list[level]; - lcl[index_yz].add_interval(interval); - - if (level > neighbour.mesh[mesh_id_t::reference].min_level()) - { - lcl_type& lclm1 = cell_list[level - 1]; - - static_nested_loop( - [&](auto stencil) - { - auto new_interval = interval >> 1; - lclm1[(index_yz >> 1) + stencil].add_interval( - {new_interval.start - config::prediction_order, new_interval.end + config::prediction_order}); - }); - } + lcl_m2[index_yz].add_interval(interval); }); } - } - this->cells()[mesh_id_t::all_cells] = {cell_list, false}; - - using box_t = Box; - - const auto& domain = this->domain(); - const auto& domain_min_indices = domain.min_indices(); - const auto& domain_max_indices = domain.max_indices(); - - const auto& subdomain = this->subdomain(); - const auto& subdomain_min_indices = subdomain.min_indices(); - const auto& subdomain_max_indices = subdomain.max_indices(); - - // add ghosts for periodicity - xt::xtensor_fixed> shift; - xt::xtensor_fixed> min_corner; - xt::xtensor_fixed> max_corner; - - shift.fill(0); - -#ifdef SAMURAI_WITH_MPI - std::size_t reference_max_level = mpi::all_reduce(world, - this->cells()[mesh_id_t::reference].max_level(), - mpi::maximum()); - std::size_t reference_min_level = mpi::all_reduce(world, - this->cells()[mesh_id_t::reference].min_level(), - mpi::minimum()); + }; - std::vector neighbourhood_extended_subdomain(this->mpi_neighbourhood().size()); - for (size_t neighbor_id = 0; neighbor_id != neighbourhood_extended_subdomain.size(); ++neighbor_id) - { - const auto& neighbor_subdomain = this->mpi_neighbourhood()[neighbor_id].mesh.subdomain(); - if (not neighbor_subdomain.empty()) + for_each_level( + this->cells()[mesh_id_t::cells], + [&](std::size_t level) { - const auto& neighbor_min_indices = neighbor_subdomain.min_indices(); - const auto& neighbor_max_indices = neighbor_subdomain.max_indices(); - for (std::size_t level = reference_min_level; level <= reference_max_level; ++level) + // own part + add_prediction_ghosts(nestedExpand(self(this->cells()[mesh_id_t::cells][level]).on(level - 2)), + intersection(nestedExpand( + self(this->cells()[mesh_id_t::cells_and_ghosts][level]).on(level - 1)), + nestedExpand(self(this->subdomain()).on(level - 1))), + level); + + // periodic part + const int delta_l = int(this->domain().level() - level); + for (const auto& d : directions) { - const std::size_t delta_l = subdomain.level() - level; - box_t box; - for (std::size_t d = 0; d < dim; ++d) - { - box.min_corner()[d] = (neighbor_min_indices[d] >> delta_l) - config::ghost_width; - box.max_corner()[d] = (neighbor_max_indices[d] >> delta_l) + config::ghost_width; - } - neighbourhood_extended_subdomain[neighbor_id][level] = {level, box}; + add_prediction_ghosts( + intersection(nestedExpand( + translate(this->cells()[mesh_id_t::cells][level], d >> delta_l).on(level - 2)), + self(this->subdomain()).on(level - 2)), + intersection(nestedExpand( + translate(this->cells()[mesh_id_t::cells_and_ghosts][level], d >> delta_l).on(level - 1)), + self(this->subdomain()).on(level - 1)), + level); } - } - } -#endif // SAMURAI_WITH_MPI - const auto& mesh_ref = this->cells()[mesh_id_t::reference]; - for (std::size_t level = 0; level <= this->max_level(); ++level) - { - const std::size_t delta_l = subdomain.level() - level; - lcl_type& lcl = cell_list[level]; - - for (std::size_t d = 0; d < dim; ++d) - { - min_corner[d] = (subdomain_min_indices[d] >> delta_l) - config::ghost_width; - max_corner[d] = (subdomain_max_indices[d] >> delta_l) + config::ghost_width; - } - lca_type lca_extended_subdomain(level, box_t(min_corner, max_corner)); - for (std::size_t d = 0; d < dim; ++d) - { - if (this->is_periodic(d)) - { - shift[d] = (domain_max_indices[d] - domain_min_indices[d]) >> delta_l; - - min_corner[d] = (domain_min_indices[d] >> delta_l) - config::ghost_width; - max_corner[d] = (domain_min_indices[d] >> delta_l) + config::ghost_width; - - lca_type lca_min(level, box_t(min_corner, max_corner)); - - min_corner[d] = (domain_max_indices[d] >> delta_l) - config::ghost_width; - max_corner[d] = (domain_max_indices[d] >> delta_l) + config::ghost_width; - - lca_type lca_max(level, box_t(min_corner, max_corner)); - - auto set1 = intersection(translate(intersection(mesh_ref[level], lca_min), shift), - intersection(lca_extended_subdomain, lca_max)); - set1( - [&](const auto& i, const auto& index_yz) - { - lcl[index_yz].add_interval(i); - }); - auto set2 = intersection(translate(intersection(mesh_ref[level], lca_max), -shift), - intersection(lca_extended_subdomain, lca_min)); - set2( - [&](const auto& i, const auto& index_yz) - { - lcl[index_yz].add_interval(i); - }); -#ifdef SAMURAI_WITH_MPI - //~ for (const auto& mpi_neighbor : this->mpi_neighbourhood()) - for (size_t neighbor_id = 0; neighbor_id != this->mpi_neighbourhood().size(); ++neighbor_id) - { - const auto& mpi_neighbor = this->mpi_neighbourhood()[neighbor_id]; - const auto& neighbor_mesh_ref = mpi_neighbor.mesh[mesh_id_t::reference]; - - auto set1_mpi = intersection(translate(intersection(neighbor_mesh_ref[level], lca_min), shift), - intersection(lca_extended_subdomain, lca_max)); - set1_mpi( - [&](const auto& i, const auto& index_yz) - { - lcl[index_yz].add_interval(i); - }); - auto set2_mpi = intersection(translate(intersection(neighbor_mesh_ref[level], lca_max), -shift), - intersection(lca_extended_subdomain, lca_min)); - set2_mpi( - [&](const auto& i, const auto& index_yz) - { - lcl[index_yz].add_interval(i); - }); - } -#endif // SAMURAI_WITH_MPI - this->cells()[mesh_id_t::all_cells][level] = {lcl}; - - /* reset variables for next iterations. */ - shift[d] = 0; - min_corner[d] = (subdomain_min_indices[d] >> delta_l) - config::ghost_width; - max_corner[d] = (subdomain_max_indices[d] >> delta_l) + config::ghost_width; - } - } - } - for (std::size_t level = 0; level < max_level; ++level) - { - lcl_type& lcl = cell_list[level + 1]; - lcl_type lcl_proj{level}; - auto expr = intersection(this->cells()[mesh_id_t::all_cells][level], this->get_union()[level]); - - expr( - [&](const auto& interval, const auto& index_yz) - { - static_nested_loop( - [&](auto s) - { - lcl[(index_yz << 1) + s].add_interval(interval << 1); - }); - lcl_proj[index_yz].add_interval(interval); - }); - this->cells()[mesh_id_t::all_cells][level + 1] = lcl; - this->cells()[mesh_id_t::proj_cells][level] = lcl_proj; - } - this->update_neighbour_subdomain(); - this->update_meshid_neighbour(mesh_id_t::all_cells); + }); + // mpi part + this->update_meshid_neighbour(mesh_id_t::cells_and_ghosts); for (auto& neighbour : this->mpi_neighbourhood()) { - for (std::size_t level = 0; level <= this->max_level(); ++level) - { - auto expr = intersection(nestedExpand(self(this->subdomain()).on(level), config::ghost_width), - neighbour.mesh[mesh_id_t::reference][level]); - expr( - [&](const auto& interval, const auto& index_yz) - { - lcl_type& lcl = cell_list[level]; - lcl[index_yz].add_interval(interval); - }); - for (std::size_t d = 0; d < dim; ++d) + for_each_level( + neighbour.mesh[mesh_id_t::cells], + [&](std::size_t level) { - if (this->is_periodic(d)) + add_prediction_ghosts( + intersection(nestedExpand(self(neighbour.mesh[mesh_id_t::cells][level]).on(level - 2)), + self(this->subdomain()).on(level - 2)), + intersection(nestedExpand( + self(neighbour.mesh[mesh_id_t::cells_and_ghosts][level]).on(level - 1)), + self(this->subdomain()).on(level - 1)), + level); + + const int delta_l = int(this->domain().level() - level); + + for (const auto& d : directions) { - auto domain_shift = get_periodic_shift(this->domain(), level, d); - auto expr_left = intersection(nestedExpand(self(this->subdomain()).on(level), config::ghost_width), - translate(neighbour.mesh[mesh_id_t::reference][level], -domain_shift)); - expr_left( - [&](const auto& interval, const auto& index_yz) - { - lcl_type& lcl = cell_list[level]; - lcl[index_yz].add_interval(interval); - }); - - auto expr_right = intersection(nestedExpand(self(this->subdomain()).on(level), config::ghost_width), - translate(neighbour.mesh[mesh_id_t::reference][level], domain_shift)); - expr_right( - [&](const auto& interval, const auto& index_yz) - { - lcl_type& lcl = cell_list[level]; - lcl[index_yz].add_interval(interval); - }); + add_prediction_ghosts( + intersection(nestedExpand( + translate(neighbour.mesh[mesh_id_t::cells][level], d >> delta_l).on(level - 2)), + self(this->subdomain()).on(level - 2)), + intersection(nestedExpand( + translate(neighbour.mesh[mesh_id_t::cells_and_ghosts][level], d >> delta_l).on(level - 1)), + self(this->subdomain()).on(level - 1)), + level); } - } - } + }); } - this->cells()[mesh_id_t::all_cells] = {cell_list, false}; - this->update_meshid_neighbour(mesh_id_t::all_cells); - } - - else - { - this->cells()[mesh_id_t::all_cells] = {cell_list, false}; - // TODO : I think we do not want to update subdomain in this case, it remains the same iteration after iteration. - this->update_neighbour_subdomain(); - this->update_meshid_neighbour(mesh_id_t::cells_and_ghosts); - this->update_meshid_neighbour(mesh_id_t::reference); + this->cells()[mesh_id_t::reference] = {cell_list, false}; + + // Some of the ghosts added by the prediction or the stencil of the scheme can be computed by projection + // as described in the following figure (1D example): + // + // other levels |||||||| + // + // level 2 |-|-|.|.|.| |-| cells + // |.| ghost cells added by the prediction + // level 1 |---|...|...| + // + // We can observe that the first ghost cell on the right of the cell at level 1 can be computed by projection + // from level 2. But we only have one ghost cell at level 2, so we need to add the other one to have enough information + // + // level 2 |-|-|.|.|.|.| |-| cells + // |.| ghost cells added by the prediction + // level 1 |---|...|...| + // + for_each_level(this->cells()[mesh_id_t::reference], + [&](auto level) + { + lcl_type& lcl = cell_list[level + 1]; + lcl_type lcl_proj{level}; + auto expr = intersection(this->cells()[mesh_id_t::reference][level], this->get_union()[level]); + + expr( + [&](const auto& interval, const auto& index_yz) + { + static_nested_loop( + [&](auto s) + { + lcl[(index_yz << 1) + s].add_interval(interval << 1); + }); + lcl_proj[index_yz].add_interval(interval); + }); + this->cells()[mesh_id_t::reference][level + 1] = lcl; + this->cells()[mesh_id_t::proj_cells][level] = lcl_proj; + }); } + this->cells()[mesh_id_t::reference] = {cell_list, false}; + this->update_meshid_neighbour(mesh_id_t::reference); + times::timers.stop("mesh construction"); } template @@ -517,8 +438,8 @@ struct fmt::formatter : formatter case samurai::MRMeshId::union_cells: name = "union cells"; break; - case samurai::MRMeshId::all_cells: - name = "all cells"; + case samurai::MRMeshId::reference: + name = "reference"; break; case samurai::MRMeshId::count: name = "count"; diff --git a/include/samurai/schemes/fv/FV_scheme.hpp b/include/samurai/schemes/fv/FV_scheme.hpp index cdfa0d2e3..ff71d85b2 100644 --- a/include/samurai/schemes/fv/FV_scheme.hpp +++ b/include/samurai/schemes/fv/FV_scheme.hpp @@ -181,8 +181,8 @@ namespace samurai { if (input_field.mesh().is_periodic()) { - std::cerr << "Error: apply_directional_bc() not implemented for non-box domains with periodic directions." << std::endl; - assert(false); + // std::cerr << "Error: apply_directional_bc() not implemented for non-box domains with periodic directions." << + // std::endl; assert(false); return; } apply_field_bc(input_field, d); diff --git a/include/samurai/static_algorithm.hpp b/include/samurai/static_algorithm.hpp index c2a3adf0a..79fa3c1cc 100644 --- a/include/samurai/static_algorithm.hpp +++ b/include/samurai/static_algorithm.hpp @@ -111,6 +111,24 @@ namespace samurai return detail::NestedExpand::run(idx, lca, width); } + template + auto nestedExpand_impl(const LCA_OR_SET& lca, std::integral_constant) + { + return nestedExpand(lca, 1); + } + + template + auto nestedExpand_impl(const LCA_OR_SET& lca, std::integral_constant) + { + return nestedExpand(nestedExpand(lca, std::integral_constant{}), 1); + } + + template + auto nestedExpand(const LCA_OR_SET& lca) + { + return nestedExpand_impl(lca, std::integral_constant{}); + } + template inline void nestedLoop(int i0, int i1, Function&& func) {