Repository: OceanParcels/parcels
Pull request #386 (2018-06-15 09:47:09)

Errorcode repeat


Patch:
diff --git a/parcels/codegenerator.py b/parcels/codegenerator.py
index 8c17011f..89c91efd 100644
--- a/parcels/codegenerator.py
+++ b/parcels/codegenerator.py
@@ -613,6 +613,30 @@ class LoopGenerator(object):
 
         ccode += [str(c.Typedef(c.GenerableStruct("", vdecl, declname=self.ptype.name)))]
 
+        args = [c.Pointer(c.Value(self.ptype.name, "particle_backup")),
+                c.Pointer(c.Value(self.ptype.name, "particle"))]
+        p_back_set_decl = c.FunctionDeclaration(c.Static(c.DeclSpecifier(c.Value("void", "set_particle_backup"),
+                                                         spec='inline')), args)
+        body = []
+        for v in self.ptype.variables:
+            if v.dtype != np.uint64 and v.name not in ['dt', 'state']:
+                body += [c.Assign(("particle_backup->%s" % v.name), ("particle->%s" % v.name))]
+        p_back_set_body = c.Block(body)
+        p_back_set = str(c.FunctionBody(p_back_set_decl, p_back_set_body))
+        ccode += [p_back_set]
+
+        args = [c.Pointer(c.Value(self.ptype.name, "particle_backup")),
+                c.Pointer(c.Value(self.ptype.name, "particle"))]
+        p_back_get_decl = c.FunctionDeclaration(c.Static(c.DeclSpecifier(c.Value("void", "get_particle_backup"),
+                                                         spec='inline')), args)
+        body = []
+        for v in self.ptype.variables:
+            if v.dtype != np.uint64 and v.name not in ['dt', 'state']:
+                body += [c.Assign(("particle->%s" % v.name), ("particle_backup->%s" % v.name))]
+        p_back_get_body = c.Block(body)
+        p_back_get = str(c.FunctionBody(p_back_get_decl, p_back_get_body))
+        ccode += [p_back_get]
+
         if c_include:
             ccode += [c_include]
 
@@ -631,24 +655,27 @@ class LoopGenerator(object):
                               + list(const_args.keys()))
         # Inner loop nest for forward runs
         sign_dt = c.Assign("sign_dt", "dt > 0 ? 1 : -1")
+        particle_backup = c.Statement("%s particle_backup" % self.ptype.name)
         sign_end_part = c.Assign("sign_end_part", "endtime - particles[p].time > 0 ? 1 : -1")
         dt_pos = c.Assign("__dt", "fmin(fabs(particles[p].dt), fabs(endtime - particles[p].time))")
         dt_0_break = c.If("particles[p].dt == 0", c.Statement("break"))
         notstarted_continue = c.If("(sign_end_part != sign_dt) && (particles[p].dt != 0)",
                                    c.Statement("continue"))
-        body = [c.Assign("res", "%s(&(particles[p]), %s)" % (funcname, fargs_str))]
+        body = [c.Statement("set_particle_backup(&particle_backup, &(particles[p]))")]
+        body += [c.Assign("res", "%s(&(particles[p]), %s)" % (funcname, fargs_str))]
         body += [c.Assign("particles[p].state", "res")]  # Store return code on particle
         body += [c.If("res == SUCCESS", c.Block([c.Statement("particles[p].time += sign_dt * __dt"),
                                                  dt_pos, dt_0_break, c.Statement("continue")]))]
-        body += [c.If("res == REPEAT", c.Block([dt_pos, c.Statement("continue")]),
-                      c.Statement("break"))]
+        body += [c.If("res == REPEAT",
+                 c.Block([c.Statement("get_particle_backup(&particle_backup, &(particles[p]))"),
+                         dt_pos, c.Statement("break")]), c.Statement("break"))]
 
         time_loop = c.While("__dt > __tol || particles[p].dt == 0", c.Block(body))
         part_loop = c.For("p = 0", "p < num_particles", "++p",
                           c.Block([sign_end_part, notstarted_continue, dt_pos, time_loop]))
         fbody = c.Block([c.Value("int", "p, sign_dt, sign_end_part"), c.Value("ErrorCode", "res"),
                          c.Value("double", "__dt, __tol"), c.Assign("__tol", "1.e-6"),
-                         sign_dt, part_loop])
+                         sign_dt, particle_backup, part_loop])
         fdecl = c.FunctionDeclaration(c.Value("void", "particle_loop"), args)
         ccode += [str(c.FunctionBody(fdecl, fbody))]
         return "\n\n".join(ccode)
diff --git a/parcels/kernel.py b/parcels/kernel.py
index 917d07e9..7c2045e6 100644
--- a/parcels/kernel.py
+++ b/parcels/kernel.py
@@ -187,7 +187,12 @@ class Kernel(object):
     def execute_python(self, pset, endtime, dt):
         """Performs the core update loop via Python"""
         sign_dt = np.sign(dt)
+
+        # back up variables in case of ErrorCode.Repeat
+        p_var_back = {}
+
         for p in pset.particles:
+            ptype = p.getPType()
             # Don't execute particles that aren't started yet
             sign_end_part = np.sign(endtime - p.time)
             if (sign_end_part != sign_dt) and (dt != 0):
@@ -196,6 +201,8 @@ class Kernel(object):
             # Compute min/max dt for first timestep
             dt_pos = min(abs(p.dt), abs(endtime - p.time))
             while dt_pos > 1e-6 or dt == 0:
+                for var in ptype.variables:
+                    p_var_back[var.name] = getattr(p, var.name)
                 try:
                     res = self.pyfunc(p, pset.fieldset, p.time, sign_dt * dt_pos)
                 except FieldSamplingError as fse:
@@ -219,8 +226,11 @@ class Kernel(object):
                     continue
                 elif res == ErrorCode.Repeat:
                     # Try again without time update
+                    for var in ptype.variables:
+                        if var.name not in ['dt', 'state']:
+                            setattr(p, var.name, p_var_back[var.name])
                     dt_pos = min(abs(p.dt), abs(endtime - p.time))
-                    continue
+                    break
                 else:
                     break  # Failure - stop time loop
 
@@ -249,15 +259,18 @@ class Kernel(object):
         # Remove all particles that signalled deletion
         remove_deleted(pset)
 
-        # Idenitify particles that threw errors
+        # Identify particles that threw errors
         error_particles = [p for p in pset.particles
-                           if p.state not in [ErrorCode.Success, ErrorCode.Repeat]]
+                           if p.state != ErrorCode.Success]
         while len(error_particles) > 0:
             # Apply recovery kernel
             for p in error_particles:
-                recovery_kernel = recovery_map[p.state]
-                p.state = ErrorCode.Success
-                recovery_kernel(p, self.fieldset, p.time, dt)
+                if p.state == ErrorCode.Repeat:
+                    p.state = ErrorCode.Success
+                else:
+                    recovery_kernel = recovery_map[p.state]
+                    p.state = ErrorCode.Success
+                    recovery_kernel(p, self.fieldset, p.time, dt)
 
             # Remove all particles that signalled deletion
             remove_deleted(pset)
@@ -269,7 +282,7 @@ class Kernel(object):
                 self.execute_python(pset, endtime, dt)
 
             error_particles = [p for p in pset.particles
-                               if p.state not in [ErrorCode.Success, ErrorCode.Repeat]]
+                               if p.state != ErrorCode.Success]
 
     def merge(self, kernel):
         funcname = self.funcname + kernel.funcname


Test patch:
diff --git a/tests/test_kernel_execution.py b/tests/test_kernel_execution.py
index ef973b76..21049b51 100644
--- a/tests/test_kernel_execution.py
+++ b/tests/test_kernel_execution.py
@@ -245,3 +245,20 @@ def test_update_kernel_in_script(fieldset, mode):
     pset.execute(pset.Kernel(MoveEast), endtime=1., dt=1.)
     pset.execute(pset.Kernel(MoveWest), endtime=2., dt=1.)
     assert np.allclose([p.lon for p in pset], 0.3, rtol=1e-5)  # should be 0.5 + 0.1 - 0.3 = 0.3
+
+
+@pytest.mark.parametrize('mode', ['scipy', 'jit'])
+def test_errorcode_repeat(fieldset, mode):
+    def simpleKernel(particle, fieldset, time, dt):
+        if particle.lon > .1 and time < 1.:
+            # if particle.lon is not re-setted before kernel repetition, it will break here
+            return ErrorCode.Error
+        particle.lon += 0.1
+        if dt > 1.49:
+            # dt is used to leave the repetition loop (dt is the only variable not re-setted)
+            return ErrorCode.Success
+        particle.dt += .1
+        return ErrorCode.Repeat
+
+    pset = ParticleSet(fieldset, pclass=ptype[mode], lon=[0.], lat=[0.])
+    pset.execute(pset.Kernel(simpleKernel), endtime=3., dt=1.)