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

fixed RK4 integration

This commit is contained in:
robin 2024-04-29 17:20:58 +02:00
parent ccac46ff40
commit 37d0593f65
3 changed files with 36 additions and 27 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

24
advection/src/RK4AdvectionKernel.cpp
View File

@ -6,30 +6,30 @@ 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);
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 + u1*.5*DT;
double lat1 = latitude + v1*.5*DT;
double lon1 = longitude + v1 * .5*DT;
double lat1 = latitude + u1*.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);
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 + u2 * 0.5 * DT;
double lat2 = latitude + v2 * 0.5 * DT;
double lon2 = longitude + v2 * 0.5 * DT;
double lat2 = latitude + u2 * 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);
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 + u3 * DT;
double lat3 = latitude + v3 * DT;
double lon3 = longitude + v3 * DT;
double lat3 = latitude + u3 * DT;
// (u4, v4) = fieldset.UV[time + particle.dt, particle.depth, lat3, lon3, particle]
auto [v4, u4] = bilinearinterpolation(*grid, time + DT, lat3, lon3);
auto [u4, v4] = 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;
double lonFinal = longitude + (v1 + 2 * v2 + 2 * v3 + v4) / 6.0 * DT;
double latFinal = latitude + (u1 + 2 * u2 + 2 * u3 + u4) / 6.0 * DT;
return {latFinal, lonFinal};
}

5
advection/src/UVGrid.cpp
View File

@ -34,6 +34,11 @@ UVGrid::UVGrid() {
}
const Vel &UVGrid::operator[](size_t timeIndex, size_t latIndex, size_t lonIndex) const {
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 * lonIndex + lonIndex;
return uvData[index];
}

34
advection/src/main.cpp
View File

@ -5,33 +5,37 @@
#include <ranges>
#include <chrono>
#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) {
static_assert(std::is_base_of<AdvectionKernel, AdvectionKernelImpl>::value, NotAKernelError);
double lat1 = latstart, lon1 = lonstart;
for(int time = 100; time <= 10000; time += DT) {
cout << "lat = " << lat1 << " lon = " << lon1 << endl;
auto [templat, templon] = kernel.advect(time, lat1, lon1);
lat1 = templat;
lon1 = templon;
}
}
int main() {
std::shared_ptr<UVGrid> uvGrid = std::make_shared<UVGrid>();
// uvGrid->streamSlice(cout, 100);
EulerAdvectionKernel kernelEuler = EulerAdvectionKernel(uvGrid);
RK4AdvectionKernel kernelRK4 = RK4AdvectionKernel(uvGrid);
double latstart = 54.860801, lonstart = 4.075492;
double latstart = 52.881770, lonstart = 3.079979;
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;
}
advectForSomeTime(*uvGrid, kernelEuler, latstart, lonstart);
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;
}
advectForSomeTime(*uvGrid, kernelRK4, latstart, lonstart);
return 0;
}