Commit f8b9e6ef authored by Alain O' Miniussi's avatar Alain O' Miniussi
Browse files

computeMedium on device

parent 289d8c25
......@@ -369,41 +369,41 @@ namespace kk = Kokkos;
/// \brief Given a set of grid coords base on inferior cells boundaries, return the coords of the centers.
/// \param rad the infs
arr1d<real>
computeMedium(arr1d<real const> rads, std::function<real (real,real)> fct) {
computeMedium(arr1d<real const> rads, int dim = -1) {
arr1d<real> radss("augmented_rads", rads.size()+1);
typedef std::pair<int,int> range;
kk::deep_copy(kk::subview(radss, range(0,rads.size())), rads);
// this one should not be needed, it's one past the end
kk::parallel_for("set_last_rad", kk::RangePolicy<>(rads.size(), radss.size()),
LBD(int i) { radss(i) = rads(i-1)+(rads(i-1)-rads(i-2)); });
arr1d<real> meds("meds", rads.size());
arr1d<real>::HostMirror hmeds = kk::create_mirror_view(meds);
int sz = rads.extent(0);
for (int i=0; i < sz-1; i++) {
real rinf = rads(i);
real rsup = rads(i+1);
hmeds(i) = fct(rinf, rsup);
if (dim<0) {
kk::parallel_for("set_meds", kk::RangePolicy<>(0, meds.size()),
LBD(int i) { meds(i) = (radss(i) + radss(i+1))/2; });
} else {
using std::pow;
kk::parallel_for("set_meds", kk::RangePolicy<>(0, meds.size()),
LBD(int i) { meds(i) = (real(dim)/(dim+1)
* (pow(radss(i+1),dim+1) - pow(radss(i),dim+1))
/ (pow(radss(i+1),dim) - pow(radss(i),dim)));
});
}
// this one should not be needed, it's one past the end
int i = sz-1;
real rinf = rads(i);
real rsup = rads(i) + (rads(i)-rads(i-1)); // assume it's mostly regular at that point
hmeds(i) = fct(rinf, rsup);
kk::deep_copy(meds, hmeds);
return meds;
}
arr1d<real>
getRadiusMeds(arr1d<real const> rads, int dim) {
return computeMedium(rads,
[=](real inf, real sup) {
return real(dim)/real(dim+1)
* (std::pow(sup,dim+1)-std::pow(inf,dim+1)) / (std::pow(sup,dim)-std::pow(inf,dim));
});
return computeMedium(rads, dim);
}
arr1d<real>
getSectorMeds(arr1d<real const> limits) {
return computeMedium(limits, [=](real inf, real sup) { return (inf+sup)/2; });
return computeMedium(limits);
}
arr1d<real>
getLayerMeds(arr1d<real const> limits) {
return computeMedium(limits, [=](real inf, real sup) { return (inf+sup)/2; });
return computeMedium(limits);
}
std::ostream&
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment