try_solve_zhao_monotonic_scalar Subroutine

public subroutine try_solve_zhao_monotonic_scalar(p, branch, x, success)

Type-B/C の定常電流式から ambient density を消去し、電位だけをbracketする。

Arguments

Type IntentOptional Attributes Name
type(zhao_params_type), intent(in) :: p
character(len=1), intent(in) :: branch
real(kind=dp), intent(out) :: x(2)
logical, intent(out) :: success

Calls

proc~~try_solve_zhao_monotonic_scalar~~CallsGraph proc~try_solve_zhao_monotonic_scalar try_solve_zhao_monotonic_scalar proc~evaluate_monotonic_stationary_phi evaluate_monotonic_stationary_phi proc~try_solve_zhao_monotonic_scalar->proc~evaluate_monotonic_stationary_phi proc~residual_norm residual_norm proc~try_solve_zhao_monotonic_scalar->proc~residual_norm proc~zhao_residuals_type_b zhao_residuals_type_b proc~try_solve_zhao_monotonic_scalar->proc~zhao_residuals_type_b proc~zhao_residuals_type_c zhao_residuals_type_c proc~try_solve_zhao_monotonic_scalar->proc~zhao_residuals_type_c proc~evaluate_monotonic_stationary_phi->proc~zhao_residuals_type_b proc~evaluate_monotonic_stationary_phi->proc~zhao_residuals_type_c proc~swe_free_current_term swe_free_current_term proc~evaluate_monotonic_stationary_phi->proc~swe_free_current_term proc~zhao_residuals_type_b->proc~swe_free_current_term proc~zhao_residuals_type_c->proc~swe_free_current_term

Called by

proc~~try_solve_zhao_monotonic_scalar~~CalledByGraph proc~try_solve_zhao_monotonic_scalar try_solve_zhao_monotonic_scalar proc~try_solve_zhao_branch_b try_solve_zhao_branch_b proc~try_solve_zhao_branch_b->proc~try_solve_zhao_monotonic_scalar proc~try_solve_zhao_branch_c try_solve_zhao_branch_c proc~try_solve_zhao_branch_c->proc~try_solve_zhao_monotonic_scalar proc~solve_zhao_branch_b solve_zhao_branch_b proc~solve_zhao_branch_b->proc~try_solve_zhao_branch_b proc~solve_zhao_branch_c solve_zhao_branch_c proc~solve_zhao_branch_c->proc~try_solve_zhao_branch_c proc~try_solve_zhao_unknowns try_solve_zhao_unknowns proc~try_solve_zhao_unknowns->proc~try_solve_zhao_branch_b proc~try_solve_zhao_unknowns->proc~try_solve_zhao_branch_c proc~evaluate_surface_current_model evaluate_surface_current_model proc~evaluate_surface_current_model->proc~try_solve_zhao_unknowns proc~solve_zhao_unknowns solve_zhao_unknowns proc~solve_zhao_unknowns->proc~try_solve_zhao_unknowns proc~evaluate_surface_closure evaluate_surface_closure proc~evaluate_surface_closure->proc~evaluate_surface_current_model

Source Code

  subroutine try_solve_zhao_monotonic_scalar(p, branch, x, success)
    type(zhao_params_type), intent(in) :: p
    character(len=1), intent(in) :: branch
    real(dp), intent(out) :: x(2)
    logical, intent(out) :: success

    integer :: iteration
    real(dp) :: phi_left, phi_right, phi_mid, residual_left, residual_right, residual_mid
    real(dp) :: density_left, density_right, density_mid, voltage_scale, phi_limit
    real(dp) :: residual(2), residual_scale
    logical :: left_ok, right_ok, mid_ok

    x = 0.0_dp
    success = .false.
    voltage_scale = max(p%t_phe_ev, p%t_swe_ev, 1.0_dp)
    residual_scale = max(p%n_swi_inf_m3, p%n_phe0_m3, 1.0_dp)
    select case (branch)
    case ('B')
      phi_left = 128.0_dp*epsilon(1.0_dp)*voltage_scale
      phi_right = voltage_scale
      phi_limit = 100.0_dp*voltage_scale
    case ('C')
      phi_left = -voltage_scale
      phi_right = -128.0_dp*epsilon(1.0_dp)*voltage_scale
      phi_limit = -100.0_dp*voltage_scale
    case default
      return
    end select

    call evaluate_monotonic_stationary_phi(p, branch, phi_left, residual_left, density_left, left_ok)
    call evaluate_monotonic_stationary_phi(p, branch, phi_right, residual_right, density_right, right_ok)
    if (.not. left_ok .or. .not. right_ok) return
    do iteration = 1, 16
      if (residual_left == 0.0_dp .or. residual_right == 0.0_dp .or. &
          sign(1.0_dp, residual_left) /= sign(1.0_dp, residual_right)) exit
      if (branch == 'B') then
        phi_right = min(phi_limit, 2.0_dp*phi_right)
        call evaluate_monotonic_stationary_phi( &
          p, branch, phi_right, residual_right, density_right, right_ok &
          )
        if (.not. right_ok .or. phi_right >= phi_limit) exit
      else
        phi_left = max(phi_limit, 2.0_dp*phi_left)
        call evaluate_monotonic_stationary_phi( &
          p, branch, phi_left, residual_left, density_left, left_ok &
          )
        if (.not. left_ok .or. phi_left <= phi_limit) exit
      end if
    end do
    if (.not. left_ok .or. .not. right_ok) return
    if (residual_left /= 0.0_dp .and. residual_right /= 0.0_dp) then
      if (sign(1.0_dp, residual_left) == sign(1.0_dp, residual_right)) return
    end if

    if (residual_left == 0.0_dp) then
      x = [phi_left, density_left]
    else if (residual_right == 0.0_dp) then
      x = [phi_right, density_right]
    else
      do iteration = 1, 160
        phi_mid = phi_left + 0.5_dp*(phi_right - phi_left)
        call evaluate_monotonic_stationary_phi( &
          p, branch, phi_mid, residual_mid, density_mid, mid_ok &
          )
        if (.not. mid_ok) return
        if (residual_mid == 0.0_dp) then
          phi_left = phi_mid
          phi_right = phi_mid
          density_left = density_mid
          density_right = density_mid
          exit
        end if
        if (sign(1.0_dp, residual_left) /= sign(1.0_dp, residual_mid)) then
          phi_right = phi_mid
          residual_right = residual_mid
          density_right = density_mid
        else
          phi_left = phi_mid
          residual_left = residual_mid
          density_left = density_mid
        end if
        if (abs(phi_right - phi_left) <= &
            256.0_dp*epsilon(1.0_dp)*max(1.0_dp, abs(phi_left), abs(phi_right))) exit
      end do
      phi_mid = phi_left + 0.5_dp*(phi_right - phi_left)
      call evaluate_monotonic_stationary_phi( &
        p, branch, phi_mid, residual_mid, density_mid, mid_ok &
        )
      if (.not. mid_ok) return
      x = [phi_mid, density_mid]
    end if

    if (branch == 'B') then
      call zhao_residuals_type_b(p, x, residual)
    else
      call zhao_residuals_type_c(p, x, residual)
    end if
    success = all(ieee_is_finite(x)) .and. all(ieee_is_finite(residual)) .and. &
              x(2) > 0.0_dp .and. residual_norm(residual) <= &
              max(nonlinear_tol, 1024.0_dp*epsilon(1.0_dp)*residual_scale)
  end subroutine try_solve_zhao_monotonic_scalar