Line data    Source code

       1             : #
       2             : # Compute relative error
       3             : #
       4             : 
       5             : const AY = 1
       6             : const ADY = 1
       7             : 
       8           0 : normcontrol(solver::Solver) = false # TODO: const trait of Solver
       9             : 
      10           0 : function scale(solver::Solver{N, T, S}, absy::SVector{N, T}) where {N, T, S}
      11           0 :     T(AY) * solver.opts.rtol * absy .+ solver.opts.atol
      12             : end
      13             : 
      14      879285 : function scale(solver::Solver{N, T, S},
      15             :                absy::SVector{N, T}, absdy::SVector{N, T},
      16             :                absh::T=one(T)) where {N, T, S}
      17      879285 :     (T(AY) * absy + T(ADY) * absh * absdy) * solver.opts.rtol .+ solver.opts.atol
      18             : end
      19             : 
      20           0 : function scalednorm(solver::Solver{N, T, S}, y::SVector{N, T}) where {N, T, S}
      21           0 :     s = scale(solver, abs.(y))
      22           0 :     scalednorm(solver, y, s)
      23             : end
      24             : 
      25      879285 : function  scalednorm(solver::Solver{N, T, S},
      26             :                      y::SVector{N, T}, scale::SVector{N, T}) where {N, T, S}
      27      879285 :     z = y ./ scale
      28             : 
      29      879285 :     if normcontrol(solver) # TODO: DECISION FIXED BY TYPE !
      30           0 :         norm(z) / N
      31             :     else
      32      879285 :         norm(z, Inf)
      33             :     end
      34             : end
      35             : 
      36      879285 : function relativeerror(solver::Solver{N, T, S},
      37             :                        yerr::SVector{N, T}) where {N, T, S}
      38      879285 :     absy = max.(abs.(solver.y), abs.(solver.ynew))
      39      879285 :     scalednorm(solver, yerr, scale(solver, absy, abs.(solver.dy), solver.absh))
      40             : end
      41             : 
      42      879285 : function setrelativeerror!(solver::Solver{N, T, S},
      43             :                            yerr::SVector{N, T}) where {N, T, S}
      44      879285 :     solver.relerr = relativeerror(solver, yerr)
      45             : end
      46             : 
      47             : #
      48             : # Compute adapted steps size
      49             : #
      50             : 
      51             : const THRESHOLD = 1.1
      52             : 
      53      879285 : isacceptable(solver::Solver{N, T, S}) where {N, T, S} =
      54             :     solver.relerr <= T(THRESHOLD)
      55             : 
      56             : const MINSCALE = 0.2
      57             : 
      58        9571 : function decreasedstepsize(solver::Solver{N, T, S}) where {N, T, S}
      59        9571 :     P = order0(solver.stepper)
      60       19142 :     scale = RHO / solver.relerr^(1 / P) # TODO: [alg] check !!!
      61           0 :     @massert scale < 1
      62        9571 :     solver.absh * max(MINSCALE, scale)
      63             : end
      64             : 
      65             : const MAXSCALE = 5
      66             : const THRESHOLD_INCREASE = 0.5
      67             : const RHO = 0.8
      68             : 
      69      513811 : function increasedstepsize(solver::Solver{N, T, S}) where {N, T, S}
      70      513811 :     if (solver.relerr < THRESHOLD_INCREASE) && (solver.absh < solver.opts.hmax)
      71             :         # Don't increase if absh was capped to hmax.
      72      454950 :         Q = order1(solver.stepper);
      73      909900 :         scale = RHO / solver.relerr ^ (1 / Q) # TODO: [alg] check !!!
      74           0 :         @massert scale/RHO > 1 # TODO: [alg] just cap?
      75      454950 :         min(solver.opts.hmax,
      76             :             solver.absh * max(min(MAXSCALE, scale), T(1)))
      77             :     else
      78       58861 :         solver.absh
      79             :     end
      80             : end
      81             : 
      82             : #
      83             : # Bisection and "final stretch"
      84             : #
      85             : 
      86             : const BISECTION_FACTOR = 0.25
      87             : 
      88             : # TODO: [alg:upstream] favor bisection factor 0.5?!!
      89             : 
      90      404910 : bisectedstepsize(solver::Solver{N, T, S}) where {N, T, S} =
      91             :     solver.absh * BISECTION_FACTOR
      92             : 
      93             : const STRETCH_FINAL = 1.1
      94             : 
      95             : """
      96             :     tnew = reachdestination(solver)
      97             : 
      98             : Ensure that `tfinal` is is hit exactly.
      99             : 
     100             : The factor `STRETCH_FINAL`
     101             : 
     102             : - allows stretching the final step if `>1`, i.e., "overshoot";
     103             :   increases `absh`
     104             : - requires a "safety margin" if `tfinal` is reached with `t+h`
     105             :   but not with `t+STETCH_FINAL*h`, i.e., don't take full step;
     106             :   decreases `hmax`
     107             : 
     108             : Returns new time `t+stepsize()` or `tfinal`; stretching changes `absh`
     109             : """
     110      869953 : function reachdestination(solver::Solver{N, T, S}) where {N, T, S}
     111           0 :     @massert STRETCH_FINAL > 0
     112             : 
     113      869953 :     absh, h = solver.absh, stepsize(solver)
     114      869953 :     t, tfinal = solver.t, solver.tfinal
     115             : 
     116      869953 :     s = max(STRETCH_FINAL, T(1))
     117             : 
     118     1727706 :     (s * absh > abs(tfinal - t)) || return t + h
     119       12200 :     (s >= 1) && return tfinal
     120           0 :     (STRETCH_FINAL * absh < abs(tfinal - t)) || return tfinal
     121             : 
     122           0 :     t + h * min(STRETCH_FINAL, T(0.5))
     123             : end