diff --git a/advection/README.md b/advection/README.md index 76cebec..0abcf14 100644 --- a/advection/README.md +++ b/advection/README.md @@ -25,6 +25,7 @@ void advectForSomeTime(const UVGrid &uvGrid, const AdvectionKernelImpl &kernel, ## Location of data The data path is hardcoded such that the following tree structure is assumed: +The current assumption is that the name of the `u`s and `v`s are flipped since this is the way the data was given to us. ``` data/ grid.h5 diff --git a/advection/src/AdvectionKernel.h b/advection/src/AdvectionKernel.h index f89caca..0170c79 100644 --- a/advection/src/AdvectionKernel.h +++ b/advection/src/AdvectionKernel.h @@ -11,7 +11,7 @@ */ class AdvectionKernel { public: - const static int DT = 100000; // Seconds + const static int DT = 60 * 15; // 60 sec/min * 15 mins /** * 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. diff --git a/advection/src/EulerAdvectionKernel.cpp b/advection/src/EulerAdvectionKernel.cpp index af3c436..01c76d9 100644 --- a/advection/src/EulerAdvectionKernel.cpp +++ b/advection/src/EulerAdvectionKernel.cpp @@ -9,5 +9,5 @@ EulerAdvectionKernel::EulerAdvectionKernel(std::shared_ptr grid): grid(g std::pair EulerAdvectionKernel::advect(int time, double latitude, double longitude) const { auto [u, v] = bilinearinterpolation(*grid, time, latitude, longitude); - return {latitude+u*DT, longitude+v*DT}; + return {latitude+metreToDegrees(u*DT), longitude+metreToDegrees(v*DT)}; } diff --git a/advection/src/EulerAdvectionKernel.h b/advection/src/EulerAdvectionKernel.h index e98297f..d9893cb 100644 --- a/advection/src/EulerAdvectionKernel.h +++ b/advection/src/EulerAdvectionKernel.h @@ -7,8 +7,8 @@ /** * Implementation of AdvectionKernel. * The basic equation is: - * new_latitude = latitude + u*DT - * new_longitude = longitude + v*DT + * new_latitude = latitude + v*DT + * new_longitude = longitude + u*DT * * Uses bilinear interpolation as implemented in interpolate.h */ diff --git a/advection/src/RK4AdvectionKernel.cpp b/advection/src/RK4AdvectionKernel.cpp index 3aa1006..ef03306 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 + metreToDegrees(v1 * 0.5*DT); - double lat1 = latitude + metreToDegrees(u1 * 0.5*DT); + double lon1 = longitude + metreToDegrees(u1 * 0.5*DT); + double lat1 = latitude + metreToDegrees(v1 * 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 + metreToDegrees(v2 * 0.5 * DT); - double lat2 = latitude + metreToDegrees(u2 * 0.5 * DT); + double lon2 = longitude + metreToDegrees(u2 * 0.5 * DT); + double lat2 = latitude + metreToDegrees(v2 * 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 + metreToDegrees(v3 * DT); - double lat3 = latitude + metreToDegrees(u3 * DT); + double lon3 = longitude + metreToDegrees(u3 * DT); + double lat3 = latitude + metreToDegrees(v3 * 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 + metreToDegrees((v1 + 2 * v2 + 2 * v3 + v4) / 6.0 * DT); - double latFinal = latitude + metreToDegrees((u1 + 2 * u2 + 2 * u3 + u4) / 6.0 * DT); + double lonFinal = longitude + metreToDegrees((u1 + 2 * u2 + 2 * u3 + u4) / 6.0 * DT); + double latFinal = latitude + metreToDegrees((v1 + 2 * v2 + 2 * v3 + v4) / 6.0 * DT); return {latFinal, lonFinal}; } diff --git a/advection/src/main.cpp b/advection/src/main.cpp index 8312e2e..eee999a 100644 --- a/advection/src/main.cpp +++ b/advection/src/main.cpp @@ -27,7 +27,7 @@ void advectForSomeTime(const UVGrid &uvGrid, const AdvectionKernelImpl &kernel, lon1 = templon; } catch (const out_of_range& e) { cerr << "broke out of loop!" << endl; - break; + time = 31536001; } } cout << setprecision(8) << latstart << "," << setprecision(8) << lonstart << ",begin" << i << "," << colour << endl; @@ -63,26 +63,32 @@ int main() { uvGrid->streamSlice(cout, 900); - RK4AdvectionKernel kernelRK4 = RK4AdvectionKernel(uvGrid); + auto kernelRK4 = RK4AdvectionKernel(uvGrid); // You can use https://maps.co/gis/ to visualise these points cout << "======= RK4 Integration =======" << endl; - 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, 53.53407391652826, 6.274975037862238, 0, "#ADD8E6"); + advectForSomeTime(*uvGrid, kernelRK4, 53.494053820479365, 5.673454142386921, 1, "#DC143C"); + advectForSomeTime(*uvGrid, kernelRK4, 53.49321966653616, 5.681867022043919, 2, "#50C878"); + advectForSomeTime(*uvGrid, kernelRK4, 53.581548701694324, 6.552600066543153, 3, "#FFEA00"); + advectForSomeTime(*uvGrid, kernelRK4, 53.431446729744124, 5.241592961691523, 4, "#663399"); + advectForSomeTime(*uvGrid, kernelRK4, 53.27913608324572, 4.82094897884165, 5, "#FFA500"); + advectForSomeTime(*uvGrid, kernelRK4, 53.18597595482688, 4.767667388308705, 6, "#008080"); + advectForSomeTime(*uvGrid, kernelRK4, 53.01592078792383, 4.6064205160882, 7, "#FFB6C1"); + advectForSomeTime(*uvGrid, kernelRK4, 52.72816940158886, 4.5853883152993635, 8, "#36454F"); // on land + advectForSomeTime(*uvGrid, kernelRK4, 52.56142091881038, 4.502661662924255, 9, "#1E90FF"); // Dodger Blue + advectForSomeTime(*uvGrid, kernelRK4, 52.23202593893584, 4.2825246383181845, 10, "#FFD700"); // Gold + advectForSomeTime(*uvGrid, kernelRK4, 52.08062567609582, 4.112864890830927, 11, "#6A5ACD"); // Slate Blue + advectForSomeTime(*uvGrid, kernelRK4, 51.89497719759734, 3.8114033568921686, 12, "#20B2AA"); // Light Sea Green + advectForSomeTime(*uvGrid, kernelRK4, 51.752848503723634, 3.664177951809339, 13, "#FF69B4"); // Hot Pink + advectForSomeTime(*uvGrid, kernelRK4, 51.64595756528835, 3.626319993352851, 14, "#800080"); // Purple + advectForSomeTime(*uvGrid, kernelRK4, 51.55140730645238, 3.4326152213887986, 15, "#FF4500"); // Orange Red + advectForSomeTime(*uvGrid, kernelRK4, 51.45679776223422, 3.4452813365018384, 16, "#A52A2A"); // Brown + advectForSomeTime(*uvGrid, kernelRK4, 51.41444662720727, 3.4648562416765363, 17, "#4682B4"); // Steel Blue + advectForSomeTime(*uvGrid, kernelRK4, 51.37421261203866, 3.2449264214689455, 18, "#FF6347"); // Tomato + advectForSomeTime(*uvGrid, kernelRK4, 51.29651848898365, 2.9547572241424773, 19, "#008000"); // Green + advectForSomeTime(*uvGrid, kernelRK4, 51.19705098468974, 2.7647654914530024, 20, "#B8860B"); // Dark Goldenrod + advectForSomeTime(*uvGrid, kernelRK4, 51.114719857442665, 2.577076679365129, 21, "#FFC0CB"); // Pink // advectForSomeTime(*uvGrid, kernelRK4, ,0); return 0; diff --git a/advection/src/readdata.cpp b/advection/src/readdata.cpp index 85a56e8..3597c5b 100644 --- a/advection/src/readdata.cpp +++ b/advection/src/readdata.cpp @@ -22,14 +22,7 @@ vector getVarVector(const NcVar &var) { } vector readHydrodynamicU() { - netCDF::NcFile data("../../../../data/hydrodynamic_U.h5", netCDF::NcFile::read); - - multimap< string, NcVar > vars = data.getVars(); - - return getVarVector(vars.find("uo")->second); -} - -vector readHydrodynamicV() { + // Vs and Us flipped cause the files are named incorrectly netCDF::NcFile data("../../../../data/hydrodynamic_V.h5", netCDF::NcFile::read); multimap< string, NcVar > vars = data.getVars(); @@ -37,6 +30,15 @@ vector readHydrodynamicV() { return getVarVector(vars.find("vo")->second); } +vector readHydrodynamicV() { + // Vs and Us flipped cause the files are named incorrectly + netCDF::NcFile data("../../../../data/hydrodynamic_U.h5", netCDF::NcFile::read); + + multimap< string, NcVar > vars = data.getVars(); + + return getVarVector(vars.find("uo")->second); +} + tuple, vector, vector> readGrid() { netCDF::NcFile data("../../../../data/grid.h5", netCDF::NcFile::read); multimap< string, NcVar > vars = data.getVars();