Merge pull request #10 from MakeNEnjoy/robin-advection

Robin advection

advection/README.md

@@ -0,0 +1,46 @@
+## What is new?
+There is one new added component: `AdvectionKernel`s which is an "interface" (i.e an abstract class).
+There are two implementations simple Euler integration called `EulerIntegrationKernel` and 
+Runge Kutta integration called `RK4AdvectionKernel`.
+
+Main function gives a good example of how to use the library. Especially the following function which prints the
+position of the particle at every time step.
+```Cpp
+template <typename AdvectionKernelImpl>
+void advectForSomeTime(const UVGrid &uvGrid, const AdvectionKernelImpl &kernel, double latstart, double lonstart) {
+
+    // Require at compile time that kernel derives from the abstract class AdvectionKernel
+    static_assert(std::is_base_of<AdvectionKernel, AdvectionKernelImpl>::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;
+    }
+}
+```
+
+
+## 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
+  hydrodynamic_U.h5
+  hydrodynamic_V.h5
+interactive-track-and-trace/
+  opening-hdf5/
+     ...
+```
+
+## Compiling
+Let the current directory be the `src` directory. Run:
+```shell
+mkdir build
+cd build
+cmake ..
+make
+```

advection/src/AdvectionKernel.h

@@ -0,0 +1,32 @@
+#ifndef ADVECTION_ADVECTIONKERNEL_H
+#define ADVECTION_ADVECTIONKERNEL_H
+
+#include <tuple>
+
+#include "Vel.h"
+
+
+/*
+ * Implement this class for every integration method.
+ */
+class AdvectionKernel {
+public:
+    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.
+     * @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) 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

advection/src/CMakeLists.txt

@@ -1,5 +1,5 @@
 cmake_minimum_required (VERSION 3.28)
-project (LinearInterpolate)
+project (Advection)
 
 set(CMAKE_CXX_STANDARD 23)
 set(CMAKE_CXX_STANDARD_REQUIRED ON)
@@ -8,7 +8,7 @@ set(CMAKE_EXPORT_COMPILE_COMMANDS ON)
 
 find_package(netCDF REQUIRED)
 
-add_executable(LinearInterpolate main.cpp
+add_executable(Advection main.cpp
         readdata.cpp
         readdata.h
         interpolate.cpp
@@ -16,7 +16,13 @@ add_executable(LinearInterpolate main.cpp
         UVGrid.cpp
         UVGrid.h
         Vel.h
-        Vel.cpp)
+        Vel.cpp
+        AdvectionKernel.h
+        EulerAdvectionKernel.cpp
+        EulerAdvectionKernel.h
+        RK4AdvectionKernel.cpp
+        RK4AdvectionKernel.h
+)
 
 execute_process(
     COMMAND nc-config --includedir
@@ -30,7 +36,7 @@ execute_process(
     OUTPUT_STRIP_TRAILING_WHITESPACE
 )
 
-target_include_directories(LinearInterpolate PUBLIC ${netCDF_INCLUDE_DIR})
+target_include_directories(Advection PUBLIC ${netCDF_INCLUDE_DIR})
 
 find_library(NETCDF_LIB NAMES netcdf-cxx4 netcdf_c++4 PATHS ${NETCDFCXX_LIB_DIR} NO_DEFAULT_PATH)
-target_link_libraries(LinearInterpolate ${NETCDF_LIB})
+target_link_libraries(Advection ${NETCDF_LIB})

advection/src/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) const {
+    auto [u, v] = bilinearinterpolate(*grid, time, latitude, longitude);
+
+    return {latitude+metreToDegrees(v*DT), longitude+metreToDegrees(u*DT)};
+}

advection/src/EulerAdvectionKernel.h

@@ -0,0 +1,25 @@
+#ifndef ADVECTION_EULERADVECTIONKERNEL_H
+#define ADVECTION_EULERADVECTIONKERNEL_H
+
+#include "AdvectionKernel.h"
+#include "UVGrid.h"
+
+/**
+ * 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) const override;
+
+};
+
+
+#endif //ADVECTION_EULERADVECTIONKERNEL_H

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 [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};
+}

advection/src/RK4AdvectionKernel.h

@@ -0,0 +1,22 @@
+#ifndef ADVECTION_RK4ADVECTIONKERNEL_H
+#define ADVECTION_RK4ADVECTIONKERNEL_H
+
+#include "AdvectionKernel.h"
+#include "UVGrid.h"
+
+/**
+ * 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) const override;
+
+};
+
+
+#endif //ADVECTION_RK4ADVECTIONKERNEL_H

advection/src/UVGrid.cpp

@@ -34,7 +34,12 @@ UVGrid::UVGrid() {
 }
 
 const Vel &UVGrid::operator[](size_t timeIndex, size_t latIndex, size_t lonIndex) const {
-    size_t index = timeIndex * (latSize * lonSize) + latIndex * lonIndex + lonIndex;
+    if(timeIndex < 0 or timeIndex >= timeSize
+        or latIndex < 0 or latIndex >= latSize
+        or lonIndex < 0 or lonIndex >= lonSize) {
+        throw std::out_of_range("Index out of bounds");
+    }
+    size_t index = timeIndex * (latSize * lonSize) + latIndex * lonSize + lonIndex;
     return uvData[index];
 }
 
@@ -58,4 +63,4 @@ void UVGrid::streamSlice(ostream &os, size_t t) {
         }
       os << endl;
     }
-}
+}

advection/src/UVGrid.h

@@ -1,5 +1,5 @@
-#ifndef LINEARINTERPOLATE_UVGRID_H
+#ifndef ADVECTION_UVGRID_H
-#define LINEARINTERPOLATE_UVGRID_H
+#define ADVECTION_UVGRID_H
 
 #include <vector>
 #include "Vel.h"
@@ -13,6 +13,9 @@ private:
 public:
     UVGrid();
 
+    /**
+     * The matrix has shape (timeSize, latSize, lonSize)
+     */
     size_t timeSize;
     size_t latSize;
     size_t lonSize;
@@ -35,17 +38,28 @@ public:
      */
     int timeStep() const;
 
+    /**
+     * times, lats, lons are vector of length timeSize, latSize, lonSize respectively.
+     * The maintain the following invariant:
+     * grid[timeIndex,latIndex,lonIndex] gives the u,v at the point with latitude at lats[latIndex],
+     * with longitude at lons[lonIndex], and with time at times[timeIndex].
+     */
     std::vector<int> times;
     std::vector<double> lats;
     std::vector<double> 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;
 
+    /**
+     * Streams a slice at timeIndex t of the matrix to the outstream given by os
+     * @param os outstream
+     * @param t index with which to slice matrix
+     */
     void streamSlice(std::ostream &os, size_t t);
 };
 
-#endif //LINEARINTERPOLATE_UVGRID_H
+#endif //ADVECTION_UVGRID_H

advection/src/Vel.h

@@ -1,5 +1,5 @@
-#ifndef LINEARINTERPOLATE_VEL_H
+#ifndef ADVECTION_VEL_H
-#define LINEARINTERPOLATE_VEL_H
+#define ADVECTION_VEL_H
 
 #include <utility>
 #include <stdexcept>
@@ -41,4 +41,4 @@ Vel operator*(Scalar scalar, const Vel& p) {
     return Vel(p.u * scalar, p.v * scalar);
 }
 
-#endif //LINEARINTERPOLATE_VEL_H
+#endif //ADVECTION_VEL_H

advection/src/interpolate.cpp

@@ -2,7 +2,7 @@
 
 using namespace std;
 
-Vel bilinearInterpolate(const UVGrid &uvGrid, int time, double lat, double lon) {
+Vel bilinearinterpolate(const UVGrid &uvGrid, int time, double lat, double lon) {
     double latStep = uvGrid.latStep();
     double lonStep = uvGrid.lonStep();
     int timeStep = uvGrid.timeStep();
@@ -36,11 +36,11 @@ Vel bilinearInterpolate(const UVGrid &uvGrid, int time, double lat, double lon)
     return point;
 }
 
-vector<Vel> bilinearInterpolate(const UVGrid &uvGrid, vector<tuple<int, double, double>> points) {
+vector<Vel> bilinearinterpolation(const UVGrid &uvGrid, vector<tuple<int, double, double>> points) {
     vector<Vel> result;
     result.reserve(points.size());
     for (auto [time, lat, lon]: points) {
-        result.push_back(bilinearInterpolate(uvGrid, time, lat, lon));
+        result.push_back(bilinearinterpolate(uvGrid, time, lat, lon));
     }
 
     return result;

advection/src/interpolate.h

@@ -1,5 +1,5 @@
-#ifndef LINEARINTERPOLATE_INTERPOLATE_H
+#ifndef ADVECTION_INTERPOLATE_H
-#define LINEARINTERPOLATE_INTERPOLATE_H
+#define ADVECTION_INTERPOLATE_H
 
 #include <vector>
 
@@ -15,7 +15,7 @@
  * @param lon longitude of point
  * @return interpolated velocity
  */
-Vel bilinearInterpolate(const UVGrid &uvGrid, int time, double lat, double lon);
+Vel bilinearinterpolate(const UVGrid &uvGrid, int time, double lat, double lon);
 
 /**
  * Helper function for bilnearly interpolating a vector of points
@@ -23,6 +23,6 @@ Vel bilinearInterpolate(const UVGrid &uvGrid, int time, double lat, double lon);
  * @param points vector of points
  * @return interpolated velocities
  */
-std::vector<Vel> bilinearInterpolate(const UVGrid &uvGrid, std::vector<std::tuple<int, double, double>> points);
+std::vector<Vel> bilinearinterpolation(const UVGrid &uvGrid, std::vector<std::tuple<int, double, double>> points);
 
-#endif //LINEARINTERPOLATE_INTERPOLATE_H
+#endif //ADVECTION_INTERPOLATE_H

advection/src/main.cpp

@@ -0,0 +1,95 @@
+#include <ranges>
+#include <iomanip>
+#include <stdexcept>
+
+#include "interpolate.h"
+#include "Vel.h"
+#include "EulerAdvectionKernel.h"
+#include "RK4AdvectionKernel.h"
+#include "interpolate.h"
+
+#define NotAKernelError "Template parameter T must derive from AdvectionKernel"
+
+using namespace std;
+
+template <typename AdvectionKernelImpl>
+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<AdvectionKernel, AdvectionKernelImpl>::value, NotAKernelError);
+
+    double lat1 = latstart, lon1 = lonstart;
+    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;
+            time = 31536001;
+        }
+    }
+    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 << bilinearinterpolate(*uvGrid, time, 53.80956379699079, -1.6496306344654406) << endl;
+    cout << bilinearinterpolate(*uvGrid, time, 55.31428895563707, -2.851581041325997) << endl;
+    cout << bilinearinterpolate(*uvGrid, time, 47.71548983067583, -1.8704054037408626) << endl;
+    cout << bilinearinterpolate(*uvGrid, time, 56.23521060314398, 8.505979324950573) << endl;
+    cout << bilinearinterpolate(*uvGrid, time, 53.135645440244375, 8.505979324950573) << endl;
+    cout << bilinearinterpolate(*uvGrid, time, 56.44761278775708, -4.140629303756164) << endl;
+    cout << bilinearinterpolate(*uvGrid, time, 52.67625153110339, 0.8978569759455872) << endl;
+    cout << bilinearinterpolate(*uvGrid, time, 52.07154079279377, 4.627951041411331) << endl;
+
+    cout << "=== ocean === (should give not 0)" << endl;
+    cout << bilinearinterpolate(*uvGrid, time, 47.43923166616274, -4.985451481829083) << endl;
+    cout << bilinearinterpolate(*uvGrid, time, 50.68943556852362, -9.306162999561733) << endl;
+    cout << bilinearinterpolate(*uvGrid, time, 53.70606799886677, -4.518347647034465) << endl;
+    cout << bilinearinterpolate(*uvGrid, time, 60.57987114267971, -12.208262973672621) << endl;
+    cout << bilinearinterpolate(*uvGrid, time, 46.532221548197285, -13.408189172582638) << endl;
+    cout << bilinearinterpolate(*uvGrid, time, 50.92725094937812, 1.3975824052375256) << endl;
+    cout << bilinearinterpolate(*uvGrid, time, 51.4028921682209, 2.4059571950925203) << endl;
+    cout << bilinearinterpolate(*uvGrid, time, 53.448445236769004, 0.7996966058017515) << endl;
+//    cout << bilinearinterpolate(*uvGrid, time, ) << endl;
+}
+
+int main() {
+    std::shared_ptr<UVGrid> uvGrid = std::make_shared<UVGrid>();
+
+    uvGrid->streamSlice(cout, 900);
+
+    auto kernelRK4 = RK4AdvectionKernel(uvGrid);
+
+//    You can use https://maps.co/gis/ to visualise these points
+    cout << "======= RK4 Integration =======" << endl;
+    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;
+}

advection/src/readdata.cpp

@@ -22,14 +22,7 @@ vector<T> getVarVector(const NcVar &var) {
 }
 
 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);
 
     multimap< string, NcVar > vars = data.getVars();
@@ -37,6 +30,15 @@ vector<double> readHydrodynamicV() {
     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() {
     netCDF::NcFile data("../../../../data/grid.h5", netCDF::NcFile::read);
     multimap< string, NcVar > vars = data.getVars();

advection/src/readdata.h

@@ -1,5 +1,5 @@
-#ifndef LINEARINTERPOLATE_READDATA_H
+#ifndef ADVECTION_READDATA_H
-#define LINEARINTERPOLATE_READDATA_H
+#define ADVECTION_READDATA_H
 
 /**
  * reads the file hydrodynamic_U.h5
@@ -19,4 +19,4 @@ std::vector<double> readHydrodynamicV();
  */
 std::tuple<std::vector<int>, std::vector<double>, std::vector<double>> readGrid();
 
-#endif //LINEARINTERPOLATE_READDATA_H
+#endif //ADVECTION_READDATA_H

linear-interpolation/README.md

@@ -1,39 +0,0 @@
-## Location of data
-The data path is hardcoded such that the following tree structure is assumed:
-```
-data/
-  grid.h5
-  hydrodynamic_U.h5
-  hydrodynamic_V.h5
-interactive-track-and-trace/
-  opening-hdf5/
-     ...
-```
-
-## Compiling
-Let the current directory be the `src` directory. Run:
-```shell
-mkdir build
-cd build
-cmake ..
-make
-```
-
-### Building with Linux
-Makes use of `mdspan` which is not supported by glibc++ at time of writing. See [compiler support](https://en.cppreference.com/w/cpp/compiler_support/23) for `mdspan`. The solution to this is to use Clang and libc++; this is configured in our CMake setup, however the default installation of the `netcdf-cxx` package on at least Arch linux (and suspectedly Debian derivatives as well) specifically builds for the glibc implementation. To get the netcdf C++ bindings functional with the libc++ implementation, one needs to build from source. On Linux, this requires a few changes to the CMake file included with the netcdf-cxx source code, which are detailed below.
-
-Step-by-step to build the program using clang++ and libc++ on linux:
- 1. Download the source code of netcdf-cxx, found at 'https://github.com/Unidata/netcdf-cxx4/releases/tag/v4.3.1' (make sure to download the release source code, as the master branch contains non-compilable code).
- 2. Edit the CMakeLists.txt file, by appending '-stdlib=libc++' to the `CMAKE_CXX_FLAGS` variable in line 430. This means line 430 should read: 
-    ```cmake 
-    SET(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -g -Wall -Wno-unused-variable -Wno-unused-parameter -stdlib=libc++")
-    ```
- 2. Build the source code with the following: 
-  ```sh 
-    mkdir build && cd build
-    cmake .. -DCMAKE_CXX_COMPILER=/usr/bin/clang++
-    make
-    ctest
-    sudo make install
-  ```
- 3. Now the code should compile through the standard steps described in the Compiling section.

linear-interpolation/src/main.cpp

@@ -1,51 +0,0 @@
-#include "interpolate.h"
-#include "Vel.h"
-#include <ranges>
-#include <chrono>
-
-using namespace std;
-
-int main() {
-    UVGrid uvGrid;
-    uvGrid.streamSlice(cout, 100);
-
-    int N = 10000000;  // Number of points
-
-    int time_start = 0;
-    int time_end = 31536000;
-
-    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});
-    }
-
-    auto start = chrono::high_resolution_clock::now();
-
-    auto x = bilinearInterpolate(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 << "Sum = " << sum << endl;
-
-    return 0;
-}