Add SSP Linear Multistep Methods with Variable Step Size

szaghi · szaghi · commit 60ac876d6a59 · 2017-03-31T16:03:34.000+02:00
Add Strong Stability Preserving Linear Multistep Methods with Variable Step Size.

Stable release, fully barckward compatible.
diff --git a/papers/announcement/errors_analysis/errors_analysis.sh b/papers/announcement/errors_analysis/errors_analysis.sh
@@ -11,6 +11,7 @@ euler=0
 lf=0
 lf_raw=0
 lmm_ssp=0
+lmm_ssp_vss=0
 ls_rk=0
 tvd_rk=0
 
@@ -20,18 +21,19 @@ function print_usage {
   echo "script automotion for FOODIE's errors analysis"
   echo "Usage: `basename $0` #solver"
   echo "Valid solver switches are:"
-  echo "  -ab      # => adams-bashforth"
-  echo "  -abm     # => adams-bashforth-moulton"
-  echo "  -am      # => adams-moulton"
-  echo "  -bdf     # => backward differentiation formula"
-  echo "  -emd-rk  # => embeded RK"
-  echo "  -euler   # => euler"
-  echo "  -lf      # => leapfrog"
-  echo "  -lf-raw  # => leapfrog raw filtered"
-  echo "  -lmm-ssp # => linear multistep methods SSP"
-  echo "  -ls-rk   # => low storage RK"
-  echo "  -tvd-rk  # => TVD RK"
-  echo "  -all     # => all solvers!"
+  echo "  -ab          # => adams-bashforth"
+  echo "  -abm         # => adams-bashforth-moulton"
+  echo "  -am          # => adams-moulton"
+  echo "  -bdf         # => backward differentiation formula"
+  echo "  -emd-rk      # => embeded RK"
+  echo "  -euler       # => euler"
+  echo "  -lf          # => leapfrog"
+  echo "  -lf-raw      # => leapfrog raw filtered"
+  echo "  -lmm-ssp     # => linear multistep methods SSP"
+  echo "  -lmm-ssp-vss # => linear multistep methods SSP with Variable Step Size"
+  echo "  -ls-rk       # => low storage RK"
+  echo "  -tvd-rk      # => TVD RK"
+  echo "  -all         # => all solvers!"
 }
 
 #parsing command line
@@ -64,6 +66,9 @@ while [ $# -gt 0 ]; do
     "-lmm-ssp")
       lmm_ssp=1
       ;;
+    "-lmm-ssp-vss")
+      lmm_ssp_vss=1
+      ;;
     "-ls-rk")
       ls_rk=1
       ;;
@@ -149,6 +154,10 @@ fi
 if [ $lmm_ssp -eq 1 ] ; then
   ../tests/oscillation -r -tf 1e6 --solver lmm-ssp --ss 3 5 -Dt 5000.d0 2500.d0 1250.d0 625.d0 320.d0 100.d0 50.d0 10.d0 --errors_analysis > analysis-lmm-ssp.log
 fi
+if [ $lmm_ssp_vss -eq 1 ] ; then
+  ../tests/oscillation -r -tf 1e6 --solver lmm-ssp-vss --ss 3 5 --order 2 -Dt 5000.d0 2500.d0 1250.d0 625.d0 320.d0 100.d0 --errors_analysis > analysis-lmm-ssp-vss.log
+  ../tests/oscillation -r -tf 1e6 --solver lmm-ssp-vss --ss 3 5 --order 3 -Dt 5000.d0 2500.d0 1250.d0 625.d0 320.d0 100.d0 --errors_analysis > analysis-lmm-ssp-vss.log
+fi
 if [ $ls_rk -eq 1 ] ; then
   ../tests/oscillation -r -tf 1e6 --solver ls-runge-kutta --ss 1 -Dt 5000.d0 2500.d0 1250.d0 625.d0 320.d0 100.d0 --errors_analysis #> analysis-leapfrog.log
   ../tests/oscillation -r -tf 1e6 --solver ls-runge-kutta --ss 5 7 -Dt 5000.d0 2500.d0 1250.d0 625.d0 320.d0 100.d0 --errors_analysis #> analysis-leapfrog.log
diff --git a/src/lib/foodie.f90 b/src/lib/foodie.f90
@@ -72,6 +72,7 @@ module foodie
 use foodie_integrator_euler_explicit, only : integrator_euler_explicit
 use foodie_integrator_leapfrog, only : integrator_leapfrog
 use foodie_integrator_lmm_ssp, only : integrator_lmm_ssp
+use foodie_integrator_lmm_ssp_vss, only : integrator_lmm_ssp_vss
 use foodie_integrator_runge_kutta_emd, only : integrator_runge_kutta_emd
 use foodie_integrator_runge_kutta_low_storage, only : integrator_runge_kutta_ls
 use foodie_integrator_runge_kutta_tvd, only : integrator_runge_kutta_tvd
@@ -87,13 +88,14 @@ module foodie
 public :: integrator_euler_explicit
 public :: integrator_leapfrog
 public :: integrator_lmm_ssp
+public :: integrator_lmm_ssp_vss
 public :: integrator_runge_kutta_emd
 public :: integrator_runge_kutta_ls
 public :: foodie_integrator
 public :: integrator_runge_kutta_tvd
 
 contains
-  function foodie_integrator(scheme, steps, stages, tolerance, nu, alpha) result(integrator)
+  function foodie_integrator(scheme, steps, stages, order, tolerance, nu, alpha) result(integrator)
   !< Return a concrete instance of [[integrator]] given a scheme selection.
   !<
   !< This is the FOODIE integrators factory.
@@ -102,6 +104,7 @@ function foodie_integrator(scheme, steps, stages, tolerance, nu, alpha) result(i
   character(*), intent(in)              :: scheme     !< Selected integrator given.
   integer(I_P), intent(in), optional    :: steps      !< Number of time steps used in multi-step schemes.
   integer(I_P), intent(in), optional    :: stages     !< Number of Runge-Kutta stages used in multi-stage schemes.
+  integer(I_P), intent(in), optional    :: order      !< Order of accuracy.
   real(R_P),    intent(in), optional    :: tolerance  !< Tolerance on the local truncation error.
   real(R_P),    intent(in), optional    :: nu         !< Williams-Robert-Asselin filter coefficient.
   real(R_P),    intent(in), optional    :: alpha      !< Robert-Asselin filter coefficient.
@@ -176,6 +179,18 @@ function foodie_integrator(scheme, steps, stages, tolerance, nu, alpha) result(i
                                     error_message='missing steps number for initializing integrator!', &
                                     is_severe=.true.)
     endif
+  case('lmm_ssp_vss')
+    allocate(integrator_lmm_ssp_vss :: integrator)
+    if ((.not.present(steps)).or.(.not.present(order))) then
+      call integrator%trigger_error(error=ERROR_MISSING_STEPS_NUMBER,                                  &
+                                    error_message='missing steps number for initializing integrator!', &
+                                    is_severe=.true.)
+    else
+      select type(integrator)
+      type is(integrator_lmm_ssp_vss)
+        call integrator%init(steps=steps, order=order)
+      endselect
+    endif
   case('runge_kutta_emd')
     allocate(integrator_runge_kutta_emd :: integrator)
     if (present(stages)) then
diff --git a/src/lib/foodie_error_codes.f90 b/src/lib/foodie_error_codes.f90
@@ -11,9 +11,11 @@ module foodie_error_codes
 public :: ERROR_BAD_STEPS_NUMBER
 public :: ERROR_MISSING_STAGES_NUMBER
 public :: ERROR_BAD_STAGES_NUMBER
+public :: ERROR_INTEGRATOR_INIT_FAIL
 
 integer(I_P), parameter :: ERROR_MISSING_STEPS_NUMBER  = 1 !< Error missing steps number (for FOODIE factory).
 integer(I_P), parameter :: ERROR_BAD_STEPS_NUMBER      = 2 !< Error bad (unsupported) steps number.
 integer(I_P), parameter :: ERROR_MISSING_STAGES_NUMBER = 3 !< Error missing stages number (for FOODIE factory).
 integer(I_P), parameter :: ERROR_BAD_STAGES_NUMBER     = 4 !< Error bad (unsupported) stages number.
+integer(I_P), parameter :: ERROR_INTEGRATOR_INIT_FAIL  = 5 !< Error integrator initialization failed.
 endmodule foodie_error_codes
diff --git a/src/lib/foodie_integrator_lmm_ssp.f90 b/src/lib/foodie_integrator_lmm_ssp.f90
@@ -115,39 +115,39 @@ subroutine init(self, steps)
   if (self%is_supported(steps)) then
     call self%destroy
     self%steps = steps
-    allocate(self%a(1:steps)) ; self%a = 0.0_R_P
-    allocate(self%b(1:steps)) ; self%b = 0.0_R_P
+    allocate(self%a(1:steps)) ; self%a = 0._R_P
+    allocate(self%b(1:steps)) ; self%b = 0._R_P
     select case(steps)
     case(3) ! LMM-SSP(3,2)
-      self%a(1) = 1.0_R_P/4.0_R_P
-      self%a(2) = 0.0_R_P
-      self%a(3) = 3.0_R_P/4.0_R_P
+      self%a(1) = 1._R_P/4._R_P
+      self%a(2) = 0._R_P
+      self%a(3) = 3._R_P/4._R_P
 
-      self%b(1) = 0.0_R_P
-      self%b(2) = 0.0_R_P
-      self%b(3) = 3.0_R_P/2.0_R_P
+      self%b(1) = 0._R_P
+      self%b(2) = 0._R_P
+      self%b(3) = 3._R_P/2._R_P
     case(4) ! LMM-SSP(4,3)
-      self%a(1) = 11.0_R_P/27.0_R_P
-      self%a(2) = 0.0_R_P
-      self%a(3) = 0.0_R_P
-      self%a(4) = 16.0_R_P/27.0_R_P
-
-      self%b(1) = 12.0_R_P/27.0_R_P
-      self%b(2) = 0.0_R_P
-      self%b(3) = 0.0_R_P
-      self%b(4) = 16.0_R_P/9.0_R_P
+      self%a(1) = 11._R_P/27._R_P
+      self%a(2) = 0._R_P
+      self%a(3) = 0._R_P
+      self%a(4) = 16._R_P/27._R_P
+
+      self%b(1) = 12._R_P/27._R_P
+      self%b(2) = 0._R_P
+      self%b(3) = 0._R_P
+      self%b(4) = 16._R_P/9._R_P
     case(5) ! LMM-SSP(5,3)
-      self%a(1) = 7.0_R_P/32.0_R_P
-      self%a(2) = 0.0_R_P
-      self%a(3) = 0.0_R_P
-      self%a(4) = 0.0_R_P
-      self%a(5) = 25.0_R_P/32.0_R_P
-
-      self%b(1) = 5.0_R_P/16.0_R_P
-      self%b(2) = 0.0_R_P
-      self%b(3) = 0.0_R_P
-      self%b(4) = 0.0_R_P
-      self%b(5) = 25.0_R_P/16.0_R_P
+      self%a(1) = 7._R_P/32._R_P
+      self%a(2) = 0._R_P
+      self%a(3) = 0._R_P
+      self%a(4) = 0._R_P
+      self%a(5) = 25._R_P/32._R_P
+
+      self%b(1) = 5._R_P/16._R_P
+      self%b(2) = 0._R_P
+      self%b(3) = 0._R_P
+      self%b(4) = 0._R_P
+      self%b(5) = 25._R_P/16._R_P
     endselect
   else
     call self%trigger_error(error=ERROR_BAD_STEPS_NUMBER,                           &
@@ -170,7 +170,9 @@ subroutine integrate(self, U, previous, Dt, t, autoupdate)
   autoupdate_ = .true. ; if (present(autoupdate)) autoupdate_ = autoupdate
   U = U * 0._R_P
   do s=1, self%steps
-    U = U + previous(s) * self%a(s) + previous(s)%t(t=t(s)) * (Dt * self%b(s))
+    if (self%a(s) /= 0._R_P) U = U + previous(s) * self%a(s)
+    if (self%b(s) /= 0._R_P) U = U + previous(s)%t(t=t(s)) * (Dt * self%b(s))
+    ! U = U + previous(s) * self%a(s) + previous(s)%t(t=t(s)) * (Dt * self%b(s))
   enddo
   if (autoupdate_) call self%update_previous(U=U, previous=previous)
   endsubroutine integrate
diff --git a/src/lib/foodie_integrator_lmm_ssp_vss.f90 b/src/lib/foodie_integrator_lmm_ssp_vss.f90
diff --git a/src/tests/accuracy/oscillation/oscillation.f90 b/src/tests/accuracy/oscillation/oscillation.f90
