From 0aa58537b112be9b133816a8215bca563ef60560 Mon Sep 17 00:00:00 2001 From: robin Date: Tue, 30 Apr 2024 12:34:10 +0200 Subject: [PATCH] fix: unit conversion and indexing --- advection/src/AdvectionKernel.h | 7 ++- advection/src/RK4AdvectionKernel.cpp | 16 +++--- advection/src/RK4AdvectionKernel.h | 1 - advection/src/UVGrid.cpp | 4 +- advection/src/UVGrid.h | 2 +- advection/src/main.cpp | 74 ++++++++++++++++++++++------ 6 files changed, 77 insertions(+), 27 deletions(-) diff --git a/advection/src/AdvectionKernel.h b/advection/src/AdvectionKernel.h index bf0f094..f89caca 100644 --- a/advection/src/AdvectionKernel.h +++ b/advection/src/AdvectionKernel.h @@ -11,7 +11,7 @@ */ class AdvectionKernel { public: - const static int DT = 50; + const static int DT = 100000; // Seconds /** * This function must take a time, latitude and longitude of a particle and must output * a new latitude and longitude after being advected once for AdvectionKernel::DT time as defined above. @@ -22,6 +22,11 @@ public: */ virtual std::pair advect(int time, double latitude, double longitude) const = 0; + // Taken from Parcels https://github.com/OceanParcels/parcels/blob/daa4b062ed8ae0b2be3d87367d6b45599d6f95db/parcels/tools/converters.py#L155 + const static double metreToDegrees(double metre) { + return metre / 1000. / 1.852 / 60.; + } + }; #endif //ADVECTION_ADVECTIONKERNEL_H diff --git a/advection/src/RK4AdvectionKernel.cpp b/advection/src/RK4AdvectionKernel.cpp index d9eb928..3aa1006 100644 --- a/advection/src/RK4AdvectionKernel.cpp +++ b/advection/src/RK4AdvectionKernel.cpp @@ -8,28 +8,28 @@ RK4AdvectionKernel::RK4AdvectionKernel(std::shared_ptr grid): grid(grid) std::pair RK4AdvectionKernel::advect(int time, double latitude, double longitude) const { auto [u1, v1] = bilinearinterpolation(*grid, time, latitude, longitude); // lon1, lat1 = (particle.lon + u1*.5*particle.dt, particle.lat + v1*.5*particle.dt); - double lon1 = longitude + v1 * .5*DT; - double lat1 = latitude + u1*.5*DT; + double lon1 = longitude + metreToDegrees(v1 * 0.5*DT); + double lat1 = latitude + metreToDegrees(u1 * 0.5*DT); // (u2, v2) = fieldset.UV[time + .5 * particle.dt, particle.depth, lat1, lon1, particle] auto [u2, v2] = bilinearinterpolation(*grid, time + 0.5*DT, lat1, lon1); // lon2, lat2 = (particle.lon + u2*.5*particle.dt, particle.lat + v2*.5*particle.dt) - double lon2 = longitude + v2 * 0.5 * DT; - double lat2 = latitude + u2 * 0.5 * DT; + double lon2 = longitude + metreToDegrees(v2 * 0.5 * DT); + double lat2 = latitude + metreToDegrees(u2 * 0.5 * DT); // (u3, v3) = fieldset.UV[time + .5 * particle.dt, particle.depth, lat2, lon2, particle] auto [u3, v3] = bilinearinterpolation(*grid, time + 0.5 * DT, lat2, lon2); // lon3, lat3 = (particle.lon + u3*particle.dt, particle.lat + v3*particle.dt) - double lon3 = longitude + v3 * DT; - double lat3 = latitude + u3 * DT; + double lon3 = longitude + metreToDegrees(v3 * DT); + double lat3 = latitude + metreToDegrees(u3 * DT); // (u4, v4) = fieldset.UV[time + particle.dt, particle.depth, lat3, lon3, particle] auto [u4, v4] = bilinearinterpolation(*grid, time + DT, lat3, lon3); - double lonFinal = longitude + (v1 + 2 * v2 + 2 * v3 + v4) / 6.0 * DT; - double latFinal = latitude + (u1 + 2 * u2 + 2 * u3 + u4) / 6.0 * DT; + double lonFinal = longitude + metreToDegrees((v1 + 2 * v2 + 2 * v3 + v4) / 6.0 * DT); + double latFinal = latitude + metreToDegrees((u1 + 2 * u2 + 2 * u3 + u4) / 6.0 * DT); return {latFinal, lonFinal}; } diff --git a/advection/src/RK4AdvectionKernel.h b/advection/src/RK4AdvectionKernel.h index 6719030..6b6c88d 100644 --- a/advection/src/RK4AdvectionKernel.h +++ b/advection/src/RK4AdvectionKernel.h @@ -1,7 +1,6 @@ #ifndef ADVECTION_RK4ADVECTIONKERNEL_H #define ADVECTION_RK4ADVECTIONKERNEL_H - #include "AdvectionKernel.h" #include "UVGrid.h" diff --git a/advection/src/UVGrid.cpp b/advection/src/UVGrid.cpp index 63f9079..1febbfe 100644 --- a/advection/src/UVGrid.cpp +++ b/advection/src/UVGrid.cpp @@ -39,7 +39,7 @@ const Vel &UVGrid::operator[](size_t timeIndex, size_t latIndex, size_t lonIndex or lonIndex < 0 or lonIndex >= lonSize) { throw std::out_of_range("Index out of bounds"); } - size_t index = timeIndex * (latSize * lonSize) + latIndex * lonIndex + lonIndex; + size_t index = timeIndex * (latSize * lonSize) + latIndex * lonSize + lonIndex; return uvData[index]; } @@ -63,4 +63,4 @@ void UVGrid::streamSlice(ostream &os, size_t t) { } os << endl; } -} +} \ No newline at end of file diff --git a/advection/src/UVGrid.h b/advection/src/UVGrid.h index 691d0df..f068ea9 100644 --- a/advection/src/UVGrid.h +++ b/advection/src/UVGrid.h @@ -49,7 +49,7 @@ public: std::vector lons; /** - * The 3D index into the data + * The 3D index into the data. The array is sized by [8761][67][116] * @return Velocity at that index */ const Vel& operator[](size_t timeIndex, size_t latIndex, size_t lonIndex) const; diff --git a/advection/src/main.cpp b/advection/src/main.cpp index 4827271..8312e2e 100644 --- a/advection/src/main.cpp +++ b/advection/src/main.cpp @@ -1,43 +1,89 @@ +#include +#include +#include + #include "interpolate.h" #include "Vel.h" #include "EulerAdvectionKernel.h" #include "RK4AdvectionKernel.h" -#include -#include +#include "interpolate.h" #define NotAKernelError "Template parameter T must derive from AdvectionKernel" using namespace std; template -void advectForSomeTime(const UVGrid &uvGrid, const AdvectionKernelImpl &kernel, double latstart, double lonstart) { +void advectForSomeTime(const UVGrid &uvGrid, const AdvectionKernelImpl &kernel, double latstart, double lonstart, int i, char colour[10]) { // Require at compile time that kernel derives from the abstract class AdvectionKernel static_assert(std::is_base_of::value, NotAKernelError); double lat1 = latstart, lon1 = lonstart; - for(int time = 100; time <= 10000; time += AdvectionKernel::DT) { - cout << "lat = " << lat1 << " lon = " << lon1 << endl; - auto [templat, templon] = kernel.advect(time, lat1, lon1); - lat1 = templat; - lon1 = templon; + for(int time = 0; time <= 31536000.; time += AdvectionKernel::DT) { +// cout << setprecision(8) << lat1 << "," << setprecision(8) << lon1 << ",end" << i << "," << colour << endl; + try { + auto [templat, templon] = kernel.advect(time, lat1, lon1); + lat1 = templat; + lon1 = templon; + } catch (const out_of_range& e) { + cerr << "broke out of loop!" << endl; + break; + } } + cout << setprecision(8) << latstart << "," << setprecision(8) << lonstart << ",begin" << i << "," << colour << endl; + cout << setprecision(8) << lat1 << "," << setprecision(8) << lon1 << ",end" << i << "," << colour << endl; } +void testGridIndexing(const UVGrid *uvGrid) { + int time = 20000; + cout << "=== land === (should all give 0)" << endl; + cout << bilinearinterpolation(*uvGrid, time, 53.80956379699079, -1.6496306344654406) << endl; + cout << bilinearinterpolation(*uvGrid, time, 55.31428895563707, -2.851581041325997) << endl; + cout << bilinearinterpolation(*uvGrid, time, 47.71548983067583, -1.8704054037408626) << endl; + cout << bilinearinterpolation(*uvGrid, time, 56.23521060314398, 8.505979324950573) << endl; + cout << bilinearinterpolation(*uvGrid, time, 53.135645440244375, 8.505979324950573) << endl; + cout << bilinearinterpolation(*uvGrid, time, 56.44761278775708, -4.140629303756164) << endl; + cout << bilinearinterpolation(*uvGrid, time, 52.67625153110339, 0.8978569759455872) << endl; + cout << bilinearinterpolation(*uvGrid, time, 52.07154079279377, 4.627951041411331) << endl; + + cout << "=== ocean === (should give not 0)" << endl; + cout << bilinearinterpolation(*uvGrid, time, 47.43923166616274, -4.985451481829083) << endl; + cout << bilinearinterpolation(*uvGrid, time, 50.68943556852362, -9.306162999561733) << endl; + cout << bilinearinterpolation(*uvGrid, time, 53.70606799886677, -4.518347647034465) << endl; + cout << bilinearinterpolation(*uvGrid, time, 60.57987114267971, -12.208262973672621) << endl; + cout << bilinearinterpolation(*uvGrid, time, 46.532221548197285, -13.408189172582638) << endl; + cout << bilinearinterpolation(*uvGrid, time, 50.92725094937812, 1.3975824052375256) << endl; + cout << bilinearinterpolation(*uvGrid, time, 51.4028921682209, 2.4059571950925203) << endl; + cout << bilinearinterpolation(*uvGrid, time, 53.448445236769004, 0.7996966058017515) << endl; +// cout << bilinearinterpolation(*uvGrid, time, ) << endl; +} int main() { std::shared_ptr uvGrid = std::make_shared(); - EulerAdvectionKernel kernelEuler = EulerAdvectionKernel(uvGrid); + uvGrid->streamSlice(cout, 900); RK4AdvectionKernel kernelRK4 = RK4AdvectionKernel(uvGrid); - double latstart = 52.881770, lonstart = 3.079979; - - cout << "======= Euler Integration =======" << endl; - advectForSomeTime(*uvGrid, kernelEuler, latstart, lonstart); +// You can use https://maps.co/gis/ to visualise these points cout << "======= RK4 Integration =======" << endl; - advectForSomeTime(*uvGrid, kernelRK4, latstart, lonstart); + advectForSomeTime(*uvGrid, kernelRK4, 54.331795276466615, 4.845871408626273, 0, "#ADD8E6"); + advectForSomeTime(*uvGrid, kernelRK4, 59.17208978388813, 0.32865481669722213, 1, "#DC143C"); + advectForSomeTime(*uvGrid, kernelRK4, 56.18661322856813, -9.521463269751877, 2, "#50C878"); + advectForSomeTime(*uvGrid, kernelRK4, 46.6048774007515, -2.8605696406405947, 3, "#FFEA00"); + advectForSomeTime(*uvGrid, kernelRK4, 51.45431808648367, 1.6682437444140332, 4, "#663399"); + advectForSomeTime(*uvGrid, kernelRK4, 51.184757012016085, -6.418923014612084, 5, "#FFA500"); + advectForSomeTime(*uvGrid, kernelRK4, 54.48325546269538, 7.167517140551792, 6, "#008080"); + advectForSomeTime(*uvGrid, kernelRK4, 55.43253322410253, -1.1712503620884716, 7, "#FFB6C1"); + advectForSomeTime(*uvGrid, kernelRK4, 48.852815074172085, 3.4294489497130516, 8, "#36454F"); // on land + advectForSomeTime(*uvGrid, kernelRK4, 58.02431905976983, 1.6892571305388995, 9, "#1E90FF"); // Dodger Blue + advectForSomeTime(*uvGrid, kernelRK4, 58.99404571805297, 3.4121513161325425, 10, "#FFD700"); // Gold + advectForSomeTime(*uvGrid, kernelRK4, 59.51039243098509, -1.6674160241521663, 11, "#6A5ACD"); // Slate Blue + advectForSomeTime(*uvGrid, kernelRK4, 60.51250220636489, 1.020893006817227, 12, "#20B2AA"); // Light Sea Green + advectForSomeTime(*uvGrid, kernelRK4, 60.38797801281417, 3.516119068711471, 13, "#FF69B4"); // Hot Pink + advectForSomeTime(*uvGrid, kernelRK4, 60.02637651315464, -2.4546004365354697, 14, "#800080"); // Purple + advectForSomeTime(*uvGrid, kernelRK4, 58.732929454411305, 3.649791893455804, 15, "#FF4500"); // Orange Red +// advectForSomeTime(*uvGrid, kernelRK4, ,0); return 0; }