PB3D [2.47]
Ideal linear high-n MHD stability in 3-D
Loading...
Searching...
No Matches
num_ops.f90
Go to the documentation of this file.
1!------------------------------------------------------------------------------!
2!> Numerical operations.
3!------------------------------------------------------------------------------!
4module num_ops
5#include <PB3D_macros.h>
6 use str_utilities
7 use messages
8 use num_vars, only: dp, iu, max_str_ln, pi
9 use num_utilities
10 use output_ops
11
12 implicit none
13 private
14 public calc_zero_hh, calc_zero_zhang
15#if ldebug
16 public debug_calc_zero
17#endif
18
19 ! global variables
20#if ldebug
21 !> \ldebug
22 logical :: debug_calc_zero = .false. !< plot debug information for calc_zero
23#endif
24
25 ! interfaces
26
27 !> \public Finds the zero of a function using Householder iteration.
28 !!
29 !! If something goes wrong, by default multiple tries can be attempted, by
30 !! backtracking on the correction by multiplying it by a relaxation factor.
31 !! This can be done \c max_nr_backtracks times.
32 !!
33 !! If still nothing is achieved, an error message is returned,
34 !! that is empty otherwise.
35 interface calc_zero_hh
36 !> \public
37 module procedure calc_zero_hh_0d
38 !> \public
39 module procedure calc_zero_hh_3d
40 end interface
41
42contains
43 !> \private 0-D version
44 !! \param[inout] fun <tt>fun(x,ord)</tt> with
45 !! - <tt>x</tt> abscissa
46 !! - <tt>ord</tt> order of derivative
47 !! - <tt>fun</tt> ordinate
48 function calc_zero_hh_0d(zero,fun,ord,guess,max_nr_backtracks,output) &
49 &result(err_msg)
50 use num_vars, only: max_it_zero, tol_zero, max_nr_backtracks_hh
51
52 ! input / output
53 real(dp), intent(inout) :: zero !< output
54 interface
55 function fun(x,ord) ! the function
56 use num_vars, only: dp
57 real(dp), intent(in) :: x
58 integer, intent(in) :: ord
59 real(dp) :: fun
60 end function fun
61 end interface
62 integer, intent(in) :: ord !< order of solution
63 real(dp), intent(in) :: guess !< first guess
64 integer, intent(in), optional :: max_nr_backtracks !< max nr. of tries with different relaxation factors
65 logical, intent(in), optional :: output !< give output on convergence
66 character(len=max_str_ln) :: err_msg !< possible error message
67
68 ! local variables
69 integer :: id, jd, kd ! counters
70 integer :: max_nr_backtracks_loc ! local max_nr_backtracks
71 logical :: output_loc ! local output
72 logical :: relaxed_enough ! whether step was relaxed enough to get a smaller new value
73 real(dp) :: zero_new ! new zero calculated
74 real(dp) :: fun_new ! function value at new zero
75 real(dp) :: corr ! correction
76 real(dp), allocatable :: fun_vals(:) ! function values
77#if ldebug
78 real(dp), allocatable :: corrs(:) ! corrections for all steps
79 real(dp), allocatable :: values(:) ! values for all steps
80#endif
81
82 ! initialize error message
83 err_msg = ''
84
85 ! set up local output
86 output_loc = .false.
87 if (present(output)) output_loc = output
88
89 ! set up zero
90 zero = guess
91
92 ! set up local max_nr_backtracks
93 max_nr_backtracks_loc = max_nr_backtracks_hh
94 if (present(max_nr_backtracks)) &
95 &max_nr_backtracks_loc = max_nr_backtracks
96
97 ! initialize function values
98 allocate(fun_vals(0:ord))
99
100#if ldebug
101 if (debug_calc_zero) then
102 ! set up corrs
103 allocate(corrs(max_it_zero))
104 allocate(values(max_it_zero))
105 end if
106#endif
107
108 ! loop over iterations
109 hh: do jd = 1,max_it_zero
110 ! calculate function values
111 do kd = 0,ord
112 fun_vals(kd) = fun(zero,kd)
113 end do
114
115 ! correction to theta_HH
116 select case (ord)
117 case (1) ! Newton-Rhapson
118 corr = -fun_vals(0)/fun_vals(1)
119 case (2) ! Halley
120 corr = -fun_vals(0)*fun_vals(1)/&
121 &(fun_vals(1)**2 - 0.5_dp*&
122 &fun_vals(0)*fun_vals(2))
123 case (3) ! 3rd order
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))
129 end select
130#if ldebug
131 if (debug_calc_zero) then
132 corrs(jd) = corr
133 values(jd) = zero
134 end if
135#endif
136 relaxed_enough = .false.
137 do id = 1,max_nr_backtracks_loc
138 zero_new = zero + corr ! propose a new zero
139
140 fun_new = fun(zero_new,0) ! calculate function value there
141
142 if (abs(fun_new)-abs(fun_vals(0)).le.0._dp) then ! error has decreased
143 relaxed_enough = .true.
144 zero = zero_new
145 exit
146 end if
147 corr = corr*0.5_dp
148 end do
149
150 ! check for relaxation
151 if (.not.relaxed_enough) then
152 err_msg = trim(i2str(max_nr_backtracks_hh))//' backtracks were &
153 &not sufficient to get a better guess for the zero'
154 zero = 0.0_dp
155 else
156 ! output
157 if (output_loc) call writo(trim(i2str(jd))//' / '&
158 &//trim(i2str(max_it_zero))//&
159 &' - relative error: '//&
160 &trim(r2str(abs(corr)))//' (tol: '//&
161 &trim(r2str(tol_zero))//')')
162
163 ! check for convergence
164 if (abs(corr).lt.tol_zero) then
165#if ldebug
166 if (debug_calc_zero) &
167 &call plot_evolution(corrs(1:jd),values(1:jd))
168#endif
169 return
170 else if (jd .eq. max_it_zero) then
171#if ldebug
172 if (debug_calc_zero) then
173 call plot_evolution(corrs,values)
174 deallocate(corrs)
175 deallocate(values)
176 end if
177#endif
178 err_msg = 'Not converged after '//trim(i2str(jd))//&
179 &' iterations, with residual '//&
180 &trim(r2strt(corr))//' and final value '//&
181 &trim(r2strt(zero))
182 zero = 0.0_dp
183 end if
184 end if
185 end do hh
186#if ldebug
187 contains
188 ! plots corrections
189 !> \private
190 subroutine plot_evolution(corrs,values)
191 ! input / output
192 real(dp), intent(in) :: corrs(:) ! corrections
193 real(dp), intent(in) :: values(:) ! values
194
195 ! local variables
196 real(dp) :: min_x, max_x ! min. and max. x for which to plot the function
197 real(dp) :: extra_x = 1._dp ! how much bigger to take the plot interval than just min and max
198 real(dp) :: delta_x ! interval between zero and guess
199 integer :: n_x = 200 ! how many x values to plot
200 real(dp), allocatable :: x_out(:,:) ! output x of plot
201 real(dp), allocatable :: y_out(:,:) ! output y of plot
202
203 ! set up min and max x
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
209
210 ! set up output of plot
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),&
216 &jd=1,n_x)]
217 y_out(:,2) = [(fun(min_x+(max_x-min_x)*(jd-1._dp)/(n_x-1),1),&
218 &jd=1,n_x)]
219
220 ! plot function
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)
224
225 ! user output
226 call writo('Last correction '//trim(r2strt(corrs(size(corrs))))&
227 &//' with tolerance '//trim(r2strt(tol_zero)))
228
229 ! plot corrections and values
230 call print_ex_2d('corrs','',corrs)
231 call print_ex_2d('values','',values)
232 end subroutine plot_evolution
233#endif
234 end function calc_zero_hh_0d
235 !> \private 3-D version
236 !! \param[inout] fun <tt>fun(dims,x,ord)</tt> with
237 !! - <tt>dims(3)</tt> dimension of abscissa and ordinate
238 !! - <tt>x(dims(1),dims(2),dims(3))</tt> abscissa
239 !! - <tt>ord</tt> order of derivative
240 !! - <tt>fun(dims(1),dims(2),dims(3))</tt> ordinate
241 function calc_zero_hh_3d(dims,zero,fun,ord,guess,max_nr_backtracks,output) &
242 &result(err_msg)
243 use num_vars, only: max_it_zero, tol_zero, max_nr_backtracks_hh
244
245 ! input / output
246 integer, intent(in) :: dims(3) !< dimensions of the problem
247 real(dp), intent(inout) :: zero(dims(1),dims(2),dims(3)) !< output
248 interface
249 function fun(dims,x,ord) ! the function
250 use num_vars, only: dp
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))
255 end function fun
256 end interface
257 integer, intent(in) :: ord !< order of solution
258 real(dp), intent(in) :: guess(dims(1),dims(2),dims(3)) !< first guess
259 integer, intent(in), optional :: max_nr_backtracks !< max nr. of backtracks
260 logical, intent(in), optional :: output !< give output on convergence
261 character(len=max_str_ln) :: err_msg !< possible error message
262
263 ! local variables
264 integer :: id, jd, kd ! counters
265 integer :: max_nr_backtracks_loc ! local max_nr_backtracks
266 integer :: mc_ind(3) ! index of maximum correction
267 logical :: output_loc ! local output
268 logical :: relaxed_enough ! relaxed enough
269 real(dp) :: zero_new(dims(1),dims(2),dims(3)) ! new zeros calculated
270 real(dp) :: fun_new(dims(1),dims(2),dims(3)) ! function values at new zero
271 real(dp) :: corr(dims(1),dims(2),dims(3)) ! correction
272 real(dp), allocatable :: fun_vals(:,:,:,:) ! function values
273#if ldebug
274 real(dp), allocatable :: corrs(:,:,:,:) ! corrections for all steps
275 real(dp), allocatable :: values(:,:,:,:) ! values for all steps
276#endif
277
278 ! initialize error message
279 err_msg = ''
280
281 ! test
282 if (ord.lt.1 .or. ord.gt.3) then
283 zero = 0.0_dp
284 err_msg = 'only orders 1 (Newton-Rhapson), 2 (Halley) and 3 &
285 &implemented, not '//trim(i2str(ord))
286 return
287 end if
288
289 ! set up local output
290 output_loc = .false.
291 if (present(output)) output_loc = output
292
293 ! set up zero
294 zero = guess
295
296 ! set up local max_nr_backtracks
297 max_nr_backtracks_loc = max_nr_backtracks_hh
298 if (present(max_nr_backtracks)) &
299 &max_nr_backtracks_loc = max_nr_backtracks
300
301 ! initialize function values
302 allocate(fun_vals(dims(1),dims(2),dims(3),0:ord))
303
304#if ldebug
305 if (debug_calc_zero) then
306 ! set up corrs
307 allocate(corrs(dims(1),dims(2),dims(3),max_it_zero))
308 allocate(values(dims(1),dims(2),dims(3),max_it_zero))
309 end if
310#endif
311
312 ! loop over iterations
313 hh: do jd = 1,max_it_zero
314 ! calculate function values
315 do kd = 0,ord
316 fun_vals(:,:,:,kd) = fun(dims,zero,kd)
317 end do
318
319 ! correction to theta_HH
320 select case (ord)
321 case (1) ! Newton-Rhapson
322 corr = -fun_vals(:,:,:,0)/fun_vals(:,:,:,1)
323 case (2) ! Halley
324 corr = -fun_vals(:,:,:,0)*fun_vals(:,:,:,1)/&
325 &(fun_vals(:,:,:,1)**2 - 0.5_dp*&
326 &fun_vals(:,:,:,0)*fun_vals(:,:,:,2))
327 case (3) ! 3rd order
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))
333 end select
334#if ldebug
335 if (debug_calc_zero) then
336 corrs(:,:,:,jd) = corr
337 values(:,:,:,jd) = zero
338 end if
339#endif
340 relaxed_enough = .false.
341 do id = 1,max_nr_backtracks_loc
342 zero_new = zero ! otherwise some values might stay unitialized
343 where (abs(corr).gt.tol_zero) zero_new = zero + corr ! propose a new zero
344
345 fun_new = fun(dims,zero_new,0) ! calculate function value there
346
347#if ldebug
348 if (debug_calc_zero) &
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)))))//&
352 &' <= 0?')
353#endif
354 if (maxval(abs(fun_new)-abs(fun_vals(:,:,:,0))).le.0._dp) then ! error has decreased
355 relaxed_enough = .true.
356 zero = zero_new
357 exit
358 end if
359 corr = corr*0.5_dp
360 end do
361
362 ! check for relaxation
363 if (.not.relaxed_enough) then
364 err_msg = trim(i2str(max_nr_backtracks_hh))//' backtracks were &
365 &not sufficient to get a better guess for the zero'
366 zero = 0.0_dp
367 exit
368 else
369 ! output
370 if (output_loc) call writo(trim(i2str(jd))//' / '&
371 &//trim(i2str(max_it_zero))//&
372 &' - maximum relative error: '//&
373 &trim(r2str(maxval(abs(corr))))//' (tol: '//&
374 &trim(r2str(tol_zero))//')')
375
376 ! check for convergence
377 if (maxval(abs(corr)).lt.tol_zero) then
378#if ldebug
379 if (debug_calc_zero) call &
380 &plot_evolution(corrs(:,:,:,1:jd),values(:,:,:,1:jd))
381#endif
382 return
383 else if (jd .eq. max_it_zero) then
384#if ldebug
385 if (debug_calc_zero) then
386 call plot_evolution(corrs,values)
387 deallocate(corrs)
388 deallocate(values)
389 end if
390#endif
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))))
397 zero = 0.0_dp
398 end if
399 end if
400 end do hh
401#if ldebug
402 contains
403 ! plots corrections
404 !> \private
405 subroutine plot_evolution(corrs,values)
406 ! input / output
407 real(dp), intent(in) :: corrs(:,:,:,:) ! corrections
408 real(dp), intent(in) :: values(:,:,:,:) ! values
409
410 ! local variables
411 integer :: n_corrs ! number of corrections
412 character(len=max_str_ln) :: var_name(1) ! names of variable
413 character(len=max_str_ln) :: file_name ! name of file
414
415 ! initialize
416 n_corrs = size(corrs,4)
417 mc_ind = maxloc(abs(corrs(:,:,:,n_corrs)))
418
419 ! user output
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)))//&
424 &' with tolerance '//trim(r2strt(tol_zero)))
425
426 ! plot corrs
427 file_name = 'corrs'
428 var_name = 'var'
429
430 call plot_hdf5(var_name,file_name,corrs,col=2,&
431 &descr='corrections')
432
433 ! plot values
434 file_name = 'values'
435 var_name = 'var'
436
437 call plot_hdf5(var_name,file_name,values,col=2,&
438 &descr='values')
439 end subroutine plot_evolution
440#endif
441 end function calc_zero_hh_3d
442
443 !> Finds the zero of a function using Zhang's method, which is simpler than
444 !! Brent's method.
445 !!
446 !! Taken from from Steven Stage's correction of Zhang's paper
447 !! \cite zhang2011improvement.
448 !!
449 !! Unlike Householder, Zhang's method needs an interval \c x_int_in to work
450 !! in, not a guess. Also, it does not require the derivative of the
451 !! function.
452 !!
453 !! The routine returns an error message if no zero is found, and which is
454 !! empty otherwise.
455 !!
456 !! \note The interval \c x_int_in needs to be so that the function values at
457 !! either end are of different value.
458 !! \param[inout] fun <tt>fun(x)</tt> with
459 !! - <tt>x</tt> abscissa
460 !! - <tt>fun</tt> ordinate
461 function calc_zero_zhang(zero,fun,x_int_in) result(err_msg)
462 use num_vars, only: max_it_zero, tol_zero
463
464 ! input / output
465 real(dp), intent(inout) :: zero !< output
466 interface
467 function fun(x) ! the function
468 use num_vars, only: dp
469 real(dp), intent(in) :: x
470 real(dp) :: fun
471 end function fun
472 end interface
473 real(dp), intent(in) :: x_int_in(2) !< interval
474 character(len=max_str_ln) :: err_msg !< possible error message
475
476 ! local variables
477 logical :: converged ! whether converged
478 integer :: jd ! counter
479 integer :: id_loc ! local id (either 1 or 2)
480 real(dp) :: x_int(2) ! points of current interval [a,b]
481 real(dp) :: x_mid ! x at midpoint [c]
482 real(dp) :: x_rule ! x found with rule [s]
483 real(dp) :: x_temp ! temporary x
484 real(dp) :: fun_int(2) ! function at x_int
485 real(dp) :: fun_mid ! function at x_mid
486 real(dp) :: fun_rule ! function at x_rule
487#if ldebug
488 real(dp), allocatable :: x_ints(:,:) ! bounds for all steps
489#endif
490
491 ! initialize error message and zero
492 err_msg = ''
493 zero = 0.0_dp
494
495 ! initialize from input
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))]
499
500 ! test whether already zero
501 if (abs(fun_int(1)).lt.tol_zero) then
502 zero = x_int(1)
503 return
504 end if
505 if (abs(fun_int(2)).lt.tol_zero) then
506 zero = x_int(2)
507 return
508 end if
509
510 ! test whether different sign
511 if (product(fun_int).gt.0) then
512 err_msg = 'Function values on starting interval have same sign'
513 return
514 end if
515
516 ! intialize converged
517 converged = .false.
518
519#if ldebug
520 if (debug_calc_zero) then
521 ! set up x_ints
522 allocate(x_ints(max_it_zero,2))
523 end if
524#endif
525
526 ! iterations
527 do jd = 1,max_it_zero
528 ! calculate midpoint and function value
529 x_mid = 0.5_dp*sum(x_int)
530 fun_mid = fun(x_mid)
531
532 if (abs(fun_int(1)-fun_mid).gt.tol_zero .and. &
533 &abs(fun_int(2)-fun_mid).gt.tol_zero) then
534 ! three different function values: use inverse quadratic
535 ! interpolation
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))
542
543 if (x_int(1).lt.x_rule .and. x_rule.lt.x_int(2)) then
544 ! found x within interval
545 else
546 ! use bisection instead
547 if (fun_int(1)*fun_mid.lt.0) then ! [a,c] contains sign change
548 id_loc = 1
549 else if (fun_int(2)*fun_mid.lt.0) then ! [b,c] contains sign change
550 id_loc = 2
551 else
552 err_msg = 'Something went wrong...'
553 converged = .true.
554 end if
555 if (.not.converged) x_rule = 0.5_dp*(x_int(id_loc)+x_mid)
556 end if
557 else
558 ! two function values overlap so use secant rule
559 if (fun_int(1)*fun_mid.lt.0) then ! [a,c] contains sign change
560 id_loc = 1
561 else if (fun_int(2)*fun_mid.lt.0) then ! [b,c] contains sign change
562 id_loc = 2
563 else
564 err_msg = 'Something went wrong...'
565 converged = .true.
566 end if
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)
570 end if
571
572 ! calculate fun(s)
573 fun_rule = fun(x_rule)
574
575 ! make sure x_mid < x_rule
576 if (x_mid.gt.x_rule) then
577 x_temp = x_rule
578 x_rule = x_mid
579 x_mid = x_temp
580 end if
581
582 ! check whether [c,s] contains root
583 if (fun_mid*fun_rule.lt.0) then
584 x_int(1) = x_mid
585 x_int(2) = x_rule
586 else
587 if (fun_int(1)*fun_mid.lt.0) then
588 x_int(2) = x_mid
589 else
590 x_int(1) = x_rule
591 end if
592 end if
593
594#if ldebug
595 if (debug_calc_zero) then
596 x_ints(jd,:) = x_int
597 end if
598#endif
599
600 ! check for convergence
601 fun_int = [fun(x_int(1)),fun(x_int(2))]
602 if (abs(fun_int(1)).lt.tol_zero) then
603 zero = x_int(1)
604 converged = .true.
605 end if
606 if (abs(fun_int(2)).lt.tol_zero) then
607 zero = x_int(2)
608 converged = .true.
609 end if
610 if ((x_int(2)-x_int(1)).lt.tol_zero) then
611 zero = 0.5_dp*sum(x_int)
612 converged = .true.
613 end if
614
615 if (converged) then
616#if ldebug
617 if (debug_calc_zero) call plot_evolution(x_ints(1:jd,:))
618#endif
619 return
620 end if
621 end do
622#if ldebug
623 contains
624 ! plots corrections
625 !> \private
626 subroutine plot_evolution(x_ints)
627 ! input / output
628 real(dp), intent(in) :: x_ints(:,:) ! x_intervals
629
630 ! local variables
631 real(dp) :: min_x, max_x ! min. and max. x for which to plot the function
632 integer :: n_x = 200 ! how many x values to plot
633 integer :: n_ints ! nr. of intervals
634 real(dp), allocatable :: x_out(:) ! output x of plot
635 real(dp), allocatable :: y_out(:) ! output y of plot
636
637 ! set up nr. of intervals
638 n_ints = size(x_ints,1)
639
640 ! set up min and max x
641 min_x = x_int_in(1)
642 max_x = x_int_in(2)
643
644 ! set up output of plot
645 allocate(x_out(n_x))
646 allocate(y_out(n_x))
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)]
649
650 ! plot function
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)))
654 call print_ex_2d('function','',y_out,x=x_out)
655
656 ! user output
657 call writo('Final interval ['//trim(r2str(x_ints(n_ints,1)))//&
658 &'..'//trim(r2str(x_ints(n_ints,2)))//'] with tolerance '//&
659 &trim(r2strt(tol_zero)))
660
661 ! plot intervals
662 call print_ex_2d(['x_intervals'],'',x_ints)
663 end subroutine plot_evolution
664#endif
665 end function calc_zero_zhang
666end module num_ops
num_ops::calc_zero_hh
Finds the zero of a function using Householder iteration.
Definition num_ops.f90:35
output_ops::plot_hdf5
Prints variables vars with names var_names in an HDF5 file with name c file_name and accompanying XDM...
Definition output_ops.f90:137
output_ops::print_ex_2d
Print 2-D output on a file.
Definition output_ops.f90:47
messages
Numerical utilities related to giving output.
Definition messages.f90:4
messages::writo
subroutine, public writo(input_str, persistent, error, warning, alert)
Write output to file identified by output_i.
Definition messages.f90:275
num_ops
Numerical operations.
Definition num_ops.f90:4
num_ops::debug_calc_zero
logical, public debug_calc_zero
plot debug information for calc_zero
Definition num_ops.f90:22
num_ops::calc_zero_zhang
character(len=max_str_ln) function, public calc_zero_zhang(zero, fun, x_int_in)
Finds the zero of a function using Zhang's method, which is simpler than Brent's method.
Definition num_ops.f90:462
num_utilities
Numerical utilities.
Definition num_utilities.f90:4
num_vars
Numerical variables used by most other modules.
Definition num_vars.f90:4
num_vars::dp
integer, parameter, public dp
double precision
Definition num_vars.f90:46
num_vars::pi
real(dp), parameter, public pi
Definition num_vars.f90:83
num_vars::iu
complex(dp), parameter, public iu
complex unit
Definition num_vars.f90:85
num_vars::max_str_ln
integer, parameter, public max_str_ln
maximum length of strings
Definition num_vars.f90:50
num_vars::max_nr_backtracks_hh
integer, public max_nr_backtracks_hh
maximum number of backtracks for Householder, relax. factors
Definition num_vars.f90:132
num_vars::max_it_zero
integer, public max_it_zero
maximum number of iterations to find zeros
Definition num_vars.f90:131
num_vars::tol_zero
real(dp), public tol_zero
tolerance for zeros
Definition num_vars.f90:133
output_ops
Operations concerning giving output, on the screen as well as in output files.
Definition output_ops.f90:5
str_utilities
Operations on strings.
Definition str_utilities.f90:4
str_utilities::i2str
elemental character(len=max_str_ln) function, public i2str(k)
Convert an integer to string.
Definition str_utilities.f90:18
str_utilities::r2str
elemental character(len=max_str_ln) function, public r2str(k)
Convert a real (double) to string.
Definition str_utilities.f90:42
str_utilities::r2strt
elemental character(len=max_str_ln) function, public r2strt(k)
Convert a real (double) to string.
Definition str_utilities.f90:54
this documentation was created with Doxygen 1.16.1 by Toon Weyens