Fixing bug in backward AnalyticalAdvection where wrong cell was used

OceanParcels · Aug 2, 2024 · 6203413 · 6203413
1 parent abbf683
commit 6203413
Show file tree

Hide file tree

Showing 2 changed files with 15 additions and 16 deletions.
diff --git a/parcels/application_kernels/advectionanalytical.h b/parcels/application_kernels/advectionanalytical.h
@@ -14,7 +14,7 @@ void compute_ds(double F0, double F1, double r, double direction, double tol,
                 double *ds, double *B, double *delta){
   double up = F0 * (1-r) + F1 * r;
   double r_target = 0;
-  if (direction * up >=0.){
+  if (up >=0.){
     r_target = 1.;
   }
   *delta = -F0;
@@ -109,26 +109,26 @@ static inline StatusCode calcAdvectionAnalytical_JIT(CField *fu, CField *fv, CFi
 
     bool updateCells = 0;
     if (fabs(xsi - 1) < tol){
-      if (dataU_2D[0][1][1] > 0){
+      if (dataU_2D[0][1][1]*direction > 0){
         *xi += 1;
         xsi = 0;
         updateCells = 1;
       }
     } else if (fabs(xsi) < tol){
-      if (dataU_2D[0][0][0] < 0){
+      if (dataU_2D[0][0][0]*direction < 0){
         *xi -= 1;
         xsi = 1;
         updateCells = 1;
       }
     }
     if (fabs(eta - 1) < tol){
-      if (dataV_2D[0][1][1] > 0){
+      if (dataV_2D[0][1][1]*direction > 0){
         *yi += 1;
         eta = 0;
         updateCells = 1;
       }
     } else if (fabs(eta) < tol){
-      if (dataV_2D[0][0][0] < 0){
+      if (dataV_2D[0][0][0]*direction < 0){
         *yi -= 1;
         eta = 1;
         updateCells = 1;
@@ -144,39 +144,39 @@ static inline StatusCode calcAdvectionAnalytical_JIT(CField *fu, CField *fv, CFi
     status = getCell3D(fw, *xi, *yi, *zi, tii, dataW_3D, first_tstep_only); CHECKSTATUS(status);
     bool updateCells = 0;
     if (fabs(xsi - 1) < tol){
-      if (dataU_3D[0][1][1][1] > 0){
+      if (dataU_3D[0][1][1][1]*direction > 0){
         *xi += 1;
         xsi = 0;
         updateCells = 1;
       }
     } else if (fabs(xsi) < tol){
-      if (dataU_3D[0][0][0][0] < 0){
+      if (dataU_3D[0][0][0][0]*direction < 0){
         *xi -= 1;
         xsi = 1;
         updateCells = 1;
       }
     }
     if (fabs(eta - 1) < tol){
-      if (dataV_3D[0][1][1][1] > 0){
+      if (dataV_3D[0][1][1][1]*direction > 0){
         *yi += 1;
         eta = 0;
         updateCells = 1;
       }
     } else if (fabs(eta) < tol){
-      if (dataV_3D[0][0][0][0] < 0){
+      if (dataV_3D[0][0][0][0]*direction < 0){
         *yi -= 1;
         eta = 1;
         updateCells = 1;
       }
     }
     if (fabs(zeta - 1) < tol){
-      if (dataW_3D[0][1][1][1] > 0){
+      if (dataW_3D[0][1][1][1]*direction > 0){
         *zi += 1;
         zeta = 0;
         updateCells = 1;
       }
     } else if (fabs(zeta) < tol){
-      if (dataW_3D[0][0][0][0] < 0){
+      if (dataW_3D[0][0][0][0]*direction < 0){
         *zi -= 1;
         zeta = 1;
         updateCells = 1;

diff --git a/tests/test_advection.py b/tests/test_advection.py
@@ -531,12 +531,11 @@ def test_uniform_analytical(mode, u, v, w, direction, tmpdir):
     fieldset.V.interp_method = 'cgrid_velocity'
 
     x0, y0, z0 = 6.1, 6.2, 20
+    runtime = 4
     pset = ParticleSet(fieldset, pclass=ptype[mode], lon=x0, lat=y0, depth=z0)
-
     outfile_path = tmpdir.join("uniformanalytical.zarr")
     outfile = pset.ParticleFile(name=outfile_path, outputdt=1, chunks=(1, 1))
-
-    pset.execute(AdvectionAnalytical, runtime=4, dt=direction,
+    pset.execute(AdvectionAnalytical, runtime=runtime, dt=direction,
                  output_file=outfile)
     assert np.abs(pset.lon - x0 - pset.time * u) < 1e-6
     assert np.abs(pset.lat - y0 - pset.time * v) < 1e-6
@@ -545,7 +544,7 @@ def test_uniform_analytical(mode, u, v, w, direction, tmpdir):
 
     ds = xr.open_zarr(outfile_path)
     times = (direction*ds['time'][:]).values.astype('timedelta64[s]')[0]
-    timeref = np.arange(5).astype('timedelta64[s]')
+    timeref = np.arange(runtime+1).astype('timedelta64[s]')
     assert np.allclose(times, timeref, atol=np.timedelta64(1, 'ms'))
     lons = ds['lon'][:].values
-    assert np.allclose(lons, x0+direction*u*np.arange(5))
+    assert np.allclose(lons, x0+direction*u*np.arange(runtime+1))