diff --git a/doc/source/programmingguide.rst b/doc/source/programmingguide.rst index 6f65742bc..d8f989df7 100644 --- a/doc/source/programmingguide.rst +++ b/doc/source/programmingguide.rst @@ -489,7 +489,36 @@ For the CSV constructor, ``nDim`` can only have the values 1 or 2. .. note:: The input CSV file is allowed to contain comments, starting with "#". Any character following - "#" is ignored by the ``LookupTable`` class. + "#" is ignored by the ``LookupTable`` class. Hence, a line starting with "#" is entirely ignored, + and so are empty lines and lines made only of white spaces. Comments and blank lines can therefore + be inserted anywhere in the file, including before the first line of the table. + +By default, the coordinates and the data are read as *lines* of the CSV file, as in the examples above. For 1D lookup +tables, they can also be read as *columns* of the CSV file, using the optional argument ``readColumns`` of the constructor: + +.. code-block:: c++ + + template + LookupTable::LookupTable(std::string filename, char delimiter, + bool errorIfOutOfBound, bool readColumns); + +With ``readColumns=true``, the 1D lookup table is read from the two first columns of the file, the first one +giving the coordinates and the second one the data: + +.. list-table:: example1D_col.csv + :widths: 25 25 + :header-rows: 0 + + * - x\ :sub:`1` + - data\ :sub:`1` + * - x\ :sub:`2` + - data\ :sub:`2` + * - x\ :sub:`3` + - data\ :sub:`3` + +.. note:: + ``readColumns`` is only available for 1D lookup tables. 2D lookup tables are always read as lines, + following the layout of ``example2D.csv`` above. Numpy constructor @@ -507,6 +536,60 @@ the constructor expects a vector of size ``nDim`` of .npy files for the 1D coord Note that the template parameter ``nDim`` should match the number of dimensions of the numpy array stored in the file ``dataSet``. +Interpolation in function space ++++++++++++++++++++++++++++++++ + +By default, ``LookupTable`` performs a multi-linear interpolation directly on the coordinates ``x`` and the data +of the table. For data spanning several orders of magnitude (cooling tables, opacity tables, etc...), it is often +more accurate to interpolate the data in *function space*, i.e. to interpolate ``func(data)`` as a function of +``func(x)``, and to transform the interpolated value back with the inverse function ``invFunc``. This is enabled with +the optional argument ``interpolateInFuncSpace`` accepted by each of the constructors: + +.. code-block:: c++ + + // 1D CSV table read as columns, interpolated in log space + LookupTable<1> table("cooling.csv", ',', true, true, true); + + // 3D numpy table, interpolated in log space + LookupTable<3> tablenpy(coordinates, "data.npy", true, true); + +When this option is enabled, ``func`` is applied to the coordinates and to the data when the table is constructed +(so that ``xinHost``, ``xinDev``, ``dataHost`` and ``dataDev`` all store transformed values), the interpolation is +performed in that space, and ``Get`` and ``GetHost`` return ``invFunc`` of the interpolated value. The coordinates +passed to ``Get`` and ``GetHost``, as well as the value they return, are therefore always expressed in the original +(untransformed) space. + +``func`` and ``invFunc`` default to the natural logarithm and the exponential (the functors ``LookupTableLog`` and +``LookupTableExp``). They can be changed when the lookup table is created, using the two optional template parameters +of ``LookupTable``, which expect functors defining ``operator()(const real)``. This operator should be decorated with +``KOKKOS_INLINE_FUNCTION``, so that the transformation can be used both on the host and on the device: + +.. code-block:: c++ + + // A user-defined transformation and its inverse + struct MyLog10 { + KOKKOS_INLINE_FUNCTION real operator() (const real x) const { + return(log10(x)); + } + }; + + struct MyPow10 { + KOKKOS_INLINE_FUNCTION real operator() (const real x) const { + return(pow(10.0,x)); + } + }; + + // A 1D CSV table interpolated in log10 space + LookupTable<1, MyLog10, MyPow10> table("cooling.csv", ',', true, true, true); + +Instances of these functors can also be given as the last two arguments of the constructors, which is useful when the +transformation carries a state (e.g. a tunable exponent). + +.. warning:: + The transformation is applied to every coordinate and to every element of the data. Since the default transformation + is the natural logarithm, the whole table should be strictly positive when the default ``func`` is used, otherwise + ``LookupTable`` triggers an error while it is being constructed. + Using the lookup table ++++++++++++++++++++++ @@ -536,6 +619,116 @@ The ``Get`` and ``GetHost`` functions expect a C array of size ``nDim`` and retu real result = csv.GetHost(y); +Accessing the neighbours used for the interpolation ++++++++++++++++++++++++++++++++++++++++++++++++++++ + +In addition to the interpolated value, ``LookupTable`` can return the elements of the table which surround the +requested coordinates, i.e. the ones the interpolation is computed from. This is useful to implement a custom +interpolation or reconstruction (e.g. a power law between two nodes), to estimate a local slope, or simply to +check which part of the table is being used. Four methods are available, ``GetNeighbours`` and ``GetNeighboursIndx`` +on the device, and ``GetNeighboursHost`` and ``GetNeighboursIndxHost`` on the host: + +.. code-block:: c++ + + // Coordinates and data of the neighbours (device and host) + void GetNeighbours(const real x[nDim], real xN[2*nDim], real dataN[1 << nDim]); + void GetNeighboursHost(const real x[nDim], real xN[2*nDim], real dataN[1 << nDim]); + + // Indices of the same neighbours (device and host) + void GetNeighboursIndx(const real x[nDim], int idx[nDim], int dataIdx[1 << nDim]); + void GetNeighboursIndxHost(const real x[nDim], int idx[nDim], int dataIdx[1 << nDim]); + +These methods do not return anything, but they fill the arrays which are given to them: + +* ``xN[2*n]`` and ``xN[2*n+1]`` are the two coordinates of the table bracketing ``x[n]`` along the dimension ``n``. +* ``dataN[v]`` is the data on the vertex ``v`` of the cell surrounding ``x``. Since a cell of a ``nDim`` table has + ``2^nDim`` vertices, ``dataN`` has ``1 << nDim`` elements. The bit ``n`` of ``v`` tells whether the vertex is on + the right (1) or on the left (0) of the dimension ``n``: in 2D, ``dataN[0]``, ``dataN[1]``, ``dataN[2]`` and + ``dataN[3]`` are therefore the data in (x\ :sub:`i`, y\ :sub:`j`), (x\ :sub:`i+1`, y\ :sub:`j`), + (x\ :sub:`i`, y\ :sub:`j+1`) and (x\ :sub:`i+1`, y\ :sub:`j+1`). +* ``idx[n]`` is the index of the left neighbour along the dimension ``n`` (the right one being ``idx[n]+1``). +* ``dataIdx[v]`` is the index of the vertex ``v`` in the data array of the table, following the same convention as + ``dataN``, so that ``dataN[v]`` and ``dataHost(dataIdx[v])`` refer to the same element. + +For instance, the 2D table ``example2D.csv`` above, interpolated in (x=2.1, y=3.5), gives: + +.. code-block:: c++ + + real x[2]; + x[0] = 2.1; + x[1] = 3.5; + + real xN[4]; // 2 coordinates for each of the 2 dimensions + real dataN[4]; // 2^2 vertices + int idx[2]; + int dataIdx[4]; + + csv.GetNeighboursHost(x, xN, dataN); + csv.GetNeighboursIndxHost(x, idx, dataIdx); + + // xN[0] and xN[1] now bracket x[0], while xN[2] and xN[3] bracket x[1], and the interpolated + // value returned by csv.GetHost(x) can be recomputed from dataN: + real dx = (x[0]-xN[0])/(xN[1]-xN[0]); + real dy = (x[1]-xN[2])/(xN[3]-xN[2]); + real value = (1-dx)*(1-dy)*dataN[0] + dx*(1-dy)*dataN[1] + + (1-dx)*dy*dataN[2] + dx*dy*dataN[3]; + +The same methods without the ``Host`` suffix can be called from within an ``idefix_for`` loop, the arrays being then +local to the loop. + +Both ``Get`` and the methods above search the table for the cell surrounding the requested coordinates. When the +interpolated value *and* its neighbours are needed for the same coordinates, this search can be performed only once, +by handing the same ``LookupTableSearchCache`` structure to each of them: whichever is called first searches the +table and stores the result in the structure, and the ones called next reuse it instead of searching the table again. +The order in which they are called is irrelevant, ``Get`` can fill the structure for ``GetNeighbours`` and +``GetNeighboursIndx``, or the other way around. + +.. code-block:: c++ + + idefix_for("loop",0, 10, KOKKOS_LAMBDA (int i) { + real x[2]; + x[0] = 2.1; + x[1] = 3.5; + + // The search performed by Get is stored in neighbours... + LookupTableSearchCache<2> neighbours; + real value = csv.Get(x, neighbours); + + // ... and is reused by the two methods below, which do not search the table again + real xN[4]; + real dataN[4]; + int idx[2]; + int dataIdx[4]; + csv.GetNeighbours(x, neighbours, xN, dataN); + csv.GetNeighboursIndx(x, neighbours, idx, dataIdx); + }); + +The very same methods are available on the host, as ``GetHost``, ``GetNeighboursHost`` and ``GetNeighboursIndxHost``. +The structure gives access to the search itself, ``neighbours.idx[n]`` being the index of the left neighbour along +the dimension ``n``, and ``neighbours.delta[n]`` the elementary ratio used to weight the two neighbours of that +dimension. + +.. note:: + The stored search is only reused when the coordinates it was computed for are the ones being requested. Calling + ``GetNeighbours`` or ``GetNeighboursIndx`` with different coordinates simply searches the table again, and updates + the structure accordingly, so that a ``LookupTableSearchCache`` can be reused from one point to the next. + +.. warning:: + The ``LookupTableSearchCache`` structure is declared *per thread*, and not by the lookup table itself: when it is used inside an + ``idefix_for`` loop, it should be declared *inside* the loop, so that each thread has its own copy. A lookup table + is shared by all of the threads of a loop, and can therefore not store anything of the sort itself. + +.. note:: + When the table is interpolated in function space (see above), ``xN`` and ``dataN`` are returned in the original + space of the table, i.e. ``invFunc`` is applied to the values which are stored in ``xin`` and ``data``. The + indices returned by ``GetNeighboursIndx`` are of course not affected by the transformation. + +.. note:: + The search of the neighbours follows exactly the same rules as ``Get``: an out of bound coordinate triggers an + error when ``errorIfOutOfBound`` is enabled, and is otherwise clamped to the closest cell of the table. If one of + the requested coordinates is a nan, ``xN`` and ``dataN`` are filled with nans, and ``idx`` and ``dataIdx`` with -1. + + .. note:: Usage examples are provided in `test/utils/lookupTable`. diff --git a/src/real_types.hpp b/src/real_types.hpp index 440c7aac9..988ca1675 100644 --- a/src/real_types.hpp +++ b/src/real_types.hpp @@ -33,6 +33,8 @@ #define TAN(x) tanf(x) #define SIN(x) sinf(x) #define COS(x) cosf(x) +#define LOG(x) logf(x) +#define EXP(x) expf(x) #define COPYSIGN(x,y) copysignf(x,y) #define ISNAN(x) isnanf(x) #define FMOD(x,y) fmodf(x,y) @@ -52,6 +54,8 @@ #define TAN(x) tan(x) #define SIN(x) sin(x) #define COS(x) cos(x) +#define LOG(x) log(x) +#define EXP(x) exp(x) #define COPYSIGN(x,y) copysign(x,y) #define ISNAN(x) isnan(x) #define FMOD(x,y) fmod(x,y) diff --git a/src/utils/lookupTable.hpp b/src/utils/lookupTable.hpp index 72d13a838..ee2e0bd15 100644 --- a/src/utils/lookupTable.hpp +++ b/src/utils/lookupTable.hpp @@ -14,18 +14,82 @@ #include "lookupTable.hpp" #include "npy.hpp" +// Default transformation (and its inverse) used when the lookup table interpolates in function +// space: the interpolation is then performed on log(x) and log(data), and the interpolated value +// is transformed back with exp. +// Any functor exposing a KOKKOS_INLINE_FUNCTION operator()(const real) can be used instead, so +// that the transformation is available both on the host and on the device. +struct LookupTableLog { + KOKKOS_INLINE_FUNCTION real operator() (const real x) const { + return(LOG(x)); + } +}; + +struct LookupTableExp { + KOKKOS_INLINE_FUNCTION real operator() (const real x) const { + return(EXP(x)); + } +}; + +// Neighbours of a point of a lookup table, as they are found by the search performed by Get: +// idx[n] is the index of the left neighbour along the dimension n, and delta[n] is the elementary +// ratio used to weight the two neighbours of that dimension. +// Get (and GetHost) fill this structure when it is given as an argument, and GetNeighbours and +// GetNeighboursIndx then reuse the search it contains instead of performing it a second time +// (which is only done when they are called with the same coordinates: they search the table again +// as usual when a different x is requested). +// This structure is kept by the caller, so that it is thread-private when it is declared inside +// an idefix_for loop. The lookup table itself is shared by all of the threads of a loop, and can +// therefore not be used to store anything of the sort. template +struct LookupTableSearchCache { + real x[kDim]; // coordinates for which the neighbours below were computed + int idx[kDim]; // index of the left neighbour along each dimension + real delta[kDim]; // elementary ratio between the two neighbours of each dimension + bool valid{false}; // whether the neighbours above have been successfully computed + + // Check whether this structure already holds the neighbours of the coordinates xIn. + // The comparison is a bit-exact "==" on purpose: this is meant to catch the case where the + // caller passes back literally the same x it just queried (e.g. Get then GetNeighbours on the + // same iteration), not to approximate "close enough" coordinates. Do not replace this with a + // tolerance-based comparison, or two distinct nearby queries could wrongly reuse a stale search. + KOKKOS_INLINE_FUNCTION + bool Matches(const real xIn[kDim]) const { + if(!valid) return(false); + for(int n = 0 ; n < kDim ; n++) { + if(x[n] != xIn[n]) return(false); + } + return(true); + } +}; + +template class LookupTable { public: LookupTable() = default; - LookupTable(std::string filename, char delimiter, bool errorIfOutOfBound = true); + // Constructor from a CSV/ASCII file. + // By default (readColumns=false), the arrays are read as *lines* of the input file: + // 1D: 1st line = coordinates (xinHost), 2nd line = data (dataHost) + // 2D: 1st line = coordinates of the 1st dimension, then each line starts with the + // coordinate of the 2nd dimension, followed by the data of that line. + // For 1D tables only, the arrays can also be read as *columns* of the input file by setting + // readColumns=true: + // 1st column = coordinates (xinHost), 2nd column = data (dataHost) + // All of the constructors accept the optional argument interpolateInFuncSpace: when it is set + // to true, the coordinates and the data are stored as func(x) and func(data), the interpolation + // is performed in that space, and Get returns invFunc of the interpolated value (see below). + LookupTable(std::string filename, char delimiter, bool errorIfOutOfBound = true, + bool readColumns = false, bool interpolateInFuncSpace = false, + TFunc func = TFunc(), TInvFunc invFunc = TInvFunc()); LookupTable(std::vector filenames, std::string dataSet, - bool errorIfOutOfBound = true); + bool errorIfOutOfBound = true, bool interpolateInFuncSpace = false, + TFunc func = TFunc(), TInvFunc invFunc = TInvFunc()); template LookupTable(Kokkos::View array, std::array,kDim>, - bool errorIfOutOfBound = true); + bool errorIfOutOfBound = true, bool interpolateInFuncSpace = false, + TFunc func = TFunc(), TInvFunc invFunc = TInvFunc()); IdefixArray1D dimensionsDev; IdefixArray1D offsetDev; // Actually sum_(n-1) (dimensions) @@ -40,20 +104,62 @@ class LookupTable { bool errorIfOutOfBound{true}; - // Generic getter for all kinds of input arrays + // When enabled, the table is interpolated in function space: xin and data store func(x) and + // func(data), and Get returns invFunc(interpolated value). func and invFunc are stored by + // value, so that they are available on the device as well as on the host. + bool interpolateInFuncSpace{false}; + TFunc func; + TInvFunc invFunc; + + // Transform the coordinates and the data of the table in function space. + // This is called on the host by the constructors, before the arrays are copied on the device. + void ToFuncSpace() { + for(int i = 0 ; i < xinHost.extent(0) ; i++) { + xinHost(i) = func(xinHost(i)); + if(!std::isfinite(xinHost(i))) { + IDEFIX_ERROR("LookupTable: the transformation of the coordinates in function space " + "produced invalid values (with the default log, the coordinates of the " + "table should all be strictly positive)"); + } + } + for(int i = 0 ; i < dataHost.extent(0) ; i++) { + dataHost(i) = func(dataHost(i)); + if(!std::isfinite(dataHost(i))) { + IDEFIX_ERROR("LookupTable: the transformation of the data in function space " + "produced invalid values (with the default log, the data of the " + "table should all be strictly positive)"); + } + } + } + + // --------------------------------------------------------------------------------------------- + // Implementation details below: these generic (Tint/Treal-templated) helpers are shared by the + // public Get/GetHost/GetNeighbours*/GetNeighboursIndx* entry points declared further down, which + // simply instantiate them with the device or host storage. They are not meant to be called + // directly, hence kept private; make this public again if another translation unit genuinely + // needs the generic (host-or-device) form directly. + // --------------------------------------------------------------------------------------------- + private: + // Search of the bracketing indices used by the interpolation, for all kinds of input arrays. + // On output, idx[n] is the index of the left neighbour along the dimension n (so that the + // interpolation is performed between xin(offset(n)+idx[n]) and xin(offset(n)+idx[n]+1)), and + // delta[n] is the elementary ratio used to weight these two neighbours. + // Returns false when the requested coordinates contain nans, and true otherwise. + // This is the only place that performs the actual search: Get, GetNeighbours and + // GetNeighboursIndx all reach it through GetIndices below, so that a bug fix here never + // needs to be duplicated elsewhere. template KOKKOS_INLINE_FUNCTION - real Get(const real x[kDim], Tint &dimensions, Tint &offset, Treal &xin, Treal &data) const { - // Fetch function that should be called inside idefix_loop - int idx[kDim]; - real delta[kDim]; - + bool GetIndices(const real x[kDim], Tint &dimensions, Tint &offset, Treal &xin, + int idx[kDim], real delta[kDim]) const { for(int n = 0 ; n < kDim ; n++) { real xstart = xin(offset(n)); real xend = xin(offset(n)+dimensions(n)-1); - real x_n = x[n]; + // When interpolating in function space, xin already stores func(x), so the coordinate + // we are looking for should be transformed accordingly + real x_n = interpolateInFuncSpace ? func(x[n]) : x[n]; - if(std::isnan(x_n)) return(NAN); + if(std::isnan(x_n)) return(false); // Compute index of closest element assuming even distribution int i; @@ -82,16 +188,20 @@ class LookupTable { // Check if resulting bounding elements are correct if(xin(offset(n) + i) > x_n || xin(offset(n) + i+1) < x_n) { // Nop, so the points are not evenly distributed - // Search for the correct index (a dicotomy would be more appropriate...) - - i = 0; - while(xin(offset(n) + i) < x_n && i < dimensions(n)-1 ) { - i++; + // Search for the correct index with a dichotomy + int ileft = 0; + int iright = dimensions(n)-1; + while(iright-ileft>1) { + int imid = (ileft + iright) / 2; + if(xin(offset(n) + imid) <= x_n) { + ileft = imid; + } else { + iright = imid; + } } - i = i-1; // i is overestimated by one + i = ileft; } } - // Store the index idx[n] = i; @@ -99,33 +209,190 @@ class LookupTable { delta[n] = (x_n - xin(offset(n) + i) ) / (xin(offset(n) + i+1) - xin(offset(n) + i)); } - // De a linear interpolation from the neightbouring points to get our value. + return(true); + } + + // Index in the data array of the vertex "vertex" of the neighbours. Each bit of "vertex" tells + // whether we are on the right (1) or on the left (0) of the corresponding dimension, so that + // there are 2^kDim vertices surrounding the point we are interpolating. + template + KOKKOS_INLINE_FUNCTION + int GetDataIndex(Tint &dimensions, const int idx[kDim], const unsigned int vertex) const { + int index = 0; + for(unsigned int m = 0 ; m < kDim ; m++) { + index = index * dimensions(m); + unsigned int myBit = 1 << m; + // If bit is set, we're doing the right vertex, otherwise we're doing the left vertex + if((vertex & myBit) > 0) { + index += idx[m]+1; + } else { + index += idx[m]; + } + } + return(index); + } + + // Fill "searchCache" with the neighbours of x, unless it already holds them (in which case the + // table is not searched again). This is what allows Get, GetNeighbours and GetNeighboursIndx to + // share a single search when they are called successively with the same coordinates. + template + KOKKOS_INLINE_FUNCTION + void SearchNeighbours(const real x[kDim], Tint &dimensions, Tint &offset, Treal &xin, + LookupTableSearchCache &searchCache) const { + // Nothing to do if the neighbours of these very coordinates are already known + if(searchCache.Matches(x)) return; + + searchCache.valid = GetIndices(x, dimensions, offset, xin, searchCache.idx, searchCache.delta); + for(int n = 0 ; n < kDim ; n++) { + searchCache.x[n] = x[n]; + } + } + + // Generic getter which stores in "searchCache" the elements of the table it used, so that a + // subsequent call to GetNeighbours or GetNeighboursIndx with the same coordinates does not + // search the table again + template + KOKKOS_INLINE_FUNCTION + real Get(const real x[kDim], Tint &dimensions, Tint &offset, Treal &xin, Treal &data, + LookupTableSearchCache &searchCache) const { + SearchNeighbours(x, dimensions, offset, xin, searchCache); + + if(!searchCache.valid) return(NAN); + + // Do a linear interpolation from the neightbouring points to get our value. real value = 0; // loop on all of the vertices of the neighbours for(unsigned int n = 0 ; n < (1 << kDim) ; n++) { - int index = 0; real weight = 1.0; for(unsigned int m = 0 ; m < kDim ; m++) { - index = index * dimensions(m); unsigned int myBit = 1 << m; // If bit is set, we're doing the right vertex, otherwise we're doing the left vertex if((n & myBit) > 0) { // We're on the right - weight = weight*delta[m]; - index += idx[m]+1; + weight = weight*searchCache.delta[m]; } else { // We're on the left - weight = weight*(1-delta[m]); - index += idx[m]; + weight = weight*(1-searchCache.delta[m]); } } - value = value + weight*data(index); + value = value + weight*data(GetDataIndex(dimensions, searchCache.idx, n)); } + // The interpolation was performed on func(data), so we transform the result back + if(interpolateInFuncSpace) value = invFunc(value); + return(value); } + // Generic getter for all kinds of input arrays. + // Implemented as a thin delegation to the neighbours-caching overload below, using a throwaway + // searchCache structure: LookupTableSearchCache is a small stack-only POD (no allocation), so + // this costs nothing extra while guaranteeing the two overloads can never disagree on the result. + template + KOKKOS_INLINE_FUNCTION + real Get(const real x[kDim], Tint &dimensions, Tint &offset, Treal &xin, Treal &data) const { + LookupTableSearchCache searchCache; + return Get(x, dimensions, offset, xin, data, searchCache); + } + + // Generic getter for the neighbours used by the interpolation, for all kinds of input arrays. + // On output, xN[2*n] and xN[2*n+1] are the coordinates bracketing x[n] along the dimension n, + // and dataN[v] is the data at the vertex v of these neighbours (see GetDataIndex for the + // convention on v). When the table is interpolated in function space, both are returned in the + // original space of the table (i.e. invFunc is applied to the values stored in the table). + // All of the outputs are set to nan when the requested coordinates contain nans. + template + KOKKOS_INLINE_FUNCTION + void GetNeighbours(const real x[kDim], Tint &dimensions, Tint &offset, Treal &xin, Treal &data, + real xN[2*kDim], real dataN[1 << kDim]) const { + LookupTableSearchCache searchCache; + GetNeighbours(x, dimensions, offset, xin, data, searchCache, xN, dataN); + } + + // Same as above, but the search is stored in (and reused from) "searchCache": the table is only + // searched again when "searchCache" does not already hold the neighbours of x, e.g. because it + // was filled by a previous call to Get with these very same coordinates. + template + KOKKOS_INLINE_FUNCTION + void GetNeighbours(const real x[kDim], Tint &dimensions, Tint &offset, Treal &xin, Treal &data, + LookupTableSearchCache &searchCache, + real xN[2*kDim], real dataN[1 << kDim]) const { + // Reuse GetNeighboursIndx for the search and the idx -> flat data offset bookkeeping, so that + // logic only lives in one place (GetDataIndex/GetNeighboursIndx). + int idx[kDim]; + int dataIdx[1 << kDim]; + GetNeighboursIndx(x, dimensions, offset, xin, searchCache, idx, dataIdx); + + if(!searchCache.valid) { + for(int n = 0 ; n < 2*kDim ; n++) xN[n] = NAN; + for(unsigned int n = 0 ; n < (1 << kDim) ; n++) dataN[n] = NAN; + return; + } + + // Coordinates of the neighbours along each dimension + for(int n = 0 ; n < kDim ; n++) { + xN[2*n] = xin(offset(n) + idx[n]); + xN[2*n+1] = xin(offset(n) + idx[n]+1); + if(interpolateInFuncSpace) { + xN[2*n] = invFunc(xN[2*n]); + xN[2*n+1] = invFunc(xN[2*n+1]); + } + } + + // Data on each vertex of the neighbours + for(unsigned int n = 0 ; n < (1 << kDim) ; n++) { + dataN[n] = data(dataIdx[n]); + if(interpolateInFuncSpace) dataN[n] = invFunc(dataN[n]); + } + } + + // Generic getter for the indices of the neighbours used by the interpolation, for all kinds of + // input arrays. On output, idx[n] is the index of the left neighbour along the dimension n + // (the right one being idx[n]+1), and dataIdx[v] is the index in the data array of the vertex v + // of these neighbours (see GetDataIndex for the convention on v). + // All of the outputs are set to -1 when the requested coordinates contain nans. + template + KOKKOS_INLINE_FUNCTION + void GetNeighboursIndx(const real x[kDim], Tint &dimensions, Tint &offset, Treal &xin, + LookupTableSearchCache &searchCache, + int idx[kDim], int dataIdx[1 << kDim]) const { + SearchNeighbours(x, dimensions, offset, xin, searchCache); + + if(!searchCache.valid) { + for(int n = 0 ; n < kDim ; n++) idx[n] = -1; + for(unsigned int n = 0 ; n < (1 << kDim) ; n++) dataIdx[n] = -1; + return; + } + + for(int n = 0 ; n < kDim ; n++) { + idx[n] = searchCache.idx[n]; + } + for(unsigned int n = 0 ; n < (1 << kDim) ; n++) { + dataIdx[n] = GetDataIndex(dimensions, searchCache.idx, n); + } + } + + +// Same as above, but the search is not stored in a caller-supplied structure: a throwaway +// LookupTableSearchCache is created and filled, and then discarded. + template + KOKKOS_INLINE_FUNCTION + void GetNeighboursIndx(const real x[kDim], Tint &dimensions, Tint &offset, Treal &xin, + int idx[kDim], int dataIdx[1 << kDim]) const { + LookupTableSearchCache searchCache; + GetNeighboursIndx(x, dimensions, offset, xin, searchCache, idx, dataIdx); + } + + + + // ----------------------------------------------------------------------------------------- + // Public interface: each of Get/GetNeighbours/GetNeighboursIndx below is a thin one-line + // dispatch to the generic (private) implementation above, selecting device or host storage, + // and optionally reusing a caller-supplied LookupTableSearchCache cache. None of them contains + // any interpolation logic of its own. + // ----------------------------------------------------------------------------------------- + public: // Getter on device KOKKOS_INLINE_FUNCTION real Get(const real x[kDim]) const { @@ -137,14 +404,86 @@ class LookupTable { real GetHost(const real x[kDim]) const { return(Get(x, dimensionsHost, offsetHost, xinHost, dataHost)); } + + // Getter for the neighbours used by the interpolation, on device + KOKKOS_INLINE_FUNCTION + void GetNeighbours(const real x[kDim], real xN[2*kDim], real dataN[1 << kDim]) const { + GetNeighbours(x, dimensionsDev, offsetDev, xinDev, dataDev, xN, dataN); + } + + // Getter for the neighbours used by the interpolation, on Host + KOKKOS_INLINE_FUNCTION + void GetNeighboursHost(const real x[kDim], real xN[2*kDim], real dataN[1 << kDim]) const { + GetNeighbours(x, dimensionsHost, offsetHost, xinHost, dataHost, xN, dataN); + } + + // Getter for the indices of the neighbours used by the interpolation, on device + KOKKOS_INLINE_FUNCTION + void GetNeighboursIndx(const real x[kDim], int idx[kDim], int dataIdx[1 << kDim]) const { + GetNeighboursIndx(x, dimensionsDev, offsetDev, xinDev, idx, dataIdx); + } + + // Getter for the indices of the neighbours used by the interpolation, on Host + KOKKOS_INLINE_FUNCTION + void GetNeighboursIndxHost(const real x[kDim], int idx[kDim], int dataIdx[1 << kDim]) const { + GetNeighboursIndx(x, dimensionsHost, offsetHost, xinHost, idx, dataIdx); + } + + // Getter on device, which stores the neighbours it used in "neighbours". Giving that same + // structure to GetNeighbours or GetNeighboursIndx below then avoids searching the table twice. + KOKKOS_INLINE_FUNCTION + real Get(const real x[kDim], LookupTableSearchCache &neighbours) const { + return(Get(x, dimensionsDev, offsetDev, xinDev, dataDev, neighbours)); + } + + // Getter on Host, which stores the neighbours it used in "neighbours" + KOKKOS_INLINE_FUNCTION + real GetHost(const real x[kDim], LookupTableSearchCache &neighbours) const { + return(Get(x, dimensionsHost, offsetHost, xinHost, dataHost, neighbours)); + } + + // Getter for the neighbours used by the interpolation, on device, reusing the search stored in + // "neighbours" when it was performed for these very same coordinates + KOKKOS_INLINE_FUNCTION + void GetNeighbours(const real x[kDim], LookupTableSearchCache &neighbours, + real xN[2*kDim], real dataN[1 << kDim]) const { + GetNeighbours(x, dimensionsDev, offsetDev, xinDev, dataDev, neighbours, xN, dataN); + } + + // Getter for the neighbours used by the interpolation, on Host, reusing the search stored in + // "neighbours" when it was performed for these very same coordinates + KOKKOS_INLINE_FUNCTION + void GetNeighboursHost(const real x[kDim], LookupTableSearchCache &neighbours, + real xN[2*kDim], real dataN[1 << kDim]) const { + GetNeighbours(x, dimensionsHost, offsetHost, xinHost, dataHost, neighbours, xN, dataN); + } + + // Getter for the indices of the neighbours, on device, reusing the search stored in + // "neighbours" when it was performed for these very same coordinates + KOKKOS_INLINE_FUNCTION + void GetNeighboursIndx(const real x[kDim], LookupTableSearchCache &searchCache, + int idx[kDim], int dataIdx[1 << kDim]) const { + GetNeighboursIndx(x, dimensionsDev, offsetDev, xinDev, searchCache, idx, dataIdx); + } + + // Getter for the indices of the neighbours, on Host, reusing the search stored in + // "neighbours" when it was performed for these very same coordinates + KOKKOS_INLINE_FUNCTION + void GetNeighboursIndxHost(const real x[kDim], LookupTableSearchCache &neighbours, + int idx[kDim], int dataIdx[1 << kDim]) const { + GetNeighboursIndx(x, dimensionsHost, offsetHost, xinHost, neighbours, idx, dataIdx); + } }; -template -LookupTable::LookupTable(std::vector filenames, +template +LookupTable::LookupTable(std::vector filenames, std::string dataSet, - bool errOOB) { + bool errOOB, bool funcSpace, TFunc funcIn, TInvFunc invFuncIn) { idfx::pushRegion("LookupTable::LookupTable"); this->errorIfOutOfBound = errOOB; + this->interpolateInFuncSpace = funcSpace; + this->func = funcIn; + this->invFunc = invFuncIn; std::vector shape; bool fortran_order; @@ -177,7 +516,7 @@ LookupTable::LookupTable(std::vector filenames, } // Allocate the required memory - //Allocate arrays so that the data fits in it + // Allocate arrays so that the data fits in it this->xinDev = IdefixArray1D ("Table_x", sizeTotal); this->dimensionsDev = IdefixArray1D ("Table_dim", kDim); this->offsetDev = IdefixArray1D ("Table_offset", kDim); @@ -233,6 +572,9 @@ LookupTable::LookupTable(std::vector filenames, } } + // Transform the table in function space if required + if(this->interpolateInFuncSpace) this->ToFuncSpace(); + // Copy to target Kokkos::deep_copy(this->xinDev ,xinHost); Kokkos::deep_copy(this->dimensionsDev, dimensionsHost); @@ -244,13 +586,23 @@ LookupTable::LookupTable(std::vector filenames, // Constructor from CSV file -template -LookupTable::LookupTable(std::string filename, char delimiter, bool errOOB) { +// The coordinates and the data are read as lines of the input file, unless readColumns is set +// to true, in which case they are read as columns of the input file (1D tables only, see the +// declaration of the class for a description of both layouts). +template +LookupTable::LookupTable(std::string filename, char delimiter, bool errOOB, + bool readColumns, bool funcSpace, TFunc funcIn, TInvFunc invFuncIn) { idfx::pushRegion("LookupTable::LookupTable"); this->errorIfOutOfBound = errOOB; + this->interpolateInFuncSpace = funcSpace; + this->func = funcIn; + this->invFunc = invFuncIn; if(kDim>2) { IDEFIX_ERROR("CSV files are only compatible with 1D and 2D tables"); } + if(kDim>1 && readColumns) { + IDEFIX_ERROR("CSV files can only be read as columns for 1D tables"); + } // Only 1 process loads the file // Size of the array int size[2]; @@ -265,21 +617,18 @@ LookupTable::LookupTable(std::string filename, char delimiter, bool errOOB if(file.is_open()) { std::string line, lineWithComments; - bool firstLine = true; - int nx = -1; + // Full content of the file, stored as a list of lines, each line being a list of values + std::vector> fileContent; while(std::getline(file, lineWithComments)) { // get rid of comments (starting with #) line = lineWithComments.substr(0, lineWithComments.find("#",0)); if (line.empty()) continue; // skip blank line - char firstChar = line.find_first_not_of(" "); + std::size_t firstChar = line.find_first_not_of(" "); if (firstChar == std::string::npos) continue; // line is all white space // Walk the line - bool firstColumn=true; - if(kDim == 1) firstColumn = false; - - std::vector dataLine; - dataLine.clear(); + std::vector lineVector; + lineVector.clear(); // make the line a string stream, and get all of the values separated by a delimiter std::stringstream str(line); std::string valueString; @@ -294,31 +643,62 @@ LookupTable::LookupTable(std::string filename, char delimiter, bool errOOB << "\" cannot be converted to real." << std::endl; IDEFIX_ERROR(errmsg); } - if(firstLine) { - xVector.push_back(value); - } else if(firstColumn) { - yVector.push_back(value); - firstColumn = false; - } else { - dataLine.push_back(value); - } + lineVector.push_back(value); } // We have finished the line - if(firstLine) { - nx = xVector.size(); - firstLine=false; - } else { - if(dataLine.size() != nx) { + fileContent.push_back(lineVector); + // When read as lines, a 1D table is fully described by the first two lines of the file, + // so we stop reading what's after them + if(kDim < 2 && !readColumns && fileContent.size() == 2) break; + } + file.close(); + // End of file reached + + if(fileContent.size() < 2) { + IDEFIX_ERROR("LookupTable: The input CSV file should contain at least two lines"); + } + + // Dispatch the content of the file in the coordinate and data containers. + // Note that dataVector is always indexed as dataVector[j][i], where i (resp. j) is the + // index along the 1st (resp. 2nd) dimension of the table. + if(readColumns) { + // (1D tables only) the coordinates are stored in the 1st column of the input file, + // while the data is stored in its 2nd column + dataVector.push_back(std::vector()); + for(int i = 0 ; i < fileContent.size() ; i++) { + if(fileContent[i].size() < 2) { + IDEFIX_ERROR("LookupTable: The input CSV file should have at least two columns " + "when the table is read as columns"); + } + xVector.push_back(fileContent[i][0]); + dataVector[0].push_back(fileContent[i][1]); + } + } else { + // (default) the arrays are stored as lines of the input file + // 1st line always gives the coordinates of the 1st dimension + xVector = fileContent[0]; + const int nx = xVector.size(); + if(kDim < 2) { + // 2nd line is the data + if(fileContent[1].size() != nx) { IDEFIX_ERROR("LookupTable: The number of columns in the input CSV " "file should be constant"); } - dataVector.push_back(dataLine); - firstLine = false; - if(kDim < 2) break; // Stop reading what's after the first two lines + dataVector.push_back(fileContent[1]); + } else { + // each of the following lines starts with the coordinate of the 2nd dimension, + // followed by the data of that line + for(int j = 1 ; j < fileContent.size() ; j++) { + if(fileContent[j].size() != nx+1) { + IDEFIX_ERROR("LookupTable: The number of columns in the input CSV " + "file should be constant"); + } + yVector.push_back(fileContent[j][0]); + dataVector.push_back(std::vector(fileContent[j].begin()+1, + fileContent[j].end())); + } } } - file.close(); - // End of file reached } else { std::stringstream errmsg; errmsg << "LookupTable: Unable to open file " << filename << std::endl; @@ -398,6 +778,9 @@ LookupTable::LookupTable(std::string filename, char delimiter, bool errOOB MPI_Bcast(dataHost.data(),dataHost.extent(0), realMPI, 0, MPI_COMM_WORLD); #endif + // Transform the table in function space if required + if(this->interpolateInFuncSpace) this->ToFuncSpace(); + // Copy to target Kokkos::deep_copy(this->xinDev ,xinHost); Kokkos::deep_copy(this->dimensionsDev, dimensionsHost); @@ -426,13 +809,16 @@ LookupTable::LookupTable(std::string filename, char delimiter, bool errOOB // Constructor from IdefixHostArray -template +template template -LookupTable::LookupTable(Kokkos::View array, +LookupTable::LookupTable(Kokkos::View array, std::array,kDim> x, - bool errOOB) { + bool errOOB, bool funcSpace, TFunc funcIn, TInvFunc invFuncIn) { idfx::pushRegion("LookupTable::LookupTable"); this->errorIfOutOfBound = errOOB; + this->interpolateInFuncSpace = funcSpace; + this->func = funcIn; + this->invFunc = invFuncIn; std::vector shape(kDim); for(int i = 0 ; i < kDim ; i++) shape[i] = x[i].extent(0); @@ -504,6 +890,9 @@ LookupTable::LookupTable(Kokkos::View array, } } + // Transform the table in function space if required + if(this->interpolateInFuncSpace) this->ToFuncSpace(); + // Copy to target Kokkos::deep_copy(this->xinDev ,xinHost); Kokkos::deep_copy(this->dimensionsDev, dimensionsHost); diff --git a/test/utils/lookupTable/main.cpp b/test/utils/lookupTable/main.cpp index 1fc660cd6..d13c47e02 100644 --- a/test/utils/lookupTable/main.cpp +++ b/test/utils/lookupTable/main.cpp @@ -10,6 +10,87 @@ // minimal skeleton to use idfx basic functions void testReduction(); +// Custom transformation (and its inverse) to check that the functions used for the interpolation +// in function space can be changed when the lookup table is created +struct MyLog10 { + KOKKOS_INLINE_FUNCTION real operator() (const real x) const { + return(log10(x)); + } +}; + +struct MyPow10 { + KOKKOS_INLINE_FUNCTION real operator() (const real x) const { + return(pow(10.0,x)); + } +}; + +// --------------------------------------------------------------------------------------------- +// Small test helpers, so that each check below states what it verifies and why it failed, instead +// of every failure printing an undifferentiated "ERROR!!" that has to be traced back to the last +// banner that was printed. +// --------------------------------------------------------------------------------------------- +namespace { + +void Banner(const std::string &what) { + idfx::cout << "--------------------------------------" << std::endl; + idfx::cout << what << std::endl; +} + +void Success() { + idfx::cout << "Success" << std::endl; +} + +// Exact equality: used for values that should reproduce bit-for-bit (e.g. requesting a table node +// directly, or comparing an int index), matching what the original tests checked with a plain !=. +template +void CheckEqual(T got, T expected, const std::string &what) { + if(got != expected) { + idfx::cerr << "ERROR!! " << what << ": got " << got << ", expected " << expected << std::endl; + exit(1); + } +} + +// Element-wise exact equality over a small fixed-size array (neighbour coordinates/data/indices). +template +void CheckArrayEqual(const T *got, const T *expected, int n, const std::string &what) { + for(int i = 0 ; i < n ; i++) { + if(got[i] != expected[i]) { + idfx::cerr << "ERROR!! " << what << " (element " << i << "): got " << got[i] + << ", expected " << expected[i] << std::endl; + exit(1); + } + } +} + +// Tolerance-based equality, for values obtained through interpolation. +void CheckClose(real got, real expected, real tol, const std::string &what) { + if(std::fabs(got - expected) > tol) { + idfx::cerr << std::scientific; + idfx::cerr << "ERROR!! " << what << ": got " << got << ", expected " << expected + << " (|diff|=" << std::fabs(got - expected) << ", tol=" << tol << ")" << std::endl; + exit(1); + } +} + +void CheckTrue(bool cond, const std::string &what) { + if(!cond) { + idfx::cerr << "ERROR!! " << what << std::endl; + exit(1); + } +} + +// Expected bracketing neighbours of the 2D CSV table (toto.csv) around (x=2.1, y=3.5): x brackets +// [2,3], y brackets [3,4], and the data on the four surrounding vertices is 5,6,6,7. This is +// checked from several angles across the file below (Get vs GetNeighbours vs GetNeighboursIndx, +// host vs device, cached vs uncached search), so it is defined once here instead of being +// retyped as magic numbers in every block. +const real kXN2D[4] = {2.0, 3.0, 3.0, 4.0}; +const real kDataN2D[4] = {5.0, 6.0, 6.0, 7.0}; +const int kIdx2D[2] = {0, 1}; +const int kDataIdx2D[4] = {1, 4, 2, 5}; +const real kValue2D = 5.6; // csv.Get({2.1, 3.5}) + +} // namespace // main function int main( int argc, char* argv[] ) @@ -35,8 +116,8 @@ int main( int argc, char* argv[] ) { idfx::initialize(); - idfx::cout << "--------------------------------------" << std::endl; - idfx::cout << "Testing 2D CSV file on device." << std::endl; + + Banner("Testing 2D CSV file on device."); IdefixArray1D arr = IdefixArray1D("Test",1); IdefixArray1D::host_mirror_type arrHost = Kokkos::create_mirror_view(arr); @@ -50,33 +131,20 @@ int main( int argc, char* argv[] ) }); Kokkos::deep_copy(arrHost , arr); - idfx::cout << "result="<1e-13) { - idfx::cerr << std::scientific; - idfx::cerr << "ERROR!!" << std::endl; - idfx::cerr << arrHost(0)-5.6; - exit(1); - } - idfx::cout << "Success" << std::endl; + CheckClose(arrHost(0), kValue2D, 1e-13, "2D CSV, device"); + Success(); - idfx::cout << "--------------------------------------" << std::endl; - idfx::cout << "Testing 2D CSV file on Host." << std::endl; + Banner("Testing 2D CSV file on Host."); real x[2]; x[0] = 2.1; x[1] = 3.5; real result = csv.GetHost(x); idfx::cout << "result="<1e-13) { - idfx::cerr << std::scientific; - idfx::cerr << "ERROR!!" << std::endl; - idfx::cerr << result-5.6; - exit(1); - } - idfx::cout << "Success" << std::endl; + CheckClose(result, kValue2D, 1e-13, "2D CSV, host"); + Success(); - idfx::cout << "--------------------------------------" << std::endl; - idfx::cout << "Testing 1D CSV file on device." << std::endl; + Banner("Testing 1D CSV file on device."); // Read 1D CSV File LookupTable<1> csv1D("toto1D.csv",','); @@ -87,30 +155,403 @@ int main( int argc, char* argv[] ) }); Kokkos::deep_copy(arrHost , arr); - idfx::cout << "result="< csv1Dcolumn("toto1Dcolumn.csv",',', true, true); + + idefix_for("loop",0, 1, KOKKOS_LAMBDA (int i) { + real x[1]; + x[0] = 2.1; + arr(i) = csv1Dcolumn.Get(x); + }); + + Kokkos::deep_copy(arrHost , arr); + idfx::cout << "result="< csv1Dlin("toto1Dlog.csv",',', true, true); + LookupTable<1> csv1Dlog("toto1Dlog.csv",',', true, true, true); + // the transformation and its inverse can also be chosen when the table is created + LookupTable<1, MyLog10, MyPow10> csv1Dlog10("toto1Dlog.csv",',', true, true, true); + + idefix_for("loop",0, 1, KOKKOS_LAMBDA (int i) { + real x[1]; + x[0] = 3.0; + arr(i) = csv1Dlog.Get(x); + }); + Kokkos::deep_copy(arrHost , arr); + idfx::cout << "result="< xNdev = IdefixArray1D("xN",4); + IdefixArray1D dataNdev = IdefixArray1D("dataN",4); + IdefixArray1D idxDev = IdefixArray1D("idx",2); + IdefixArray1D dataIdxDev = IdefixArray1D("dataIdx",4); + + idefix_for("neighbours",0, 1, KOKKOS_LAMBDA (int i) { + real xq[2]; + xq[0] = 2.1; + xq[1] = 3.5; + real xN[4]; + real dataN[4]; + int idx[2]; + int dataIdx[4]; + csv.GetNeighbours(xq, xN, dataN); + csv.GetNeighboursIndx(xq, idx, dataIdx); + for(int n = 0 ; n < 4 ; n++) { + xNdev(n) = xN[n]; + dataNdev(n) = dataN[n]; + dataIdxDev(n) = dataIdx[n]; + } + idxDev(0) = idx[0]; + idxDev(1) = idx[1]; + }); + + auto xNHost = Kokkos::create_mirror_view(xNdev); + auto dataNHost = Kokkos::create_mirror_view(dataNdev); + auto idxHost = Kokkos::create_mirror_view(idxDev); + auto dataIdxHost = Kokkos::create_mirror_view(dataIdxDev); + Kokkos::deep_copy(xNHost, xNdev); + Kokkos::deep_copy(dataNHost, dataNdev); + Kokkos::deep_copy(idxHost, idxDev); + Kokkos::deep_copy(dataIdxHost, dataIdxDev); + + idfx::cout << "2D: x=" << xNHost(0) << "," << xNHost(1) << " y=" << xNHost(2) << "," + << xNHost(3) << " data=" << dataNHost(0) << "," << dataNHost(1) << "," + << dataNHost(2) << "," << dataNHost(3) << " idx=" << idxHost(0) << "," + << idxHost(1) << std::endl; + real xN[4], dataN[4]; + int idx[2], dataIdx[4]; + for(int n = 0 ; n < 4 ; n++) { xN[n] = xNHost(n); dataN[n] = dataNHost(n); } + idx[0] = idxHost(0); idx[1] = idxHost(1); + for(int n = 0 ; n < 4 ; n++) dataIdx[n] = dataIdxHost(n); + CheckArrayEqual(xN, kXN2D, 4, "2D neighbours (device), coordinates"); + CheckArrayEqual(dataN, kDataN2D, 4, "2D neighbours (device), data"); + CheckArrayEqual(idx, kIdx2D, 2, "2D neighbours (device), idx"); + CheckArrayEqual(dataIdx, kDataIdx2D, 4, "2D neighbours (device), dataIdx"); + Success(); + } + + Banner("Testing the 1D neighbours on the edges of the table."); + { + // toto1D.csv holds x=1,2,3 and data=2,4,6. Whatever the requested value, we expect the two + // neighbours bracketing it, i.e. the last two nodes when we sit on the upper edge + real xq[1]; + real xN[2]; + real dataN[2]; + real expected[3][2] = {{1.0,2.0}, {2.0,3.0}, {2.0,3.0}}; + real xRequest[3] = {1.0, 2.0, 3.0}; + for(int n = 0 ; n < 3 ; n++) { + xq[0] = xRequest[n]; + csv1D.GetNeighboursHost(xq, xN, dataN); + idfx::cout << "x=" << xq[0] << " -> neighbours " << xN[0] << "," << xN[1] + << " (data " << dataN[0] << "," << dataN[1] << ")" << std::endl; + CheckArrayEqual(xN, expected[n], 2, "1D edge neighbours, coordinates"); + CheckTrue(xN[0] <= xq[0] && xN[1] >= xq[0], + "1D edge neighbours do not bracket the requested value"); + } + Success(); + } + + Banner("Testing the reuse of the search between Get and GetNeighbours on Host."); + { + real xq[2]; + xq[0] = 2.1; + xq[1] = 3.5; + // The neighbours found by Get are stored in nb... + LookupTableSearchCache<2> nb; + result = csv.GetHost(xq, nb); + CheckClose(result, kValue2D, 1e-13, "Get-then-GetNeighbours (host), value"); + CheckTrue(nb.valid, "Get-then-GetNeighbours (host), search validity"); + CheckArrayEqual(nb.idx, kIdx2D, 2, "Get-then-GetNeighbours (host), cached idx"); + + // ... and are reused (not computed again) by the getters below + real xN[4]; + real dataN[4]; + int idx[2]; + int dataIdx[4]; + csv.GetNeighboursHost(xq, nb, xN, dataN); + csv.GetNeighboursIndxHost(xq, nb, idx, dataIdx); + CheckArrayEqual(xN, kXN2D, 4, "Get-then-GetNeighbours (host), coordinates"); + CheckArrayEqual(dataN, kDataN2D, 4, "Get-then-GetNeighbours (host), data"); + CheckArrayEqual(idx, kIdx2D, 2, "Get-then-GetNeighbours (host), idx"); + CheckArrayEqual(dataIdx, kDataIdx2D, 4, "Get-then-GetNeighbours (host), dataIdx"); + + // Check that the table is really not searched again for the same coordinates: we corrupt + // the stored search, and check that the corrupted result is the one which is used + nb.idx[0] = 1; + csv.GetNeighboursIndxHost(xq, nb, idx, dataIdx); + CheckEqual(idx[0], 1, "Get-then-GetNeighbours (host): the stored search was not reused"); + + // ... while different coordinates trigger a new search, as usual + real xq2[2]; + xq2[0] = 2.9; + xq2[1] = 2.5; + csv.GetNeighboursIndxHost(xq2, nb, idx, dataIdx); + CheckTrue(idx[0] == 0 && idx[1] == 0, + "Get-then-GetNeighbours (host): different coordinates did not trigger a new search"); + Success(); + } + + Banner("Testing the search performed by GetNeighbours first, and reused by Get, on Host."); + { + real xq[2]; + xq[0] = 2.1; + xq[1] = 3.5; + + // No call to Get yet: GetNeighbours searches the table as usual, and stores the result + LookupTableSearchCache<2> nb; + real xN[4]; + real dataN[4]; + csv.GetNeighboursHost(xq, nb, xN, dataN); + CheckTrue(nb.valid, "GetNeighbours-then-Get (host), search validity"); + CheckArrayEqual(nb.idx, kIdx2D, 2, "GetNeighbours-then-Get (host), cached idx"); + CheckArrayEqual(xN, kXN2D, 4, "GetNeighbours-then-Get (host), coordinates"); + CheckArrayEqual(dataN, kDataN2D, 4, "GetNeighbours-then-Get (host), data"); + + // Get now reuses that search for the same coordinates + result = csv.GetHost(xq, nb); + CheckClose(result, kValue2D, 1e-13, "GetNeighbours-then-Get (host), value"); + + // Same check as before, the other way around: we corrupt the ratio stored by + // GetNeighbours, and check that Get uses it instead of computing it again. + // With delta[0]=0.5 instead of 0.1, the interpolation gives (5+6+6+7)/4=6 + nb.delta[0] = 0.5; + result = csv.GetHost(xq, nb); + CheckClose(result, 6.0, 1e-13, + "GetNeighbours-then-Get (host): the stored search was not reused by Get"); + + // ... while different coordinates make Get search the table again + real xq2[2]; + xq2[0] = 2.9; + xq2[1] = 2.5; + result = csv.GetHost(xq2, nb); + CheckClose(result, csv.GetHost(xq2), 1e-13, + "GetNeighbours-then-Get (host): different coordinates did not trigger a new search"); + CheckEqual(nb.idx[1], 0, + "GetNeighbours-then-Get (host): different coordinates did not trigger a new search"); + + // The same holds when the search comes from GetNeighboursIndx, here on the 1D table + LookupTableSearchCache<1> nb1D; + real xq1[1]; + xq1[0] = 2.1; + int idx1[1]; + int dataIdx1[2]; + csv1D.GetNeighboursIndxHost(xq1, nb1D, idx1, dataIdx1); + CheckTrue(nb1D.valid, "GetNeighboursIndx-then-Get (host, 1D), search validity"); + CheckEqual(idx1[0], 1, "GetNeighboursIndx-then-Get (host, 1D), idx"); + nb1D.delta[0] = 0.0; // x is now on the left neighbour, so we expect its data + result = csv1D.GetHost(xq1, nb1D); + CheckClose(result, 4.0, 1e-13, + "GetNeighboursIndx-then-Get (host, 1D): the stored search was not reused by Get"); + Success(); } - idfx::cout << "Success" << std::endl; - idfx::cout << "--------------------------------------" << std::endl; - idfx::cout << "Testing 3D npy file on device." << std::endl; + Banner("Testing the reuse of the search between Get and GetNeighbours on device."); + { + IdefixArray1D xNdev = IdefixArray1D("xN",4); + IdefixArray1D dataNdev = IdefixArray1D("dataN",4); + IdefixArray1D idxDev = IdefixArray1D("idx",2); + IdefixArray1D dataIdxDev = IdefixArray1D("dataIdx",4); + + idefix_for("neighbours",0, 1, KOKKOS_LAMBDA (int i) { + real xq[2]; + xq[0] = 2.1; + xq[1] = 3.5; + // the structure is local to the loop, and is therefore private to each thread + LookupTableSearchCache<2> nb; + real xN[4]; + real dataN[4]; + int idx[2]; + int dataIdx[4]; + arr(i) = csv.Get(xq, nb); + csv.GetNeighbours(xq, nb, xN, dataN); + csv.GetNeighboursIndx(xq, nb, idx, dataIdx); + for(int n = 0 ; n < 4 ; n++) { + xNdev(n) = xN[n]; + dataNdev(n) = dataN[n]; + dataIdxDev(n) = dataIdx[n]; + } + idxDev(0) = idx[0]; + idxDev(1) = idx[1]; + }); + + Kokkos::deep_copy(arrHost , arr); + auto xNHost = Kokkos::create_mirror_view(xNdev); + auto dataNHost = Kokkos::create_mirror_view(dataNdev); + auto idxHost = Kokkos::create_mirror_view(idxDev); + auto dataIdxHost = Kokkos::create_mirror_view(dataIdxDev); + Kokkos::deep_copy(xNHost, xNdev); + Kokkos::deep_copy(dataNHost, dataNdev); + Kokkos::deep_copy(idxHost, idxDev); + Kokkos::deep_copy(dataIdxHost, dataIdxDev); + + // Same test on device, but with GetNeighbours called first and Get reusing its search. + // The ratio stored by GetNeighbours is corrupted in between, so that a value of 6 (instead + // of 5.6) proves that Get did not search the table again + IdefixArray1D valuesDev = IdefixArray1D("values",2); + idefix_for("neighboursFirst",0, 1, KOKKOS_LAMBDA (int i) { + real xq[2]; + xq[0] = 2.1; + xq[1] = 3.5; + LookupTableSearchCache<2> nb; + real xN[4]; + real dataN[4]; + // no call to Get yet: the search is performed here for the first time + csv.GetNeighbours(xq, nb, xN, dataN); + valuesDev(0) = csv.Get(xq, nb); + nb.delta[0] = 0.5; + valuesDev(1) = csv.Get(xq, nb); + }); + auto valuesHost = Kokkos::create_mirror_view(valuesDev); + Kokkos::deep_copy(valuesHost, valuesDev); + idfx::cout << "neighbours first: value=" << valuesHost(0) + << " (with a corrupted stored search)=" << valuesHost(1) << std::endl; + CheckClose(valuesHost(0), kValue2D, 1e-13, + "GetNeighbours-then-Get (device): the stored search was not reused by Get"); + CheckClose(valuesHost(1), 6.0, 1e-13, + "GetNeighbours-then-Get (device): the stored search was not reused by Get"); + + idfx::cout << "value=" << arrHost(0) << " x=" << xNHost(0) << "," << xNHost(1) + << " y=" << xNHost(2) << "," << xNHost(3) + << " data=" << dataNHost(0) << "," << dataNHost(1) << "," << dataNHost(2) + << "," << dataNHost(3) << " idx=" << idxHost(0) << "," << idxHost(1) << std::endl; + CheckClose(arrHost(0), kValue2D, 1e-13, "Get-then-GetNeighbours (device), value"); + real xN[4], dataN[4]; + int idx[2], dataIdx[4]; + for(int n = 0 ; n < 4 ; n++) { xN[n] = xNHost(n); dataN[n] = dataNHost(n); } + idx[0] = idxHost(0); idx[1] = idxHost(1); + for(int n = 0 ; n < 4 ; n++) dataIdx[n] = dataIdxHost(n); + CheckArrayEqual(xN, kXN2D, 4, "Get-then-GetNeighbours (device), coordinates"); + CheckArrayEqual(dataN, kDataN2D, 4, "Get-then-GetNeighbours (device), data"); + CheckArrayEqual(idx, kIdx2D, 2, "Get-then-GetNeighbours (device), idx"); + CheckArrayEqual(dataIdx, kDataIdx2D, 4, "Get-then-GetNeighbours (device), dataIdx"); + Success(); + } + + Banner("Testing 3D npy file on device."); // Read npy File std::vector coords({"x.npy","y.npy","z.npy"}); @@ -125,36 +566,21 @@ int main( int argc, char* argv[] ) }); Kokkos::deep_copy(arrHost , arr); - idfx::cout << "result="<1e-13) { - idfx::cerr << std::scientific; - idfx::cerr << "ERROR!!" << std::endl; - idfx::cerr << arrHost(0)-13.6; - exit(1); - } - idfx::cout << "Success" << std::endl; - - idfx::cout << "--------------------------------------" << std::endl; - idfx::cout << "Testing 3D npy file on host." << std::endl; + CheckClose(arrHost(0), 13.6, 1e-13, "3D npy, device"); + Success(); + Banner("Testing 3D npy file on host."); real y[3]; y[0] = 2.7; y[1] = 7.4; y[2] = 3.9; result = csvnpy.GetHost(y); - idfx::cout << "result="<< result << std::endl; - if(std::fabs(result - 13.6)>1e-13) { - idfx::cerr << std::scientific; - idfx::cerr << "ERROR!!" << std::endl; - idfx::cerr << result-13.6; - exit(1); - } - idfx::cout << "Success" << std::endl; - idfx::cout << "--------------------------------------" << std::endl; - idfx::cout << "Done." << std::endl; + CheckClose(result, 13.6, 1e-13, "3D npy, host"); + Success(); + Banner("Done."); } Kokkos::finalize(); #ifdef WITH_MPI diff --git a/test/utils/lookupTable/testmelib.py b/test/utils/lookupTable/testmelib.py index b41266147..5a5b942af 100644 --- a/test/utils/lookupTable/testmelib.py +++ b/test/utils/lookupTable/testmelib.py @@ -1,11 +1,15 @@ import numpy as np +from scipy.interpolate import RegularGridInterpolator def MakeNumpyFile(): - x = np.arange(1, 10, 1.0) - y = np.arange(5, 10, 1.0) - z = np.arange(2, 5, 1.0) + x = 2 ** np.arange(np.log2(1), np.log2(10), 0.2) + y = 2 ** np.arange(np.log2(5), np.log2(10), 0.2) + z = 2 ** np.arange(np.log2(2), np.log2(5), 0.2) + print(x) + print(y) + print(z) xp, yp, zp = np.meshgrid(x, y, z, indexing="ij") data = xp + 2 * yp - zp @@ -14,6 +18,10 @@ def MakeNumpyFile(): np.save("y.npy", y) np.save("z.npy", z) np.save("data.npy", data) - # show the expected result - # f=RegularGridInterpolator((x, y, z), data) - # print(f([2.7,7.4,3.9])) + f = RegularGridInterpolator((x, y, z), data) + return f([2.7, 7.4, 3.9]) + + +if __name__ == "__main__": + f = MakeNumpyFile() + print("expected result:%f", f) diff --git a/test/utils/lookupTable/toto1Dcolumn.csv b/test/utils/lookupTable/toto1Dcolumn.csv new file mode 100644 index 000000000..291aa31cb --- /dev/null +++ b/test/utils/lookupTable/toto1Dcolumn.csv @@ -0,0 +1,5 @@ +# test csv file, same 1D table as toto1D.csv, but stored as columns + +1.0, 2.0 +2.0, 4.0 +3.0, 6.0 diff --git a/test/utils/lookupTable/toto1Dlog.csv b/test/utils/lookupTable/toto1Dlog.csv new file mode 100644 index 000000000..2cb5d5d51 --- /dev/null +++ b/test/utils/lookupTable/toto1Dlog.csv @@ -0,0 +1,6 @@ +# test csv file, power law data = x^2 stored as columns +# used to test the interpolation in function space + +1.0, 1.0 +2.0, 4.0 +4.0, 16.0