Skip to content

Commit

Permalink
Back and forth MLS
Browse files Browse the repository at this point in the history
  • Loading branch information
@mrlag31
mrlag31 committed Aug 24, 2023
1 parent 220407c commit dcdedd3
Showing 1 changed file with 184 additions and 108 deletions.
292 changes: 184 additions & 108 deletions examples/moving_least_squares/moving_least_squares.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,19 +9,11 @@
* SPDX-License-Identifier: BSD-3-Clause *
****************************************************************************/

// Example taken from DataTransferKit
// (https://github.com/ORNL-CEES/DataTransferKit)
// with MLS resolution from
// (http://dx.doi.org/10.1016/j.jcp.2015.11.055)
// and
// (A conservative mesh-free approach for fluid-structure interface problems)

#include <ArborX.hpp>

#include <Kokkos_Core.hpp>

#include <iostream>
#include <limits>
#include <tuple>
#include <vector>

Expand All @@ -33,153 +25,237 @@
using ExecutionSpace = Kokkos::DefaultExecutionSpace;
using MemorySpace = ExecutionSpace::memory_space;

using HostExecutionSpace = Kokkos::DefaultHostExecutionSpace;
using HostMemorySpace = HostExecutionSpace::memory_space;

// Function to approximate
KOKKOS_INLINE_FUNCTION float manufactured_solution(ArborX::Point const &p)
struct Step
{
return Kokkos::cos(5 * p[2]) * p[0] + p[1] + 1;
}
KOKKOS_INLINE_FUNCTION static float eval(ArborX::Point const &p)
{
return !Kokkos::signbit(p[0]) * 1.f;
}

int main(int argc, char *argv[])
{
MPI_Init(&argc, &argv);
Kokkos::ScopeGuard guard(argc, argv);
template <class... Properties>
static Kokkos::View<float *, Properties...>
map(ExecutionSpace const &space,
Kokkos::View<ArborX::Point *, Properties...> const &ps)
{
Kokkos::View<float *, Properties...> evals("Example::evals", ps.extent(0));
Kokkos::parallel_for(
"Example::evaluation",
Kokkos::RangePolicy<ExecutionSpace>(space, 0, ps.extent(0)),
KOKKOS_LAMBDA(int const i) { evals(i) = eval(ps(i)); });
return evals;
}
};

ExecutionSpace space{};
MPI_Comm mpi_comm = MPI_COMM_WORLD;
Kokkos::Array<Kokkos::View<ArborX::Point *, MemorySpace>, 2>
createPointClouds(HostExecutionSpace const &hspace, ExecutionSpace const &space,
MPI_Comm comm, std::size_t points_num)
{
int mpi_size, mpi_rank;
MPI_Comm_size(mpi_comm, &mpi_size);
MPI_Comm_rank(mpi_comm, &mpi_rank);
MPI_Comm_size(comm, &mpi_size);
MPI_Comm_rank(comm, &mpi_rank);

Kokkos::Array<Kokkos::View<ArborX::Point *, HostMemorySpace>, 2>
point_clouds_host{Kokkos::View<ArborX::Point *, HostMemorySpace>(
Kokkos::view_alloc(Kokkos::WithoutInitializing,
"Example::points_cloud_0"),
points_num),
Kokkos::View<ArborX::Point *, HostMemorySpace>(
Kokkos::view_alloc(Kokkos::WithoutInitializing,
"Example::points_cloud_1"),
points_num)};

static constexpr std::size_t total_source_points = 1024 * 512;
std::size_t local_source_points_num = total_source_points / mpi_size;
static constexpr std::size_t total_target_points = 1024;
std::size_t local_target_points_num = total_target_points / mpi_size;
static constexpr double cube_side = 5;

Kokkos::View<ArborX::Point *, MemorySpace> source_points(
Kokkos::view_alloc(Kokkos::WithoutInitializing, "Example::source_points"),
local_source_points_num);
auto source_points_host = Kokkos::create_mirror_view(source_points);
Kokkos::View<ArborX::Point *, MemorySpace> target_points(
Kokkos::view_alloc(Kokkos::WithoutInitializing, "Example::target_points"),
local_source_points_num);
auto target_points_host = Kokkos::create_mirror_view(target_points);

// source and target points are within a 5x5x5 cube
if (mpi_rank == 0)
{
Kokkos::View<ArborX::Point *, Kokkos::HostSpace> all_source_points(
Kokkos::view_alloc(Kokkos::WithoutInitializing,
"Example::all_source_points"),
total_source_points);
filledBoxCloud(cube_side / 2, all_source_points);
MPI_Scatter(
all_source_points.data(), local_source_points_num * 3 * sizeof(float),
MPI_BYTE, source_points_host.data(),
local_source_points_num * 3 * sizeof(float), MPI_BYTE, 0, mpi_comm);

Kokkos::View<ArborX::Point *, Kokkos::HostSpace> all_target_points(
Kokkos::view_alloc(Kokkos::WithoutInitializing,
"Example::all_target_points"),
total_target_points);
filledBoxCloud(cube_side / 2, all_target_points);
MPI_Scatter(
all_target_points.data(), local_target_points_num * 3 * sizeof(float),
MPI_BYTE, target_points_host.data(),
local_target_points_num * 3 * sizeof(float), MPI_BYTE, 0, mpi_comm);
Kokkos::Array<Kokkos::View<ArborX::Point *, HostMemorySpace>, 2>
all_point_clouds{Kokkos::View<ArborX::Point *, HostMemorySpace>(
Kokkos::view_alloc(Kokkos::WithoutInitializing,
"Example::all_points_cloud_0"),
points_num * mpi_size),
Kokkos::View<ArborX::Point *, HostMemorySpace>(
Kokkos::view_alloc(Kokkos::WithoutInitializing,
"Example::all_points_cloud_1"),
points_num * mpi_size)};

filledBoxCloud(.5, all_point_clouds[0]);
filledBoxCloud(.5, all_point_clouds[1]);

MPI_Scatter(all_point_clouds[0].data(), points_num * 3 * sizeof(float),
MPI_BYTE, point_clouds_host[0].data(),
points_num * 3 * sizeof(float), MPI_BYTE, 0, comm);
MPI_Scatter(all_point_clouds[1].data(), points_num * 3 * sizeof(float),
MPI_BYTE, point_clouds_host[1].data(),
points_num * 3 * sizeof(float), MPI_BYTE, 0, comm);
}
else
{
MPI_Scatter(nullptr, local_source_points_num * 3 * sizeof(float), MPI_BYTE,
source_points_host.data(),
local_source_points_num * 3 * sizeof(float), MPI_BYTE, 0,
mpi_comm);

MPI_Scatter(nullptr, local_target_points_num * 3 * sizeof(float), MPI_BYTE,
target_points_host.data(),
local_target_points_num * 3 * sizeof(float), MPI_BYTE, 0,
mpi_comm);
MPI_Scatter(nullptr, 0, MPI_BYTE, point_clouds_host[0].data(),
points_num * 3 * sizeof(float), MPI_BYTE, 0, comm);
MPI_Scatter(nullptr, 0, MPI_BYTE, point_clouds_host[1].data(),
points_num * 3 * sizeof(float), MPI_BYTE, 0, comm);
}

Kokkos::deep_copy(space, source_points, source_points_host);
Kokkos::deep_copy(space, target_points, target_points_host);

// Create the transform from a point cloud to another
MovingLeastSquares<MemorySpace, float> mls(mpi_comm, space, source_points,
target_points, Details::degree<2>,
Details::wendland<0>);

// Compute source values
Kokkos::View<float *, MemorySpace> source_values("Example::source_values",
local_source_points_num);
Kokkos::parallel_for(
"Example::source_evaluation",
Kokkos::RangePolicy<ExecutionSpace>(space, 0, local_source_points_num),
"Example::flatten_points",
Kokkos::RangePolicy<HostExecutionSpace>(hspace, 0, points_num),
KOKKOS_LAMBDA(int const i) {
source_values(i) = manufactured_solution(source_points(i));
point_clouds_host[0](i)[2] = 0;
point_clouds_host[1](i)[2] = 0;
});

// Compute target values from source ones
auto target_values = mls.apply(space, source_values);
auto target_values_host = Kokkos::create_mirror_view(target_values);
Kokkos::deep_copy(space, target_values_host, target_values);
Kokkos::Array<Kokkos::View<ArborX::Point *, MemorySpace>, 2> point_clouds{
Kokkos::View<ArborX::Point *, MemorySpace>(
Kokkos::view_alloc(Kokkos::WithoutInitializing,
"Example::points_cloud_0"),
points_num),
Kokkos::View<ArborX::Point *, MemorySpace>(
Kokkos::view_alloc(Kokkos::WithoutInitializing,
"Example::points_cloud_1"),
points_num)};
Kokkos::deep_copy(space, point_clouds[0], point_clouds_host[0]);
Kokkos::deep_copy(space, point_clouds[1], point_clouds_host[1]);

// Compute target values via evaluation
Kokkos::View<float *, MemorySpace> target_values_exact(
"Example::target_values_exact", local_target_points_num);
Kokkos::parallel_for(
"Example::target_evaluation",
Kokkos::RangePolicy<ExecutionSpace>(space, 0, local_target_points_num),
KOKKOS_LAMBDA(int const i) {
target_values_exact(i) = manufactured_solution(target_points(i));
});
return point_clouds;
}

template <typename Deg, typename RBF>
Kokkos::Array<MovingLeastSquares<MemorySpace, float>, 2> createMLSObjects(
MPI_Comm comm, ExecutionSpace const &space,
Kokkos::View<ArborX::Point *, MemorySpace> const &point_clouds_0,
Kokkos::View<ArborX::Point *, MemorySpace> const &point_clouds_1,
Deg const &deg, RBF const &rbf)
{
return {MovingLeastSquares<MemorySpace, float>(comm, space, point_clouds_0,
point_clouds_1, deg, rbf),
MovingLeastSquares<MemorySpace, float>(comm, space, point_clouds_1,
point_clouds_0, deg, rbf)};
}

void doError(MPI_Comm comm, ExecutionSpace const &space,
Kokkos::View<ArborX::Point *, MemorySpace> const &points,
Kokkos::View<float *, MemorySpace> const &approx,
Kokkos::View<float *, MemorySpace> const &values)
{
int mpi_size, mpi_rank;
MPI_Comm_size(comm, &mpi_size);
MPI_Comm_rank(comm, &mpi_rank);

// Compute local error
static constexpr float epsilon = std::numeric_limits<float>::epsilon();
using ErrType = typename Kokkos::MaxLoc<float, std::size_t>::value_type;
ErrType error{0, 0};
float error_sum = 0;
Kokkos::parallel_reduce(
"Example::error_computation",
Kokkos::RangePolicy<ExecutionSpace>(space, 0, local_target_points_num),
KOKKOS_LAMBDA(int const i, ErrType &loc_error) {
float abs_error =
Kokkos::abs(target_values(i) - target_values_exact(i));
float abs_value = Kokkos::abs(target_values_exact(i)) +
epsilon;

if (loc_error.val < abs_error / abs_value)
Kokkos::RangePolicy<ExecutionSpace>(space, 0, approx.extent(0)),
KOKKOS_LAMBDA(int const i, ErrType &loc_error, float &loc_error_sum) {
float abs_error = Kokkos::abs(approx(i) - values(i));

loc_error_sum += abs_error;
if (loc_error.val < abs_error)
{
loc_error.val = abs_error / abs_value;
loc_error.val = abs_error;
loc_error.loc = i;
}
},
Kokkos::MaxLoc<float, std::size_t>(error));
Kokkos::MaxLoc<float, std::size_t>(error), Kokkos::Sum<float>(error_sum));

auto approx_host = Kokkos::create_mirror_view(approx);
auto points_host = Kokkos::create_mirror_view(points);
Kokkos::deep_copy(space, approx_host, approx);
Kokkos::deep_copy(space, points_host, points);

std::tuple<float, ArborX::Point, float> error_obj{
error.val, target_points_host(error.loc), target_values_host(error.loc)};
error.val, points_host(error.loc), approx_host(error.loc)};

// Compute global error
if (mpi_rank == 0)
{
float error_sum_global;
std::vector<decltype(error_obj)> all_error_obj(mpi_size);
MPI_Gather(&error_obj, sizeof(decltype(error_obj)), MPI_BYTE,
all_error_obj.data(), sizeof(decltype(error_obj)), MPI_BYTE, 0,
mpi_comm);
comm);
MPI_Reduce(&error_sum, &error_sum_global, 1, MPI_FLOAT, MPI_SUM, 0, comm);

for (int i = 0; i < mpi_size; i++)
if (std::get<0>(error_obj) < std::get<0>(all_error_obj[i]))
error_obj = all_error_obj[i];

float error = std::get<0>(error_obj), approx = std::get<2>(error_obj);
auto point = std::get<1>(error_obj);
std::cout << "Maximum error: " << error << " at point " << point[0] << ", "
<< point[1] << ", " << point[2]
<< "\nTrue value: " << manufactured_solution(point)
<< "\nComputed: " << approx << std::endl;
std::cout << "Mean error: "
<< error_sum_global / (points.extent(0) * mpi_size)
<< "\nMaximum error: " << error << " at point " << point[0]
<< ", " << point[1] << "\n True value: " << Step::eval(point)
<< "\n Computed: " << approx << std::endl;
}
else
{
MPI_Gather(&error_obj, sizeof(decltype(error_obj)), MPI_BYTE, nullptr,
sizeof(decltype(error_obj)), MPI_BYTE, 0, mpi_comm);
sizeof(decltype(error_obj)), MPI_BYTE, 0, comm);
MPI_Reduce(&error_sum, nullptr, 1, MPI_FLOAT, MPI_SUM, 0, comm);
}
}

Kokkos::View<float *, MemorySpace>
doOne(MPI_Comm comm, ExecutionSpace const &space,
Kokkos::View<ArborX::Point *, MemorySpace> const &tgt,
Kokkos::View<float *, MemorySpace> const &values,
Kokkos::View<float *, MemorySpace> const &true_values,
MovingLeastSquares<MemorySpace, float> &mls)
{
auto tgt_values = mls.apply(space, values);
doError(comm, space, tgt, tgt_values, true_values);
return tgt_values;
}

int main(int argc, char *argv[])
{
static constexpr std::size_t total_points = 1024 * 128;
static constexpr std::size_t num_back_forth = 50;
static constexpr auto deg = Details::degree<4>;
static constexpr auto rbf = Details::wu<2>;

MPI_Init(&argc, &argv);
Kokkos::ScopeGuard guard(argc, argv);

ExecutionSpace space{};
HostExecutionSpace host_space{};
MPI_Comm mpi_comm = MPI_COMM_WORLD;
int mpi_size, mpi_rank;
MPI_Comm_size(mpi_comm, &mpi_size);
MPI_Comm_rank(mpi_comm, &mpi_rank);

auto point_clouds =
createPointClouds(host_space, space, mpi_comm, total_points / mpi_size);

// Create the transform from a point cloud to another
auto mlss = createMLSObjects(mpi_comm, space, point_clouds[0],
point_clouds[1], deg, rbf);

Kokkos::Array<Kokkos::View<float *, MemorySpace>, 2> true_values{
Step::map(space, point_clouds[0]), Step::map(space, point_clouds[1])};

Kokkos::View<float *, MemorySpace> source_values = true_values[0];
for (int i = 0; i < num_back_forth * 2; i++)
{
if (mpi_rank == 0)
std::cout << "=== TURN " << i + 1 << std::endl;

Kokkos::View<ArborX::Point *, MemorySpace> target =
point_clouds[1 - (i % 2)];
Kokkos::View<float *, MemorySpace> tgt_true_values =
true_values[1 - (i % 2)];
MovingLeastSquares<MemorySpace, float> &mls = mlss[i % 2];

source_values =
doOne(mpi_comm, space, target, source_values, tgt_true_values, mls);

if (mpi_rank == 0)
std::cout << "===\n" << std::endl;
}

MPI_Finalize();
Expand Down

0 comments on commit dcdedd3

Please sign in to comment.