Context Navigation

wave_equations_integrator.F90 @ 585

Revision 585, 23.7 KB checked in by micael, 2 years ago (diff)
More work to transform the wavefunctions into spinors. Derivatives of Dirac's wavefunctions are now computed using the wave-equation, instead of using numerical methods.
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	!Set the parameters needed to compute the derivatives of the functions
348	call set_derivs_params(qn, e, wave_eq, potential)
349
350	!Integrate the equation from r0 to ri
351	ierr = ode_solve(ri, integrator%stp, integrator%evl, integrator%ctrl, (r0+ri)/M_TWO, &
352	integrator%hmin, integrator%nstepmax, nstep, r, g, f, derivs)
353	if (ierr /= 0) then
354	if (in_debug_mode) then
355	call wave_equations_integrator_debug(integrator, r0, ri, nstep, r, g, f)
356	call wave_equations_derivs_debug()
357	end if
358	message(1) = "Error in subtoutine outward_integration. Error message:"
359	select case (ierr)
360	case (100)
361	message(2) = "stepsize underflow"
362	case (101)
363	message(2) = "integration diverged"
364	case default
365	call gsl_strerror(ierr, message(2))
366	end select
367	call write_fatal(2)
368	end if
369
370	!Copy the functions to new arrays
371	allocate(r_out(nstep), wf_out(nstep, wf_dim), wfp_out(nstep, wf_dim))
372	r_out = r(1:nstep)
373	select case (wave_eq)
374	case (SCHRODINGER)
375	wf_out(1:nstep,1) = g(1:nstep)
376	wfp_out(1:nstep,1) = f(1:nstep)
377	case (SCALAR_REL)
378	wf_out(1:nstep,1) = g(1:nstep)
379	do i = 1, nstep
380	wfp_out(i,1) = (M_TWOM_C + (e - v(potential, r_out(i), qn))/M_C)f(i)
381	end do
382	case (DIRAC)
383	wf_out(1:nstep,1) = g(1:nstep)
384	wf_out(1:nstep,2) = f(1:nstep)
385	do i = 1, nstep
386	call derivs(r_out(i), wf_out(i,:), wfp_out(i,:))
387	end do
388	end select
389
390	!Deallocate arrays
391	deallocate(r, f, g)
392
393	!Unset the derivatives parameters
394	call unset_derivs_params()
395
396	call pop_sub()
397	end subroutine outward_integration
398
399	subroutine inward_integration(qn, e, wave_eq, wf_dim, potential, integrator, nstep, r_in, wf_in, wfp_in)
400	!-----------------------------------------------------------------------!
401	! Integrates the radial wave-equation inward from a point at infinity !
402	! to the classical turning point. !
403	! !
404	! qn - set of quantum numbers !
405	! e - energy !
406	! wave_eq - wave-equation to integrate !
407	! wf_dim - dimension of the wavefunction spinor !
408	! potential - potential object !
409	! integrator - integrator object !
410	! nstep - number of steps taken by the solver !
411	! r_in - mesh used by the ODE solver !
412	! wf_in - wavefunction !
413	! wfp_in - wavefunction derivative !
414	!-----------------------------------------------------------------------!
415	type(qn_t), intent(in) :: qn
416	real(R8), intent(in) :: e
417	integer, intent(in) :: wave_eq, wf_dim
418	type(potential_t), intent(in) :: potential
419	type(integrator_t), intent(inout) :: integrator
420	integer, intent(out) :: nstep
421	real(R8), pointer :: r_in(:), wf_in(:,:), wfp_in(:,:)
422
423	integer :: lpp, ierr, i
424	real(R8) :: vinf, kpp, rinf, ri, minf
425	real(R8), allocatable :: r(:), g(:), f(:)
426
427	call push_sub("inward_integration")
428
429	!Allocate work arrays
430	allocate(r(integrator%nstepmax), g(integrator%nstepmax), f(integrator%nstepmax))
431
432	!Get the initial and final points
433	ri = classical_turning_point(potential, e, qn)
434	rinf = practical_infinity(e, wave_eq, integrator%tol)
435
436	!Values of the wavefunction and its derivative at infinity
437	r(1) = rinf
438	select case (wave_eq)
439	case (SCHRODINGER)
440	g(1) = SCHRODINGER_FINF
441	f(1) = -sqrt(-M_TWOe)g(1)
442	case (DIRAC)
443	!Values of the wavefunctions at infinity
444	lpp = int(qn%j + M_HALF)*(qn%j - M_HALF)
445	if (lpp /= 0) lpp = lpp/qn%l
446	kpp = sqrt(-e)*sqrt(M_TWO + e/M_C2)
447	vinf = v(potential, rinf, qn)
448	g(1) = DIRAC_FINF
449
450	f(1) = (e - vinf)/M_C* &
451	(gsl_sf_bessel_knu_scaled(kpp*rinf, lpp + M_HALF) + (real(qn%k - lpp, R8) - M_ONE)/rinf)/ &
452	(kppgsl_sf_bessel_knu_scaled(kpprinf, lpp + 1.5_r8))*g(1)
453	case (SCALAR_REL)
454	g(1) = SCHRODINGER_FINF
455	vinf = v(potential, rinf, qn)
456	minf = M_ONE + (e - vinf)/M_TWO/M_C2
457	f(1) = -sqrt(-M_TWOminfe)/(M_TWOminfM_C)*g(1)
458	end select
459
460	!Set the parameters needed to compute the derivatives of the functions
461	call set_derivs_params(qn, e, wave_eq, potential)
462
463	!Integrate the equation from rinf to ri
464	ierr = ode_solve(ri, integrator%stp, integrator%evl, integrator%ctrl, (ri+rinf)/M_TWO, &
465	integrator%hmin, integrator%nstepmax, nstep, r, g, f, derivs)
466	if (ierr /= 0) then
467	if (in_debug_mode) then
468	call wave_equations_integrator_debug(integrator, rinf, ri, nstep, r, g, f)
469	call wave_equations_derivs_debug()
470	end if
471	message(1) = "Error in subtoutine inward_integration. Error message:"
472	select case (ierr)
473	case (100)
474	message(2) = "stepsize underflow"
475	case (101)
476	message(2) = "integration diverged"
477	case default
478	call gsl_strerror(ierr, message(2))
479	end select
480	call write_fatal(2)
481	end if
482
483	!Copy the functions to new arrays
484	allocate(r_in(nstep), wf_in(nstep, wf_dim), wfp_in(nstep, wf_dim))
485	r_in = r(1:nstep)
486	select case (wave_eq)
487	case (SCHRODINGER)
488	wf_in(1:nstep,1) = g(1:nstep)
489	wfp_in(1:nstep,1) = f(1:nstep)
490	case (SCALAR_REL)
491	wf_in(1:nstep,1) = g(1:nstep)
492	do i = 1, nstep
493	wfp_in(i,1) = (M_TWOM_C + (e - v(potential, r_in(i), qn))/M_C)f(i)
494	end do
495	case (DIRAC)
496	wf_in(1:nstep,1) = g(1:nstep)
497	wf_in(1:nstep,2) = f(1:nstep)
498	do i = 1, nstep
499	call derivs(r_in(i), wf_in(i,:), wfp_in(i,:))
500	end do
501	end select
502
503	!Deallocate arrays
504	deallocate(r, g, f)
505
506	!Unset the derivatives parameters
507	call unset_derivs_params()
508
509	call pop_sub()
510	end subroutine inward_integration
511
512	function practical_infinity(e, wave_eq, tol)
513	!-----------------------------------------------------------------------!
514	! Returns the value of the practical infinity. The practical infinity !
515	! is such that wavefunction at the practical infinity is of the order !
516	! of magnitude of SCHRODINGER_FINF or DIRAC_FINF. !
517	! !
518	! e - energy !
519	! wave_eq - wave-equation to use !
520	! tol - tolerance !
521	!-----------------------------------------------------------------------!
522	real(R8), intent(in) :: e
523	integer, intent(in) :: wave_eq
524	real(R8), intent(in) :: tol
525	real(R8) :: practical_infinity
526
527	real(R8) :: k, x, f, xm, d
528
529	call push_sub("practical_infinity")
530
531	select case (wave_eq)
532	case (SCHRODINGER)
533	if (e == M_ZERO) then
534	practical_infinity = -log(SCHRODINGER_FINF)/sqrt(M_TWO*1.0e-12)
535	else
536	practical_infinity = -log(SCHRODINGER_FINF)/sqrt(-M_TWO*e)
537	end if
538	case (DIRAC)
539	k = sqrt(-e)*sqrt(M_TWO + e/M_C2)
540	x = M_ZERO
541	f = -M_HALFlog(x) - kx - log(DIRAC_FINF)
542	d = -log(DIRAC_FINF)/k - x
543	do
544	d = d*M_HALF
545	xm = x + d
546	if (f(-M_HALFlog(xm) - k*xm - log(DIRAC_FINF)) > M_ZERO) x = xm
547	if (-M_HALFlog(xm) - kxm - log(DIRAC_FINF) == M_ZERO .or. d <= tol) exit
548	end do
549	practical_infinity = xm
550	case (SCALAR_REL)
551	if (e == M_ZERO) then
552	practical_infinity = -log(SCHRODINGER_FINF)/sqrt(M_TWO(M_ONE - 1.0e-12/M_TWO/M_C2)1.0e-12)
553	else
554	practical_infinity = -log(SCHRODINGER_FINF)/sqrt(-M_TWO(M_ONE + e/M_TWO/M_C2)e)
555	end if
556
557	end select
558
559	call pop_sub()
560	end function practical_infinity
561
562	subroutine wave_equations_integrator_debug(integrator, ri, rf, nstep, r, g, f)
563	!-----------------------------------------------------------------------!
564	! Prints debug information to the "debug_info/integrator" file. !
565	! !
566	! integrator - integrator object !
567	! ri - initial integration point !
568	! rf - final integration point !
569	! nstep - number of steps taken by the solver !
570	! r - mesh used by the ODE solver !
571	! g - wavefunction !
572	! f - wavefunction !
573	!-----------------------------------------------------------------------!
574	type(integrator_t), intent(in) :: integrator
575	real(R8), intent(in) :: ri, rf
576	integer, intent(in) :: nstep
577	real(R8), intent(in) :: r(integrator%nstepmax), g(integrator%nstepmax), f(integrator%nstepmax)
578
579	integer :: unit, i
580
581	call push_sub("wave_equations_integrator_debug")
582
583	call io_open(unit, file='debug_info/integrator')
584	write(unit,'("ODE Integrator tolerance: ",ES10.3E2)') integrator%tol
585	write(unit,'("ODE Integrator maximum number of steps: ",I6)') integrator%nstepmax
586	write(unit,'("ODE Integrator minimum step size: ",ES10.3E2)') integrator%hmin
587	close(unit)
588
589	call io_open(unit, file='debug_info/wavefunctions')
590	write(unit,'("# Integration Starting Point: ",ES10.3E2)') ri
591	write(unit,'("# Integration Ending Point: ",ES10.3E2)') rf
592	do i = 1, nstep
593	write(unit,'(ES10.3E2,1X,ES10.3E2,1X,ES10.3E2)') r(i), g(i), f(i)
594	end do
595	close(unit)
596
597	call pop_sub()
598	end subroutine wave_equations_integrator_debug
599
600	subroutine wave_equations_integrator_end(integrator_sp, integrator_dp)
601	!-----------------------------------------------------------------------!
602	! Frees all the memory associated with the GSL objects used by the !
603	! ODE solver. !
604	!-----------------------------------------------------------------------!
605	type(integrator_t), intent(inout) :: integrator_sp, integrator_dp
606
607	call push_sub("wave_equations_integrator_end")
608
609	!Free GSL objects for simple precision integrations
610	call gsl_odeiv_step_free(integrator_sp%stp)
611	call gsl_odeiv_evolve_free(integrator_sp%evl)
612	call gsl_odeiv_control_free(integrator_sp%ctrl)
613
614	!Free GSL objects for double precision integrations
615	call gsl_odeiv_step_free(integrator_dp%stp)
616	call gsl_odeiv_evolve_free(integrator_dp%evl)
617	call gsl_odeiv_control_free(integrator_dp%ctrl)
618
619	call pop_sub()
620	end subroutine wave_equations_integrator_end
621
622	end module wave_equations_integrator_m

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

Download in other formats:

Original Format