SteadyState solver robustness updates

[sjdaines]

New options for SteadyState.solve_ptc (ie apply to each outer timestep):

• SteadyState.RestartSmallValuesCallback that provides a callback function with signature

  rsvc(state, tmodel, deltat, model, modeldata)
  

  that modifies state to reset variables with values < modify_threshold to value modify_val, hence restarting the Newton solve from modify_val on the next iteration.

  Add to step_callbacks list eg

  restart_small_values_callback = PALEOmodel.SteadyState.RestartSmallValuesCallback(
      modeldata.solver_view_all.stateexplicit;
      exclude_variables=["atm.N_mr", "atm.CH21_mr"], # exclude variables that do actually reach values < threshold
      modify_threshold = mr_threshold,
      modify_val = 1e-30,
  )
  
  PALEOmodel.SteadyState.steadystate_ptcForwardDiff(
      ...
      step_callbacks = [restart_small_values_callback, ],
      ...
  )
  

• CheckValuesCallback (provides additional conditions to accept a timestep), provides a callback function with signature

  cvc(state, tmodel, deltat, model, modeldata)::Bool
  

  that checks variables have check_min_threshold < values < check_max_threshold and rejects timestep if not.

  Add to ``check_callbacks` list, eg

  check_small_values_callback = PALEOmodel.SteadyState.CheckValuesCallback(
    modeldata.solver_view_all.stateexplicit;
    exclude_variables=["atm.N_mr"], # exclude variables that do actually reach values < threshold
    check_min_value = 2*mr_threshold,
    verbose=false,
  )  
  
  PALEOmodel.SteadyState.steadystate_ptcForwardDiff(
      ...
      check_callbacks = [check_small_values_callback, ],    
      ...
  )
  

New options for use with NLsolve.jl ie apply to each Newton iteration (add to solvekwargs argument):

• SolverFunctions.StepClampMultAll! callback (restricts Newton step size to multiple or fraction mr_maxstep of current value, as well as projecting new value to bounded region mr_min - mr_max)
  Requires https://github.com/PALEOtoolkit/NLsolve.jl project_region branch, with new option apply_step! to NLsolve:

  PALEOmodel.SteadyState.steadystate_ptcForwardDiff(
      ...
      solvekwargs = (apply_step! =  SolverFunctions.StepClampMultAll!(mr_threshold, mr_max, mr_maxstep), ...)
      ...
  )
  

• Add sparse linear solvers adapted to NLsolve interface (linsolve argument, supplied via solvekwargs)

   PALEOmodel.SteadyState.steadystate_ptcForwardDiff(
      ...
      solvekwargs = (linsolve = PALEOmodel.SolverFunctions.SparseLinsolveUMFPACK(), ...) # UMFPACK with iterative refinement reenabled
      ...
  )
  
  PALEOmodel.SteadyState.steadystate_ptcForwardDiff(
      ...
      solvekwargs = (linsolve = PALEOmodel.SolverFunctions.SparseLinsolveSparspak64x2() , ...) # high accuracy quad precision lu factorisation
      ...
  )