djairoh
/
interactive-track-and-trace
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

fix: us and vs, tested better

This commit is contained in:
robin 2024-05-01 12:51:20 +02:00
parent 0aa58537b1
commit 84674c36de
7 changed files with 47 additions and 38 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
advection/README.md
View File

@ -25,6 +25,7 @@ void advectForSomeTime(const UVGrid &uvGrid, const AdvectionKernelImpl &kernel,
## Location of data ## Location of data
The data path is hardcoded such that the following tree structure is assumed: 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/ data/
grid.h5 grid.h5

2
advection/src/AdvectionKernel.h
View File

@ -11,7 +11,7 @@
*/ */
class AdvectionKernel { class AdvectionKernel {
public: 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 * 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. * a new latitude and longitude after being advected once for AdvectionKernel::DT time as defined above.

2
advection/src/EulerAdvectionKernel.cpp
View File

@ -9,5 +9,5 @@ EulerAdvectionKernel::EulerAdvectionKernel(std::shared_ptr<UVGrid> grid): grid(g
std::pair<double, double> EulerAdvectionKernel::advect(int time, double latitude, double longitude) const { std::pair<double, double> EulerAdvectionKernel::advect(int time, double latitude, double longitude) const {
auto [u, v] = bilinearinterpolation(*grid, time, latitude, longitude); auto [u, v] = bilinearinterpolation(*grid, time, latitude, longitude);
return {latitude+u*DT, longitude+v*DT}; return {latitude+metreToDegrees(u*DT), longitude+metreToDegrees(v*DT)};
} }

4
advection/src/EulerAdvectionKernel.h
View File

@ -7,8 +7,8 @@
/** /**
* Implementation of AdvectionKernel. * Implementation of AdvectionKernel.
* The basic equation is: * The basic equation is:
* new_latitude = latitude + u*DT * new_latitude = latitude + v*DT
* new_longitude = longitude + v*DT * new_longitude = longitude + u*DT
* *
* Uses bilinear interpolation as implemented in interpolate.h * Uses bilinear interpolation as implemented in interpolate.h
*/ */

16
advection/src/RK4AdvectionKernel.cpp
View File

@ -8,28 +8,28 @@ RK4AdvectionKernel::RK4AdvectionKernel(std::shared_ptr<UVGrid> grid): grid(grid)
std::pair<double, double> RK4AdvectionKernel::advect(int time, double latitude, double longitude) const { std::pair<double, double> RK4AdvectionKernel::advect(int time, double latitude, double longitude) const {
auto [u1, v1] = bilinearinterpolation(*grid, time, latitude, longitude); auto [u1, v1] = bilinearinterpolation(*grid, time, latitude, longitude);
// lon1, lat1 = (particle.lon + u1*.5*particle.dt, particle.lat + v1*.5*particle.dt); // lon1, lat1 = (particle.lon + u1*.5*particle.dt, particle.lat + v1*.5*particle.dt);
double lon1 = longitude + metreToDegrees(v1 * 0.5*DT); double lon1 = longitude + metreToDegrees(u1 * 0.5*DT);
double lat1 = latitude + 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] // (u2, v2) = fieldset.UV[time + .5 * particle.dt, particle.depth, lat1, lon1, particle]
auto [u2, v2] = bilinearinterpolation(*grid, time + 0.5*DT, lat1, lon1); 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) // lon2, lat2 = (particle.lon + u2*.5*particle.dt, particle.lat + v2*.5*particle.dt)
double lon2 = longitude + metreToDegrees(v2 * 0.5 * DT); double lon2 = longitude + metreToDegrees(u2 * 0.5 * DT);
double lat2 = latitude + 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] // (u3, v3) = fieldset.UV[time + .5 * particle.dt, particle.depth, lat2, lon2, particle]
auto [u3, v3] = bilinearinterpolation(*grid, time + 0.5 * DT, lat2, lon2); auto [u3, v3] = bilinearinterpolation(*grid, time + 0.5 * DT, lat2, lon2);
// lon3, lat3 = (particle.lon + u3*particle.dt, particle.lat + v3*particle.dt) // lon3, lat3 = (particle.lon + u3*particle.dt, particle.lat + v3*particle.dt)
double lon3 = longitude + metreToDegrees(v3 * DT); double lon3 = longitude + metreToDegrees(u3 * DT);
double lat3 = latitude + metreToDegrees(u3 * DT); double lat3 = latitude + metreToDegrees(v3 * DT);
// (u4, v4) = fieldset.UV[time + particle.dt, particle.depth, lat3, lon3, particle] // (u4, v4) = fieldset.UV[time + particle.dt, particle.depth, lat3, lon3, particle]
auto [u4, v4] = bilinearinterpolation(*grid, time + DT, lat3, lon3); auto [u4, v4] = bilinearinterpolation(*grid, time + DT, lat3, lon3);
double lonFinal = longitude + metreToDegrees((v1 + 2 * v2 + 2 * v3 + v4) / 6.0 * DT); double lonFinal = longitude + metreToDegrees((u1 + 2 * u2 + 2 * u3 + u4) / 6.0 * DT);
double latFinal = latitude + 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}; return {latFinal, lonFinal};
} }

42
advection/src/main.cpp
View File

@ -27,7 +27,7 @@ void advectForSomeTime(const UVGrid &uvGrid, const AdvectionKernelImpl &kernel,
lon1 = templon; lon1 = templon;
} catch (const out_of_range& e) { } catch (const out_of_range& e) {
cerr << "broke out of loop!" << endl; cerr << "broke out of loop!" << endl;
break; time = 31536001;
} }
} }
cout << setprecision(8) << latstart << "," << setprecision(8) << lonstart << ",begin" << i << "," << colour << endl; cout << setprecision(8) << latstart << "," << setprecision(8) << lonstart << ",begin" << i << "," << colour << endl;
@ -63,26 +63,32 @@ int main() {
uvGrid->streamSlice(cout, 900); uvGrid->streamSlice(cout, 900);
RK4AdvectionKernel kernelRK4 = RK4AdvectionKernel(uvGrid); auto kernelRK4 = RK4AdvectionKernel(uvGrid);
// You can use https://maps.co/gis/ to visualise these points // You can use https://maps.co/gis/ to visualise these points
cout << "======= RK4 Integration =======" << endl; cout << "======= RK4 Integration =======" << endl;
advectForSomeTime(*uvGrid, kernelRK4, 54.331795276466615, 4.845871408626273, 0, "#ADD8E6"); advectForSomeTime(*uvGrid, kernelRK4, 53.53407391652826, 6.274975037862238, 0, "#ADD8E6");
advectForSomeTime(*uvGrid, kernelRK4, 59.17208978388813, 0.32865481669722213, 1, "#DC143C"); advectForSomeTime(*uvGrid, kernelRK4, 53.494053820479365, 5.673454142386921, 1, "#DC143C");
advectForSomeTime(*uvGrid, kernelRK4, 56.18661322856813, -9.521463269751877, 2, "#50C878"); advectForSomeTime(*uvGrid, kernelRK4, 53.49321966653616, 5.681867022043919, 2, "#50C878");
advectForSomeTime(*uvGrid, kernelRK4, 46.6048774007515, -2.8605696406405947, 3, "#FFEA00"); advectForSomeTime(*uvGrid, kernelRK4, 53.581548701694324, 6.552600066543153, 3, "#FFEA00");
advectForSomeTime(*uvGrid, kernelRK4, 51.45431808648367, 1.6682437444140332, 4, "#663399"); advectForSomeTime(*uvGrid, kernelRK4, 53.431446729744124, 5.241592961691523, 4, "#663399");
advectForSomeTime(*uvGrid, kernelRK4, 51.184757012016085, -6.418923014612084, 5, "#FFA500"); advectForSomeTime(*uvGrid, kernelRK4, 53.27913608324572, 4.82094897884165, 5, "#FFA500");
advectForSomeTime(*uvGrid, kernelRK4, 54.48325546269538, 7.167517140551792, 6, "#008080"); advectForSomeTime(*uvGrid, kernelRK4, 53.18597595482688, 4.767667388308705, 6, "#008080");
advectForSomeTime(*uvGrid, kernelRK4, 55.43253322410253, -1.1712503620884716, 7, "#FFB6C1"); advectForSomeTime(*uvGrid, kernelRK4, 53.01592078792383, 4.6064205160882, 7, "#FFB6C1");
advectForSomeTime(*uvGrid, kernelRK4, 48.852815074172085, 3.4294489497130516, 8, "#36454F"); // on land advectForSomeTime(*uvGrid, kernelRK4, 52.72816940158886, 4.5853883152993635, 8, "#36454F"); // on land
advectForSomeTime(*uvGrid, kernelRK4, 58.02431905976983, 1.6892571305388995, 9, "#1E90FF"); // Dodger Blue advectForSomeTime(*uvGrid, kernelRK4, 52.56142091881038, 4.502661662924255, 9, "#1E90FF"); // Dodger Blue
advectForSomeTime(*uvGrid, kernelRK4, 58.99404571805297, 3.4121513161325425, 10, "#FFD700"); // Gold advectForSomeTime(*uvGrid, kernelRK4, 52.23202593893584, 4.2825246383181845, 10, "#FFD700"); // Gold
advectForSomeTime(*uvGrid, kernelRK4, 59.51039243098509, -1.6674160241521663, 11, "#6A5ACD"); // Slate Blue advectForSomeTime(*uvGrid, kernelRK4, 52.08062567609582, 4.112864890830927, 11, "#6A5ACD"); // Slate Blue
advectForSomeTime(*uvGrid, kernelRK4, 60.51250220636489, 1.020893006817227, 12, "#20B2AA"); // Light Sea Green advectForSomeTime(*uvGrid, kernelRK4, 51.89497719759734, 3.8114033568921686, 12, "#20B2AA"); // Light Sea Green
advectForSomeTime(*uvGrid, kernelRK4, 60.38797801281417, 3.516119068711471, 13, "#FF69B4"); // Hot Pink advectForSomeTime(*uvGrid, kernelRK4, 51.752848503723634, 3.664177951809339, 13, "#FF69B4"); // Hot Pink
advectForSomeTime(*uvGrid, kernelRK4, 60.02637651315464, -2.4546004365354697, 14, "#800080"); // Purple advectForSomeTime(*uvGrid, kernelRK4, 51.64595756528835, 3.626319993352851, 14, "#800080"); // Purple
advectForSomeTime(*uvGrid, kernelRK4, 58.732929454411305, 3.649791893455804, 15, "#FF4500"); // Orange Red 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); // advectForSomeTime(*uvGrid, kernelRK4, ,0);
return 0; return 0;

18
advection/src/readdata.cpp
View File

@ -22,14 +22,7 @@ vector<T> getVarVector(const NcVar &var) {
} }
vector<double> readHydrodynamicU() { vector<double> readHydrodynamicU() {
netCDF::NcFile data("../../../../data/hydrodynamic_U.h5", netCDF::NcFile::read); // Vs and Us flipped cause the files are named incorrectly
multimap< string, NcVar > vars = data.getVars();
return getVarVector<double>(vars.find("uo")->second);
}
vector<double> readHydrodynamicV() {
netCDF::NcFile data("../../../../data/hydrodynamic_V.h5", netCDF::NcFile::read); netCDF::NcFile data("../../../../data/hydrodynamic_V.h5", netCDF::NcFile::read);
multimap< string, NcVar > vars = data.getVars(); multimap< string, NcVar > vars = data.getVars();
@ -37,6 +30,15 @@ vector<double> readHydrodynamicV() {
return getVarVector<double>(vars.find("vo")->second); return getVarVector<double>(vars.find("vo")->second);
} }
vector<double> 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<double>(vars.find("uo")->second);
}
tuple<vector<int>, vector<double>, vector<double>> readGrid() { tuple<vector<int>, vector<double>, vector<double>> readGrid() {
netCDF::NcFile data("../../../../data/grid.h5", netCDF::NcFile::read); netCDF::NcFile data("../../../../data/grid.h5", netCDF::NcFile::read);
multimap< string, NcVar > vars = data.getVars(); multimap< string, NcVar > vars = data.getVars();