euler int

advection/src/AdvectionKernel.h

@@ -5,7 +5,7 @@
 
 #include "Vel.h"
 
-#define DT 4000
+#define DT 50
 
 class AdvectionKernel {
 public:

advection/src/CMakeLists.txt

@@ -20,6 +20,8 @@ add_executable(Advection main.cpp
         AdvectionKernel.h
         EulerAdvectionKernel.cpp
         EulerAdvectionKernel.h
+        RK4AdvectionKernel.cpp
+        RK4AdvectionKernel.h
 )
 
 execute_process(

advection/src/EulerAdvectionKernel.cpp

@@ -1,11 +1,13 @@
 
 #include "EulerAdvectionKernel.h"
+#include "interpolate.h"
 
 using namespace std;
 
-EulerAdvectionKernel::EulerAdvectionKernel(std::shared_ptr<UVGrid> grid) { }
+EulerAdvectionKernel::EulerAdvectionKernel(std::shared_ptr<UVGrid> grid): grid(grid) { }
 
 std::pair<double, double> EulerAdvectionKernel::advect(int time, double latitude, double longitude) const {
-    auto [u,v] = (*grid)[time, latitude, longitude];
+    auto [u, v] = bilinearinterpolation(*grid, time, latitude, longitude);
+
     return {latitude+u*DT, longitude+v*DT};
 }

advection/src/RK4AdvectionKernel.cpp

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

advection/src/RK4AdvectionKernel.h

@@ -0,0 +1,18 @@
+#ifndef ADVECTION_RK4ADVECTIONKERNEL_H
+#define ADVECTION_RK4ADVECTIONKERNEL_H
+
+
+#include "AdvectionKernel.h"
+#include "UVGrid.h"
+
+class RK4AdvectionKernel: public AdvectionKernel {
+private:
+    std::shared_ptr<UVGrid> grid;
+public:
+    explicit RK4AdvectionKernel(std::shared_ptr<UVGrid> grid);
+    std::pair<double, double> advect(int time, double latitude, double longitude) const override;
+
+};
+
+
+#endif //ADVECTION_RK4ADVECTIONKERNEL_H

advection/src/main.cpp

@@ -1,52 +1,37 @@
 #include "interpolate.h"
 #include "Vel.h"
+#include "EulerAdvectionKernel.h"
+#include "RK4AdvectionKernel.h"
 #include <ranges>
 #include <chrono>
 
 using namespace std;
 
 int main() {
-//    UVGrid uvGrid;
     std::shared_ptr<UVGrid> uvGrid = std::make_shared<UVGrid>();
-    uvGrid->streamSlice(cout, 100);
+//    uvGrid->streamSlice(cout, 100);
 
-    int N = 10000000;  // Number of points
+    EulerAdvectionKernel kernelEuler = EulerAdvectionKernel(uvGrid);
+    RK4AdvectionKernel kernelRK4 = RK4AdvectionKernel(uvGrid);
 
-    int time_start = 0;
-    int time_end = 31536000;
+    double latstart = 54.860801, lonstart = 4.075492;
 
-    double lat_start = 46.125;
-    double lat_end = 62.625;
-
-    double lon_start = -15.875;
-    double lon_end = 12.875;
-
-    double lat_step = (lat_end - lat_start) / (N - 1);
-    double lon_step = (lon_end - lon_start) / (N - 1);
-    int time_step = (time_end - time_start) / (N - 1);
-
-    vector<tuple<int, double, double>> points;
-
-    for(int i = 0; i < N; i++) {
-        points.push_back({time_start+time_step*i, lat_start+lat_step*i, lon_start+lon_step*i});
+    double lat1 = latstart, lon1 = lonstart;
+    cout << "======= Euler Integration =======" << endl;
+    for(int time = 100; time <= 10000; time += DT) {
+        cout << "lat = " << lat1 << " lon = " << lon1 << endl;
+        auto [templat, templon] = kernelEuler.advect(time, lat1, lon1);
+        lat1 = templat;
+        lon1 = templon;
     }
-
-    auto start = chrono::high_resolution_clock::now();
-
-    auto x = bilinearinterpolation(*uvGrid, points);
-
-    auto stop = chrono::high_resolution_clock::now();
-
-    auto duration = chrono::duration_cast<std::chrono::milliseconds >(stop - start);
-
-    cout << "Time taken for " << N << " points was " << duration.count()/1000. << " seconds\n";
-
-    // Do something with result in case of optimisation
-    double sum = 0;
-    for (auto [u,v]: x) {
-        sum += u + v;
+    cout << "======= RK4 Integration =======" << endl;
+    lat1 = latstart, lon1 = lonstart;
+    for(int time = 100; time <= 10000; time += DT) {
+        cout << "lat = " << lat1 << " lon = " << lon1 << endl;
+        auto [templat, templon] = kernelRK4.advect(time, lat1, lon1);
+        lat1 = templat;
+        lon1 = templon;
     }
-    cout << "Sum = " << sum << endl;
 
     return 0;
 }