5 #include <PB3D_macros.h>
37 module procedure calc_zero_hh_0d
39 module procedure calc_zero_hh_3d
48 function calc_zero_hh_0d(zero,fun,ord,guess,max_nr_backtracks,output) &
53 real(dp),
intent(inout) :: zero
57 real(
dp),
intent(in) :: x
58 integer,
intent(in) :: ord
62 integer,
intent(in) :: ord
63 real(dp),
intent(in) :: guess
64 integer,
intent(in),
optional :: max_nr_backtracks
65 logical,
intent(in),
optional :: output
66 character(len=max_str_ln) :: err_msg
70 integer :: max_nr_backtracks_loc
72 logical :: relaxed_enough
76 real(dp),
allocatable :: fun_vals(:)
78 real(dp),
allocatable :: corrs(:)
79 real(dp),
allocatable :: values(:)
87 if (
present(output)) output_loc = output
94 if (
present(max_nr_backtracks)) &
95 &max_nr_backtracks_loc = max_nr_backtracks
98 allocate(fun_vals(0:ord))
112 fun_vals(kd) = fun(zero,kd)
118 corr = -fun_vals(0)/fun_vals(1)
120 corr = -fun_vals(0)*fun_vals(1)/&
121 &(fun_vals(1)**2 - 0.5_dp*&
122 &fun_vals(0)*fun_vals(2))
124 corr = -(6*fun_vals(0)*fun_vals(1)**2-&
125 &3*fun_vals(0)**2*fun_vals(2))/&
126 &(6*fun_vals(1)**3 - 6*fun_vals(0)*&
127 &fun_vals(1)*fun_vals(2)+&
128 &fun_vals(0)**2*fun_vals(3))
136 relaxed_enough = .false.
137 do id = 1,max_nr_backtracks_loc
138 zero_new = zero + corr
140 fun_new = fun(zero_new,0)
142 if (abs(fun_new)-abs(fun_vals(0)).le.0._dp)
then
143 relaxed_enough = .true.
151 if (.not.relaxed_enough)
then
153 ¬ sufficient to get a better guess for the zero'
157 if (output_loc)
call writo(trim(
i2str(jd))//
' / '&
159 &
' - relative error: '//&
160 &trim(
r2str(abs(corr)))//
' (tol: '//&
167 &
call plot_evolution(corrs(1:jd),values(1:jd))
173 call plot_evolution(corrs,values)
178 err_msg =
'Not converged after '//trim(
i2str(jd))//&
179 &
' iterations, with residual '//&
180 &trim(
r2strt(corr))//
' and final value '//&
190 subroutine plot_evolution(corrs,values)
192 real(dp),
intent(in) :: corrs(:)
193 real(dp),
intent(in) :: values(:)
196 real(dp) :: min_x, max_x
197 real(dp) :: extra_x = 1._dp
200 real(dp),
allocatable :: x_out(:,:)
201 real(dp),
allocatable :: y_out(:,:)
204 min_x = min(zero,guess)
205 max_x = max(zero,guess)
206 delta_x = max_x-min_x
207 min_x = min_x - delta_x*extra_x
208 max_x = max_x + delta_x*extra_x
211 allocate(x_out(n_x,2))
212 allocate(y_out(n_x,2))
213 x_out(:,1) = [(min_x+(max_x-min_x)*(jd-1._dp)/(n_x-1),jd=1,n_x)]
214 x_out(:,2) = [(min_x+(max_x-min_x)*(jd-1._dp)/(n_x-1),jd=1,n_x)]
215 y_out(:,1) = [(fun(min_x+(max_x-min_x)*(jd-1._dp)/(n_x-1),0),&
217 y_out(:,2) = [(fun(min_x+(max_x-min_x)*(jd-1._dp)/(n_x-1),1),&
221 call writo(
'Initial guess: '//trim(
r2str(guess)))
222 call writo(
'Resulting zero: '//trim(
r2str(zero)))
223 call print_ex_2d([
'function ',
'derivative'],
'',y_out,x=x_out)
226 call writo(
'Last correction '//trim(
r2strt(corrs(
size(corrs))))&
232 end subroutine plot_evolution
234 end function calc_zero_hh_0d
241 function calc_zero_hh_3d(dims,zero,fun,ord,guess,max_nr_backtracks,output) &
246 integer,
intent(in) :: dims(3)
247 real(dp),
intent(inout) :: zero(dims(1),dims(2),dims(3))
249 function fun(dims,x,ord)
251 integer,
intent(in) :: dims(3)
252 real(
dp),
intent(in) :: x(dims(1),dims(2),dims(3))
253 integer,
intent(in) :: ord
254 real(
dp) :: fun(dims(1),dims(2),dims(3))
257 integer,
intent(in) :: ord
258 real(dp),
intent(in) :: guess(dims(1),dims(2),dims(3))
259 integer,
intent(in),
optional :: max_nr_backtracks
260 logical,
intent(in),
optional :: output
261 character(len=max_str_ln) :: err_msg
264 integer :: id, jd, kd
265 integer :: max_nr_backtracks_loc
267 logical :: output_loc
268 logical :: relaxed_enough
269 real(dp) :: zero_new(dims(1),dims(2),dims(3))
270 real(dp) :: fun_new(dims(1),dims(2),dims(3))
271 real(dp) :: corr(dims(1),dims(2),dims(3))
272 real(dp),
allocatable :: fun_vals(:,:,:,:)
274 real(dp),
allocatable :: corrs(:,:,:,:)
275 real(dp),
allocatable :: values(:,:,:,:)
282 if (ord.lt.1 .or. ord.gt.3)
then
284 err_msg =
'only orders 1 (Newton-Rhapson), 2 (Halley) and 3 &
285 &implemented, not '//trim(
i2str(ord))
291 if (
present(output)) output_loc = output
298 if (
present(max_nr_backtracks)) &
299 &max_nr_backtracks_loc = max_nr_backtracks
302 allocate(fun_vals(dims(1),dims(2),dims(3),0:ord))
307 allocate(corrs(dims(1),dims(2),dims(3),
max_it_zero))
308 allocate(values(dims(1),dims(2),dims(3),
max_it_zero))
316 fun_vals(:,:,:,kd) = fun(dims,zero,kd)
322 corr = -fun_vals(:,:,:,0)/fun_vals(:,:,:,1)
324 corr = -fun_vals(:,:,:,0)*fun_vals(:,:,:,1)/&
325 &(fun_vals(:,:,:,1)**2 - 0.5_dp*&
326 &fun_vals(:,:,:,0)*fun_vals(:,:,:,2))
328 corr = -(6*fun_vals(:,:,:,0)*fun_vals(:,:,:,1)**2-&
329 &3*fun_vals(:,:,:,0)**2*fun_vals(:,:,:,2))/&
330 &(6*fun_vals(:,:,:,1)**3 - 6*fun_vals(:,:,:,0)*&
331 &fun_vals(:,:,:,1)*fun_vals(:,:,:,2)+&
332 &fun_vals(:,:,:,0)**2*fun_vals(:,:,:,3))
336 corrs(:,:,:,jd) = corr
337 values(:,:,:,jd) = zero
340 relaxed_enough = .false.
341 do id = 1,max_nr_backtracks_loc
343 where (abs(corr).gt.
tol_zero) zero_new = zero + corr
345 fun_new = fun(dims,zero_new,0)
349 &
call writo(
'backtrack '//trim(
i2str(id))//
' of max. '//&
350 &trim(
i2str(max_nr_backtracks_loc))//
' - criterion: '//&
351 &trim(
r2str(maxval(abs(fun_new)-abs(fun_vals(:,:,:,0)))))//&
354 if (maxval(abs(fun_new)-abs(fun_vals(:,:,:,0))).le.0._dp)
then
355 relaxed_enough = .true.
363 if (.not.relaxed_enough)
then
365 ¬ sufficient to get a better guess for the zero'
370 if (output_loc)
call writo(trim(
i2str(jd))//
' / '&
372 &
' - maximum relative error: '//&
373 &trim(
r2str(maxval(abs(corr))))//
' (tol: '//&
377 if (maxval(abs(corr)).lt.
tol_zero)
then
380 &plot_evolution(corrs(:,:,:,1:jd),values(:,:,:,1:jd))
386 call plot_evolution(corrs,values)
391 mc_ind = maxloc(abs(corr))
392 err_msg =
'Not converged after '//trim(
i2str(jd))//&
393 &
' iterations, with maximum residual '//&
394 &trim(
r2strt(corr(mc_ind(1),mc_ind(2),mc_ind(3))))&
395 &//
' and final value '//trim(
r2strt(&
396 &zero(mc_ind(1),mc_ind(2),mc_ind(3))))
405 subroutine plot_evolution(corrs,values)
407 real(dp),
intent(in) :: corrs(:,:,:,:)
408 real(dp),
intent(in) :: values(:,:,:,:)
412 character(len=max_str_ln) :: var_name(1)
413 character(len=max_str_ln) :: file_name
416 n_corrs =
size(corrs,4)
417 mc_ind = maxloc(abs(corrs(:,:,:,n_corrs)))
420 call writo(
'Last maximum correction '//&
421 &trim(
r2strt(corrs(mc_ind(1),mc_ind(2),mc_ind(3),n_corrs)))&
422 &//
' and final value '//trim(
r2strt(&
423 &values(mc_ind(1),mc_ind(2),mc_ind(3),n_corrs)))//&
430 call plot_hdf5(var_name,file_name,corrs,col=2,&
431 &descr=
'corrections')
437 call plot_hdf5(var_name,file_name,values,col=2,&
439 end subroutine plot_evolution
441 end function calc_zero_hh_3d
465 real(dp),
intent(inout) :: zero
469 real(
dp),
intent(in) :: x
473 real(dp),
intent(in) :: x_int_in(2)
474 character(len=max_str_ln) :: err_msg
484 real(dp) :: fun_int(2)
488 real(dp),
allocatable :: x_ints(:,:)
496 x_int(1) = minval(x_int_in)
497 x_int(2) = maxval(x_int_in)
498 fun_int = [fun(x_int(1)),fun(x_int(2))]
501 if (abs(fun_int(1)).lt.
tol_zero)
then
505 if (abs(fun_int(2)).lt.
tol_zero)
then
511 if (product(fun_int).gt.0)
then
512 err_msg =
'Function values on starting interval have same sign'
529 x_mid = 0.5_dp*sum(x_int)
532 if (abs(fun_int(1)-fun_mid).gt.
tol_zero .and. &
533 &abs(fun_int(2)-fun_mid).gt.
tol_zero)
then
536 x_rule = fun_mid*fun_int(2)*x_int(1)/&
537 &((fun_int(1)-fun_mid)*(fun_int(1)-fun_int(2))) + &
538 &fun_int(1)*fun_int(2)*x_mid/&
539 &((fun_mid-fun_int(1))*(fun_mid-fun_int(2))) + &
540 &fun_int(1)*fun_mid*x_int(2)/&
541 &((fun_int(2)-fun_mid)*(fun_int(2)-fun_mid))
543 if (x_int(1).lt.x_rule .and. x_rule.lt.x_int(2))
then
547 if (fun_int(1)*fun_mid.lt.0)
then
549 else if (fun_int(2)*fun_mid.lt.0)
then
552 err_msg =
'Something went wrong...'
555 if (.not.converged) x_rule = 0.5_dp*(x_int(id_loc)+x_mid)
559 if (fun_int(1)*fun_mid.lt.0)
then
561 else if (fun_int(2)*fun_mid.lt.0)
then
564 err_msg =
'Something went wrong...'
567 if (.not.converged) x_rule = x_int(id_loc) - &
568 &fun_int(id_loc)*(x_int(id_loc)-x_mid)/&
569 &(fun_int(id_loc)-fun_mid)
573 fun_rule = fun(x_rule)
576 if (x_mid.gt.x_rule)
then
583 if (fun_mid*fun_rule.lt.0)
then
587 if (fun_int(1)*fun_mid.lt.0)
then
601 fun_int = [fun(x_int(1)),fun(x_int(2))]
602 if (abs(fun_int(1)).lt.
tol_zero)
then
606 if (abs(fun_int(2)).lt.
tol_zero)
then
610 if ((x_int(2)-x_int(1)).lt.
tol_zero)
then
611 zero = 0.5_dp*sum(x_int)
626 subroutine plot_evolution(x_ints)
628 real(dp),
intent(in) :: x_ints(:,:)
631 real(dp) :: min_x, max_x
634 real(dp),
allocatable :: x_out(:)
635 real(dp),
allocatable :: y_out(:)
638 n_ints =
size(x_ints,1)
647 x_out = [(min_x+(max_x-min_x)*(jd-1._dp)/(n_x-1),jd=1,n_x)]
648 y_out = [(fun(min_x+(max_x-min_x)*(jd-1._dp)/(n_x-1)),jd=1,n_x)]
651 call writo(
'Initial interval: ['//trim(
r2str(x_int_in(1)))//&
652 &
'..'//trim(
r2str(x_int_in(2)))//
']')
653 call writo(
'Resulting zero: '//trim(
r2str(zero)))
657 call writo(
'Final interval ['//trim(
r2str(x_ints(n_ints,1)))//&
658 &
'..'//trim(
r2str(x_ints(n_ints,2)))//
'] with tolerance '//&
663 end subroutine plot_evolution