diff --git a/doc/source/programmingguide.rst b/doc/source/programmingguide.rst index 6f65742bc..9d05e7c69 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 ``LookupTableNeighbours`` 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... + LookupTableNeighbours<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 ``LookupTableNeighbours`` can be reused from one point to the next. + +.. warning:: + The structure is declared by the caller, 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..3892fc319 100644 --- a/src/utils/lookupTable.hpp +++ b/src/utils/lookupTable.hpp @@ -14,18 +14,78 @@ #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 LookupTableNeighbours { + 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 + 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,6 +100,119 @@ class LookupTable { bool errorIfOutOfBound{true}; + // 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)"); + } + } + } + + // Generic search of the neighbours used for 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. + template + KOKKOS_INLINE_FUNCTION + 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); + // 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(false); + + // Compute index of closest element assuming even distribution + int i; + + // Check that we're within bounds + if(x_n < xstart) { + if(errorIfOutOfBound) { + Kokkos::abort("LookupTable:: ERROR! Attempt to interpolate below your lower bound."); + } else { + x_n = xstart; + i = 0; + } + } else if( x_n > xend) { + if(errorIfOutOfBound) { + Kokkos::abort("LookupTable:: ERROR! Attempt to interpolate above your upper bound."); + } else { + // We set x_n=xend, and we do the interpolation between xin(dim-2) and xin(dim-1), + // so i= dim-2 + i = dimensions(n)-2; + x_n = xend; + } + } else { + // Bounds are fine, + i = static_cast ( (x_n - xstart) / (xend - xstart) * (dimensions(n)-1)); + if(i == dimensions(n)-1) i = dimensions(n)-2; + // 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++; + } + i = i-1; // i is overestimated by one + } + } + + // Store the index + idx[n] = i; + + // Store the elementary ratio + delta[n] = (x_n - xin(offset(n) + i) ) / (xin(offset(n) + i+1) - xin(offset(n) + i)); + } + + 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); + } + // Generic getter for all kinds of input arrays template KOKKOS_INLINE_FUNCTION @@ -51,7 +224,9 @@ class LookupTable { 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); @@ -99,7 +274,7 @@ 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. + // Do a linear interpolation from the neightbouring points to get our value. real value = 0; // loop on all of the vertices of the neighbours @@ -123,9 +298,148 @@ class LookupTable { value = value + weight*data(index); } + // The interpolation was performed on func(data), so we transform the result back + if(interpolateInFuncSpace) value = invFunc(value); + return(value); } + // Fill "neighbours" 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, + LookupTableNeighbours &neighbours) const { + // Nothing to do if the neighbours of these very coordinates are already known + if(neighbours.Matches(x)) return; + + neighbours.valid = GetIndices(x, dimensions, offset, xin, neighbours.idx, neighbours.delta); + for(int n = 0 ; n < kDim ; n++) { + neighbours.x[n] = x[n]; + } + } + + // Generic getter which stores in "neighbours" 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, + LookupTableNeighbours &neighbours) const { + SearchNeighbours(x, dimensions, offset, xin, neighbours); + + if(!neighbours.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++) { + real weight = 1.0; + for(unsigned int m = 0 ; m < kDim ; 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*neighbours.delta[m]; + } else { + // We're on the left + weight = weight*(1-neighbours.delta[m]); + } + } + value = value + weight*data(GetDataIndex(dimensions, neighbours.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 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 { + LookupTableNeighbours neighbours; + GetNeighbours(x, dimensions, offset, xin, data, neighbours, xN, dataN); + } + + // Same as above, but the search is stored in (and reused from) "neighbours": the table is only + // searched again when "neighbours" 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, + LookupTableNeighbours &neighbours, + real xN[2*kDim], real dataN[1 << kDim]) const { + SearchNeighbours(x, dimensions, offset, xin, neighbours); + + if(!neighbours.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) + neighbours.idx[n]); + xN[2*n+1] = xin(offset(n) + neighbours.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(GetDataIndex(dimensions, neighbours.idx, 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, + int idx[kDim], int dataIdx[1 << kDim]) const { + LookupTableNeighbours neighbours; + GetNeighboursIndx(x, dimensions, offset, xin, neighbours, idx, dataIdx); + } + + // Same as above, but the search is stored in (and reused from) "neighbours", exactly like the + // GetNeighbours variant above + template + KOKKOS_INLINE_FUNCTION + void GetNeighboursIndx(const real x[kDim], Tint &dimensions, Tint &offset, Treal &xin, + LookupTableNeighbours &neighbours, + int idx[kDim], int dataIdx[1 << kDim]) const { + SearchNeighbours(x, dimensions, offset, xin, neighbours); + + if(!neighbours.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] = neighbours.idx[n]; + } + for(unsigned int n = 0 ; n < (1 << kDim) ; n++) { + dataIdx[n] = GetDataIndex(dimensions, neighbours.idx, n); + } + } + // Getter on device KOKKOS_INLINE_FUNCTION real Get(const real x[kDim]) const { @@ -137,14 +451,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], LookupTableNeighbours &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], LookupTableNeighbours &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], LookupTableNeighbours &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], LookupTableNeighbours &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], LookupTableNeighbours &neighbours, + int idx[kDim], int dataIdx[1 << kDim]) const { + GetNeighboursIndx(x, dimensionsDev, offsetDev, xinDev, neighbours, 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], LookupTableNeighbours &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 +563,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 +619,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 +633,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 +664,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 +690,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 +825,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 +856,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 +937,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..fdec245f4 100644 --- a/test/utils/lookupTable/main.cpp +++ b/test/utils/lookupTable/main.cpp @@ -10,6 +10,20 @@ // 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)); + } +}; + // main function int main( int argc, char* argv[] ) @@ -109,6 +123,497 @@ int main( int argc, char* argv[] ) } idfx::cout << "Success" << std::endl; + idfx::cout << "--------------------------------------" << std::endl; + idfx::cout << "Testing 1D CSV file read as columns on device." << std::endl; + // Read the same 1D table, but stored as columns of the CSV file + LookupTable<1> 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="<1e-13) { + idfx::cerr << std::scientific; + idfx::cerr << "ERROR!!" << std::endl; + idfx::cerr << arrHost(0)-9.0; + exit(1); + } + idfx::cout << "Success" << std::endl; + + idefix_for("loop",0, 1, KOKKOS_LAMBDA (int i) { + real x[1]; + x[0] = 3.0; + arr(i) = csv1Dlog10.Get(x); + }); + + Kokkos::deep_copy(arrHost , arr); + + idfx::cout << "result (with log10)="<1e-13) { + idfx::cerr << std::scientific; + idfx::cerr << "ERROR!!" << std::endl; + idfx::cerr << arrHost(0)-9.0; + exit(1); + } + idfx::cout << "Success" << std::endl; + + idfx::cout << "--------------------------------------" << std::endl; + idfx::cout << "Testing interpolation in function space on Host." << std::endl; + + x[0] = 3.0; + result = csv1Dlog.GetHost(x); + idfx::cout << "result="<1e-13) { + idfx::cerr << std::scientific; + idfx::cerr << "ERROR!!" << std::endl; + idfx::cerr << result-9.0; + exit(1); + } + result = csv1Dlog10.GetHost(x); + idfx::cout << "result (with log10)="<1e-13) { + idfx::cerr << std::scientific; + idfx::cerr << "ERROR!!" << std::endl; + idfx::cerr << result-9.0; + exit(1); + } + // the same table, interpolated linearly (default behaviour) + result = csv1Dlin.GetHost(x); + idfx::cout << "result (linear interpolation)="<1e-13) { + idfx::cerr << std::scientific; + idfx::cerr << "ERROR!!" << std::endl; + idfx::cerr << result-10.0; + exit(1); + } + idfx::cout << "Success" << std::endl; + + idfx::cout << "--------------------------------------" << std::endl; + idfx::cout << "Testing the neighbours used for the interpolation on Host." << std::endl; + { + // 1D table (x=1,2,3 and data=2,4,6) interpolated in x=2.1 + real xq[1]; + xq[0] = 2.1; + real xN[2]; + real dataN[2]; + int idx[1]; + int dataIdx[2]; + csv1D.GetNeighboursHost(xq, xN, dataN); + csv1D.GetNeighboursIndxHost(xq, idx, dataIdx); + idfx::cout << "1D: x=" << xN[0] << "," << xN[1] << " data=" << dataN[0] << "," << dataN[1] + << " idx=" << idx[0] << " dataIdx=" << dataIdx[0] << "," << dataIdx[1] + << std::endl; + if(xN[0] != 2.0 || xN[1] != 3.0 || dataN[0] != 4.0 || dataN[1] != 6.0 + || idx[0] != 1 || dataIdx[0] != 1 || dataIdx[1] != 2) { + idfx::cerr << "ERROR!!" << std::endl; + exit(1); + } + + // 2D table interpolated in (2.1,3.5). The vertices are ordered so that the bit n of the + // vertex index tells whether we are on the right (1) or on the left (0) of the dimension n + real xq2[2]; + xq2[0] = 2.1; + xq2[1] = 3.5; + real xN2[4]; + real dataN2[4]; + int idx2[2]; + int dataIdx2[4]; + csv.GetNeighboursHost(xq2, xN2, dataN2); + csv.GetNeighboursIndxHost(xq2, idx2, dataIdx2); + idfx::cout << "2D: x=" << xN2[0] << "," << xN2[1] << " y=" << xN2[2] << "," << xN2[3] + << " data=" << dataN2[0] << "," << dataN2[1] << "," << dataN2[2] << "," + << dataN2[3] << " idx=" << idx2[0] << "," << idx2[1] << std::endl; + if(xN2[0] != 2.0 || xN2[1] != 3.0 || xN2[2] != 3.0 || xN2[3] != 4.0) { + idfx::cerr << "ERROR!!" << std::endl; + exit(1); + } + if(dataN2[0] != 5.0 || dataN2[1] != 6.0 || dataN2[2] != 6.0 || dataN2[3] != 7.0) { + idfx::cerr << "ERROR!!" << std::endl; + exit(1); + } + if(idx2[0] != 0 || idx2[1] != 1 + || dataIdx2[0] != 1 || dataIdx2[1] != 4 || dataIdx2[2] != 2 || dataIdx2[3] != 5) { + idfx::cerr << "ERROR!!" << std::endl; + exit(1); + } + + // The neighbours should reproduce the value returned by GetHost + real dx = (xq2[0]-xN2[0])/(xN2[1]-xN2[0]); + real dy = (xq2[1]-xN2[2])/(xN2[3]-xN2[2]); + real interpolated = (1-dx)*(1-dy)*dataN2[0] + dx*(1-dy)*dataN2[1] + + (1-dx)*dy*dataN2[2] + dx*dy*dataN2[3]; + idfx::cout << "2D: interpolation from the neighbours=" << interpolated << std::endl; + if(std::fabs(interpolated - csv.GetHost(xq2))>1e-13) { + idfx::cerr << "ERROR!!" << std::endl; + exit(1); + } + + // When the table is interpolated in function space, the neighbours are returned in the + // original space of the table (x=1,2,4 and data=1,4,16 here) + real xq3[1]; + xq3[0] = 3.0; + real xN3[2]; + real dataN3[2]; + int idx3[1]; + int dataIdx3[2]; + csv1Dlog.GetNeighboursHost(xq3, xN3, dataN3); + csv1Dlog.GetNeighboursIndxHost(xq3, idx3, dataIdx3); + idfx::cout << "1D (function space): x=" << xN3[0] << "," << xN3[1] + << " data=" << dataN3[0] << "," << dataN3[1] << " idx=" << idx3[0] << std::endl; + if(std::fabs(xN3[0]-2.0)>1e-13 || std::fabs(xN3[1]-4.0)>1e-13 + || std::fabs(dataN3[0]-4.0)>1e-13 || std::fabs(dataN3[1]-16.0)>1e-13 + || idx3[0] != 1 || dataIdx3[0] != 1 || dataIdx3[1] != 2) { + idfx::cerr << "ERROR!!" << std::endl; + exit(1); + } + idfx::cout << "Success" << std::endl; + } + + idfx::cout << "--------------------------------------" << std::endl; + idfx::cout << "Testing the neighbours used for the interpolation on device." << std::endl; + { + 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; + 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; + if(xNHost(0) != 2.0 || xNHost(1) != 3.0 || xNHost(2) != 3.0 || xNHost(3) != 4.0) { + idfx::cerr << "ERROR!!" << std::endl; + exit(1); + } + if(dataNHost(0) != 5.0 || dataNHost(1) != 6.0 || dataNHost(2) != 6.0 + || dataNHost(3) != 7.0) { + idfx::cerr << "ERROR!!" << std::endl; + exit(1); + } + if(idxHost(0) != 0 || idxHost(1) != 1 || dataIdxHost(0) != 1 || dataIdxHost(1) != 4 + || dataIdxHost(2) != 2 || dataIdxHost(3) != 5) { + idfx::cerr << "ERROR!!" << std::endl; + exit(1); + } + idfx::cout << "Success" << std::endl; + } + + idfx::cout << "--------------------------------------" << std::endl; + idfx::cout << "Testing the 1D neighbours on the edges of the table." << std::endl; + { + // 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; + if(xN[0] != expected[n][0] || xN[1] != expected[n][1]) { + idfx::cerr << "ERROR!!" << std::endl; + exit(1); + } + if(xN[0] > xq[0] || xN[1] < xq[0]) { + idfx::cerr << "ERROR!! the neighbours do not bracket the requested value" << std::endl; + exit(1); + } + } + idfx::cout << "Success" << std::endl; + } + + idfx::cout << "--------------------------------------" << std::endl; + idfx::cout << "Testing the reuse of the search between Get and GetNeighbours on Host." + << std::endl; + { + real xq[2]; + xq[0] = 2.1; + xq[1] = 3.5; + // The neighbours found by Get are stored in nb... + LookupTableNeighbours<2> nb; + result = csv.GetHost(xq, nb); + if(std::fabs(result - 5.6)>1e-13 || !nb.valid || nb.idx[0] != 0 || nb.idx[1] != 1) { + idfx::cerr << "ERROR!!" << std::endl; + exit(1); + } + + // ... 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); + if(xN[0] != 2.0 || xN[1] != 3.0 || xN[2] != 3.0 || xN[3] != 4.0) { + idfx::cerr << "ERROR!!" << std::endl; + exit(1); + } + if(dataN[0] != 5.0 || dataN[1] != 6.0 || dataN[2] != 6.0 || dataN[3] != 7.0) { + idfx::cerr << "ERROR!!" << std::endl; + exit(1); + } + if(idx[0] != 0 || idx[1] != 1 + || dataIdx[0] != 1 || dataIdx[1] != 4 || dataIdx[2] != 2 || dataIdx[3] != 5) { + idfx::cerr << "ERROR!!" << std::endl; + exit(1); + } + + // 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); + if(idx[0] != 1) { + idfx::cerr << "ERROR!! the stored search was not reused" << std::endl; + exit(1); + } + + // ... 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); + if(idx[0] != 0 || idx[1] != 0) { + idfx::cerr << "ERROR!! the search was not performed again" << std::endl; + exit(1); + } + idfx::cout << "Success" << std::endl; + } + + idfx::cout << "--------------------------------------" << std::endl; + idfx::cout << "Testing the search performed by GetNeighbours first, and reused by Get, " + "on Host." << std::endl; + { + 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 + LookupTableNeighbours<2> nb; + real xN[4]; + real dataN[4]; + csv.GetNeighboursHost(xq, nb, xN, dataN); + if(!nb.valid || nb.idx[0] != 0 || nb.idx[1] != 1) { + idfx::cerr << "ERROR!! the search was not performed" << std::endl; + exit(1); + } + if(xN[0] != 2.0 || xN[1] != 3.0 || dataN[0] != 5.0 || dataN[3] != 7.0) { + idfx::cerr << "ERROR!!" << std::endl; + exit(1); + } + + // Get now reuses that search for the same coordinates + result = csv.GetHost(xq, nb); + if(std::fabs(result - 5.6)>1e-13) { + idfx::cerr << "ERROR!!" << std::endl; + exit(1); + } + + // 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); + if(std::fabs(result - 6.0)>1e-13) { + idfx::cerr << "ERROR!! the stored search was not reused by Get" << std::endl; + exit(1); + } + + // ... 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); + if(std::fabs(result - csv.GetHost(xq2))>1e-13 || nb.idx[1] != 0) { + idfx::cerr << "ERROR!! the search was not performed again" << std::endl; + exit(1); + } + + // The same holds when the search comes from GetNeighboursIndx, here on the 1D table + LookupTableNeighbours<1> nb1D; + real xq1[1]; + xq1[0] = 2.1; + int idx1[1]; + int dataIdx1[2]; + csv1D.GetNeighboursIndxHost(xq1, nb1D, idx1, dataIdx1); + if(!nb1D.valid || idx1[0] != 1) { + idfx::cerr << "ERROR!!" << std::endl; + exit(1); + } + nb1D.delta[0] = 0.0; // x is now on the left neighbour, so we expect its data + result = csv1D.GetHost(xq1, nb1D); + if(std::fabs(result - 4.0)>1e-13) { + idfx::cerr << "ERROR!! the stored search was not reused by Get" << std::endl; + exit(1); + } + idfx::cout << "Success" << std::endl; + } + + idfx::cout << "--------------------------------------" << std::endl; + idfx::cout << "Testing the reuse of the search between Get and GetNeighbours on device." + << std::endl; + { + 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 + LookupTableNeighbours<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; + LookupTableNeighbours<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; + if(std::fabs(valuesHost(0) - 5.6)>1e-13 || std::fabs(valuesHost(1) - 6.0)>1e-13) { + idfx::cerr << "ERROR!! the stored search was not reused by Get" << std::endl; + exit(1); + } + + 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; + if(std::fabs(arrHost(0) - 5.6)>1e-13) { + idfx::cerr << "ERROR!!" << std::endl; + exit(1); + } + if(xNHost(0) != 2.0 || xNHost(1) != 3.0 || xNHost(2) != 3.0 || xNHost(3) != 4.0) { + idfx::cerr << "ERROR!!" << std::endl; + exit(1); + } + if(dataNHost(0) != 5.0 || dataNHost(1) != 6.0 || dataNHost(2) != 6.0 + || dataNHost(3) != 7.0) { + idfx::cerr << "ERROR!!" << std::endl; + exit(1); + } + if(idxHost(0) != 0 || idxHost(1) != 1 || dataIdxHost(0) != 1 || dataIdxHost(1) != 4 + || dataIdxHost(2) != 2 || dataIdxHost(3) != 5) { + idfx::cerr << "ERROR!!" << std::endl; + exit(1); + } + idfx::cout << "Success" << std::endl; + } + idfx::cout << "--------------------------------------" << std::endl; idfx::cout << "Testing 3D npy file on device." << std::endl; // Read npy File 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