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

Make CO2GasPvt GPU compatible #4218

Draft
wants to merge 2 commits into
base: master
Choose a base branch
from
Draft

Make CO2GasPvt GPU compatible #4218

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
e8a2973
first attempt at making tabulated2dfunction gpu compatible
Sep 20, 2024
a77fbde
correct template template syntax
Sep 23, 2024
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
4 changes: 2 additions & 2 deletions opm/material/common/UniformXTabulated2DFunction.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -28,8 +28,8 @@

namespace Opm {

template<class Scalar>
void UniformXTabulated2DFunction<Scalar>::print(std::ostream& os) const
template <typename Scalar, template <class...> class ContainerT>
void UniformXTabulated2DFunction<Scalar, ContainerT>::print(std::ostream& os) const
{
Scalar x0 = xMin();
Scalar x1 = xMax();
Expand Down
101 changes: 66 additions & 35 deletions opm/material/common/UniformXTabulated2DFunction.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -33,6 +33,8 @@
#include <opm/material/common/Valgrind.hpp>
#include <opm/material/common/MathToolbox.hpp>

#include <opm/common/utility/gpuDecorators.hpp>

#include <cassert>
#include <cmath>
#include <iosfwd>
Expand All @@ -49,8 +51,11 @@ namespace Opm {
* "Uniform on the X-axis" means that all Y sampling points must be located along a line
* for this value. This class can be used when the sampling points are calculated at run
* time.
*
* By default the tabulated data is stored in std::vectors, opting out of this is usually done
* with GPUBuffers and Views
*/
template <class Scalar>
template <class Scalar, template<class...> class ContainerT = std::vector>
class UniformXTabulated2DFunction
{
public:
Expand All @@ -72,13 +77,13 @@ class UniformXTabulated2DFunction
Vertical
};

explicit UniformXTabulated2DFunction(const InterpolationPolicy interpolationGuide = Vertical)
OPM_HOST_DEVICE explicit UniformXTabulated2DFunction(const InterpolationPolicy interpolationGuide = Vertical)
: interpolationGuide_(interpolationGuide)
{ }

UniformXTabulated2DFunction(const std::vector<Scalar>& xPos,
const std::vector<Scalar>& yPos,
const std::vector<std::vector<SamplePoint>>& samples,
OPM_HOST_DEVICE UniformXTabulated2DFunction(const ContainerT<Scalar>& xPos,
const ContainerT<Scalar>& yPos,
const ContainerT<ContainerT<SamplePoint>>& samples,
InterpolationPolicy interpolationGuide)
: samples_(samples)
, xPos_(xPos)
Expand All @@ -89,91 +94,91 @@ class UniformXTabulated2DFunction
/*!
* \brief Returns the minimum of the X coordinate of the sampling points.
*/
Scalar xMin() const
OPM_HOST_DEVICE Scalar xMin() const
{ return xPos_.front(); }

/*!
* \brief Returns the maximum of the X coordinate of the sampling points.
*/
Scalar xMax() const
OPM_HOST_DEVICE Scalar xMax() const
{ return xPos_.back(); }

/*!
* \brief Returns the value of the X coordinate of the sampling points.
*/
Scalar xAt(size_t i) const
OPM_HOST_DEVICE Scalar xAt(size_t i) const
{ return xPos_[i]; }

/*!
* \brief Returns the value of the Y coordinate of a sampling point.
*/
Scalar yAt(size_t i, size_t j) const
OPM_HOST_DEVICE Scalar yAt(size_t i, size_t j) const
{ return std::get<1>(samples_[i][j]); }

/*!
* \brief Returns the value of a sampling point.
*/
Scalar valueAt(size_t i, size_t j) const
OPM_HOST_DEVICE Scalar valueAt(size_t i, size_t j) const
{ return std::get<2>(samples_[i][j]); }

/*!
* \brief Returns the number of sampling points in X direction.
*/
size_t numX() const
OPM_HOST_DEVICE size_t numX() const
{ return xPos_.size(); }

/*!
* \brief Returns the minimum of the Y coordinate of the sampling points for a given column.
*/
Scalar yMin(unsigned i) const
OPM_HOST_DEVICE Scalar yMin(unsigned i) const
{ return std::get<1>(samples_.at(i).front()); }

/*!
* \brief Returns the maximum of the Y coordinate of the sampling points for a given column.
*/
Scalar yMax(unsigned i) const
OPM_HOST_DEVICE Scalar yMax(unsigned i) const
{ return std::get<1>(samples_.at(i).back()); }

/*!
* \brief Returns the number of sampling points in Y direction a given column.
*/
size_t numY(unsigned i) const
OPM_HOST_DEVICE size_t numY(unsigned i) const
{ return samples_.at(i).size(); }

/*!
* \brief Return the position on the x-axis of the i-th interval.
*/
Scalar iToX(unsigned i) const
OPM_HOST_DEVICE Scalar iToX(unsigned i) const
{
assert(i < numX());

return xPos_.at(i);
}

const std::vector<std::vector<SamplePoint>>& samples() const
OPM_HOST_DEVICE const ContainerT<ContainerT<SamplePoint>>& samples() const
{
return samples_;
}

const std::vector<Scalar>& xPos() const
OPM_HOST_DEVICE const ContainerT<Scalar>& xPos() const
{
return xPos_;
}

const std::vector<Scalar>& yPos() const
OPM_HOST_DEVICE const ContainerT<Scalar>& yPos() const
{
return yPos_;
}

InterpolationPolicy interpolationGuide() const
OPM_HOST_DEVICE InterpolationPolicy interpolationGuide() const
{
return interpolationGuide_;
}

/*!
* \brief Return the position on the y-axis of the j-th interval.
*/
Scalar jToY(unsigned i, unsigned j) const
OPM_HOST_DEVICE Scalar jToY(unsigned i, unsigned j) const
{
assert(i < numX());
assert(size_t(j) < samples_[i].size());
Expand All @@ -185,7 +190,7 @@ class UniformXTabulated2DFunction
* \brief Return the interval index of a given position on the x-axis.
*/
template <class Evaluation>
unsigned xSegmentIndex(const Evaluation& x,
OPM_HOST_DEVICE unsigned xSegmentIndex(const Evaluation& x,
[[maybe_unused]] bool extrapolate = false) const
{
assert(extrapolate || (xMin() <= x && x <= xMax()));
Expand Down Expand Up @@ -222,7 +227,7 @@ class UniformXTabulated2DFunction
* the range of the segment. In particular this happens for the extrapolation case.
*/
template <class Evaluation>
Evaluation xToAlpha(const Evaluation& x, unsigned segmentIdx) const
OPM_HOST_DEVICE Evaluation xToAlpha(const Evaluation& x, unsigned segmentIdx) const
{
Scalar x1 = xPos_[segmentIdx];
Scalar x2 = xPos_[segmentIdx + 1];
Expand All @@ -233,7 +238,7 @@ class UniformXTabulated2DFunction
* \brief Return the interval index of a given position on the y-axis.
*/
template <class Evaluation>
unsigned ySegmentIndex(const Evaluation& y, unsigned xSampleIdx,
OPM_HOST_DEVICE unsigned ySegmentIndex(const Evaluation& y, unsigned xSampleIdx,
[[maybe_unused]] bool extrapolate = false) const
{
assert(xSampleIdx < numX());
Expand Down Expand Up @@ -271,7 +276,7 @@ class UniformXTabulated2DFunction
* the range of the segment. In particular this happens for the extrapolation case.
*/
template <class Evaluation>
Evaluation yToBeta(const Evaluation& y, unsigned xSampleIdx, unsigned ySegmentIdx) const
OPM_HOST_DEVICE Evaluation yToBeta(const Evaluation& y, unsigned xSampleIdx, unsigned ySegmentIdx) const
{
assert(xSampleIdx < numX());
assert(ySegmentIdx < numY(xSampleIdx) - 1);
Expand All @@ -288,7 +293,7 @@ class UniformXTabulated2DFunction
* \brief Returns true iff a coordinate lies in the tabulated range
*/
template <class Evaluation>
bool applies(const Evaluation& x, const Evaluation& y) const
OPM_HOST_DEVICE bool applies(const Evaluation& x, const Evaluation& y) const
{
if (x < xMin() || xMax() < x)
return false;
Expand Down Expand Up @@ -316,7 +321,7 @@ class UniformXTabulated2DFunction
* range, a \c Opm::NumericalIssue exception is thrown.
*/
template <class Evaluation>
Evaluation eval(const Evaluation& x, const Evaluation& y, bool extrapolate=false) const
OPM_HOST_DEVICE Evaluation eval(const Evaluation& x, const Evaluation& y, bool extrapolate=false) const
{
Evaluation alpha, beta1, beta2;
unsigned i, j1, j2;
Expand All @@ -325,7 +330,7 @@ class UniformXTabulated2DFunction
}

template <class Evaluation>
void findPoints(unsigned& i,
OPM_HOST_DEVICE void findPoints(unsigned& i,
unsigned& j1,
unsigned& j2,
Evaluation& alpha,
Expand All @@ -336,6 +341,7 @@ class UniformXTabulated2DFunction
bool extrapolate) const
{
#ifndef NDEBUG
#if !OPM_IS_INSIDE_DEVICE_FUNCTION
if (!extrapolate && !applies(x, y)) {
if constexpr (std::is_floating_point_v<Evaluation>) {
throw NumericalProblem("Attempt to get undefined table value (" +
Expand All @@ -347,6 +353,7 @@ class UniformXTabulated2DFunction
std::to_string(y.value()) + ")");
}
};
#endif
#endif

// bi-linear interpolation: first, calculate the x and y indices in the lookup
Expand Down Expand Up @@ -393,18 +400,21 @@ class UniformXTabulated2DFunction
}

template <class Evaluation>
Evaluation eval(const unsigned& i, const unsigned& j1, const unsigned& j2, const Evaluation& alpha,const Evaluation& beta1,const Evaluation& beta2) const
OPM_HOST_DEVICE Evaluation eval(const unsigned& i, const unsigned& j1, const unsigned& j2, const Evaluation& alpha,const Evaluation& beta1,const Evaluation& beta2) const
{
// evaluate the two function values for the same y value ...
const Evaluation& s1 = valueAt(i, j1)*(1.0 - beta1) + valueAt(i, j1 + 1)*beta1;
const Evaluation& s2 = valueAt(i + 1, j2)*(1.0 - beta2) + valueAt(i + 1, j2 + 1)*beta2;

#if OPM_IS_INSIDE_DEVICE_FUNCTION
Valgrind::CheckDefined(s1);
Valgrind::CheckDefined(s2);

#endif
// ... and combine them using the x position
const Evaluation& result = s1*(1.0 - alpha) + s2*alpha;
#if OPM_IS_INSIDE_DEVICE_FUNCTION
Valgrind::CheckDefined(result);
#endif

return result;
}
Expand All @@ -426,7 +436,7 @@ class UniformXTabulated2DFunction
// this is slow, but so what?
xPos_.insert(xPos_.begin(), nextX);
yPos_.insert(yPos_.begin(), std::numeric_limits<Scalar>::lowest() / 2);
samples_.insert(samples_.begin(), std::vector<SamplePoint>());
samples_.insert(samples_.begin(), ContainerT<SamplePoint>());
return 0;
}
throw std::invalid_argument("Sampling points should be specified either monotonically "
Expand All @@ -437,6 +447,7 @@ class UniformXTabulated2DFunction
* \brief Append a sample point.
*
* Returns the i index of the new point within its line.
* Avoid running this in the GPU code
*/
size_t appendSamplePoint(size_t i, Scalar y, Scalar value)
{
Expand All @@ -457,7 +468,6 @@ class UniformXTabulated2DFunction
}
return 0;
}

throw std::invalid_argument("Sampling points must be specified in either monotonically "
"ascending or descending order.");
}
Expand All @@ -470,7 +480,7 @@ class UniformXTabulated2DFunction
*/
void print(std::ostream& os) const;

bool operator==(const UniformXTabulated2DFunction<Scalar>& data) const {
OPM_HOST_DEVICE bool operator==(const UniformXTabulated2DFunction<Scalar, ContainerT>& data) const {
return this->xPos() == data.xPos() &&
this->yPos() == data.yPos() &&
this->samples() == data.samples() &&
Expand All @@ -481,14 +491,35 @@ class UniformXTabulated2DFunction
// the vector which contains the values of the sample points
// f(x_i, y_j). don't use this directly, use getSamplePoint(i,j)
// instead!
std::vector<std::vector<SamplePoint> > samples_;
ContainerT<ContainerT<SamplePoint> > samples_;

// the position of each vertical line on the x-axis
std::vector<Scalar> xPos_;
ContainerT<Scalar> xPos_;
// the position on the y-axis of the guide point
std::vector<Scalar> yPos_;
ContainerT<Scalar> yPos_;
InterpolationPolicy interpolationGuide_;
};
} // namespace Opm

namespace Opm::gpuistl {
// Take in a CO2 type that exists on the CPU and move the data into GPU buffers
// We make the CO2 object on the CPU and just move it to avoid implementing vector::insert on GPU
template<class Scalar, class GpuCO2Type, template<class> class BufferType>
GpuCO2Type moveToGpu(UniformXTabulated2DFunction<Scalar> cpuCO2){

using InterpolationPolicy = typename UniformXTabulated2DFunction<Scalar>::InterpolationPolicy;
using SamplePoint = typename UniformXTabulated2DFunction<Scalar>::SamplePoint;

// TODO cast the vectors into GPUBuffers
BufferType<BufferType<SamplePoint>> samples(cpuCO2.samples());
BufferType<Scalar> xPos(cpuCO2.xpos());
BufferType<Scalar> yPos(cpuCO2.ypos());
InterpolationPolicy interpolationGuide = cpuCO2.interpolationGuide();

return UniformXTabulated2DFunction<Scalar, BufferType>(samples, xPos, yPos, interpolationGuide);
}

// TODO: Create a view so this CO2 object can be cheaply copied
}

#endif