Context Navigation

wave_equations_integrator.F90 @ 594

Revision 594, 24.0 KB checked in by micael, 2 years ago (diff)
The extra MGGA term in the Hamiltonian was not correct. Changed the way how this term is applied.
Property svn:eol-style set to `native` Property svn:keywords set to `Id`

Line
1	!! Copyright (C) 2004-2011 M. Oliveira, F. Nogueira
2	!!
3	!! This program is free software; you can redistribute it and/or modify
4	!! it under the terms of the GNU General Public License as published by
5	!! the Free Software Foundation; either version 2, or (at your option)
6	!! any later version.
7	!!
8	!! This program is distributed in the hope that it will be useful,
9	!! but WITHOUT ANY WARRANTY; without even the implied warranty of
10	!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
11	!! GNU General Public License for more details.
12	!!
13	!! You should have received a copy of the GNU General Public License
14	!! along with this program; if not, write to the Free Software
15	!! Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
16	!! 02111-1307, USA.
17	!!
18	!! $Id$
19
20	#include "global.h"
21
22	module wave_equations_integrator_m
23	use global_m
24	use messages_m
25	use io_m
26	use gsl_interface_m
27	use oct_parser_m
28	use quantum_numbers_m
29	use potentials_m
30	use wave_equations_derivs_m
31	implicit none
32
33
34	!---Interfaces---!
35
36	interface
37	function ode_solve(x2, stp, evl, ctrl, h1, hmin, nstepmax, nstep, x, y1, y2, myfunc)
38	use global_m
39	implicit none
40	integer(POINTER_SIZE), intent(inout) :: stp, evl, ctrl
41	integer, intent(in) :: nstepmax
42	integer, intent(out) :: nstep
43	real(R8), intent(in) :: x2, h1, hmin
44	real(R8), intent(inout) :: x(nstepmax), y1(nstepmax), y2(nstepmax)
45	integer :: ode_solve
46	interface
47	subroutine myfunc(t, y, f)
48	use global_m
49	real(R8), intent(in) :: t
50	real(R8), intent(in) :: y(2)
51	real(R8), intent(out) :: f(2)
52	end subroutine myfunc
53	end interface
54	end function ode_solve
55	end interface
56
57
58	!---Derived Data Types---!
59
60	type integrator_t
61	private
62	integer :: nstepmax
63	integer(POINTER_SIZE) :: evl, stp, ctrl
64	real(R8) :: hmin, tol
65	end type integrator_t
66
67
68	!---Global Variables---!
69
70	integer, parameter :: SCHRODINGER = EQ1, &
71	DIRAC = EQ2, &
72	SCALAR_REL = EQ3
73
74	integer, parameter :: M_RK2 = 1, &
75	M_RK4 = 2, &
76	M_RKF4 = 3, &
77	M_RKCK4 = 4, &
78	M_RKPD8 = 5
79
80
81	!Wavefunctions values at infinity
82	real(R8), parameter :: SCHRODINGER_FINF = 1.0E-20_r8
83	real(R8), parameter :: DIRAC_FINF = 1.0E-20_r8
84
85
86	!---Public/Private Statements---!
87
88	private
89	public :: integrator_t, &
90	wave_equations_integrator_init, &
91	outward_integration, &
92	inward_integration, &
93	practical_infinity, &
94	wave_equations_integrator_end, &
95	SCHRODINGER, &
96	DIRAC, &
97	SCALAR_REL, &
98	SCHRODINGER_FINF, &
99	DIRAC_FINF, &
100	M_RK2, &
101	M_RK4, &
102	M_RKF4, &
103	M_RKCK4, &
104	M_RKPD8
105
106	contains
107
108	subroutine wave_equations_integrator_init(integrator_sp, integrator_dp)
109	!-----------------------------------------------------------------------!
110	! Initializes integrator objects by reading some parameters from the !
111	! input file and by initializing GSL ODE solver objects. !
112	! Two objects are initialized: one for single precision integrations !
113	! and a second one for double precision integrations. !
114	! !
115	! integrator_sp - single-precision integrator object !
116	! integrator_dp - double-precision integrator object !
117	!-----------------------------------------------------------------------!
118	type(integrator_t), intent(out) :: integrator_sp, integrator_dp
119
120	integer :: stepping_func
121	real(R8) :: tol
122
123	call push_sub("wave_equations_integrator_init")
124
125	message(1) = ""
126	message(2) = "Initializing Wave-Equations Integrator"
127	call write_info(2)
128
129	!Which stepping function will we use
130	call oct_parse_int('SteppingFunction', 5, stepping_func)
131	select case (stepping_func)
132	case (M_RK2)
133	message(1) = " Stepping function: Embedded 2nd order Runge-Kutta"
134	message(2) = " with 3rd order error estimate"
135	call write_info(2,20)
136	case (M_RK4)
137	message(1) = " Stepping function: 4th order (classical) Runge-Kutta method"
138	call write_info(1,20)
139	case (M_RKF4)
140	message(1) = " Stepping function: Embedded 4th order Runge-Kutta-Fehlberg "
141	message(2) = " method with 5th order error estimate"
142	call write_info(2,20)
143	case (M_RKCK4)
144	message(1) = " Stepping function: Embedded 4th order Runge-Kutta Cash-Karp"
145	message(2) = " method with 5th order error estimate"
146	call write_info(2,20)
147	case (M_RKPD8)
148	message(1) = " Stepping function: Embedded 8th order Runge-Kutta Prince-Dormand"
149	message(2) = " method with 9th order error estimate"
150	call write_info(2,20)
151	! case (6)
152	! message(1) = " Stepping function: Implicit 4th order Runge-Kutta at Gaussian points"
153	! call write_info(1,20)
154	! case (7)
155	! message(1) = " Stepping function: M=2 implicit Gear method"
156	! call write_info(1,20)
157	case default
158	message(1) = "SteppingFunction must be between 1 and 5."
159	call write_fatal(1)
160	end select
161
162	!Read the ODE integrator tolerance
163	call oct_parse_float('ODEIntTolerance', 1.0E-6_r8, tol)
164	if (tol <= M_ZERO) then
165	message(1) = "ODEIntTolerance must be greater than zero."
166	call write_fatal(1)
167	end if
168	write(message(1),'(2X,"ODE Integrator tolerance: ",ES10.3E2)') tol
169	call write_info(1,20)
170
171	!Read the maximum number of steps the ODE integrator is allowed to take
172	call oct_parse_int('MaxSteps', 50000, integrator_sp%nstepmax)
173	if (integrator_sp%nstepmax <= M_ZERO) then
174	message(1) = "MaxSteps must be greater than zero."
175	call write_fatal(1)
176	end if
177	write(message(1),'(2X,"ODE Integrator maximum number of steps: ",I6)') integrator_sp%nstepmax
178	call write_info(1,20)
179
180	!What is the minimum step size allowed
181	call oct_parse_float('MinimumStep', 1.0E-15_r8, integrator_sp%hmin)
182	if (integrator_sp%hmin <= M_ZERO) then
183	message(1) = "MinimumStep must be greater than zero."
184	call write_fatal(1)
185	end if
186	write(message(1),'(2X,"ODE Integrator minimum step size: ",ES10.3E2)') integrator_sp%hmin
187	call write_info(1,20)
188
189	integrator_dp = integrator_sp
190	integrator_sp%tol = sqrt(tol)
191	integrator_dp%tol = tol**2
192
193	!Initialize GSL objects
194	call gsl_odeiv_step_alloc(stepping_func, 2, integrator_sp%stp)
195	call gsl_odeiv_evolve_alloc(2, integrator_sp%evl)
196	call gsl_odeiv_control_standart_new(integrator_sp%ctrl, 1.0E-22_r8, integrator_sp%tol, M_ONE, M_ONE)
197
198	call gsl_odeiv_step_alloc(stepping_func, 2, integrator_dp%stp)
199	call gsl_odeiv_evolve_alloc(2, integrator_dp%evl)
200	call gsl_odeiv_control_standart_new(integrator_dp%ctrl, 1.0E-22_r8, integrator_dp%tol, M_ONE, M_ONE)
201
202	call pop_sub()
203	end subroutine wave_equations_integrator_init
204
205	subroutine outward_integration(qn, e, wave_eq, wf_dim, potential, integrator, nstep, r_out, wf_out, wfp_out, rc)
206	!-----------------------------------------------------------------------!
207	! Integrates the radial wave-equation outward from a point near the !
208	! origin to the classical turning point. !
209	! !
210	! qn - set of quantum numbers !
211	! e - energy !
212	! wave_eq - wave-equation to integrate !
213	! wf_dim - dimension of the wavefunction spinor !
214	! potential - potential object !
215	! integrator - integrator object !
216	! nstep - number of steps taken by the solver !
217	! r_out - mesh used by the ODE solver !
218	! wf_out - wavefunction !
219	! wfp_out - wavefunction derivative !
220	! rc - if present the integration will stop at rc instead of !
221	! stopping at the classical turning point !
222	!-----------------------------------------------------------------------!
223	type(qn_t), intent(in) :: qn
224	real(R8), intent(in) :: e
225	integer, intent(in) :: wave_eq, wf_dim
226	type(potential_t), intent(in) :: potential
227	type(integrator_t), intent(inout) :: integrator
228	integer, intent(out) :: nstep
229	real(R8), pointer :: r_out(:), wf_out(:,:), wfp_out(:,:)
230	real(R8), intent(in), optional :: rc
231
232	integer :: ierr, i
233	real(R8) :: z, vp0, r0, ri, s, a0, a1, a2, a3, b0, b1, b2, w
234	real(R8), allocatable :: r(:), g(:), f(:)
235
236	call push_sub("outward_integration")
237
238	!Allocate work arrays
239	allocate(r(integrator%nstepmax), g(integrator%nstepmax), f(integrator%nstepmax))
240
241	!Get the initial and final points
242	r0 = smallest_safe_r(potential)
243	if (present(rc)) then
244	ri = rc
245	else
246	ri = classical_turning_point(potential, e, qn)
247	end if
248
249	! For a potential v(r) = vp(r) - Z/r, the general solutions of the equations
250	! when r->0 are of the form:
251	!
252	! g(r) = r(s-1) (a0 + a1 r + a2 r2 + ...)
253	! f(r) = r(s-1) (b0 + b1 r + b2 r2 + ...)
254	!
255	! Schrodinger equation:
256	! s = l
257	! a0 = 0
258	! a2 = - Z/(l + 1) a1
259	! a3 = [-Z a2 + (vp - e) a1]/(2l + 3)
260	! b0 = l a1
261	! b1 = (l + 1) a2
262	! b2 = (l + 2) a3
263	!
264	! Scalar-relativistic equation:
265	! s = sqrt(l(l + 1) + 1 - (Z/c)**2)
266	! a1 = [(e + 2c2 - vp) b0 - Z/c (2e + 2c2 - vp) a0]/[(2s + 1)c]
267	! a2 = [(e + 2c2 - vp)(vp - e)/c a0 + (e + 2c2 - vp) b1 - Z/c (2e + 2c**2 - vp) a1]/[2(2s + 2)c]
268	! b0 = c/Z (s - 1) a0
269	! b1 = c/Z [s a1 - (e + 2c**2 - vp)/c b0]
270	! b2 = c/Z [(s + 1) a2 - (e + 2c**2 - vp)/c b1]
271	!
272	! Dirac equation:
273	! s = sqrt(k2 - (Z/c)2)
274	! a0 = Z/c b0/(s + k) v a0 = - c/Z (s - k) b0
275	! a1 = [(s + 1 - k)(e + 2c**2 - vp) b0 - Z/c (e - vp) a0]/(2s + 1)/c
276	! a2 = [(s + 2 - k)(e + 2c**2 - vp) b1 - Z/c (e - vp) a1]/(2s + 2)/(2c)
277	! b0 = -Z/c a0/(s - k) v b0 = c/Z (s + k) a0
278	! b1 = [-(e - vp)/c a0 - Z/c a1]/(s + 1 - k)
279	! b2 = [-(e - vp)/c a1 - Z/c a2]/(s + 2 - k)
280	!
281	! NOTE: for the scalar-relativistic equation and for Dirac's
282	! equation Z=0 is a special case.
283	r(1) = r0
284	z = potential_nuclear_charge(potential)
285	vp0 = v(potential, r0, qn) + z/r0
286	select case (wave_eq)
287	case (SCHRODINGER)
288	s = qn%l
289	a0 = M_ZERO
290	a1 = M_ONE
291	a2 = -z/(s + M_ONE)*a1
292	a3 = (-za2 + (vp0 - e)a1)/(M_TWO*s + M_THREE)
293	b0 = s*a1
294	b1 = (s + 1)*a2
295	b2 = (s + 2)*a3
296	case (DIRAC)
297	if (z /= M_ZERO) then
298	a0 = M_ONE
299	s = sqrt(qn%k2 - (z/M_C)2)
300	b0 = M_C/z(s + qn%k)a0
301	else
302	if (qn%k < 0) then
303	s = -qn%k
304	a0 = M_ONE
305	b0 = M_ZERO
306	else
307	s = qn%k
308	a0 = M_ZERO
309	b0 = M_ONE
310	end if
311	end if
312	w = e + M_TWOM_C*2 - vp0
313	a1 = (w(s + M_ONE - qn%k)b0 - z/M_C(e - vp0)a0)/(M_TWO*s + M_ONE)/M_C
314	b1 = -((e - vp0)a0 + za1)/(s + M_ONE - qn%k)/M_C
315	a2 = (w(s + M_TWO - qn%k)b1 - z/M_C(e - vp0)a1)/(M_TWOs + M_TWO)/(M_TWOM_C)
316	b2 = -((e - vp0)a1 + za2)/(s + M_TWO - qn%k)/M_C
317	a3 = M_ZERO
318	case (SCALAR_REL)
319	w = e + M_TWOM_C*2 - vp0
320	if (z /= M_ZERO) then
321	s = sqrt(qn%l(qn%l + M_ONE) + M_ONE - (z/M_C)*2)
322	a0 = M_ONE
323	b0 = M_C/z(s - M_ONE)a0
324	a1 = (wb0 - z/M_CM_TWO(e + M_C - vp0)a0)/(M_TWO*s + M_ONE)/M_C
325	b1 = M_C/z(sa1 - w*b0/M_C)
326	a2 = (w(vp0 - e)a0/M_C + wb1 - z/M_CM_TWO(e + M_C - vp0)a1)/(M_FOUR*s + M_FOUR)/M_C
327	b2 = M_C/z((s + M_ONE)a2 - w*b1/M_C)
328	a3 = M_ZERO
329	else
330	if (qn%l == 0) then
331	s = M_ONE
332	a0 = M_ONE
333	b1 = (vp0 - e)*a0/M_THREE
334	a2 = w/(M_TWOM_C)b1
335	b0 = M_ZERO ; a1 = M_ZERO ; b2 = M_ZERO ; a3 = M_ZERO
336	else
337	s = qn%l
338	b0 = M_ONE
339	a1 = w/(sM_C)b0
340	a0 = M_ZERO ; b1 = M_ZERO ; a2 = M_ZERO ; b2 = M_ZERO ; a3 = M_ZERO
341	end if
342	end if
343	end select
344	g(1) = r0*(s - M_ONE)(a0 + a1r0 + a2r0*2 + a3r0**3)
345	f(1) = r0*(s - M_ONE)(b0 + b1r0 + b2r0**2)
346
347	! MGGA term
348	if (wave_eq == SCHRODINGER) then
349	f(1) = f(1)(M_ONE + M_TWOvtau(potential, r0, qn))
350	end if
351
352	!Set the parameters needed to compute the derivatives of the functions
353	call set_derivs_params(qn, e, wave_eq, potential)
354
355	!Integrate the equation from r0 to ri
356	ierr = ode_solve(ri, integrator%stp, integrator%evl, integrator%ctrl, (r0+ri)/M_TWO, &
357	integrator%hmin, integrator%nstepmax, nstep, r, g, f, derivs)
358	if (ierr /= 0) then
359	if (in_debug_mode) then
360	call wave_equations_integrator_debug(integrator, r0, ri, nstep, r, g, f)
361	call wave_equations_derivs_debug()
362	end if
363	message(1) = "Error in subtoutine outward_integration. Error message:"
364	select case (ierr)
365	case (100)
366	message(2) = "stepsize underflow"
367	case (101)
368	message(2) = "integration diverged"
369	case default
370	call gsl_strerror(ierr, message(2))
371	end select
372	call write_fatal(2)
373	end if
374
375	!Copy the functions to new arrays
376	allocate(r_out(nstep), wf_out(nstep, wf_dim), wfp_out(nstep, wf_dim))
377	r_out = r(1:nstep)
378	select case (wave_eq)
379	case (SCHRODINGER)
380	wf_out(1:nstep,1) = g(1:nstep)
381	do i = 1, nstep
382	wfp_out(i,1) = f(i)/(M_ONE + M_TWO*vtau(potential, r_out(i), qn))
383	end do
384	case (SCALAR_REL)
385	wf_out(1:nstep,1) = g(1:nstep)
386	do i = 1, nstep
387	wfp_out(i,1) = (M_TWOM_C + (e - v(potential, r_out(i), qn))/M_C)f(i)
388	end do
389	case (DIRAC)
390	wf_out(1:nstep,1) = g(1:nstep)
391	wf_out(1:nstep,2) = f(1:nstep)
392	do i = 1, nstep
393	call derivs(r_out(i), wf_out(i,:), wfp_out(i,:))
394	end do
395	end select
396
397	!Deallocate arrays
398	deallocate(r, f, g)
399
400	!Unset the derivatives parameters
401	call unset_derivs_params()
402
403	call pop_sub()
404	end subroutine outward_integration
405
406	subroutine inward_integration(qn, e, wave_eq, wf_dim, potential, integrator, nstep, r_in, wf_in, wfp_in)
407	!-----------------------------------------------------------------------!
408	! Integrates the radial wave-equation inward from a point at infinity !
409	! to the classical turning point. !
410	! !
411	! qn - set of quantum numbers !
412	! e - energy !
413	! wave_eq - wave-equation to integrate !
414	! wf_dim - dimension of the wavefunction spinor !
415	! potential - potential object !
416	! integrator - integrator object !
417	! nstep - number of steps taken by the solver !
418	! r_in - mesh used by the ODE solver !
419	! wf_in - wavefunction !
420	! wfp_in - wavefunction derivative !
421	!-----------------------------------------------------------------------!
422	type(qn_t), intent(in) :: qn
423	real(R8), intent(in) :: e
424	integer, intent(in) :: wave_eq, wf_dim
425	type(potential_t), intent(in) :: potential
426	type(integrator_t), intent(inout) :: integrator
427	integer, intent(out) :: nstep
428	real(R8), pointer :: r_in(:), wf_in(:,:), wfp_in(:,:)
429
430	integer :: lpp, ierr, i
431	real(R8) :: vinf, kpp, rinf, ri, minf
432	real(R8), allocatable :: r(:), g(:), f(:)
433
434	call push_sub("inward_integration")
435
436	!Allocate work arrays
437	allocate(r(integrator%nstepmax), g(integrator%nstepmax), f(integrator%nstepmax))
438
439	!Get the initial and final points
440	ri = classical_turning_point(potential, e, qn)
441	rinf = practical_infinity(e, wave_eq, integrator%tol)
442
443	!Values of the wavefunction and its derivative at infinity
444	r(1) = rinf
445	select case (wave_eq)
446	case (SCHRODINGER)
447	g(1) = SCHRODINGER_FINF
448	f(1) = -sqrt(-M_TWOe)g(1)(M_ONE + M_TWOvtau(potential, rinf, qn))
449	case (DIRAC)
450	!Values of the wavefunctions at infinity
451	lpp = int(qn%j + M_HALF)*(qn%j - M_HALF)
452	if (lpp /= 0) lpp = lpp/qn%l
453	kpp = sqrt(-e)*sqrt(M_TWO + e/M_C2)
454	vinf = v(potential, rinf, qn)
455	g(1) = DIRAC_FINF
456
457	f(1) = (e - vinf)/M_C* &
458	(gsl_sf_bessel_knu_scaled(kpp*rinf, lpp + M_HALF) + (real(qn%k - lpp, R8) - M_ONE)/rinf)/ &
459	(kppgsl_sf_bessel_knu_scaled(kpprinf, lpp + 1.5_r8))*g(1)
460	case (SCALAR_REL)
461	g(1) = SCHRODINGER_FINF
462	vinf = v(potential, rinf, qn)
463	minf = M_ONE + (e - vinf)/M_TWO/M_C2
464	f(1) = -sqrt(-M_TWOminfe)/(M_TWOminfM_C)*g(1)
465	end select
466
467	!Set the parameters needed to compute the derivatives of the functions
468	call set_derivs_params(qn, e, wave_eq, potential)
469
470	!Integrate the equation from rinf to ri
471	ierr = ode_solve(ri, integrator%stp, integrator%evl, integrator%ctrl, (ri+rinf)/M_TWO, &
472	integrator%hmin, integrator%nstepmax, nstep, r, g, f, derivs)
473	if (ierr /= 0) then
474	if (in_debug_mode) then
475	call wave_equations_integrator_debug(integrator, rinf, ri, nstep, r, g, f)
476	call wave_equations_derivs_debug()
477	end if
478	message(1) = "Error in subtoutine inward_integration. Error message:"
479	select case (ierr)
480	case (100)
481	message(2) = "stepsize underflow"
482	case (101)
483	message(2) = "integration diverged"
484	case default
485	call gsl_strerror(ierr, message(2))
486	end select
487	call write_fatal(2)
488	end if
489
490	!Copy the functions to new arrays
491	allocate(r_in(nstep), wf_in(nstep, wf_dim), wfp_in(nstep, wf_dim))
492	r_in = r(1:nstep)
493	select case (wave_eq)
494	case (SCHRODINGER)
495	wf_in(1:nstep,1) = g(1:nstep)
496	do i = 1, nstep
497	wfp_in(i,1) = f(i)/(M_ONE + M_TWO*vtau(potential, r_in(i), qn))
498	end do
499	case (SCALAR_REL)
500	wf_in(1:nstep,1) = g(1:nstep)
501	do i = 1, nstep
502	wfp_in(i,1) = (M_TWOM_C + (e - v(potential, r_in(i), qn))/M_C)f(i)
503	end do
504	case (DIRAC)
505	wf_in(1:nstep,1) = g(1:nstep)
506	wf_in(1:nstep,2) = f(1:nstep)
507	do i = 1, nstep
508	call derivs(r_in(i), wf_in(i,:), wfp_in(i,:))
509	end do
510	end select
511
512	!Deallocate arrays
513	deallocate(r, g, f)
514
515	!Unset the derivatives parameters
516	call unset_derivs_params()
517
518	call pop_sub()
519	end subroutine inward_integration
520
521	function practical_infinity(e, wave_eq, tol)
522	!-----------------------------------------------------------------------!
523	! Returns the value of the practical infinity. The practical infinity !
524	! is such that wavefunction at the practical infinity is of the order !
525	! of magnitude of SCHRODINGER_FINF or DIRAC_FINF. !
526	! !
527	! e - energy !
528	! wave_eq - wave-equation to use !
529	! tol - tolerance !
530	!-----------------------------------------------------------------------!
531	real(R8), intent(in) :: e
532	integer, intent(in) :: wave_eq
533	real(R8), intent(in) :: tol
534	real(R8) :: practical_infinity
535
536	real(R8) :: k, x, f, xm, d
537
538	call push_sub("practical_infinity")
539
540	select case (wave_eq)
541	case (SCHRODINGER)
542	if (e == M_ZERO) then
543	practical_infinity = -log(SCHRODINGER_FINF)/sqrt(M_TWO*1.0e-12)
544	else
545	practical_infinity = -log(SCHRODINGER_FINF)/sqrt(-M_TWO*e)
546	end if
547	case (DIRAC)
548	k = sqrt(-e)*sqrt(M_TWO + e/M_C2)
549	x = M_ZERO
550	f = -M_HALFlog(x) - kx - log(DIRAC_FINF)
551	d = -log(DIRAC_FINF)/k - x
552	do
553	d = d*M_HALF
554	xm = x + d
555	if (f(-M_HALFlog(xm) - k*xm - log(DIRAC_FINF)) > M_ZERO) x = xm
556	if (-M_HALFlog(xm) - kxm - log(DIRAC_FINF) == M_ZERO .or. d <= tol) exit
557	end do
558	practical_infinity = xm
559	case (SCALAR_REL)
560	if (e == M_ZERO) then
561	practical_infinity = -log(SCHRODINGER_FINF)/sqrt(M_TWO(M_ONE - 1.0e-12/M_TWO/M_C2)1.0e-12)
562	else
563	practical_infinity = -log(SCHRODINGER_FINF)/sqrt(-M_TWO(M_ONE + e/M_TWO/M_C2)e)
564	end if
565
566	end select
567
568	call pop_sub()
569	end function practical_infinity
570
571	subroutine wave_equations_integrator_debug(integrator, ri, rf, nstep, r, g, f)
572	!-----------------------------------------------------------------------!
573	! Prints debug information to the "debug_info/integrator" file. !
574	! !
575	! integrator - integrator object !
576	! ri - initial integration point !
577	! rf - final integration point !
578	! nstep - number of steps taken by the solver !
579	! r - mesh used by the ODE solver !
580	! g - wavefunction !
581	! f - wavefunction !
582	!-----------------------------------------------------------------------!
583	type(integrator_t), intent(in) :: integrator
584	real(R8), intent(in) :: ri, rf
585	integer, intent(in) :: nstep
586	real(R8), intent(in) :: r(integrator%nstepmax), g(integrator%nstepmax), f(integrator%nstepmax)
587
588	integer :: unit, i
589
590	call push_sub("wave_equations_integrator_debug")
591
592	call io_open(unit, file='debug_info/integrator')
593	write(unit,'("ODE Integrator tolerance: ",ES10.3E2)') integrator%tol
594	write(unit,'("ODE Integrator maximum number of steps: ",I6)') integrator%nstepmax
595	write(unit,'("ODE Integrator minimum step size: ",ES10.3E2)') integrator%hmin
596	close(unit)
597
598	call io_open(unit, file='debug_info/wavefunctions')
599	write(unit,'("# Integration Starting Point: ",ES10.3E2)') ri
600	write(unit,'("# Integration Ending Point: ",ES10.3E2)') rf
601	do i = 1, nstep
602	write(unit,'(ES10.3E2,1X,ES10.3E2,1X,ES10.3E2)') r(i), g(i), f(i)
603	end do
604	close(unit)
605
606	call pop_sub()
607	end subroutine wave_equations_integrator_debug
608
609	subroutine wave_equations_integrator_end(integrator_sp, integrator_dp)
610	!-----------------------------------------------------------------------!
611	! Frees all the memory associated with the GSL objects used by the !
612	! ODE solver. !
613	!-----------------------------------------------------------------------!
614	type(integrator_t), intent(inout) :: integrator_sp, integrator_dp
615
616	call push_sub("wave_equations_integrator_end")
617
618	!Free GSL objects for simple precision integrations
619	call gsl_odeiv_step_free(integrator_sp%stp)
620	call gsl_odeiv_evolve_free(integrator_sp%evl)
621	call gsl_odeiv_control_free(integrator_sp%ctrl)
622
623	!Free GSL objects for double precision integrations
624	call gsl_odeiv_step_free(integrator_dp%stp)
625	call gsl_odeiv_evolve_free(integrator_dp%evl)
626	call gsl_odeiv_control_free(integrator_dp%ctrl)
627
628	call pop_sub()
629	end subroutine wave_equations_integrator_end
630
631	end module wave_equations_integrator_m

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format