Merge pull request clawpack#532 from hadjimy/lmmvss_conditions

Revise 3rd-order LMM's step size condition for starting values.

Author: ketch

src/pyclaw/sharpclaw/solver.py

@@ -191,6 +191,7 @@ def __init__(self,riemann_solver=None,claw_package=None):
         self.prev_dt_values = []
         self.prev_dtFE_values = []
         self.check_lmm_cond = False
+        self.lmm_cond = True
 
         # Call general initialization function
         super(SharpClawSolver,self).__init__(riemann_solver,claw_package)
@@ -290,7 +291,7 @@ def step(self,solution,take_one_step,tstart,tend):
             self.dq_dt = self.dq(state) / self.dt
 
         if 'LMM' in self.time_integrator:
-            self.update_saved_values(state,step_index)
+            step_index = self.update_saved_values(state,step_index)
 
         self.get_dt(solution.t,tstart,tend,take_one_step)
 
@@ -448,20 +449,21 @@ def update_saved_values(self,state,step_index):
                 dtFE = self.dt / cfl * self.cfl_max / self.sspcoeff
 
                 if self.time_integrator == 'SSPLMMk3' and self.check_lmm_cond:
-                    self.accept_step = self.check_3rd_ord_cond(state,step_index,dtFE)
-                    if not self.accept_step:
+                    self.lmm_cond = self.check_3rd_ord_cond(state,step_index,dtFE)
+                    if not self.lmm_cond:
+                        self.accept_step = False
                         state.q[:] = self._registers[-1].q
+                        self.dq_dt = self.prev_dq_dt_values[-1]
                         state.t = self._registers[-1].t
                         self.status['numsteps'] -= 1
-                        return
+                        return self.status['numsteps'] + 1
 
                 if step_index <= len(self._registers):  # Startup
                     if self.time_integrator == 'SSPLMMk3':
                         self.prev_dq_dt_values.append(self.dq_dt)
                     self.prev_dt_values.append(self.dt_old)
                     self.prev_dtFE_values.append(dtFE)
                 else:
-                    self.sspcoeff0 = self.sspcoeff 
                     if self.time_integrator == 'SSPLMMk3':
                         self.prev_dq_dt_values = self.prev_dq_dt_values[1:] + self.prev_dq_dt_values[:1]
                         self.prev_dq_dt_values[-1] = self.dq_dt
@@ -483,6 +485,8 @@ def update_saved_values(self,state,step_index):
             self._registers[-1].q[:] = state.q
             self._registers[-1].t = state.t
 
+        return step_index
+
 
     def check_3rd_ord_cond(self,state,step_index,dtFE):
         r"""
@@ -493,19 +497,19 @@ def check_3rd_ord_cond(self,state,step_index,dtFE):
         If the conditions are violated we muct retrieve the previous solution and discard
         that step; otherwise the step is accepted.
         """
-        accept_step = True
+        lmm_cond = True
 
         if step_index <= len(self._registers):
             rho = 0.6 if len(self._registers) == 4 else 0.57
             if self.dt > rho * dtFE:
-                accept_step = False
+                lmm_cond = False
 
         rhoFE = 0.9 if len(self._registers) == 4 else 0.962
         dtFEm1 = self.prev_dtFE_values[-1]
         if rhoFE * dtFEm1 > dtFE or dtFE > dtFEm1 / rhoFE:
-            accept_step = False
+            lmm_cond = False
 
-        return accept_step
+        return lmm_cond
 
 
     def _set_mthlim(self):
@@ -652,18 +656,20 @@ def get_dt_new(self):
                 # Step-size update of starting methods
                 sspcoeff_ratio = self.sspcoeff0/self.sspcoeff
                 self.dt = sspcoeff_ratio * self.dt * self.cfl_desired / cfl
-                if self.time_integrator == 'SSPLMMk3' and self.check_lmm_cond:
+                if self.time_integrator == 'SSPLMMk3' and self.check_lmm_cond and not self.lmm_cond:
                     rho = 0.6 if len(self._registers)== 4 else 0.57
                     self.dt = rho * self.dt
             else:
+                # Step size selection guarantees CFL condition is satified.
+                # Only need to check 3rd-order LMM's condition
                 if self.accept_step:
                     s = len(self._registers)
                     p = int(self.time_integrator[-1])
                     mu = min([self.prev_dtFE_values[i] for i in range(s)])
                     H = sum(self.prev_dt_values[1:])
                     self.dt = H * mu / (H + (p-1)*mu)
                 elif self.time_integrator == 'SSPLMMk3' and self.check_lmm_cond:
-                    self.dt = 0.5*self.dt
+                    self.dt = 0.5 * self.dt
 
         # Step-size update for RK methods
         else: