Hello,
I am trying to use lammps c++ api. Lammps has been compiled with Kokkos/CUDA in mind.
My test case is a simple simulation box creation and creating some atoms using the lammps_create_atoms C API.
Unfortunately, as soon as I do a lmp->input->one("run 0") to ask for the neighbour list to be rebuilt after the creation of some atoms, I get an CUDA Illegal memory access.
I tracked down the issue using gdb:
LAMMPS_NS::AtomKokkos::map_set_device() at lammps-29Aug2024/src/KOKKOS/atom_map_kokkos.cpp:258
lammps_create_atoms() at lammps-29Aug2024/src/library.cpp:5556
It seems there is an issue with the tag deep_copy.
My code is the following:
//LAMMPS includes
#include "lammps.h"
#include "input.h"
#include "atom_masks.h"
#include "atom.h"
#include "library.h"
#include "universe.h"
#include "fix.h"
#include "domain.h"
#include "irregular.h"
#include "error.h"
#include "info.h"
#include "platform.h"
#include "pair.h"
#include "force.h"
#include "version.h"
#include "kokkos.h"
/*#include "accelerator_kokkos.h"
#include "neighbor_kokkos.h"
#include "neigh_list_kokkos.h"
#include "kokkos_type.h"
#include "Kokkos_Core.hpp"
#include "Kokkos_StdAlgorithms.hpp"
*/
#include <mpi.h>
#include <iostream>
#include <algorithm>
#include <ranges>
#include <random>
//--------------------------------------------------------------------------------------------------
/**
* @brief Flatten a vector, using ranges and concept!
* @link https://stackoverflow.com/a/71680819
*/
template<typename R>
concept nested_range = std::ranges::input_range<R> && std::ranges::range<std::ranges::range_reference_t<R>>;
struct flatten_t
{
template<nested_range R>
auto operator()(R && r) const
{
return std::forward<R>(r) | std::ranges::views::transform(*this) | std::ranges::views::join;
}
template<typename T>
auto operator()(T && t) const
{
return std::forward<T>(t);
}
};
template<typename T>
auto operator |(T && t, flatten_t f)
{
return f(std::forward<T>(t));
}
constexpr flatten_t flatten;
//--------------------------------------------------------------------------------------------------
/**
* @brief Initialize the lammps object
* @param lammps handle
*/
void Initialize(LAMMPS_NS::LAMMPS *lmp){
lmp->input->one("clear");
lmp->input->one("units metal");
lmp->input->one("dimension 3");
lmp->input->one("boundary p p f");
lmp->input->one("atom_style atomic");
lmp->input->one("atom_modify map array");
lmp->input->one("region box block 0 100 0 100 0 300");
lmp->input->one("create_box 1 box");
lmp->input->one("pair_style tersoff");
lmp->input->one("pair_coeff * * ~/lammps-29Aug2024/potentials/Si.tersoff Si");
lmp->input->one("mass 1 28.0855"); //Silicium
lmp->input->one("compute 1 all pe/atom");
lmp->input->one("thermo 10");
lmp->input->one("thermo_style custom step temp ke pe etotal");
return void();
}
//--------------------------------------------------------------------------------------------------
/**
* @brief test lammps_create_atoms function
* @param lammps handle
* @param k step
*/
void test_create_atoms(LAMMPS_NS::LAMMPS *lmp, const unsigned long long int& k){
//Create a grid
std::vector<std::array<double,3>> PES;
for(int i=0;i<100; ++i){
for(int j=0;j<100; ++j)
PES.push_back({0.5+i, 0.5+j,1.0*k});
}
std::vector<std::array<double,3>> part_PES;
std::ranges::sample(PES, std::back_inserter(part_PES), 50, std::mt19937_64 { std::random_device{}()});
//view::common trick from htpps://stackoverfow.com/a/63116423
auto view_flattened = part_PES | flatten | std::ranges::views::common;
std::vector<double> out(view_flattened.begin(), view_flattened.end());
std::vector<int> types(part_PES.size(),1);
lammps_create_atoms(lmp, part_PES.size(), nullptr, types.data(), out.data(), nullptr, nullptr, 1);
return void();
}
//--------------------------------------------------------------------------------------------------
//--------------------------------------------------------------------------------------------------
//--------------------------------------------------------------------------------------------------
int main(int argc, char **argv){
setenv("OMP_PROC_BIND", "spread", true);
setenv("OMP_PLACES", "threads", true);
const char* lmpargv[] {"liblammps", "-nocite", "-echo", "none", "-log", "none", "-k", "on", "g", "1", "t", "16", "-sf", "kk", "-pk", "kokkos", "newton", "on"
, "neigh", "half", "comm", "device"};
constexpr int lmpargc = sizeof(lmpargv)/sizeof(const char *);
MPI_Init(&argc, &argv);
LAMMPS_NS::LAMMPS *lmp;
lmp = new LAMMPS_NS::LAMMPS(lmpargc, (char** )lmpargv, MPI_COMM_WORLD);
std::cout << "LAMMPS version id = " << lmp->num_ver << std::endl;
Initialize(lmp);
lmp->input->one("run 0");
lmp->input->one("write_dump all custom/gz ~/test_lmp/test_0.atom.gz type id x y z");
unsigned long long int k(0);
for(;k<235;++k){
test_create_atoms(lmp,k);
lmp->input->one("run 0");
lmp->input->one("write_dump all custom/gz ~/test_lmp/test_" + std::to_string(k) + ".atom.gz type id x y z");
}
lmp->input->one("write_dump all custom/gz ~/test_lmp/test_" + std::to_string(k) + ".atom.gz type id x y z");
delete lmp;
Kokkos::finalize();
MPI_Barrier(MPI_COMM_WORLD);
MPI_Finalize();
return 0;
}
From the code above, if I remove the run 0 from the for loop, the program runs correctly.
My question is, is this a bug or I am just not using correctly the API?
Thanks in advance,
Regards