From a10c7945f9681e103e98237fdd089bf80b51f3f9 Mon Sep 17 00:00:00 2001 From: Aaron Ridley Date: Thu, 13 Mar 2025 06:55:10 -0400 Subject: [PATCH 01/10] FEAT: pass latitudes from one grid to the next --- src/grid_match.cpp | 168 ++++++++++++++++++++++++++++++++++++++++++--- 1 file changed, 159 insertions(+), 9 deletions(-) diff --git a/src/grid_match.cpp b/src/grid_match.cpp index 25da6188..39639dab 100644 --- a/src/grid_match.cpp +++ b/src/grid_match.cpp @@ -3,6 +3,65 @@ #include "aether.h" +// ----------------------------------------------------------------------------- +// Send arrays of variables to other processors on the given grid. +// ----------------------------------------------------------------------------- + +bool exchange_information(int64_t *nPointsToPass, + std::vector varToSend, + int64_t *nPointsToReceive, + std::vector varToReceive) { + + int64_t jNode, iPt, iTag, iProcTo, iProcFrom; + std::vector requests(nGrids); + + // Here we send the message into the wind: + // - if it is the same processor, just copy the information + // - if it is a different processor, send the data + for (jNode = 0; jNode < nGrids ; jNode++) { + if (jNode == iGrid) { + for (iPt = 0; iPt < nPointsToPass[jNode]; iPt ++) { + varToReceive[jNode][iPt] = varToSend[jNode][iPt]; + } + } else { + iProcTo = iMember * nGrids + jNode; + // iTag is a unique id allowing all processors to + // communicate asynchronously + iTag = iProc * 10000 + iProcTo; + MPI_Isend(varToSend[jNode], + nPointsToPass[jNode] * sizeof(precision_t), + MPI_BYTE, + iProcTo, + iTag, + aether_comm, + &requests[jNode]); + } + } + + // Wait for everyone to get the information that was sent: + for (jNode = 0; jNode < nGrids ; jNode++) + if (jNode != iGrid) + MPI_Wait(&requests[jNode], MPI_STATUS_IGNORE); + + // Receive it into the receiving array: + for (jNode = 0; jNode < nGrids ; jNode++) + if (jNode != iGrid) { + iProcFrom = iMember * nGrids + jNode; + // Rebuid the unique id: + iTag = iProcFrom * 10000 + iProc; + MPI_Recv(varToReceive[jNode], + nPointsToReceive[jNode] * sizeof(precision_t), + MPI_BYTE, + jNode, + iTag, + aether_comm, + MPI_STATUS_IGNORE); + } + + MPI_Barrier(aether_comm); + return true; +} + bool grid_match(Grid gGrid, Grid mGrid, Quadtree gQuadtree, @@ -17,8 +76,16 @@ bool grid_match(Grid gGrid, precision_t lon, lat; precision_t normX, normY, normZ; arma_vec norms(3); - int64_t iNode; + int64_t jNode, kNode; + int64_t *nPointsToPass = static_cast(malloc(nGrids * sizeof(int64_t))); + int64_t *nPointsToReceive = static_cast(malloc(nGrids * sizeof(int64_t))); + int64_t *nPointsDummy = static_cast(malloc(nGrids * sizeof(int64_t))); + + for (jNode = 0; jNode < nGrids ; jNode++) + nPointsToPass[jNode] = 0; + // This is not the most efficient way to do this, but the first pass, let's + // just count how many points we need to send to the other processors: for (iX = mGCs; iX < mnX - mGCs; iX++) { for (iY = mGCs; iY < mnY - mGCs; iY++) { for (iZ = mGCs; iZ < mnZ - mGCs; iZ++) { @@ -28,19 +95,102 @@ bool grid_match(Grid gGrid, norms(0) = lon / cPI; norms(1) = lat / cPI; norms(2) = 0.0; - iNode = gQuadtree.find_point(norms); + jNode = gQuadtree.find_point(norms); } else { norms = sphere_to_cube(lon, lat); - iNode = gQuadtree.find_point(norms); + jNode = gQuadtree.find_point(norms); + } + if (jNode < 0 || jNode >= nGrids) { + std::cout << "out of bounds!!! " << jNode << "\n"; } - std::cout << "lon, lat, node: " << lon*cRtoD << " " - << lat*cRtoD << " " - << norms(0) << " " - << norms(1) << " " - << norms(2) << " " - << iNode << "\n"; + nPointsToPass[jNode] = nPointsToPass[jNode]+1; + /* std::cout << "lon, lat, node: " << lon*cRtoD << " " + << lat*cRtoD << " " + << norms(0) << " " + << norms(1) << " " + << norms(2) << " " + << jNode << " " + << iProc << " " + << nPoints[jNode] << "\n"; */ } } } + std::cout << "made it here: " << iProc << "\n"; + MPI_Barrier(aether_comm); + + for (jNode = 0; jNode < nGrids ; jNode++) + std::cout << "nPtsToPass : " << iProc << " " << nPointsToPass[jNode] << "\n"; + + std::cout << "sending number of points :\n"; + + // This section sends the number of points that need to be transfered to each processor. + // Then the processor saves the number of points, so it can be remembered, and both the + // sender and receiver will have the information. + for (jNode = 0; jNode < nGrids ; jNode++) { + if (jNode == iGrid) { + for (kNode = 0; kNode < nGrids ; kNode++) + nPointsDummy[kNode] = nPointsToPass[kNode]; + } + MPI_Bcast(nPointsDummy, nGrids, MPI_INT64_T, jNode, aether_comm); + nPointsToReceive[jNode] = nPointsDummy[iGrid]; + } + + MPI_Barrier(aether_comm); + + for (jNode = 0; jNode < nGrids ; jNode++) { + std::cout << "nPtsToReceive : " << iProc << " " << jNode << " " << nPointsToReceive[jNode] << "\n"; + MPI_Barrier(aether_comm); + } + + // Now we need to create an array of send points and an array of receive points. + std::vector latsToPass(nGrids); + std::vector lonsToPass(nGrids); + std::vector altsToPass(nGrids); + for (jNode = 0; jNode < nGrids ; jNode++) { + latsToPass[jNode] = static_cast(malloc(nPointsToPass[jNode] * sizeof(precision_t))); + lonsToPass[jNode] = static_cast(malloc(nPointsToPass[jNode] * sizeof(precision_t))); + altsToPass[jNode] = static_cast(malloc(nPointsToPass[jNode] * sizeof(precision_t))); + } + + std::vector latsToInterTo(nGrids); + std::vector lonsToInterTo(nGrids); + std::vector altsToInterTo(nGrids); + for (jNode = 0; jNode < nGrids ; jNode++) { + latsToInterTo[jNode] = static_cast(malloc(nPointsToReceive[jNode] * sizeof(precision_t))); + lonsToInterTo[jNode] = static_cast(malloc(nPointsToReceive[jNode] * sizeof(precision_t))); + altsToInterTo[jNode] = static_cast(malloc(nPointsToReceive[jNode] * sizeof(precision_t))); + } + + // now, the second pass, let's store the information so we can pass it: + for (jNode = 0; jNode < nGrids ; jNode++) + nPointsToPass[jNode] = 0; + for (iX = mGCs; iX < mnX - mGCs; iX++) { + for (iY = mGCs; iY < mnY - mGCs; iY++) { + for (iZ = mGCs; iZ < mnZ - mGCs; iZ++) { + lon = mGrid.geoLon_scgc(iX, iY, iZ); + lat = mGrid.geoLat_scgc(iX, iY, iZ); + if (gGrid.iGridShape_ == gGrid.iSphere_) { + norms(0) = lon / cPI; + norms(1) = lat / cPI; + norms(2) = 0.0; + jNode = gQuadtree.find_point(norms); + } else { + norms = sphere_to_cube(lon, lat); + jNode = gQuadtree.find_point(norms); + } + latsToPass[jNode][nPointsToPass[jNode]] = lat; + lonsToPass[jNode][nPointsToPass[jNode]] = lon; + altsToPass[jNode][nPointsToPass[jNode]] = mGrid.geoAlt_scgc(iX, iY, iZ); + nPointsToPass[jNode] = nPointsToPass[jNode]+1; + } + } + } + bool didWork; + didWork = exchange_information(nPointsToPass, + latsToPass, + nPointsToReceive, + latsToInterTo); + + return true; } From 243b38debb42cc2dde9a9dab9a8a515f9bece2ff Mon Sep 17 00:00:00 2001 From: Aaron Ridley Date: Tue, 20 May 2025 19:57:08 -0400 Subject: [PATCH 02/10] FEAT: get grid to grid message passing to work --- src/grid_match.cpp | 164 +++++++++++++++++++++++++++++++++++++++------ 1 file changed, 145 insertions(+), 19 deletions(-) diff --git a/src/grid_match.cpp b/src/grid_match.cpp index 39639dab..ecb62c92 100644 --- a/src/grid_match.cpp +++ b/src/grid_match.cpp @@ -62,11 +62,27 @@ bool exchange_information(int64_t *nPointsToPass, return true; } -bool grid_match(Grid gGrid, - Grid mGrid, +// ----------------------------------------------------------------------------- +// This function: +// on the requesting information side: +// - figures out which processor each point of the other grid is on +// - counts the points for each processor +// - exchanges how many points to pass for each processor +// - makes lists of coordinates to send to each processor +// - sends those lists +// on the interpolator side: +// - builds interpolators for the requested information +// ----------------------------------------------------------------------------- + +bool grid_match(Grid &gGrid, + Grid &mGrid, Quadtree gQuadtree, Quadtree mQuadtree) { + std::string function = "grid_match"; + static int iFunction = -1; + report.enter(function, iFunction); + // Let's do magnetic to geographic first: int64_t iX, mnX = mGrid.get_nX(); @@ -86,6 +102,7 @@ bool grid_match(Grid gGrid, // This is not the most efficient way to do this, but the first pass, let's // just count how many points we need to send to the other processors: + mGrid.gridToGridMap.set_size(mnX, mnY, mnZ); for (iX = mGCs; iX < mnX - mGCs; iX++) { for (iY = mGCs; iY < mnY - mGCs; iY++) { for (iZ = mGCs; iZ < mnZ - mGCs; iZ++) { @@ -103,6 +120,7 @@ bool grid_match(Grid gGrid, if (jNode < 0 || jNode >= nGrids) { std::cout << "out of bounds!!! " << jNode << "\n"; } + mGrid.gridToGridMap(iX, iY, iZ) = jNode; nPointsToPass[jNode] = nPointsToPass[jNode]+1; /* std::cout << "lon, lat, node: " << lon*cRtoD << " " << lat*cRtoD << " " @@ -115,13 +133,13 @@ bool grid_match(Grid gGrid, } } } - std::cout << "made it here: " << iProc << "\n"; MPI_Barrier(aether_comm); - for (jNode = 0; jNode < nGrids ; jNode++) - std::cout << "nPtsToPass : " << iProc << " " << nPointsToPass[jNode] << "\n"; - - std::cout << "sending number of points :\n"; + if (report.test_verbose(3)) { + for (jNode = 0; jNode < nGrids ; jNode++) + std::cout << "nPtsToPass : " << iProc << " " << nPointsToPass[jNode] << "\n"; + std::cout << "sending number of points :\n"; + } // This section sends the number of points that need to be transfered to each processor. // Then the processor saves the number of points, so it can be remembered, and both the @@ -135,27 +153,23 @@ bool grid_match(Grid gGrid, nPointsToReceive[jNode] = nPointsDummy[iGrid]; } - MPI_Barrier(aether_comm); - - for (jNode = 0; jNode < nGrids ; jNode++) { - std::cout << "nPtsToReceive : " << iProc << " " << jNode << " " << nPointsToReceive[jNode] << "\n"; - MPI_Barrier(aether_comm); + if (report.test_verbose(3)) { + for (jNode = 0; jNode < nGrids ; jNode++) { + std::cout << "nPtsToReceive : " << iProc << " " << jNode << " " << nPointsToReceive[jNode] << "\n"; + } } // Now we need to create an array of send points and an array of receive points. std::vector latsToPass(nGrids); std::vector lonsToPass(nGrids); std::vector altsToPass(nGrids); - for (jNode = 0; jNode < nGrids ; jNode++) { - latsToPass[jNode] = static_cast(malloc(nPointsToPass[jNode] * sizeof(precision_t))); - lonsToPass[jNode] = static_cast(malloc(nPointsToPass[jNode] * sizeof(precision_t))); - altsToPass[jNode] = static_cast(malloc(nPointsToPass[jNode] * sizeof(precision_t))); - } - std::vector latsToInterTo(nGrids); std::vector lonsToInterTo(nGrids); std::vector altsToInterTo(nGrids); for (jNode = 0; jNode < nGrids ; jNode++) { + latsToPass[jNode] = static_cast(malloc(nPointsToPass[jNode] * sizeof(precision_t))); + lonsToPass[jNode] = static_cast(malloc(nPointsToPass[jNode] * sizeof(precision_t))); + altsToPass[jNode] = static_cast(malloc(nPointsToPass[jNode] * sizeof(precision_t))); latsToInterTo[jNode] = static_cast(malloc(nPointsToReceive[jNode] * sizeof(precision_t))); lonsToInterTo[jNode] = static_cast(malloc(nPointsToReceive[jNode] * sizeof(precision_t))); altsToInterTo[jNode] = static_cast(malloc(nPointsToReceive[jNode] * sizeof(precision_t))); @@ -186,11 +200,123 @@ bool grid_match(Grid gGrid, } } bool didWork; + // Pass first coordinate (lons) + didWork = exchange_information(nPointsToPass, + lonsToPass, + nPointsToReceive, + lonsToInterTo); + // Pass second coordinate (lats) didWork = exchange_information(nPointsToPass, latsToPass, nPointsToReceive, latsToInterTo); + // Pass third coordinate (alts): + didWork = exchange_information(nPointsToPass, + altsToPass, + nPointsToReceive, + altsToInterTo); + if (report.test_verbose(2)) { + for (jNode = 0; jNode < nGrids ; jNode++) { + std::cout << "Received the following points from iGrid = " << jNode << "\n"; + std::cout << " -> points received : " << nPointsToReceive[jNode] << "\n"; + for (int64_t iPt = 0; iPt < nPointsToReceive[jNode]; iPt++) + std::cout << " -> " << iPt << " " + << lonsToInterTo[jNode][iPt] << " " + << latsToInterTo[jNode][iPt] << " " + << altsToInterTo[jNode][iPt] << "\n"; + } + } - return true; + struct grid_to_grid_t oneGrid; + + int64_t nPts; + for (jNode = 0; jNode < nGrids ; jNode++) { + // These are backwards now, since we will switch sender and reciever: + oneGrid.nPts = nPointsToReceive[jNode]; + oneGrid.nPtsReceive = nPointsToPass[jNode]; + oneGrid.iProcTo = iMember * nGrids + jNode; + if (report.test_verbose(2)) + std::cout << "Making interpolation coefficients for : " << jNode + << "; points : " << oneGrid.nPts << "\n"; + if (oneGrid.nPts > 0) { + // Interpolation function takes vectors, + // so transfer these arrays to vectors: + std::vector Lons(oneGrid.nPts); + std::vector Lats(oneGrid.nPts); + std::vector Alts(oneGrid.nPts); + for (int64_t iPt = 0; iPt < oneGrid.nPts; iPt++) { + Lons[iPt] = lonsToInterTo[jNode][iPt]; + Lats[iPt] = latsToInterTo[jNode][iPt]; + Alts[iPt] = altsToInterTo[jNode][iPt]; + } + oneGrid.interpCoefs = gGrid.get_interpolation_coefs(Lons, Lats, Alts); + } + gGrid.gridToGridCoefs.push_back(oneGrid); + } + + report.exit(function); + return didWork; } + +bool get_data_from_other_grid(Grid &gGrid, + Grid &mGrid, + arma_cube &gData, + arma_cube &mData) { + + std::string function = "get_data_from_other_grid"; + static int iFunction = -1; + report.enter(function, iFunction); + + int64_t jNode, iPt; + std::vector dataToSend(nGrids); + std::vector dataToReceive(nGrids); + int64_t *nPointsToSend = static_cast(malloc(nGrids * sizeof(int64_t))); + int64_t *nPointsToReceive = static_cast(malloc(nGrids * sizeof(int64_t))); + + for (jNode = 0; jNode < nGrids ; jNode++) { + if (report.test_verbose(2)) + std::cout << "nPts : " << jNode << " " << gGrid.gridToGridCoefs[jNode].nPts << "\n"; + nPointsToSend[jNode] = gGrid.gridToGridCoefs[jNode].nPts; + nPointsToReceive[jNode] = gGrid.gridToGridCoefs[jNode].nPtsReceive; + dataToSend[jNode] = static_cast(malloc(gGrid.gridToGridCoefs[jNode].nPts * sizeof(precision_t))); + dataToReceive[jNode] = static_cast(malloc(gGrid.gridToGridCoefs[jNode].nPtsReceive * sizeof(precision_t))); + std::vector values = gGrid.get_interpolation_values(gData, gGrid.gridToGridCoefs[jNode].interpCoefs); + + for (iPt = 0; iPt < gGrid.gridToGridCoefs[jNode].nPts; iPt++) { + dataToSend[jNode][iPt] = values[iPt]; + if (report.test_verbose(2)) + std::cout << "datatosend : " << iPt << " " << dataToSend[jNode][iPt] << "\n"; + } + } + bool didWork = exchange_information(nPointsToSend, + dataToSend, + nPointsToReceive, + dataToReceive); + int64_t iX, mnX = mGrid.get_nX(); + int64_t iY, mnY = mGrid.get_nY(); + int64_t iZ, mnZ = mGrid.get_nZ(); + int64_t mGCs = mGrid.get_nGCs(); + std::vector iCounter(nGrids); + for (jNode = 0; jNode < nGrids ; jNode++) + iCounter[jNode] = 0; + + for (iX = mGCs; iX < mnX - mGCs; iX++) { + for (iY = mGCs; iY < mnY - mGCs; iY++) { + for (iZ = mGCs; iZ < mnZ - mGCs; iZ++) { + jNode = mGrid.gridToGridMap(iX, iY, iZ); + if (report.test_verbose(2)) { + std::cout << "unpacking point : " << iX << " " << iY << " " << iZ << " " << jNode << " " + << iCounter[jNode] << " " << dataToReceive[jNode][iCounter[jNode]] << "\n"; + } + + mData(iX, iY, iZ) = dataToReceive[jNode][iCounter[jNode]]; + iCounter[jNode] = iCounter[jNode]+1; + } + } + } + + report.exit(function); + return true; + +} \ No newline at end of file From fc6f18b497d14aaefb209ba379a923d005811690 Mon Sep 17 00:00:00 2001 From: Aaron Ridley Date: Tue, 20 May 2025 19:57:50 -0400 Subject: [PATCH 03/10] FEAT: needed for grid to grid MP --- include/calc_grid_derived.h | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/include/calc_grid_derived.h b/include/calc_grid_derived.h index 28d5237f..046fbd46 100644 --- a/include/calc_grid_derived.h +++ b/include/calc_grid_derived.h @@ -19,9 +19,14 @@ arma_vec calc_bin_widths(arma_vec centers); // ---------------------------------------------------------------------------- // A helper function for mapping grids // ---------------------------------------------------------------------------- -bool grid_match(Grid gGrid, - Grid mGrid, +bool grid_match(Grid &gGrid, + Grid &mGrid, Quadtree gQuadtree, Quadtree mQuadtree); +bool get_data_from_other_grid(Grid &gGrid, + Grid &mGrid, + arma_cube &gData, + arma_cube &mData); + #endif // INCLUDE_CALC_GRID_DERIVED_H_ From 406642dfaa91ca37d0a031a9e3296da51a03aef7 Mon Sep 17 00:00:00 2001 From: Aaron Ridley Date: Tue, 20 May 2025 19:58:44 -0400 Subject: [PATCH 04/10] FEAT: needed for grid to grid MP --- include/grid.h | 80 +++++++++++++++++++++++++++++++++++++------------- 1 file changed, 59 insertions(+), 21 deletions(-) diff --git a/include/grid.h b/include/grid.h index fed807e9..b56fead2 100644 --- a/include/grid.h +++ b/include/grid.h @@ -11,6 +11,33 @@ // Grid class // ---------------------------------------------------------------------------- + struct interp_coef_t + { + // The point is inside the cube of [iRow, iRow+1], [iCol, iCol+1], [iAlt, iAlt+1] + uint64_t iRow; + uint64_t iCol; + uint64_t iAlt; + // The coefficients along row, column and altitude + precision_t rRow; + precision_t rCol; + precision_t rAlt; + // Whether the point is within this grid or not + bool in_grid; + // If this is set to true: + bool above_grid, below_grid; + // do interpolation in lat and lon, but extrapolate in altitude + }; + + struct grid_to_grid_t { + int64_t iProcTo; + int64_t nPts; + int64_t nPtsReceive; + std::vector interpCoefs; + std::vector valueToSend; + std::vector valueToReceive; + }; + + class Grid { @@ -20,13 +47,24 @@ class Grid const int iDipole_ = 3; int iGridShape_ = -1; - // Armidillo Cube Versions: + // The index and coefficient used for interpolation + // Each point is processed by the function set_interpolation_coefs and stored + // in the form of this structure. + // If the point is out of the grid, in_grid = false and all other members are undefined + + std::vector gridToGridCoefs; + arma::Cube gridToGridMap; + + // Armadillo Cube Versions: // Cell Center Coordinates arma_cube geoLon_scgc, geoX_scgc; arma_cube geoLat_scgc, geoY_scgc; arma_cube geoAlt_scgc, geoZ_scgc; arma_cube geoLocalTime_scgc; + // This is an array for testing things: + arma_cube test_scgc; + // Reference coordinate arma_cube refx_scgc, refy_scgc; @@ -434,6 +472,22 @@ class Grid bool set_interpolation_coefs(const std::vector &Lons, const std::vector &Lats, const std::vector &Alts); + + /** + * \brief Set the interpolation coefficients + * \param Lons The longitude of points + * \param Lats The latitude of points + * \param Alts The altitude of points + * \pre This instance is an geo grid + * \pre Lons, Lats and Alts have the same size + * \return list of interpolation coefficients + */ + + std::vector get_interpolation_coefs( + const std::vector &Lons, + const std::vector &Lats, + const std::vector &Alts); + /** * \brief Create a map of geographic locations to data and do the interpolation * \param data The value at the positions of geoLon, geoLat, and geoAlt @@ -443,6 +497,8 @@ class Grid * an empty vector if the data is not the same size as the geo grid. */ std::vector get_interpolation_values(const arma_cube &data) const; + std::vector get_interpolation_values(arma_cube data, + std::vector coefArray); private: bool IsGeoGrid; @@ -513,24 +569,6 @@ class Grid bool col_max_exclusive; }; - // The index and coefficient used for interpolation - // Each point is processed by the function set_interpolation_coefs and stored - // in the form of this structure. - // If the point is out of the grid, in_grid = false and all other members are undefined - struct interp_coef_t - { - // The point is inside the cube of [iRow, iRow+1], [iCol, iCol+1], [iAlt, iAlt+1] - uint64_t iRow; - uint64_t iCol; - uint64_t iAlt; - // The coefficients along row, column and altitude - precision_t rRow; - precision_t rCol; - precision_t rAlt; - // Whether the point is within this grid or not - bool in_grid; - }; - // Return the index of the last element that has altitude smaller than or euqal to the input uint64_t search_altitude(const precision_t alt_in) const; @@ -540,11 +578,11 @@ class Grid void get_cubesphere_grid_range(struct cubesphere_range &cr) const; // Helper function for set_interpolation_coefs - void set_interp_coef_sphere(const sphere_range &sr, + struct interp_coef_t get_interp_coef_sphere(const sphere_range &sr, const precision_t lon_in, const precision_t lat_in, const precision_t alt_in); - void set_interp_coef_cubesphere(const cubesphere_range &cr, + struct interp_coef_t get_interp_coef_cubesphere(const cubesphere_range &cr, const precision_t lon_in, const precision_t lat_in, const precision_t alt_in); From b365ed045fd52c6482f79156b4a3f13812209848 Mon Sep 17 00:00:00 2001 From: Aaron Ridley Date: Tue, 20 May 2025 19:59:30 -0400 Subject: [PATCH 05/10] FEAT: testing grid to grid MP --- src/advance.cpp | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/src/advance.cpp b/src/advance.cpp index 023a73d5..c6252c40 100644 --- a/src/advance.cpp +++ b/src/advance.cpp @@ -46,6 +46,11 @@ bool advance(Planets &planet, didWork = neutralsMag.check_for_nonfinites("Top of Advance - ion grid"); } + // here we are going to grab stuff from the neutral grid and put it on the + // ion grid + didWork = get_data_from_other_grid(gGrid, mGrid, neutrals.temperature_scgc, mGrid.test_scgc); + + gGrid.calc_sza(planet, time); mGrid.calc_sza(planet, time); From e9c01dea0f36fd321b6977f02f0b20955a3674d8 Mon Sep 17 00:00:00 2001 From: Aaron Ridley Date: Tue, 20 May 2025 20:00:02 -0400 Subject: [PATCH 06/10] FEAT: testing grid to grid MP --- src/grid.cpp | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/grid.cpp b/src/grid.cpp index 9ed462c3..eb9f8798 100644 --- a/src/grid.cpp +++ b/src/grid.cpp @@ -80,6 +80,9 @@ Grid::Grid(std::string gridtype) { geoAlt_scgc.set_size(nX, nY, nZ); geoLocalTime_scgc.set_size(nX, nY, nZ); + test_scgc.set_size(nX, nY, nZ); + test_scgc.zeros(); + refx_scgc.set_size(nX, nY, nZ); refy_scgc.set_size(nX, nY, nZ); refx_angle.set_size(nX, nY, nZ); From 7bde45e5c566d58ef7efeda5f2c8be5e2d3c3d17 Mon Sep 17 00:00:00 2001 From: Aaron Ridley Date: Tue, 20 May 2025 20:00:43 -0400 Subject: [PATCH 07/10] STY: remove space --- src/main/main.cpp | 1 - 1 file changed, 1 deletion(-) diff --git a/src/main/main.cpp b/src/main/main.cpp index 9d0a0a3f..946e293e 100644 --- a/src/main/main.cpp +++ b/src/main/main.cpp @@ -100,7 +100,6 @@ int main() { didWork = mGrid.init_geo_grid(quadtree, planet); mGrid.set_IsGeoGrid(false); } - didWork = grid_match(gGrid, mGrid, quadtree, quadtree_ion); // Initialize Neutrals on geographic grid: From 21c934dd277c30ee0881d9c381225606b3d18295 Mon Sep 17 00:00:00 2001 From: Aaron Ridley Date: Tue, 20 May 2025 20:01:27 -0400 Subject: [PATCH 08/10] FEAT: make an output file for test results --- src/output.cpp | 11 ++++++++++- 1 file changed, 10 insertions(+), 1 deletion(-) diff --git a/src/output.cpp b/src/output.cpp index 2a16a8c0..97a85a1c 100644 --- a/src/output.cpp +++ b/src/output.cpp @@ -41,6 +41,9 @@ std::string get_filename_from_type(std::string type_output) { if (type_output == "therm") filename = "3DTH"; + if (type_output == "test") + filename = "3DTE"; + return filename; } @@ -153,7 +156,6 @@ bool output(const Neutrals &neutrals, store_variable("density_" + neutrals.species[iSpecies].cName, neutrals.density_unit, neutrals.species[iSpecies].density_scgc); - // Neutral Temperature: if (type_output == "neutrals" || type_output == "states") @@ -322,6 +324,13 @@ bool output(const Neutrals &neutrals, grid.cent_acc_vcgc[2]); } + // Neutral Temperature: + if (type_output == "test") + AllOutputContainers[iOutput]. + store_variable("test_grid", + "none", + grid.test_scgc); + // ------------------------------------------------------------ // Set output file names From b7b2449d6ec2f8e394d939101906e54981334f03 Mon Sep 17 00:00:00 2001 From: Aaron Ridley Date: Tue, 20 May 2025 20:05:06 -0400 Subject: [PATCH 09/10] FEAT: allow above and below grid interpolation --- src/solver_grid_interpolation.cpp | 183 +++++++++++++++++++++++++----- 1 file changed, 155 insertions(+), 28 deletions(-) diff --git a/src/solver_grid_interpolation.cpp b/src/solver_grid_interpolation.cpp index 39bb2dc3..29379f34 100644 --- a/src/solver_grid_interpolation.cpp +++ b/src/solver_grid_interpolation.cpp @@ -184,7 +184,8 @@ void Grid::get_cubesphere_grid_range(struct cubesphere_range &cr) const { // Almost the copy of interp_sphere_linear_helper // -------------------------------------------------------------------------- -void Grid::set_interp_coef_sphere(const sphere_range &sr, +struct interp_coef_t Grid::get_interp_coef_sphere( + const sphere_range &sr, const precision_t lon_in, const precision_t lat_in, const precision_t alt_in) { @@ -200,13 +201,15 @@ void Grid::set_interp_coef_sphere(const sphere_range &sr, // Determine whether the point is inside this grid // Treat north pole specially because latitude is inclusive for both -cPI/2 and cPI/2 + // Dpn't check for altitude here! if (lon_in < sr.lon_min || lon_in >= sr.lon_max || lat_in < sr.lat_min - || lat_in > sr.lat_max || (lat_in == sr.lat_max && sr.lat_max != cPI / 2) - || alt_in < sr.alt_min || alt_in > sr.alt_max) { - interp_coefs.push_back(coef); - return; + || lat_in > sr.lat_max || (lat_in == sr.lat_max && sr.lat_max != cPI / 2)) { + return coef; } + // This point is in the grid! + coef.in_grid = true; + // ASSUMPTION: LONGITUDE AND LATITUDE ARE LINEARLY SPACED, nGCs >= 1 // For the cell containing it, directly calculate its x and y index // Find its z index using binary search @@ -228,13 +231,27 @@ void Grid::set_interp_coef_sphere(const sphere_range &sr, // The altitude may not be linearly spaced, so use binary search to find // the first element smaller than or equal to the altitude of the give point // Implemented in search_altitude - coef.iAlt = search_altitude(alt_in); - coef.rAlt = (alt_in - geoAlt_scgc(0, 0, coef.iAlt)) - / (geoAlt_scgc(0, 0, coef.iAlt + 1) - geoAlt_scgc(0, 0, coef.iAlt)); - // Put the coefficient into the vector - coef.in_grid = true; - interp_coefs.push_back(coef); + if (alt_in < sr.alt_min) { + coef.iAlt = nGCs; + coef.rAlt = alt_in - sr.alt_min; + coef.below_grid = true; + coef.above_grid = false; + } else { + if (alt_in > sr.alt_max) { + coef.iAlt = nAlts - nGCs; + coef.rAlt = alt_in - sr.alt_max; + coef.below_grid = false; + coef.above_grid = true; + } else { + coef.iAlt = search_altitude(alt_in); + coef.rAlt = (alt_in - geoAlt_scgc(0, 0, coef.iAlt)) + / (geoAlt_scgc(0, 0, coef.iAlt + 1) - geoAlt_scgc(0, 0, coef.iAlt)); + coef.below_grid = false; + coef.above_grid = false; + } + } + return coef; } // -------------------------------------------------------------------------- @@ -242,7 +259,7 @@ void Grid::set_interp_coef_sphere(const sphere_range &sr, // Almost the copy of interp_cubesphere_linear_helper // -------------------------------------------------------------------------- -void Grid::set_interp_coef_cubesphere(const cubesphere_range &cr, +struct interp_coef_t Grid::get_interp_coef_cubesphere(const cubesphere_range &cr, const precision_t lon_in, const precision_t lat_in, const precision_t alt_in) { @@ -259,8 +276,7 @@ void Grid::set_interp_coef_cubesphere(const cubesphere_range &cr, // Determine whether the projection point is on the surface of the grid if (surface_in != cr.surface_number) { - interp_coefs.push_back(coef); - return; + return coef; } // Calculate the theoretical fractional row index and column index @@ -280,12 +296,13 @@ void Grid::set_interp_coef_cubesphere(const cubesphere_range &cr, || row_frac_index > row_index_max || (row_frac_index == row_index_max && cr.row_max_exclusive) || col_frac_index > col_index_max || (col_frac_index == col_index_max && - cr.col_max_exclusive) - || alt_in < cr.alt_min || alt_in > cr.alt_max) { - interp_coefs.push_back(coef); - return; + cr.col_max_exclusive)) { + return coef; } + // This point is in the grid! + coef.in_grid = true; + // Get the real integer index and the interpolation coefficient uint64_t row_index, col_index, alt_index; precision_t rRow, rCol, rAlt; @@ -303,13 +320,27 @@ void Grid::set_interp_coef_cubesphere(const cubesphere_range &cr, coef.rCol = col_frac_index - coef.iCol; coef.iCol += nGCs - 1; // Use binary search to find the index for altitude - coef.iAlt = search_altitude(alt_in); - coef.rAlt = (alt_in - geoAlt_scgc(0, 0, coef.iAlt)) - / (geoAlt_scgc(0, 0, coef.iAlt + 1) - geoAlt_scgc(0, 0, coef.iAlt)); - // Put the coefficient into the vector - coef.in_grid = true; - interp_coefs.push_back(coef); + if (alt_in < cr.alt_min) { + coef.iAlt = nGCs; + coef.rAlt = 0.0; + coef.below_grid = true; + coef.above_grid = false; + } else { + if (alt_in > cr.alt_max) { + coef.iAlt = nAlts - nGCs; + coef.rAlt = 0.0; + coef.below_grid = false; + coef.above_grid = true; + } else { + coef.iAlt = search_altitude(alt_in); + coef.rAlt = (alt_in - geoAlt_scgc(0, 0, coef.iAlt)) + / (geoAlt_scgc(0, 0, coef.iAlt + 1) - geoAlt_scgc(0, 0, coef.iAlt)); + coef.below_grid = false; + coef.above_grid = false; + } + } + return coef; } // -------------------------------------------------------------------------- @@ -319,6 +350,8 @@ void Grid::set_interp_coef_cubesphere(const cubesphere_range &cr, bool Grid::set_interpolation_coefs(const std::vector &Lons, const std::vector &Lats, const std::vector &Alts) { + + struct interp_coef_t coef; // If this is not a geo grid, return false if (!IsGeoGrid) return false; @@ -337,21 +370,84 @@ bool Grid::set_interpolation_coefs(const std::vector &Lons, get_cubesphere_grid_range(cr); // Calculate the index and coefficients for each point - for (size_t i = 0; i < Lons.size(); ++i) - set_interp_coef_cubesphere(cr, Lons[i], Lats[i], Alts[i]); + for (size_t i = 0; i < Lons.size(); ++i) { + coef = get_interp_coef_cubesphere(cr, Lons[i], Lats[i], Alts[i]); + interp_coefs.push_back(coef); + } + } else { // Calculate the range of the grid struct sphere_range sr; get_sphere_grid_range(sr); // Calculate the index and coefficients for each point - for (size_t i = 0; i < Lons.size(); ++i) - set_interp_coef_sphere(sr, Lons[i], Lats[i], Alts[i]); + for (size_t i = 0; i < Lons.size(); ++i) { + coef = get_interp_coef_sphere(sr, Lons[i], Lats[i], Alts[i]); + interp_coefs.push_back(coef); + } } return true; } +// -------------------------------------------------------------------------- +// Set the interpolation coefficients +// (v2 - return a list of interpolation coefficients) +// -------------------------------------------------------------------------- + +std::vector Grid::get_interpolation_coefs( + const std::vector &Lons, + const std::vector &Lats, + const std::vector &Alts) { + + int64_t nPts = Lons.size(), iPt; + std::vector listOfCoefs; + struct interp_coef_t singleCoef; + bool isBad = false; + + // If this is not a geo grid, return false + if (!IsGeoGrid) + isBad = true; + + // If the size of Lons, Lats and Alts are not the same, return false + if (Lons.size() != Lats.size() || Lats.size() != Alts.size()) + isBad = true; + + if (isBad) { + for (iPt = 0; iPt < nPts; ++iPt) { + // Put the coefficient into the vector + singleCoef.in_grid = false; + listOfCoefs.push_back(singleCoef); + } + return listOfCoefs; + } + + // Handle according to whether it is cubesphere or not + if (IsCubeSphereGrid) { + // Calculate the range of the grid + struct cubesphere_range cr; + get_cubesphere_grid_range(cr); + + // Calculate the index and coefficients for each point + for (iPt = 0; iPt < nPts; ++iPt) { + singleCoef = get_interp_coef_cubesphere(cr, Lons[iPt], Lats[iPt], Alts[iPt]); + listOfCoefs.push_back(singleCoef); + } + } else { + // Calculate the range of the grid + struct sphere_range sr; + get_sphere_grid_range(sr); + + // Calculate the index and coefficients for each point + for (iPt = 0; iPt < nPts; ++iPt) { + singleCoef = get_interp_coef_sphere(sr, Lons[iPt], Lats[iPt], Alts[iPt]); + listOfCoefs.push_back(singleCoef); + } + } + + return listOfCoefs; +} + // -------------------------------------------------------------------------- // Do the interpolation based on the coefficients stored in interp_coefs // -------------------------------------------------------------------------- @@ -381,3 +477,34 @@ std::vector Grid::get_interpolation_values( return ans; } + +// -------------------------------------------------------------------------- +// Do the interpolation based on the coefficients passed in +// -------------------------------------------------------------------------- + +std::vector Grid::get_interpolation_values(arma_cube data, + std::vector coefArray ) { + std::vector ans; + + // If the size of data is not the same as the size of grid, return an empty vector + if (data.n_rows != nLons || data.n_cols != nLats || data.n_slices != nAlts) + return ans; + + for (auto &it : coefArray) { + // Do interpolation if in_grid = true. Push cNinf otherwise + if (it.in_grid) { + ans.push_back(interpolate_unit_cube( + data.subcube(it.iRow, it.iCol, it.iAlt, unit_cube_size), + it.rRow, + it.rCol, + it.rAlt + )); + // Add std::cout if needed here + // std::cout << "iProc = " << iProc << " interpolates the point successfully\n"; + } else + ans.push_back(cNinf); + } + + return ans; +} + From 5cd781caf258f92c02bf432341d7eda69dc92183 Mon Sep 17 00:00:00 2001 From: Aaron Ridley Date: Thu, 27 Nov 2025 13:24:00 -0500 Subject: [PATCH 10/10] BUG: proper merge between the Aarons codes --- include/grid.h | 101 +++++++------- src/grid_match.cpp | 114 +++++++++------ src/solver_grid_interpolation.cpp | 225 +++++++++++++++++++++--------- 3 files changed, 282 insertions(+), 158 deletions(-) diff --git a/include/grid.h b/include/grid.h index a5f4c0b1..0f7a79c7 100644 --- a/include/grid.h +++ b/include/grid.h @@ -45,30 +45,31 @@ struct cubesphere_chars { // Grid class // ---------------------------------------------------------------------------- - struct interp_coef_t { - // The point is inside the cube of [iRow, iRow+1], [iCol, iCol+1], [iAlt, iAlt+1] - uint64_t iRow; - uint64_t iCol; - uint64_t iAlt; - // The coefficients along row, column and altitude - precision_t rRow; - precision_t rCol; - precision_t rAlt; - // Whether the point is within this grid or not - bool in_grid; - // If this is set to true: - bool above_grid, below_grid; - // do interpolation in lat and lon, but extrapolate in altitude - }; +struct interp_coef_t { + // The point is inside the cube of: + // [iRow, iRow+1], [iCol, iCol+1] [iAlt, iAlt+1] + uint64_t iRow; + uint64_t iCol; + uint64_t iAlt; + // The coefficients along row, column and altitude + precision_t rRow; + precision_t rCol; + precision_t rAlt; + // Whether the point is within this grid or not + bool in_grid; + // If this is set to true: + bool above_grid, below_grid; + // do interpolation in lat and lon, but extrapolate in altitude +}; - struct grid_to_grid_t { - int64_t iProcTo; - int64_t nPts; - int64_t nPtsReceive; - std::vector interpCoefs; - std::vector valueToSend; - std::vector valueToReceive; - }; +struct grid_to_grid_t { + int64_t iProcTo; + int64_t nPts; + int64_t nPtsReceive; + std::vector interpCoefs; + std::vector valueToSend; + std::vector valueToReceive; +}; class Grid { public: @@ -552,12 +553,12 @@ class Grid { * \pre Lons, Lats and Alts have the same size * \return list of interpolation coefficients */ - + std::vector get_interpolation_coefs( - const std::vector &Lons, - const std::vector &Lats, - const std::vector &Alts); - + const std::vector &Lons, + const std::vector &Lats, + const std::vector &Alts); + /** * \brief Set the interpolation coefficients for the dipole grid * \param Lons The longitude of points @@ -666,18 +667,18 @@ class Grid { // Each point is processed by the function set_interpolation_coefs and stored // in the form of this structure. // If the point is out of the grid, in_grid = false and all other members are undefined - struct interp_coef_t { - // The point is inside the cube of [iRow, iRow+1], [iCol, iCol+1], [iAlt, iAlt+1] - uint64_t iRow; - uint64_t iCol; - uint64_t iAlt; - // The coefficients along row, column and altitude - precision_t rRow; - precision_t rCol; - precision_t rAlt; - // Whether the point is within this grid or not - bool in_grid; - }; + //struct interp_coef_t { + // // The point is inside the cube of [iRow, iRow+1], [iCol, iCol+1], [iAlt, iAlt+1] + // uint64_t iRow; + // uint64_t iCol; + // uint64_t iAlt; + // // The coefficients along row, column and altitude + // precision_t rRow; + // precision_t rCol; + // precision_t rAlt; + // // Whether the point is within this grid or not + // bool in_grid; + //}; // Calculate the range of a spherical grid void get_sphere_grid_range(struct sphere_range &sr) const; @@ -688,18 +689,18 @@ class Grid { // Helper function for set_interpolation_coefs struct interp_coef_t get_interp_coef_sphere(const sphere_range &sr, - const precision_t lon_in, - const precision_t lat_in, - const precision_t alt_in); + const precision_t lon_in, + const precision_t lat_in, + const precision_t alt_in); struct interp_coef_t get_interp_coef_cubesphere(const cubesphere_range &cr, - const precision_t lon_in, - const precision_t lat_in, - const precision_t alt_in); + const precision_t lon_in, + const precision_t lat_in, + const precision_t alt_in); // (note these are magnetic coordinates) - void set_interp_coef_dipole(const dipole_range &dr, - const precision_t lon_in, - const precision_t lat_in, - const precision_t alt_in); + struct interp_coef_t get_interp_coef_dipole(const dipole_range &dr, + const precision_t lon_in, + const precision_t lat_in, + const precision_t alt_in); // Processed interpolation coefficients std::vector interp_coefs; diff --git a/src/grid_match.cpp b/src/grid_match.cpp index 908baf8b..28d1c2b1 100644 --- a/src/grid_match.cpp +++ b/src/grid_match.cpp @@ -20,9 +20,8 @@ bool exchange_information(int64_t *nPointsToPass, // - if it is a different processor, send the data for (jNode = 0; jNode < nGrids ; jNode++) { if (jNode == iGrid) { - for (iPt = 0; iPt < nPointsToPass[jNode]; iPt ++) { + for (iPt = 0; iPt < nPointsToPass[jNode]; iPt ++) varToReceive[jNode][iPt] = varToSend[jNode][iPt]; - } } else { iProcTo = iMember * nGrids + jNode; // iTag is a unique id allowing all processors to @@ -74,8 +73,8 @@ bool exchange_information(int64_t *nPointsToPass, // - builds interpolators for the requested information // ----------------------------------------------------------------------------- -bool grid_match(Grid &gGrid, - Grid &mGrid, +bool grid_match(Grid &gGrid, + Grid &mGrid, Quadtree gQuadtree, Quadtree mQuadtree) { @@ -93,8 +92,10 @@ bool grid_match(Grid &gGrid, precision_t normX, normY, normZ; arma_vec norms(3); int64_t jNode, kNode; - int64_t *nPointsToPass = static_cast(malloc(nGrids * sizeof(int64_t))); - int64_t *nPointsToReceive = static_cast(malloc(nGrids * sizeof(int64_t))); + int64_t *nPointsToPass = static_cast(malloc(nGrids * sizeof( + int64_t))); + int64_t *nPointsToReceive = static_cast(malloc(nGrids * sizeof( + int64_t))); int64_t *nPointsDummy = static_cast(malloc(nGrids * sizeof(int64_t))); for (jNode = 0; jNode < nGrids ; jNode++) @@ -103,6 +104,7 @@ bool grid_match(Grid &gGrid, // This is not the most efficient way to do this, but the first pass, let's // just count how many points we need to send to the other processors: mGrid.gridToGridMap.set_size(mnX, mnY, mnZ); + for (iX = mGCs; iX < mnX - mGCs; iX++) { for (iY = mGCs; iY < mnY - mGCs; iY++) { for (iZ = mGCs; iZ < mnZ - mGCs; iZ++) { @@ -118,11 +120,12 @@ bool grid_match(Grid &gGrid, norms = sphere_to_cube(lon, lat); jNode = gQuadtree.find_point(norms); } - if (jNode < 0 || jNode >= nGrids) { - std::cout << "out of bounds!!! " << jNode << "\n"; - } + + if (jNode < 0 || jNode >= nGrids) + std::cout << "out of bounds!!! " << jNode << "\n"; + mGrid.gridToGridMap(iX, iY, iZ) = jNode; - nPointsToPass[jNode] = nPointsToPass[jNode]+1; + nPointsToPass[jNode] = nPointsToPass[jNode] + 1; /* std::cout << "lon, lat, node: " << lon*cRtoD << " " << lat*cRtoD << " " << norms(0) << " " @@ -134,11 +137,13 @@ bool grid_match(Grid &gGrid, } } } + MPI_Barrier(aether_comm); if (report.test_verbose(3)) { for (jNode = 0; jNode < nGrids ; jNode++) std::cout << "nPtsToPass : " << iProc << " " << nPointsToPass[jNode] << "\n"; + std::cout << "sending number of points :\n"; } @@ -150,14 +155,15 @@ bool grid_match(Grid &gGrid, for (kNode = 0; kNode < nGrids ; kNode++) nPointsDummy[kNode] = nPointsToPass[kNode]; } + MPI_Bcast(nPointsDummy, nGrids, MPI_INT64_T, jNode, aether_comm); nPointsToReceive[jNode] = nPointsDummy[iGrid]; } if (report.test_verbose(3)) { - for (jNode = 0; jNode < nGrids ; jNode++) { - std::cout << "nPtsToReceive : " << iProc << " " << jNode << " " << nPointsToReceive[jNode] << "\n"; - } + for (jNode = 0; jNode < nGrids ; jNode++) + std::cout << "nPtsToReceive : " << iProc << " " << jNode << " " << + nPointsToReceive[jNode] << "\n"; } // Now we need to create an array of send points and an array of receive points. @@ -167,24 +173,33 @@ bool grid_match(Grid &gGrid, std::vector latsToInterTo(nGrids); std::vector lonsToInterTo(nGrids); std::vector altsToInterTo(nGrids); + for (jNode = 0; jNode < nGrids ; jNode++) { - latsToPass[jNode] = static_cast(malloc(nPointsToPass[jNode] * sizeof(precision_t))); - lonsToPass[jNode] = static_cast(malloc(nPointsToPass[jNode] * sizeof(precision_t))); - altsToPass[jNode] = static_cast(malloc(nPointsToPass[jNode] * sizeof(precision_t))); - latsToInterTo[jNode] = static_cast(malloc(nPointsToReceive[jNode] * sizeof(precision_t))); - lonsToInterTo[jNode] = static_cast(malloc(nPointsToReceive[jNode] * sizeof(precision_t))); - altsToInterTo[jNode] = static_cast(malloc(nPointsToReceive[jNode] * sizeof(precision_t))); + latsToPass[jNode] = static_cast(malloc(nPointsToPass[jNode] * + sizeof(precision_t))); + lonsToPass[jNode] = static_cast(malloc(nPointsToPass[jNode] * + sizeof(precision_t))); + altsToPass[jNode] = static_cast(malloc(nPointsToPass[jNode] * + sizeof(precision_t))); + latsToInterTo[jNode] = static_cast(malloc( + nPointsToReceive[jNode] * sizeof(precision_t))); + lonsToInterTo[jNode] = static_cast(malloc( + nPointsToReceive[jNode] * sizeof(precision_t))); + altsToInterTo[jNode] = static_cast(malloc( + nPointsToReceive[jNode] * sizeof(precision_t))); } // now, the second pass, let's store the information so we can pass it: for (jNode = 0; jNode < nGrids ; jNode++) nPointsToPass[jNode] = 0; + for (iX = mGCs; iX < mnX - mGCs; iX++) { for (iY = mGCs; iY < mnY - mGCs; iY++) { for (iZ = mGCs; iZ < mnZ - mGCs; iZ++) { lon = mGrid.geoLon_scgc(iX, iY, iZ); lat = mGrid.geoLat_scgc(iX, iY, iZ); - if (gGrid.iGridShape_ == gGrid.iSphere_) { + + if (gGrid.iGridShape_ == iSphere_) { norms(0) = lon / cPI; norms(1) = lat / cPI; norms(2) = 0.0; @@ -193,13 +208,15 @@ bool grid_match(Grid &gGrid, norms = sphere_to_cube(lon, lat); jNode = gQuadtree.find_point(norms); } + latsToPass[jNode][nPointsToPass[jNode]] = lat; lonsToPass[jNode][nPointsToPass[jNode]] = lon; altsToPass[jNode][nPointsToPass[jNode]] = mGrid.geoAlt_scgc(iX, iY, iZ); - nPointsToPass[jNode] = nPointsToPass[jNode]+1; + nPointsToPass[jNode] = nPointsToPass[jNode] + 1; } } } + bool didWork; // Pass first coordinate (lons) didWork = exchange_information(nPointsToPass, @@ -221,38 +238,45 @@ bool grid_match(Grid &gGrid, for (jNode = 0; jNode < nGrids ; jNode++) { std::cout << "Received the following points from iGrid = " << jNode << "\n"; std::cout << " -> points received : " << nPointsToReceive[jNode] << "\n"; - for (int64_t iPt = 0; iPt < nPointsToReceive[jNode]; iPt++) + + for (int64_t iPt = 0; iPt < nPointsToReceive[jNode]; iPt++) std::cout << " -> " << iPt << " " - << lonsToInterTo[jNode][iPt] << " " - << latsToInterTo[jNode][iPt] << " " - << altsToInterTo[jNode][iPt] << "\n"; + << lonsToInterTo[jNode][iPt] << " " + << latsToInterTo[jNode][iPt] << " " + << altsToInterTo[jNode][iPt] << "\n"; } } struct grid_to_grid_t oneGrid; int64_t nPts; + for (jNode = 0; jNode < nGrids ; jNode++) { // These are backwards now, since we will switch sender and reciever: oneGrid.nPts = nPointsToReceive[jNode]; oneGrid.nPtsReceive = nPointsToPass[jNode]; oneGrid.iProcTo = iMember * nGrids + jNode; + if (report.test_verbose(2)) - std::cout << "Making interpolation coefficients for : " << jNode - << "; points : " << oneGrid.nPts << "\n"; + std::cout << "Making interpolation coefficients for : " << jNode + << "; points : " << oneGrid.nPts << "\n"; + if (oneGrid.nPts > 0) { - // Interpolation function takes vectors, + // Interpolation function takes vectors, // so transfer these arrays to vectors: std::vector Lons(oneGrid.nPts); std::vector Lats(oneGrid.nPts); std::vector Alts(oneGrid.nPts); + for (int64_t iPt = 0; iPt < oneGrid.nPts; iPt++) { Lons[iPt] = lonsToInterTo[jNode][iPt]; Lats[iPt] = latsToInterTo[jNode][iPt]; Alts[iPt] = altsToInterTo[jNode][iPt]; } + oneGrid.interpCoefs = gGrid.get_interpolation_coefs(Lons, Lats, Alts); } + gGrid.gridToGridCoefs.push_back(oneGrid); } @@ -260,7 +284,7 @@ bool grid_match(Grid &gGrid, return didWork; } -bool get_data_from_other_grid(Grid &gGrid, +bool get_data_from_other_grid(Grid &gGrid, Grid &mGrid, arma_cube &gData, arma_cube &mData) { @@ -272,24 +296,33 @@ bool get_data_from_other_grid(Grid &gGrid, int64_t jNode, iPt; std::vector dataToSend(nGrids); std::vector dataToReceive(nGrids); - int64_t *nPointsToSend = static_cast(malloc(nGrids * sizeof(int64_t))); - int64_t *nPointsToReceive = static_cast(malloc(nGrids * sizeof(int64_t))); + int64_t *nPointsToSend = static_cast(malloc(nGrids * sizeof( + int64_t))); + int64_t *nPointsToReceive = static_cast(malloc(nGrids * sizeof( + int64_t))); for (jNode = 0; jNode < nGrids ; jNode++) { if (report.test_verbose(2)) - std::cout << "nPts : " << jNode << " " << gGrid.gridToGridCoefs[jNode].nPts << "\n"; + std::cout << "nPts : " << jNode << " " << gGrid.gridToGridCoefs[jNode].nPts << + "\n"; + nPointsToSend[jNode] = gGrid.gridToGridCoefs[jNode].nPts; nPointsToReceive[jNode] = gGrid.gridToGridCoefs[jNode].nPtsReceive; - dataToSend[jNode] = static_cast(malloc(gGrid.gridToGridCoefs[jNode].nPts * sizeof(precision_t))); - dataToReceive[jNode] = static_cast(malloc(gGrid.gridToGridCoefs[jNode].nPtsReceive * sizeof(precision_t))); - std::vector values = gGrid.get_interpolation_values(gData, gGrid.gridToGridCoefs[jNode].interpCoefs); - + dataToSend[jNode] = static_cast(malloc( + gGrid.gridToGridCoefs[jNode].nPts * sizeof(precision_t))); + dataToReceive[jNode] = static_cast(malloc( + gGrid.gridToGridCoefs[jNode].nPtsReceive * sizeof(precision_t))); + std::vector values = gGrid.get_interpolation_values(gData, + gGrid.gridToGridCoefs[jNode].interpCoefs); + for (iPt = 0; iPt < gGrid.gridToGridCoefs[jNode].nPts; iPt++) { dataToSend[jNode][iPt] = values[iPt]; + if (report.test_verbose(2)) std::cout << "datatosend : " << iPt << " " << dataToSend[jNode][iPt] << "\n"; } } + bool didWork = exchange_information(nPointsToSend, dataToSend, nPointsToReceive, @@ -299,6 +332,7 @@ bool get_data_from_other_grid(Grid &gGrid, int64_t iZ, mnZ = mGrid.get_nZ(); int64_t mGCs = mGrid.get_nGCs(); std::vector iCounter(nGrids); + for (jNode = 0; jNode < nGrids ; jNode++) iCounter[jNode] = 0; @@ -306,13 +340,15 @@ bool get_data_from_other_grid(Grid &gGrid, for (iY = mGCs; iY < mnY - mGCs; iY++) { for (iZ = mGCs; iZ < mnZ - mGCs; iZ++) { jNode = mGrid.gridToGridMap(iX, iY, iZ); + if (report.test_verbose(2)) { - std::cout << "unpacking point : " << iX << " " << iY << " " << iZ << " " << jNode << " " - << iCounter[jNode] << " " << dataToReceive[jNode][iCounter[jNode]] << "\n"; + std::cout << "unpacking point : " << iX << " " << iY << " " << iZ << " " << + jNode << " " + << iCounter[jNode] << " " << dataToReceive[jNode][iCounter[jNode]] << "\n"; } mData(iX, iY, iZ) = dataToReceive[jNode][iCounter[jNode]]; - iCounter[jNode] = iCounter[jNode]+1; + iCounter[jNode] = iCounter[jNode] + 1; } } } diff --git a/src/solver_grid_interpolation.cpp b/src/solver_grid_interpolation.cpp index 0ec38266..e0cba3ca 100644 --- a/src/solver_grid_interpolation.cpp +++ b/src/solver_grid_interpolation.cpp @@ -3,7 +3,9 @@ #include "aether.h" -// Hepler varialbes / function begins. These are only used inside this cpp file and neither declared nor visible in any other file +// Hepler variables / function begins. +// These are only used inside this cpp file and neither declared +// nor visible in any other file // The size of a 2*2*2 arma cube const arma::SizeCube unit_cube_size = arma::size(2, 2, 2); @@ -69,17 +71,48 @@ int64_t get_cube_surface_number(const arma_vec &point_in) { } } -// Helper variables / function ends. The following are all member functions of Grid class +// Helper variables / function ends. The following are all member +// functions of Grid class // -------------------------------------------------------------------------- -// Return the index of the last element that has a value smaller than or equal to the input +// Return the index of the last element that has altitude smaller than +// or equal to the input +// -------------------------------------------------------------------------- + +uint64_t Grid::search_altitude(const precision_t alt_in) const { + // Copy from std::upper_bound. Can't directly use it + // mainly because geoAlt_scgc(0, 0, *) can't be formed as an iterator + uint64_t first, last, len; + first = nGCs; + last = nAlts - nGCs; + len = last - first; + + while (len > 0) { + uint64_t half = len >> 1; + uint64_t mid = first + half; + + if (geoAlt_scgc(0, 0, mid) > alt_in) + len = half; + + else { + first = mid + 1; + len = len - half - 1; + } + } + + return first - 1; +} + +// -------------------------------------------------------------------------- +// Return the index of the last element that has a value smaller than +// or equal to the input // - Optional argument (nGCs=0) since we cannot see grid info. // -------------------------------------------------------------------------- +// this replaces the above + uint64_t bisect_search_array(precision_t val_in, arma_vec ref_arr, int64_t nGCs = 0) { - // Copy from std::upper_bound. Can't directly use it - // mainly because geoAlt_scgc(0, 0, *) can't be formed as an iterator uint64_t first, last, len; first = nGCs; last = ref_arr.size(); @@ -181,6 +214,7 @@ void Grid::get_cubesphere_grid_range(struct cubesphere_range &cr) const { } } + // -------------------------------------------------------------------------- // Get the range of a Dipole grid // -------------------------------------------------------------------------- @@ -208,24 +242,27 @@ void Grid::get_dipole_grid_range(struct dipole_range &dr) const { // Almost the copy of interp_sphere_linear_helper // -------------------------------------------------------------------------- -struct interp_coef_t Grid::get_interp_coef_sphere( - const sphere_range &sr, - const precision_t lon_in, - const precision_t lat_in, - const precision_t alt_in) { + +struct interp_coef_t Grid::get_interp_coef_sphere(const sphere_range &sr, + const precision_t lon_in, + const precision_t lat_in, + const precision_t alt_in) { + // WARNING: IF WE ARE DEALING WITH LESS THAN THE WHOLE EARTH, THEN ALL THE POINTS WITH // LONGITUDE = geo_grid_input.lon_max = settings["GeoGrid"]["MaxLon"] // OR LATITUDE = geo_grid_input.lat_max = settings["GeoGrid"]["MaxLat"] // ARE EXCLUDED. // TO FIX IT, EACH GRID SHOULD BE ABLE TO ACCESS THE MaxLon and MaxLat - // The structure which will be put into the interp_coefs. Initialize in_grid to be false + // The structure which will be put into the interp_coefs. + // Initialize in_grid to be false struct interp_coef_t coef; coef.in_grid = false; // Determine whether the point is inside this grid - // Treat north pole specially because latitude is inclusive for both -cPI/2 and cPI/2 - // Dpn't check for altitude here! + // Treat north pole specially because latitude is inclusive for + // both -cPI/2 and cPI/2 + // Don't check for altitude here! if (lon_in < sr.lon_min || lon_in >= sr.lon_max || lat_in < sr.lat_min || lat_in > sr.lat_max || (lat_in == sr.lat_max && sr.lat_max != cPI / 2)) { return coef; @@ -268,9 +305,13 @@ struct interp_coef_t Grid::get_interp_coef_sphere( coef.below_grid = false; coef.above_grid = true; } else { - coef.iAlt = search_altitude(alt_in); - coef.rAlt = (alt_in - geoAlt_scgc(0, 0, coef.iAlt)) - / (geoAlt_scgc(0, 0, coef.iAlt + 1) - geoAlt_scgc(0, 0, coef.iAlt)); + coef.iAlt = bisect_search_array(alt_in, + geoAlt_scgc.tube(coef.iRow, coef.iCol), + nGCs); + coef.rAlt = + (alt_in - geoAlt_scgc(coef.iRow, coef.iCol, coef.iAlt)) + / (geoAlt_scgc(coef.iRow, coef.iCol, coef.iAlt + 1) - + geoAlt_scgc(coef.iRow, coef.iCol, coef.iAlt)); coef.below_grid = false; coef.above_grid = false; } @@ -343,36 +384,63 @@ struct interp_coef_t Grid::get_interp_coef_cubesphere(const cubesphere_range &cr coef.iCol = static_cast(col_frac_index); coef.rCol = col_frac_index - coef.iCol; coef.iCol += nGCs - 1; - // Use binary search to find the index for altitude (handles oblate planets) - coef.iAlt = bisect_search_array(alt_in, geoAlt_scgc.tube(coef.iRow, coef.iCol), - nGCs); - coef.rAlt = (alt_in - geoAlt_scgc(coef.iRow, coef.iCol, coef.iAlt)) - / (geoAlt_scgc(coef.iRow, coef.iCol, coef.iAlt + 1) - geoAlt_scgc(coef.iRow, - coef.iCol, coef.iAlt)); - // Put the coefficient into the vector - coef.in_grid = true; - interp_coefs.push_back(coef); + + // The altitude may not be linearly spaced, so use binary search to find + // the first element smaller than or equal to the altitude of the give point + // Implemented in search_altitude + + if (alt_in < cr.alt_min) { + coef.iAlt = nGCs; + coef.rAlt = alt_in - cr.alt_min; + coef.below_grid = true; + coef.above_grid = false; + } else { + if (alt_in > cr.alt_max) { + coef.iAlt = nAlts - nGCs; + coef.rAlt = alt_in - cr.alt_max; + coef.below_grid = false; + coef.above_grid = true; + } else { + coef.iAlt = bisect_search_array(alt_in, + geoAlt_scgc.tube(coef.iRow, coef.iCol), + nGCs); + coef.rAlt = + (alt_in - geoAlt_scgc(coef.iRow, coef.iCol, coef.iAlt)) + / (geoAlt_scgc(coef.iRow, coef.iCol, coef.iAlt + 1) - + geoAlt_scgc(coef.iRow, coef.iCol, coef.iAlt)); + coef.below_grid = false; + coef.above_grid = false; + } + } + return coef; } -void Grid::set_interp_coef_dipole(const dipole_range &dr, - const precision_t lon_in, - const precision_t lat_in, - const precision_t alt_in) { - // The structure which will be put into the interp_coefs. Initialize in_grid to be false +struct interp_coef_t Grid::get_interp_coef_dipole(const dipole_range &dr, + const precision_t lon_in, + const precision_t lat_in, + const precision_t alt_in) { + + // The structure which will be put into the interp_coefs. Initialize + // in_grid to be false struct interp_coef_t coef; coef.in_grid = false; - // Determine whether the point is inside this grid - // Treat north pole specially because latitude is inclusive for both -cPI/2 and cPI/2 - if (lon_in < dr.lon_min || lon_in >= dr.lon_max || lat_in < dr.lat_min - || lat_in > dr.lat_max || (lat_in == dr.lat_max && dr.lat_max != cPI / 2) + // Determine whether the point is inside this grid Treat north pole + // specially because latitude is inclusive for both -cPI/2 and cPI/2 + if (lon_in < dr.lon_min || + lon_in >= dr.lon_max || + lat_in < dr.lat_min || + lat_in > dr.lat_max || + (lat_in == dr.lat_max && dr.lat_max != cPI / 2) || alt_in < dr.alt_min || alt_in > dr.alt_max) { - interp_coefs.push_back(coef); - return; + return coef; } + // Put the coefficient into the vector + coef.in_grid = true; + // ASSUMPTION: LONGITUDE IS LINEARLY SPACED, nGCs >= 1 // For the cell containing it, directly calculate its x index // Find y & z indices using a bisecting search @@ -383,10 +451,10 @@ void Grid::set_interp_coef_dipole(const dipole_range &dr, coef.iRow = static_cast(coef.rRow); // Calculate the fractional part, which is the ratio for Longitude coef.rRow -= coef.iRow; - // The actual x-axis index of the bottom-left of the cube used for interpolation + // The actual x-axis index of the bottom-left of the cube used for + // interpolation coef.iRow += nGCs - 1; - // Different from the sphere, latitude & altitude are not evenly spaced. // Use the bisect search function for both. @@ -395,42 +463,44 @@ void Grid::set_interp_coef_dipole(const dipole_range &dr, coef.iCol = bisect_search_array(abs(lat_in), abs(j_center_scgc.tube(coef.iRow, coef.iCol)), nGCs); - // need alt index to find lat coef - coef.iAlt = bisect_search_array(alt_in, k_center_scgc.tube(coef.iRow, - coef.iCol), - nGCs); - - // then we can do the ratios: - coef.rCol = (lat_in - magLat_scgc(coef.iRow, coef.iCol, coef.iAlt)) - / (magLat_scgc(coef.iRow, coef.iCol + 1, coef.iAlt) - - magLat_scgc(coef.iRow, coef.iCol, coef.iAlt)); - - coef.rAlt = (alt_in - geoAlt_scgc(coef.iRow, coef.iCol, coef.iAlt)) - / (geoAlt_scgc(coef.iRow, coef.iCol, coef.iAlt + 1) - - geoAlt_scgc(coef.iRow, coef.iCol, coef.iAlt)); - - if (alt_in < cr.alt_min) { + // Use binary search to find the index for altitude + if (alt_in < dr.alt_min) { coef.iAlt = nGCs; coef.rAlt = 0.0; coef.below_grid = true; coef.above_grid = false; } else { - if (alt_in > cr.alt_max) { + if (alt_in > dr.alt_max) { coef.iAlt = nAlts - nGCs; coef.rAlt = 0.0; coef.below_grid = false; coef.above_grid = true; } else { - coef.iAlt = search_altitude(alt_in); - coef.rAlt = (alt_in - geoAlt_scgc(0, 0, coef.iAlt)) - / (geoAlt_scgc(0, 0, coef.iAlt + 1) - geoAlt_scgc(0, 0, coef.iAlt)); + // Use binary search to find the index for altitude (handles + // oblate planets) + + // need alt index to find lat coef + coef.iAlt = bisect_search_array(alt_in, + k_center_scgc.tube(coef.iRow, coef.iCol), + nGCs); + // then we can do the ratios: + coef.rCol = + (lat_in - magLat_scgc(coef.iRow, coef.iCol, coef.iAlt)) + / (magLat_scgc(coef.iRow, coef.iCol + 1, coef.iAlt) + - magLat_scgc(coef.iRow, coef.iCol, coef.iAlt)); + coef.rAlt = + (alt_in - geoAlt_scgc(coef.iRow, coef.iCol, coef.iAlt)) + / (geoAlt_scgc(coef.iRow, coef.iCol, coef.iAlt + 1) - + geoAlt_scgc(coef.iRow, coef.iCol, coef.iAlt)); coef.below_grid = false; coef.above_grid = false; } } + return coef; } + // -------------------------------------------------------------------------- // Set the interpolation coefficients // -------------------------------------------------------------------------- @@ -458,9 +528,12 @@ bool Grid::set_interpolation_coefs(const std::vector &i_coords, report.enter(function, iFunction); report.print(1, "interpolation gridtype : " + gridType); + + struct interp_coef_t coef; // If the size of Lons, Lats and Alts are not the same, return false - if (i_coords.size() != j_coords.size() || j_coords.size() != k_coords.size()) { + if (i_coords.size() != j_coords.size() || + j_coords.size() != k_coords.size()) { report.error("Length of i,j,k vectors do not match!"); return false; } @@ -468,16 +541,22 @@ bool Grid::set_interpolation_coefs(const std::vector &i_coords, // Clear the previous interpolation coefficients interp_coefs.clear(); - if (iGridShape_ == iCubesphere_) { - report.print(1, "interpolation grid is cubesphere"); - + // --------------------------------------------------- + // Cubesphere + if (IsCubeSphereGrid) { // Calculate the range of the grid struct cubesphere_range cr; get_cubesphere_grid_range(cr); // Calculate the index and coefficients for each point - for (size_t i = 0; i < i_coords.size(); ++i) - set_interp_coef_cubesphere(cr, i_coords[i], j_coords[i], k_coords[i]); + for (size_t i = 0; i < i_coords.size(); ++i) { + coef = get_interp_coef_cubesphere(cr, + i_coords[i], + j_coords[i], + k_coords[i]); + interp_coefs.push_back(coef); + } + } if (iGridShape_ == iSphere_) { @@ -487,8 +566,13 @@ bool Grid::set_interpolation_coefs(const std::vector &i_coords, get_sphere_grid_range(sr); // Calculate the index and coefficients for each point - for (size_t i = 0; i < i_coords.size(); ++i) - set_interp_coef_sphere(sr, i_coords[i], j_coords[i], k_coords[i]); + for (size_t i = 0; i < i_coords.size(); ++i) { + coef = get_interp_coef_sphere(sr, + i_coords[i], + j_coords[i], + k_coords[i]); + interp_coefs.push_back(coef); + } } if (iGridShape_ == iDipole_) { // IsDipole @@ -504,8 +588,9 @@ bool Grid::set_interpolation_coefs(const std::vector &i_coords, int64_t iLoc, nPts = i_coords.size(); std::vector mlon(nPts), p_coord(nPts), q_coord(nPts), dipijk(3); - // these are the magnetic coordinates. A temporary step! - // this is a vector of cubes with shape (nPts, 1, 1) - avoids having to overload things + // these are the magnetic coordinates. A temporary step! this is + // a vector of cubes with shape (nPts, 1, 1) - avoids having to + // overload things std::vector magCoords; if (areLocsGeo) { @@ -561,8 +646,10 @@ bool Grid::set_interpolation_coefs(const std::vector &i_coords, } // Calculate the index and coefficients for each point - for (size_t i = 0; i < i_coords.size(); ++i) - set_interp_coef_dipole(dr, mlon[i], p_coord[i], q_coord[i]); + for (size_t i = 0; i < i_coords.size(); ++i) { + coef = get_interp_coef_dipole(dr, mlon[i], p_coord[i], q_coord[i]); + interp_coefs.push_back(coef); + } } report.exit(function);