|
! ================================================================================================================================== |
|
! |
|
! [reference] |
|
! https://www.slideshare.net/HideoHirose/homotopy-67988638 |
|
! |
|
! ================================================================================================================================== |
|
|
|
module mod_solver |
|
|
|
! <module>s to import |
|
use, intrinsic :: iso_fortran_env |
|
use, intrinsic :: ieee_arithmetic |
|
use, non_intrinsic :: mod_utility |
|
use, non_intrinsic :: mod_eqn_func |
|
|
|
|
|
! require all variables to be explicitly declared |
|
implicit none |
|
|
|
|
|
! <type>s for this <module> |
|
|
|
type typ_num_trial |
|
integer(INT32) :: max |
|
integer(INT32) :: end |
|
end type typ_num_trial |
|
|
|
type typ_equation |
|
real(REAL64) :: self |
|
real(REAL64) :: drv |
|
end type typ_equation |
|
|
|
type typ_solution |
|
real(REAL64) :: self |
|
real(REAL64) :: error |
|
end type typ_solution |
|
|
|
type typ_trial |
|
type(typ_equation) :: equation |
|
type(typ_solution) :: solution |
|
end type typ_trial |
|
|
|
type typ_result |
|
type(typ_num_trial) , private :: num_trial |
|
type(typ_trial) , private , allocatable :: trial(:) |
|
end type typ_result |
|
|
|
|
|
! constant(s) for this <module> |
|
integer(INT32) , parameter , private :: num_trial_begin = 0_INT32 |
|
|
|
|
|
! accessibility of the <subroutine>s and <function>s in this <module> |
|
|
|
! kind: function |
|
private :: eval_convergence |
|
private :: get_num_trial_end |
|
|
|
! kind: subroutine |
|
private :: eval_allowable_error |
|
private :: eval_num_tiral_max |
|
private :: eval_solution_init |
|
private :: initialize_instance |
|
private :: save_result |
|
public :: solver_newton |
|
private :: update_solution_newton |
|
|
|
|
|
! contained <subroutine>s and <function>s are below |
|
contains |
|
|
|
|
|
pure function eval_convergence (trial_next, trial_pre, allowable_error) result(stat) |
|
|
|
! argument(s) for this <subroutine> |
|
type(typ_trial) , intent(in) :: trial_next |
|
type(typ_trial) , intent(in) :: trial_pre |
|
real(REAL64) , intent(in) :: allowable_error |
|
|
|
! return value of this <function> |
|
logical :: stat |
|
|
|
stat = ieee_is_nan( trial_next%solution%self ) |
|
|
|
if (stat) then |
|
return |
|
else |
|
stat = ( trial_next%solution%error .le. allowable_error * abs(trial_pre%solution%self) ) |
|
end if |
|
|
|
return |
|
|
|
end function eval_convergence |
|
|
|
|
|
pure function get_num_trial_end (instance) result(num_trial_end) |
|
|
|
! argument(s) for this <function> |
|
type(typ_result) , intent(in) :: instance |
|
|
|
! return value of this <function> |
|
integer(INT32) :: num_trial_end |
|
|
|
num_trial_end = instance%num_trial%end |
|
|
|
end function |
|
|
|
|
|
subroutine eval_allowable_error (val) |
|
|
|
! argument(s) for this <subroutine> |
|
real(REAL64) , intent(in) :: val |
|
|
|
if ( .not. ieee_is_finite(val) ) then |
|
write(unit= ERROR_UNIT, fmt='(A)') 'The allowable error of the solution must be a finite IEEE value.' |
|
stop |
|
else if ( val .ge. 1.0_REAL64 ) then |
|
write(unit= ERROR_UNIT, fmt='(A)') 'The allowable error of the solution must be less then 1.' |
|
stop |
|
end if |
|
|
|
return |
|
|
|
end subroutine eval_allowable_error |
|
|
|
|
|
subroutine eval_num_tiral_max (val) |
|
|
|
! argument(s) for this <subroutine> |
|
integer(INT32) , intent(in) :: val |
|
|
|
select case (val) |
|
case(:0) |
|
write(unit= ERROR_UNIT, fmt='(A)') 'The number of steps to update the solution must be greater than zero.' |
|
stop |
|
case default |
|
return |
|
end select |
|
|
|
end subroutine eval_num_tiral_max |
|
|
|
|
|
subroutine eval_solution_init (val) |
|
|
|
! argument(s) for this <subroutine> |
|
real(REAL64) , intent(in) :: val |
|
|
|
if ( .not. ieee_is_finite(val) ) then |
|
write(unit= ERROR_UNIT, fmt='(A)') 'The initial condition of the solution must be a finite IEEE value.' |
|
stop |
|
end if |
|
|
|
return |
|
|
|
end subroutine eval_solution_init |
|
|
|
|
|
subroutine initialize_instance (instance, num_trial_max, solution_init) |
|
|
|
! argument(s) for this <subroutine> |
|
type(typ_result) , intent(inout) :: instance |
|
integer(INT32) , intent(in) :: num_trial_max |
|
real(REAL64) , intent(in) :: solution_init |
|
|
|
|
|
! variable(s) for this <subroutine> |
|
integer :: buffer_stat |
|
type(typ_equation) :: init_equation |
|
type(typ_solution) :: init_solution |
|
type(typ_trial) :: init_trial |
|
|
|
|
|
! STEP.01 |
|
! prepare support instances to initialize the target instance |
|
init_equation % self = ieee_value(x= solution_init, class= ieee_quiet_nan) |
|
init_equation % drv = ieee_value(x= solution_init, class= ieee_quiet_nan) |
|
init_solution % self = ieee_value(x= solution_init, class= ieee_quiet_nan) |
|
init_solution % error = ieee_value(x= solution_init, class= ieee_quiet_nan) |
|
init_trial % equation = init_equation |
|
init_trial % solution = init_solution |
|
|
|
|
|
! STEP.02 |
|
! try to allocate the member arrays |
|
allocate( instance%trial(num_trial_begin:num_trial_max), stat= buffer_stat, source= init_trial ) |
|
call eval_stat_allocate(stat= buffer_stat, name= 'typ_result%trial(:)' ) |
|
|
|
|
|
! STEP.03 |
|
! initialize the member variables |
|
instance % num_trial % max = num_trial_max |
|
instance % num_trial % end = num_trial_max |
|
instance % trial(num_trial_begin) % solution % self = solution_init |
|
instance % trial(num_trial_begin) % equation % self = eqn_f ( solution_init ) |
|
instance % trial(num_trial_begin) % equation % drv = eqn_df ( solution_init ) |
|
|
|
|
|
! STEP.END |
|
return |
|
|
|
end subroutine initialize_instance |
|
|
|
|
|
subroutine save_result (unit, instance) |
|
|
|
! argument(s) for this <subroutine> |
|
integer , intent(in) :: unit |
|
type(typ_result) , intent(in) :: instance |
|
|
|
! variable(s) for this <subroutine> |
|
integer(INT32) :: iter_trial |
|
|
|
do iter_trial = num_trial_begin , get_num_trial_end(instance) , 1_INT32 |
|
write(unit= unit, fmt='(I3,4(1X,ES24.16E3))') &! |
|
iter_trial , &! |
|
instance%trial(iter_trial)%solution%self , &! |
|
instance%trial(iter_trial)%solution%error , &! |
|
instance%trial(iter_trial)%equation%self , &! |
|
instance%trial(iter_trial)%equation%drv |
|
end do |
|
|
|
end subroutine save_result |
|
|
|
|
|
subroutine solver_newton (file, solution_init, num_trial_max, allowable_error) |
|
|
|
! argument(s) for this <subroutine> |
|
character(len=*) , intent(in) :: file ! file path to save the solution |
|
real(REAL64) , intent(in) :: solution_init ! initial condition of the newton method |
|
integer(INT32) , intent(in) :: num_trial_max ! number of trials to update the numerical solution |
|
real(REAL64) , intent(in) :: allowable_error ! allowable error of the numerical solution |
|
|
|
|
|
! variable(s) (scalar) for this <subroutine> |
|
integer :: buffer_stat ! unit number to save the result |
|
integer :: save_unit ! unit number to save the result |
|
integer(INT32) :: iter_trial ! iterator |
|
integer(INT32) :: indx_trial_pre ! index to store the pre step' one |
|
|
|
|
|
! variable(s) (vector) for this <subroutine> |
|
type(typ_result) :: result |
|
|
|
|
|
! STEP.01 |
|
! check given arguments |
|
call eval_allowable_error (val= allowable_error) |
|
call eval_num_tiral_max (val= num_trial_max) |
|
call eval_solution_init (val= solution_init) |
|
|
|
|
|
! STEP.02 |
|
! initialize the instance to store the result |
|
call initialize_instance (instance= result, num_trial_max= num_trial_max, solution_init= solution_init) |
|
|
|
|
|
! STEP.03 |
|
! open a file to save the result |
|
open(&! |
|
action = 'write' , &! |
|
file = file , &! |
|
form = 'formatted' , &! |
|
iostat = buffer_stat , &! |
|
newunit = save_unit , &! |
|
status = 'replace' &! |
|
) |
|
|
|
call eval_iostat_open(iostat= buffer_stat, file= file) |
|
|
|
|
|
! STEP.04 |
|
! compute the numerical solution using Newton's method |
|
do iter_trial = (num_trial_begin + 1_INT32) , num_trial_max , 1_INT32 |
|
|
|
indx_trial_pre = iter_trial - 1_INT32 |
|
|
|
call update_solution_newton( trial_next= result%trial(iter_trial), trial_pre= result%trial(indx_trial_pre) ) |
|
|
|
if ( eval_convergence(trial_next= result%trial(iter_trial), trial_pre= result%trial(indx_trial_pre), allowable_error= allowable_error) ) then |
|
result%num_trial%end = iter_trial |
|
exit |
|
end if |
|
|
|
end do |
|
|
|
|
|
! STEP.05 |
|
! save the computed result |
|
call save_result(unit= save_unit, instance= result) |
|
|
|
|
|
! STEP.06 |
|
! close the file to save the result |
|
close(unit= save_unit, status='keep') |
|
|
|
|
|
! STEP.END |
|
return |
|
|
|
end subroutine solver_newton |
|
|
|
|
|
subroutine update_solution_newton (trial_next, trial_pre) |
|
|
|
! argument(s) for this <subroutine> |
|
type(typ_trial) , intent(inout) :: trial_next |
|
type(typ_trial) , intent(in) :: trial_pre |
|
|
|
trial_next%solution%error = trial_pre%equation%self / trial_pre%equation%drv |
|
trial_next%solution%self = trial_pre%solution%self - trial_next%solution%error |
|
trial_next%solution%error = abs( trial_next%solution%error ) |
|
trial_next%equation%self = eqn_f ( trial_next%solution%self ) |
|
trial_next%equation%drv = eqn_df ( trial_next%solution%self ) |
|
|
|
end subroutine |
|
|
|
end module mod_solver |
|
|
|
! ================================================================================================================================== |
|
! EOF |
|
! ================================================================================================================================== |
