solving adapt dt rkck

Author: Diego Vaca

src/common/m_constants.fpp

@@ -66,8 +66,7 @@ module m_constants
     ! Adaptive rkck constants
     real(kind(0d0)), parameter :: verysmall_dt = 1.d-14 !< Very small dt, stop stepper
     real(kind(0d0)), parameter :: SAFETY = 0.9d0 !< Next dt will be maximum 0.9*dt if truncation error is above tolerance.
-    real(kind(0d0)), parameter :: RNDDEC = 1.0d05 !< Need to round the relative truncation error (5th decimal) to avoid the inclusion of random decimals
-    ! after dividing calculated tolerance by a very small number (rkck_tolerance)
+    real(kind(0d0)), parameter :: RNDDEC = 1.0d8 !< Round calculated dt (16th digit) to avoid the inclusion of random decimals
     real(kind(0d0)), parameter :: PSHRNK = -0.25d0 !< Factor to reduce dt when truncation error above tolerance
     real(kind(0d0)), parameter :: SHRNKDT = 0.5d0 !< Factor to reduce dt due to negative bubble radius
     real(kind(0d0)), parameter :: ERRCON = 1.89d-4 !< Limit to slightly increase dt when truncation error is between ERRCON and 1

src/simulation/m_bubbles_EL.fpp

@@ -1354,24 +1354,32 @@ contains
 
             if ((rkck_errmax > 1d0)) then   ! Truncation error too large, reduce dt and restart time step
                 restart_rkck_step = .true.
-                htemp = SAFETY*dt*((floor(rkck_errmax*RNDDEC)/RNDDEC)**PSHRNK)
+                htemp = SAFETY*dt*(rkck_errmax**PSHRNK)
                 dt = sign(max(abs(htemp), (1d0 - SAFETY)*abs(dt)), dt)  ! No more than a factor of 10.
-                if (proc_rank == 0) print *, '>>>>> WARNING: Truncation error found. Reducing dt and restaring time step, now dt: ', dt
             else                            ! Step succeeded. Compute size of next step.
                 if (rkck_errmax > ERRCON) then
-                    dt = SAFETY*dt*((floor(rkck_errmax*RNDDEC)/RNDDEC)**PGROW) ! No more than a factor of 5 increase.
+                    dt = SAFETY*dt*(rkck_errmax**PGROW) ! No more than a factor of 5 increase.
                 else
                     dt = (1d0/SHRNKDT)*dt ! Truncation error too small (< 1.89e-4), increase time step
                 end if
             end if
 
+            !dt precision accuracy is 16 digits
+            dt = (ceiling(dt*RNDDEC)*RNDDEC + ceiling(dt*(RNDDEC**2d0) - ceiling(dt*RNDDEC)*RNDDEC))/(RNDDEC**2d0)
             dt = min(dt, dt_max)
+
+            if (dt < 0) call s_mpi_abort('dt must not be negative')
             if (num_procs > 1) then
                 call s_mpi_allreduce_min(dt, aux_glb)
                 dt = aux_glb
             end if
             !$acc update device(dt)
 
+            if (restart_rkck_step) then
+                if (proc_rank == 0) print '("WARNING: Truncation error found. Restaring time step, and now dt = "ES16.6"")', &
+                    dt
+            end if
+
         end if
 
     end subroutine s_compute_rkck_dt