From 41034352afb40bf80470b9231484f38e2a54e023 Mon Sep 17 00:00:00 2001 From: dutta-alankar Date: Wed, 19 Aug 2026 11:58:51 +0200 Subject: [PATCH 1/5] Allow 1D lookup tables to be read as columns of a CSV file The CSV constructor of LookupTable now accepts an optional readColumns argument. When it is set to true (1D tables only), xinHost is read from the first column of the ASCII file and dataHost from its second column, instead of the first and second lines of the file. Both are then shared on the device as before. The default behaviour (arrays read as lines) and the parsing of 2D tables are left unchanged. Also document the handling of comments and blank lines in CSV files, and add a test of the new column layout in test/utils/lookupTable. --- doc/source/programmingguide.rst | 31 ++++++- src/utils/lookupTable.hpp | 102 +++++++++++++++++------- test/utils/lookupTable/main.cpp | 34 ++++++++ test/utils/lookupTable/toto1Dcolumn.csv | 5 ++ 4 files changed, 142 insertions(+), 30 deletions(-) create mode 100644 test/utils/lookupTable/toto1Dcolumn.csv diff --git a/doc/source/programmingguide.rst b/doc/source/programmingguide.rst index 6f65742bc..ad2a70c4f 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 diff --git a/src/utils/lookupTable.hpp b/src/utils/lookupTable.hpp index 72d13a838..63566ae5a 100644 --- a/src/utils/lookupTable.hpp +++ b/src/utils/lookupTable.hpp @@ -18,7 +18,16 @@ 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) + LookupTable(std::string filename, char delimiter, bool errorIfOutOfBound = true, + bool readColumns = false); LookupTable(std::vector filenames, std::string dataSet, bool errorIfOutOfBound = true); @@ -177,7 +186,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); @@ -244,13 +253,20 @@ LookupTable::LookupTable(std::vector filenames, // Constructor from CSV file +// 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) { +LookupTable::LookupTable(std::string filename, char delimiter, bool errOOB, + bool readColumns) { idfx::pushRegion("LookupTable::LookupTable"); this->errorIfOutOfBound = errOOB; 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 +281,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 +307,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; diff --git a/test/utils/lookupTable/main.cpp b/test/utils/lookupTable/main.cpp index 1fc660cd6..b07672236 100644 --- a/test/utils/lookupTable/main.cpp +++ b/test/utils/lookupTable/main.cpp @@ -109,6 +109,40 @@ 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="< Date: Wed, 19 Aug 2026 13:43:04 +0200 Subject: [PATCH 2/5] LookupTable interpolation on function space --- doc/source/programmingguide.rst | 54 ++++++++++++++ src/real_types.hpp | 4 ++ src/utils/lookupTable.hpp | 104 +++++++++++++++++++++++---- test/utils/lookupTable/main.cpp | 88 +++++++++++++++++++++++ test/utils/lookupTable/toto1Dlog.csv | 6 ++ 5 files changed, 241 insertions(+), 15 deletions(-) create mode 100644 test/utils/lookupTable/toto1Dlog.csv diff --git a/doc/source/programmingguide.rst b/doc/source/programmingguide.rst index ad2a70c4f..536b92102 100644 --- a/doc/source/programmingguide.rst +++ b/doc/source/programmingguide.rst @@ -536,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 ++++++++++++++++++++++ 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 63566ae5a..99e52e466 100644 --- a/src/utils/lookupTable.hpp +++ b/src/utils/lookupTable.hpp @@ -14,7 +14,24 @@ #include "lookupTable.hpp" #include "npy.hpp" -template +// 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)); + } +}; + +template class LookupTable { public: LookupTable() = default; @@ -26,15 +43,21 @@ class LookupTable { // 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 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) @@ -49,6 +72,34 @@ 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 getter for all kinds of input arrays template KOKKOS_INLINE_FUNCTION @@ -60,7 +111,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); @@ -108,7 +161,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 @@ -132,6 +185,9 @@ 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); } @@ -148,12 +204,15 @@ class LookupTable { } }; -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; @@ -242,6 +301,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); @@ -256,11 +318,14 @@ LookupTable::LookupTable(std::vector filenames, // 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) { +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"); } @@ -442,6 +507,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); @@ -470,13 +538,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); @@ -548,6 +619,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 b07672236..d46933932 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[] ) @@ -143,6 +157,80 @@ int main( int argc, char* argv[] ) } idfx::cout << "Success" << std::endl; + idfx::cout << "--------------------------------------" << std::endl; + idfx::cout << "Testing interpolation in function space on device." << std::endl; + // data = x^2 sampled in x=1,2,4. A linear interpolation in x=3 gives 10, while an + // interpolation performed on log(x) and log(data) gives the exact result 9 + LookupTable<1> 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 3D npy file on device." << std::endl; // Read npy File 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 From e58ff2d9af0d0b6f053c0bdfa2005f4daf126711 Mon Sep 17 00:00:00 2001 From: dutta-alankar Date: Wed, 19 Aug 2026 15:08:56 +0200 Subject: [PATCH 3/5] Add accessors to the neighbours used by the LookupTable interpolation GetNeighbours and GetNeighboursIndx (and their Host counterparts) return the elements of the table surrounding the requested coordinates: the two coordinates bracketing the input along each dimension, the data on the 2^nDim vertices of the surrounding cell, and the corresponding indices in the table. When the table is interpolated in function space, the returned coordinates and data are transformed back with invFunc, so that they are expressed in the original space of the table. Get is left unchanged, and the neighbour search is shared by the new methods through the GetIndices and GetDataIndex helpers. --- doc/source/programmingguide.rst | 68 +++++++++++++ src/utils/lookupTable.hpp | 167 ++++++++++++++++++++++++++++++++ test/utils/lookupTable/main.cpp | 138 ++++++++++++++++++++++++++ 3 files changed, 373 insertions(+) diff --git a/doc/source/programmingguide.rst b/doc/source/programmingguide.rst index 536b92102..714ef63e4 100644 --- a/doc/source/programmingguide.rst +++ b/doc/source/programmingguide.rst @@ -619,6 +619,74 @@ 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. + +.. 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/utils/lookupTable.hpp b/src/utils/lookupTable.hpp index 99e52e466..bf8a26d0e 100644 --- a/src/utils/lookupTable.hpp +++ b/src/utils/lookupTable.hpp @@ -100,6 +100,91 @@ class LookupTable { } } + // 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 @@ -191,6 +276,64 @@ class LookupTable { 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 { + int idx[kDim]; + real delta[kDim]; + + if(!GetIndices(x, dimensions, offset, xin, idx, delta)) { + 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(GetDataIndex(dimensions, 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 { + real delta[kDim]; + + if(!GetIndices(x, dimensions, offset, xin, idx, delta)) { + for(int n = 0 ; n < kDim ; n++) idx[n] = -1; + for(unsigned int n = 0 ; n < (1 << kDim) ; n++) dataIdx[n] = -1; + return; + } + + for(unsigned int n = 0 ; n < (1 << kDim) ; n++) { + dataIdx[n] = GetDataIndex(dimensions, idx, n); + } + } + // Getter on device KOKKOS_INLINE_FUNCTION real Get(const real x[kDim]) const { @@ -202,6 +345,30 @@ 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); + } }; template diff --git a/test/utils/lookupTable/main.cpp b/test/utils/lookupTable/main.cpp index d46933932..26e1ac829 100644 --- a/test/utils/lookupTable/main.cpp +++ b/test/utils/lookupTable/main.cpp @@ -231,6 +231,144 @@ int main( int argc, char* argv[] ) } 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 3D npy file on device." << std::endl; // Read npy File From 7109232cf973708d22e542faa815d74f83a0a6b1 Mon Sep 17 00:00:00 2001 From: dutta-alankar Date: Wed, 19 Aug 2026 15:26:05 +0200 Subject: [PATCH 4/5] Reuse the neighbour search between Get and the neighbour getters Get (and GetHost) now accept a LookupTableNeighbours structure in which they store the search they performed. Giving that same structure to GetNeighbours or GetNeighboursIndx makes them reuse this search instead of walking the table a second time, and they search it again as usual as soon as different coordinates are requested. The structure is held by the caller rather than by the lookup table, so that it is thread-private when it is declared inside an idefix_for: a lookup table is shared by all of the threads of a loop, and can hence not cache anything itself. The existing Get is left untouched. --- doc/source/programmingguide.rst | 40 ++++++++ src/utils/lookupTable.hpp | 169 ++++++++++++++++++++++++++++++-- test/utils/lookupTable/main.cpp | 150 ++++++++++++++++++++++++++++ 3 files changed, 350 insertions(+), 9 deletions(-) diff --git a/doc/source/programmingguide.rst b/doc/source/programmingguide.rst index 714ef63e4..994a94649 100644 --- a/doc/source/programmingguide.rst +++ b/doc/source/programmingguide.rst @@ -676,6 +676,46 @@ For instance, the 2D table ``example2D.csv`` above, interpolated in (x=2.1, y=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 a ``LookupTableNeighbours`` structure to ``Get``: the search it stores is then reused by +``GetNeighbours`` and ``GetNeighboursIndx`` instead of being done a second time. + +.. 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 diff --git a/src/utils/lookupTable.hpp b/src/utils/lookupTable.hpp index bf8a26d0e..3892fc319 100644 --- a/src/utils/lookupTable.hpp +++ b/src/utils/lookupTable.hpp @@ -31,6 +31,34 @@ struct LookupTableExp { } }; +// 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: @@ -276,6 +304,59 @@ class LookupTable { 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 @@ -286,10 +367,21 @@ class LookupTable { 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 { - int idx[kDim]; - real delta[kDim]; + 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(!GetIndices(x, dimensions, offset, xin, idx, delta)) { + 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; @@ -297,8 +389,8 @@ class LookupTable { // 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); + 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]); @@ -307,7 +399,7 @@ class LookupTable { // Data on each vertex of the neighbours for(unsigned int n = 0 ; n < (1 << kDim) ; n++) { - dataN[n] = data(GetDataIndex(dimensions, idx, n)); + dataN[n] = data(GetDataIndex(dimensions, neighbours.idx, n)); if(interpolateInFuncSpace) dataN[n] = invFunc(dataN[n]); } } @@ -321,16 +413,30 @@ class LookupTable { KOKKOS_INLINE_FUNCTION void GetNeighboursIndx(const real x[kDim], Tint &dimensions, Tint &offset, Treal &xin, int idx[kDim], int dataIdx[1 << kDim]) const { - real delta[kDim]; + 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(!GetIndices(x, dimensions, offset, xin, idx, delta)) { + 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, idx, n); + dataIdx[n] = GetDataIndex(dimensions, neighbours.idx, n); } } @@ -369,6 +475,51 @@ class LookupTable { 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 diff --git a/test/utils/lookupTable/main.cpp b/test/utils/lookupTable/main.cpp index 26e1ac829..64cebd7ce 100644 --- a/test/utils/lookupTable/main.cpp +++ b/test/utils/lookupTable/main.cpp @@ -369,6 +369,156 @@ int main( int argc, char* argv[] ) 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 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); + + 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 From 955350aca3d45bf3801ba6f0e458176e9ff22cbb Mon Sep 17 00:00:00 2001 From: dutta-alankar Date: Wed, 19 Aug 2026 15:42:32 +0200 Subject: [PATCH 5/5] Cover both orders of the LookupTable neighbour search reuse GetNeighbours and GetNeighboursIndx already search the table when they are given a structure which does not hold the requested coordinates yet, and Get reuses whatever search that structure holds, so the two orders are equivalent. Add the tests which check it, on the host and on the device, and document that the order of the calls is irrelevant. --- doc/source/programmingguide.rst | 6 ++- test/utils/lookupTable/main.cpp | 95 +++++++++++++++++++++++++++++++++ 2 files changed, 99 insertions(+), 2 deletions(-) diff --git a/doc/source/programmingguide.rst b/doc/source/programmingguide.rst index 994a94649..9d05e7c69 100644 --- a/doc/source/programmingguide.rst +++ b/doc/source/programmingguide.rst @@ -678,8 +678,10 @@ 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 a ``LookupTableNeighbours`` structure to ``Get``: the search it stores is then reused by -``GetNeighbours`` and ``GetNeighboursIndx`` instead of being done a second time. +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++ diff --git a/test/utils/lookupTable/main.cpp b/test/utils/lookupTable/main.cpp index 64cebd7ce..fdec245f4 100644 --- a/test/utils/lookupTable/main.cpp +++ b/test/utils/lookupTable/main.cpp @@ -453,6 +453,75 @@ int main( int argc, char* argv[] ) 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; @@ -494,6 +563,32 @@ int main( int argc, char* argv[] ) 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)