advance_nonlinear_term_beuler Subroutine

  subroutine advance_nonlinear_term_beuler(g_state, phinew, aparnew, bparnew, dt, &
       fields_func, source_func)
    use theta_grid, only: ntgrid
    use gs2_time, only: code_dt_min
    use mp, only: mp_abort, max_allreduce, proc0
    use run_parameters, only: immediate_reset
    use gs2_layouts, only: g_lo
    use array_utils, only: copy, gs2_max_abs, zero_array
    use warning_helpers, only: exactly_equal
    implicit none
    complex, dimension (-ntgrid:, :, g_lo%llim_proc:), intent(in out) :: g_state
    complex, dimension (-ntgrid:,:,:), intent (in out) :: phinew, aparnew, bparnew
    complex, dimension(:, :, :), allocatable :: g_cur, g_local, gnl_1, gnl_2
    real, intent(in) :: dt
    real :: internal_time, time_step, max_vel
    integer :: counter, iglo
    integer, parameter :: iteration_limit = 20
    real, parameter :: time_adjust_factor = 2
    logical, parameter :: debug = .false.
    real :: inverse_error, error, actual_dt, half_time, cfl_time_step
    real, dimension(3) :: errors
    real, save :: time_step_last = -1
    procedure(invert_field_func) :: fields_func
    procedure(source_term_func) :: source_func
    logical :: reject_step, cfl_limited
    complex, dimension(:, :, :), allocatable :: error_array

    ! Initialise internal variables
    if (time_step_last < 0) time_step_last = dt
    internal_time = 0.0
    time_step = time_step_last

    allocate(error_array(-ntgrid:ntgrid, 2, g_lo%llim_proc:g_lo%ulim_alloc))
    allocate(g_cur(-ntgrid:ntgrid, 2, g_lo%llim_proc:g_lo%ulim_alloc))
    allocate(g_local(-ntgrid:ntgrid, 2, g_lo%llim_proc:g_lo%ulim_alloc))
    allocate(gnl_1(-ntgrid:ntgrid, 2, g_lo%llim_proc:g_lo%ulim_alloc))
    allocate(gnl_2(-ntgrid:ntgrid, 2, g_lo%llim_proc:g_lo%ulim_alloc))
    call copy(g_state, g_cur)
    call zero_array(error_array) ; call zero_array(g_local)
    call zero_array(gnl_1); call zero_array(gnl_2)

    ! Keep on advancing until we have advanced by dt
    do while (internal_time < dt)
       ! Reset counter and error
       counter = 0
       error = 1.0
       errors = 0.0

       ! Limit the actual time step that we take to avoid overshooting (i.e.
       ! to keep internal_time <= dt).
       actual_dt = min(time_step, dt - internal_time)
       half_time = 0.5 * actual_dt

       ! Calculate the current NL source term, for use in error
       ! estimates later.
       call source_func(g_state, gnl_2, phinew, aparnew, bparnew, &
            max_vel, need_to_adjust = .not. advance_nonadiabatic_dfn, &
            calculate_cfl_limit = .false.)
       nrhs_calls = nrhs_calls + 1

       ! Try to take a step of size actual_dt using fixed point iteration with
       ! backwards Euler.
       do while (error > convergence_tolerance .and. counter < iteration_limit)

          ! Get the current nonlinear source term and cfl limit -- note we don't use the
          ! cfl limit here
          call source_func(g_state, gnl_1, phinew, aparnew, bparnew, &
               max_vel, need_to_adjust = .not. advance_nonadiabatic_dfn, &
               calculate_cfl_limit = .false.)
          nrhs_calls = nrhs_calls + 1

          ! Calculate the estimate of the new g
          !$OMP PARALLEL DO DEFAULT(none) &
          !$OMP PRIVATE(iglo) &
          !$OMP SHARED(g_lo, g_local, g_cur, half_time, gnl_1) &
          !$OMP SCHEDULE(static)
          do iglo = g_lo%llim_proc, g_lo%ulim_proc
             g_local(:, :, iglo) = g_cur(:, :, iglo) + half_time * gnl_1(:, :, iglo)
          end do
          !$OMP END PARALLEL DO

          ! Calculate the maximum absolute error and the maximum value of the previous solution

          ! Note we _could_ merge the following loop witth the gs2_max_abs kernel
          ! so that we only loop over the error_array once (and actually don't
          ! need to explicitly store the full array)
          !$OMP PARALLEL DO DEFAULT(none) &
          !$OMP PRIVATE(iglo) &
          !$OMP SHARED(g_lo, error_array, g_state, g_local) &
          !$OMP SCHEDULE(static)
          do iglo = g_lo%llim_proc, g_lo%ulim_proc
             error_array(:, :, iglo) = g_state(:, :, iglo) - g_local(:, :, iglo)
          end do
          !$OMP END PARALLEL DO
          errors(1) = gs2_max_abs(error_array)
          errors(2) = gs2_max_abs(g_local)
          call max_allreduce(errors)

          ! Calculate an approximate relative error
          error = errors(1) / errors(2)

          ! Update internals -- increment counter and update the 'previous' solution with the new one.
          counter = counter + 1
          call copy(g_local, g_state)

          ! Update the fields consistently with the current solution
          call fields_func(g_state, phinew, aparnew, bparnew, local_only = .true.)
       end do

       ! Now we need to decide if the last step was successful or not. There are two exit conditions, either the
       ! error was small enough (good step) or the iteration limit was exceeded (bad step).
       if (counter >= iteration_limit) then
          ! If it was a bad step then throw away these changes and reduce the time step
          if (proc0 .and. debug) print*,'Too many iterations, reducing step from',time_step,'to',time_step/time_adjust_factor,'with internal_time = ',internal_time

          ! Reduce the time step
          time_step = time_step / time_adjust_factor
          ! If the time step is now too small then abort
          if (time_step < code_dt_min) call mp_abort("Time step has become too small in advance_nonlinear_term.", .true.)

          ! If we've taken too many iterations then we want to drop the time step and try
          ! again. To try again we must ensure that gnew has been reset to g, and the fields
          ! are consistent with this.
          reject_step = .true.
       else
          ! If we have met the convergence condition then we should estimate the error
          ! on the solution and see if we need to reject the step or how we should adjust
          ! the step size.
          ! To estimate the error we compare the backward Euler update with an approximate
          ! Crank-Nicolson update.
          ! g_new_beuler = g + 0.5 * actual_dt * gnl_1(g_new_beuler)
          ! g_new_cn = g + 0.5 * actual_dt * (gnl_1(g_new_beuler) + gnl_1(g))/2
          ! abs_error  = |g_new_beuler - g_new_cn|
          !      = 0.5 * actual_dt * 0.5 * (gnl_1(g_new_beuler) - gnl_1(g)
          ! Where gnl_1(h) is the source func evaluated with gnew=h.

          ! Here we calculate gnl_1(g_new_beuler) exactly. We could probably
          ! rely on using the existing value of gnl_1(g_local) as, whilst this is
          ! strictly from the previous iteration of the fixed-point scheme, we know
          ! it is within the convergence tolerance of the true value.
          call source_func(g_state, gnl_1, phinew, aparnew, bparnew, &
               max_vel, need_to_adjust = .not. advance_nonadiabatic_dfn, &
               calculate_cfl_limit = .true.)
          nrhs_calls = nrhs_calls + 1

          ! Note we _could_ merge the following loop witth the gs2_max_abs kernel
          ! so that we only loop over the error_array once (and actually don't
          ! need to explicitly store the full array)
          !$OMP PARALLEL DO DEFAULT(none) &
          !$OMP PRIVATE(iglo) &
          !$OMP SHARED(g_lo, error_array, gnl_1, gnl_2) &
          !$OMP SCHEDULE(static)
          do iglo = g_lo%llim_proc, g_lo%ulim_proc
             error_array(:, :, iglo) = gnl_1(:, :, iglo) - gnl_2(:, :, iglo)
          end do
          !$OMP END PARALLEL DO
          errors(1) = gs2_max_abs(error_array) * (0.5 * half_time)
          errors(3) = max_vel
          call max_allreduce(errors)
          inverse_error = get_inverse_error(errors)
          cfl_time_step = 1.0 / errors(3)

          !Calculate new time step based on error
          time_step = min(actual_dt * sqrt(inverse_error), cfl_time_step)
          cfl_limited = time_step < actual_dt .and. exactly_equal(time_step, cfl_time_step)
          time_step = time_step * time_step_safety_factor

          ! If error is too large reject the step
          if ( inverse_error < 1.0 .or. (cfl_limited .and. immediate_reset)) then
             if (proc0 .and. debug) print*,'Attempted beuler step failed due to error or cfl, adapting time step'
             ! Reset as we're retrying the step.
             reject_step = .true.
          else
             reject_step = .false.
          end if
       end if

       if (reject_step) then
          ! If the time step is now too small then abort
          if (time_step < code_dt_min) call mp_abort("Time step has become too small in advance_nonlinear_term.", .true.)
          call copy(g_cur, g_state)
          ! If we store the initial field then we could perhaps just copy here instead.
          ! This could be faster as it avoids communications involved in the velocity integration.
          call calculate_fields(fields_func, g_state, phinew, aparnew, bparnew)
          nfailed_steps = nfailed_steps + 1
       else
          ! If we have met the error condition and not taken too many steps
          ! then g_state represents our new solution so make sure we update g to reflect this
          internal_time = internal_time + actual_dt
          call copy(g_state, g_cur)
          nsuccessful_steps = nsuccessful_steps + 1
       end if

    end do

    ! Store the final time step size so we can start from this in the next callg
    time_step_last = time_step

    if(proc0.and.debug) print*,'Explicit done in ',counter,'steps with final error',inverse_error,'and step',time_step,'final time',internal_time

  end subroutine advance_nonlinear_term_beuler