feat: boundary conditions

particle-track-and-trace/src/advection/kernel/AdvectionKernel.h

@@ -0,0 +1,31 @@
+#ifndef ADVECTIONKERNEL_H
+#define ADVECTIONKERNEL_H
+
+#include <tuple>
+
+#include "../Vel.h"
+
+/*
+ * Implement this class for every integration method.
+ */
+class AdvectionKernel {
+public:
+    /**
+     * 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.
+     * @param time Time since the beginning of the data
+     * @param latitude Latitude of particle
+     * @param longitude Longitude of particle
+     * @return A pair of latitude and longitude of particle.
+     */
+    virtual std::pair<double, double> advect(int time, double latitude, double longitude, int dt) 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.;
+    }
+
+    virtual ~AdvectionKernel() = default; // Apparently I need this, idk why
+};
+
+#endif //ADVECTIONKERNEL_H

particle-track-and-trace/src/advection/kernel/EulerAdvectionKernel.cpp

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

particle-track-and-trace/src/advection/kernel/EulerAdvectionKernel.h

@@ -0,0 +1,26 @@
+#ifndef EULERADVECTIONKERNEL_H
+#define EULERADVECTIONKERNEL_H
+
+#include "AdvectionKernel.h"
+#include "../UVGrid.h"
+#include <memory>
+
+/**
+ * Implementation of AdvectionKernel.
+ * The basic equation is:
+ * new_latitude = latitude + v*DT
+ * new_longitude = longitude + u*DT
+ *
+ * Uses bilinear interpolation as implemented in interpolate.h
+ */
+class EulerAdvectionKernel: public AdvectionKernel {
+private:
+    std::shared_ptr<UVGrid> grid;
+public:
+    explicit EulerAdvectionKernel(std::shared_ptr<UVGrid> grid);
+    std::pair<double, double> advect(int time, double latitude, double longitude, int dt) const override;
+
+};
+
+
+#endif //EULERADVECTIONKERNEL_H

particle-track-and-trace/src/advection/kernel/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, int dt) const {
+    auto [u1, v1] = bilinearinterpolate(*grid, time, latitude, longitude);
+//    lon1, lat1 = (particle.lon + u1*.5*particle.dt, particle.lat + v1*.5*particle.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] = bilinearinterpolate(*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(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] = bilinearinterpolate(*grid, time + 0.5 * dt, lat2, lon2);
+
+//    lon3, lat3 = (particle.lon + u3*particle.dt, particle.lat + v3*particle.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] = bilinearinterpolate(*grid, time + dt, lat3, lon3);
+
+    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};
+}

particle-track-and-trace/src/advection/kernel/RK4AdvectionKernel.h

@@ -0,0 +1,23 @@
+#ifndef RK4ADVECTIONKERNEL_H
+#define RK4ADVECTIONKERNEL_H
+
+#include "AdvectionKernel.h"
+#include "../UVGrid.h"
+#include <memory>
+
+/**
+ * Implementation of Advection kernel using RK4 integration
+ * See https://en.wikipedia.org/wiki/Runge%E2%80%93Kutta_methods for more details.
+ * Uses bilinear interpolation as implemented in interpolate.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, int dt) const override;
+
+};
+
+
+#endif //RK4ADVECTIONKERNEL_H

particle-track-and-trace/src/advection/kernel/SnapBoundaryConditionKernel.cpp

@@ -0,0 +1,26 @@
+#include "SnapBoundaryConditionKernel.h"
+
+SnapBoundaryConditionKernel::SnapBoundaryConditionKernel(std::unique_ptr<AdvectionKernel> kernel,
+                                                         std::shared_ptr<UVGrid> uvGrid) :
+                                                         kernel(std::move(kernel)),
+                                                         uvGrid(uvGrid) { }
+std::pair<double, double> SnapBoundaryConditionKernel::advect(int time, double latitude, double longitude, int dt) const {
+  auto [newLat, newLon] = kernel->advect(time, latitude, longitude, dt);
+  double minLat = uvGrid->lats.front();
+  double maxLat = uvGrid->lats.back();
+  double minLon = uvGrid->lons.front();
+  double maxLon = uvGrid->lons.back();
+  if (newLat < minLat) {
+    newLat = minLat;
+  }
+  if (newLat > maxLat) {
+    newLat = maxLat;
+  }
+  if (newLon < minLon) {
+    newLon = minLon;
+  }
+  if (newLon > maxLon) {
+    newLon = maxLon;
+  }
+  return {newLat, newLon};
+}

particle-track-and-trace/src/advection/kernel/SnapBoundaryConditionKernel.h

@@ -0,0 +1,22 @@
+#ifndef SNAPBOUNDARYCONDITIONKERNEL_H
+#define SNAPBOUNDARYCONDITIONKERNEL_H
+
+#include <memory>
+#include "AdvectionKernel.h"
+#include "../UVGrid.h"
+
+/**
+ * When
+ */
+class SnapBoundaryConditionKernel: public AdvectionKernel {
+  std::unique_ptr<AdvectionKernel> kernel;
+  std::shared_ptr<UVGrid> uvGrid;
+public:
+  SnapBoundaryConditionKernel(std::unique_ptr<AdvectionKernel> kernel, std::shared_ptr<UVGrid> uvGrid);
+
+private:
+  std::pair<double, double> advect(int time, double latitude, double longitude, int dt) const override;
+};
+
+
+#endif //SNAPBOUNDARYCONDITIONKERNEL_H