V3FIT
integration_path.f
Go to the documentation of this file.
1 !*******************************************************************************
4 !
5 ! Note separating the Doxygen comment block here so detailed decription is
6 ! found in the Module not the file.
7 !
10 !*******************************************************************************
11 
12  MODULE integration_path
13  USE stel_kinds
14  USE stel_constants, ONLY : pi
15  USE mpi_inc
16 
17  IMPLICIT NONE
18 
19 !*******************************************************************************
20 ! integration_path module parameters
21 !*******************************************************************************
23  INTEGER, PARAMETER :: path_none_type = -1
25  INTEGER, PARAMETER :: path_add_type = 0
27  INTEGER, PARAMETER :: path_gleg_type = 1
29  INTEGER, PARAMETER :: path_hp_glep_type = 2
30 
32  REAL (rprec), PARAMETER :: path_default_dx = 0.0025
33 
34 !*******************************************************************************
35 ! DERIVED-TYPE DECLARATIONS
36 ! 1) vertex
37 !
38 !*******************************************************************************
39 !-------------------------------------------------------------------------------
41 !-------------------------------------------------------------------------------
42  TYPE path_int_class
44  INTEGER :: method = path_none_type
45 
47  REAL (rprec) :: dx
49  REAL (rprec) :: length
50 
52  REAL (rprec), DIMENSION(:), POINTER :: weights
54  REAL (rprec), DIMENSION(:), POINTER :: absc
55  END TYPE
56 
57 !-------------------------------------------------------------------------------
60 !-------------------------------------------------------------------------------
61  TYPE vertex
63  REAL (rprec), DIMENSION(3) :: position
65  TYPE (vertex), POINTER :: next => null()
66  END TYPE
67 
68 !*******************************************************************************
69 ! INTERFACE BLOCKS
70 !*******************************************************************************
71 !-------------------------------------------------------------------------------
74 !-------------------------------------------------------------------------------
75  INTERFACE path_construct
76  MODULE PROCEDURE path_construct_int, &
77  & path_construct_vertex
78  END INTERFACE
79 
80 !-------------------------------------------------------------------------------
83 !-------------------------------------------------------------------------------
84  INTERFACE path_destruct
85  MODULE PROCEDURE path_destruct_int, &
86  & path_destruct_vertex
87  END INTERFACE
88 
89 !-------------------------------------------------------------------------------
91 !-------------------------------------------------------------------------------
92  INTERFACE check
93  MODULE PROCEDURE check_log, check_real, check_int
94  END INTERFACE
95 
96  PRIVATE :: check, check_real, check_log, check_int, &
97  & integrate, integrate_add, integrate_gleg, &
98  & integrate_hp_gleg, search, test_function
99 
100  CONTAINS
101 
102 !*******************************************************************************
103 ! CONSTRUCTION SUBROUTINES
104 !*******************************************************************************
105 !-------------------------------------------------------------------------------
114 !-------------------------------------------------------------------------------
115  FUNCTION path_construct_int(method, npoints, length)
116 
117  IMPLICIT NONE
118 
119 ! Declare Arguments
120  TYPE (path_int_class), POINTER :: path_construct_int
121  CHARACTER (len=*), INTENT(in) :: method
122  INTEGER, INTENT(in) :: npoints
123  REAL (rprec), INTENT(in) :: length
124 
125 ! Start of executable code
126  ALLOCATE(path_construct_int)
127 
128  path_construct_int%dx = path_default_dx
129 
130  SELECT CASE (method)
131 
132  CASE ('add')
133  path_construct_int%method = path_add_type
134  path_construct_int%weights => null()
135  path_construct_int%absc => null()
136  RETURN
137 
138  CASE ('gleg')
139  path_construct_int%method = path_gleg_type
140 
141  CASE ('hp_gleg')
142  path_construct_int%method = path_hp_glep_type
143  path_construct_int%length = length
144 
145  CASE DEFAULT
146  path_construct_int%method = path_none_type
147 
148  END SELECT
149 
150  ALLOCATE(path_construct_int%weights(npoints))
151  ALLOCATE(path_construct_int%absc(npoints))
152  CALL path_get_gaussqad_weights(0.0_rprec, 1.0_rprec, &
153  & path_construct_int%absc, &
154  & path_construct_int%weights)
155 
156  END FUNCTION
157 
158 !-------------------------------------------------------------------------------
165 !-------------------------------------------------------------------------------
166  FUNCTION path_construct_vertex(position)
167 
168  IMPLICIT NONE
169 
170 ! Declare Arguments
171  TYPE (vertex), POINTER :: path_construct_vertex
172  REAL (rprec), DIMENSION(3), INTENT(in) :: position
173 
174 ! Start of executable code
175  ALLOCATE(path_construct_vertex)
176 
177  path_construct_vertex%position = position
178 
179  END FUNCTION
180 
181 !*******************************************************************************
182 ! DESTRUCTION SUBROUTINES
183 !*******************************************************************************
184 !-------------------------------------------------------------------------------
190 !-------------------------------------------------------------------------------
191  SUBROUTINE path_destruct_int(this)
192 
193  IMPLICIT NONE
194 
195 ! Declare Arguments
196  TYPE (path_int_class), POINTER :: this
197 
198 ! Start of executable code
199  IF (ASSOCIATED(this%weights)) THEN
200  DEALLOCATE(this%weights)
201  this%weights => null()
202  END IF
203 
204  IF (ASSOCIATED(this%absc)) THEN
205  DEALLOCATE(this%absc)
206  this%absc => null()
207  END IF
208 
209  DEALLOCATE(this)
210 
211  END SUBROUTINE
212 
213 !-------------------------------------------------------------------------------
220 !-------------------------------------------------------------------------------
221  RECURSIVE SUBROUTINE path_destruct_vertex(this)
222 
223  IMPLICIT NONE
224 
225 ! Declare Arguments
226  TYPE (vertex), POINTER :: this
227 
228 ! Start of executable code
229  IF (ASSOCIATED(this%next)) THEN
230  CALL path_destruct(this%next)
231  this%next => null()
232  END IF
233 
234  DEALLOCATE(this)
235  this => null()
236 
237  END SUBROUTINE
238 
239 !*******************************************************************************
240 ! UTILITY SUBROUTINES
241 !*******************************************************************************
242 !-------------------------------------------------------------------------------
251 !-------------------------------------------------------------------------------
252  RECURSIVE SUBROUTINE path_append_vertex(this, position)
253 
254  IMPLICIT NONE
255 
256 ! Declare Arguments
257  TYPE (vertex), POINTER :: this
258  REAL (rprec), DIMENSION(3), INTENT(in) :: position
259 
260 ! Start of executable code
261  IF (ASSOCIATED(this)) THEN
262  IF (ASSOCIATED(this%next)) THEN
263  CALL path_append_vertex(this%next, position)
264  ELSE
265  this%next => path_construct(position)
266  END IF
267  ELSE
268  this => path_construct(position)
269  END IF
270 
271  END SUBROUTINE
272 
273 !-------------------------------------------------------------------------------
289 !-------------------------------------------------------------------------------
290  RECURSIVE FUNCTION path_integrate(this, path, &
291  & integration_function, context) &
292  & result(total)
293 
294  IMPLICIT NONE
295 
296 ! Declare Arguments
297  REAL (rprec) :: total
298  TYPE (path_int_class), INTENT(in) :: this
299  TYPE (vertex), INTENT(in) :: path
300  CHARACTER (len=1), INTENT(in) :: context(:)
301  INTERFACE
302  FUNCTION integration_function(context, xcart, dxcart, &
303  & length, dx)
304  USE stel_kinds
305  REAL (rprec) :: integration_function
306  CHARACTER (len=1), INTENT(in) :: context(:)
307  REAL (rprec), DIMENSION(3), INTENT(in) :: xcart
308  REAL (rprec), DIMENSION(3), INTENT(in) :: dxcart
309  REAL (rprec), INTENT(in) :: length
310  REAL (rprec), INTENT(in) :: dx
311  END FUNCTION
312  END INTERFACE
313 
314 ! Start of executable code
315  If (ASSOCIATED(path%next)) THEN
316  total = path_integrate(this, path%next, integration_function, &
317  & context)
318  total = total + integrate(this, context, path, path%next, &
319  & integration_function)
320  ELSE
321  total = 0.0
322  END IF
323 
324  END FUNCTION
325 
326 !-------------------------------------------------------------------------------
341 !-------------------------------------------------------------------------------
342  RECURSIVE FUNCTION path_search(path, search_function, context, &
343  & found) RESULT(xcart)
344 
345  IMPLICIT NONE
346 
347  REAL (rprec), DIMENSION(3) :: xcart
348  TYPE (vertex), INTENT(in) :: path
349  CHARACTER (len=1), INTENT(in) :: context(:)
350  LOGICAL, INTENT(out) :: found
351  INTERFACE
352  FUNCTION search_function(context, xcart1, xcart2)
353  USE stel_kinds
354  LOGICAL :: search_function
355  CHARACTER (len=1), INTENT(in) :: context(:)
356  REAL (rprec), DIMENSION(3), INTENT(in) :: xcart1
357  REAL (rprec), DIMENSION(3), INTENT(in) :: xcart2
358  END FUNCTION
359  END INTERFACE
360 
361 ! Start of executable code
362  found = .false.
363 
364  IF (ASSOCIATED(path%next)) THEN
365  found = search(context, path, path%next, search_function, &
366  & xcart)
367  IF (.not.found) THEN
368  xcart = path_search(path%next, search_function, context, &
369  & found)
370  END IF
371  END IF
372 
373  END FUNCTION
374 
375 !*******************************************************************************
376 ! PRIVATE
377 !*******************************************************************************
378 !-------------------------------------------------------------------------------
393 !-------------------------------------------------------------------------------
394  FUNCTION integrate(this, context, vertex1, vertex2, &
395  & integration_function)
396 
397  IMPLICIT NONE
398 
399 ! Declare Arguments
400  REAL (rprec) :: integrate
401  TYPE (path_int_class), INTENT(in) :: this
402  CHARACTER(len=1), INTENT(in) :: context(:)
403  TYPE (vertex), INTENT(in) :: vertex1
404  TYPE (vertex), INTENT(in) :: vertex2
405  INTERFACE
406  FUNCTION integration_function(context, xcart, dxcart, &
407  & length, dx)
408  USE stel_kinds
409  REAL (rprec) :: integration_function
410  CHARACTER (len=1), INTENT(in) :: context(:)
411  REAL (rprec), DIMENSION(3), INTENT(in) :: xcart
412  REAL (rprec), DIMENSION(3), INTENT(in) :: dxcart
413  REAL (rprec), INTENT(in) :: length
414  REAL (rprec), INTENT(in) :: dx
415  END FUNCTION
416  END INTERFACE
417 
418 ! Start of executable code
419  SELECT CASE (this%method)
420 
421  CASE (path_add_type)
422  integrate = integrate_add(this, context, vertex1, vertex2, &
423  & integration_function)
424 
425  CASE (path_gleg_type)
426  integrate = integrate_gleg(this, context, vertex1, &
427  & vertex2, integration_function)
428 
429  CASE (path_hp_glep_type)
430  integrate = &
431  & integrate_hp_gleg(this, context, vertex1, vertex2, &
432  & integration_function)
433 
434  CASE DEFAULT
435  integrate = 0.0
436 
437  END SELECT
438 
439  END FUNCTION
440 
441 !-------------------------------------------------------------------------------
457 !-------------------------------------------------------------------------------
458  FUNCTION integrate_add(this, context, vertex1, vertex2, &
459  & integration_function)
460 
461  IMPLICIT NONE
462 
463 ! Declare Arguments
464  REAL (rprec) :: integrate_add
465  TYPE (path_int_class), INTENT(in) :: this
466  CHARACTER(len=1), INTENT(in) :: context(:)
467  TYPE (vertex), INTENT(in) :: vertex1
468  TYPE (vertex), INTENT(in) :: vertex2
469  interface
470  FUNCTION integration_function(context, xcart, dxcart, &
471  & length, dx)
472  USE stel_kinds
473  REAL (rprec) :: integration_function
474  CHARACTER (len=1), INTENT(in) :: context(:)
475  REAL (rprec), DIMENSION(3), INTENT(in) :: xcart, dxcart
476  REAL (rprec), INTENT(in) :: length, dx
477  END FUNCTION
478  END INTERFACE
479 
480 ! local variables
481  REAL (rprec), DIMENSION(3) :: xcart
482  REAL (rprec), DIMENSION(3) :: dxcart
483  REAL (rprec) :: length
484  REAL (rprec) :: dx
485  INTEGER :: i, nsteps
486 
487 ! Start of executable code
488 ! Determine the number of integration steps to take by dividing the path length
489 ! by the step length and rounding to the nearest integer.
490  dxcart = vertex2%position - vertex1%position
491  length = sqrt(dot_product(dxcart, dxcart))
492  nsteps = int(length/this%dx)
493 ! Choose the actual step size.
494  dxcart = dxcart/nsteps
495  dx = sqrt(dot_product(dxcart, dxcart))
496 
497 
498 ! Integrate the length in addition.
499  integrate_add = 0.0
500  length = 0.0
501  xcart = vertex1%position
502  DO i = 1, nsteps
503  xcart = xcart + dxcart
504  length = length + dx
505  integrate_add = integrate_add &
506  & + integration_function(context, xcart, dxcart, &
507  & length, dx)
508  END DO
509 
510  END FUNCTION
511 
512 !-------------------------------------------------------------------------------
528 !-------------------------------------------------------------------------------
529  FUNCTION integrate_gleg(this, context, vertex1, vertex2, &
530  & integration_function)
531 
532  IMPLICIT NONE
533 
534 ! Declare Arguments
535  REAL (rprec) :: integrate_gleg
536  TYPE (path_int_class), INTENT(in) :: this
537  CHARACTER(len=1), INTENT(in) :: context(:)
538  TYPE (vertex), INTENT(in) :: vertex1
539  TYPE (vertex), INTENT(in) :: vertex2
540  INTERFACE
541  FUNCTION integration_function(context, xcart, dxcart, &
542  & length, dx)
543  USE stel_kinds
544  REAL (rprec) :: integration_function
545  CHARACTER (len=1), INTENT(in) :: context(:)
546  REAL (rprec), DIMENSION(3), INTENT(in) :: xcart
547  REAL (rprec), DIMENSION(3), INTENT(in) :: dxcart
548  REAL (rprec), INTENT(in) :: length
549  REAL (rprec), INTENT(in) :: dx
550  END FUNCTION
551  END INTERFACE
552 
553 ! local variables
554  REAL (rprec), DIMENSION(3) :: xcart
555  REAL (rprec), DIMENSION(3) :: xcart1
556  REAL (rprec), DIMENSION(3) :: dxcart
557  REAL (rprec) :: length
558  REAL (rprec) :: temp_length
559  REAL (rprec) :: dx
560  INTEGER :: i
561 
562 ! Start of executable code
563 
564 ! Determine the number of integration steps to take by dividing the path length
565 ! by the step length and rounding to the nearest integer.
566  dxcart = vertex2%position - vertex1%position
567  length = sqrt(dot_product(dxcart, dxcart))
568  dxcart(:) = dxcart(:)/length
569 
570 ! Integrate the length using Gauss-Legendre quadrature
571  integrate_gleg = 0.0
572  xcart = vertex1%position
573  DO i = 1, SIZE(this%weights)
574  xcart = vertex1%position + this%absc(i)*dxcart*length
575  temp_length = this%absc(i)*length
576 
577  integrate_gleg = integrate_gleg + &
578  & this%weights(i)*integration_function(context, xcart, &
579  & dxcart, temp_length, &
580  & 1.0_rprec)
581  END DO
582  integrate_gleg = integrate_gleg*length
583 
584  END FUNCTION
585 
586 !-------------------------------------------------------------------------------
604 !-------------------------------------------------------------------------------
605  FUNCTION integrate_hp_gleg(this, context, vertex1, vertex2, &
606  & integration_function)
607 
608  IMPLICIT NONE
609 
610 ! Declare Arguments
611  REAL (rprec) :: integrate_hp_gleg
612  TYPE (path_int_class), INTENT(in) :: this
613  CHARACTER(len=1), INTENT(in) :: context(:)
614  TYPE (vertex), INTENT(in) :: vertex1
615  TYPE (vertex), INTENT(in) :: vertex2
616  interface
617  function integration_function(context, xcart, dxcart, &
618  & length, dx)
619  USE stel_kinds
620  REAL (rprec) :: integration_function
621  CHARACTER (len=1), INTENT(in) :: context(:)
622  REAL (rprec), DIMENSION(3), INTENT(in) :: xcart
623  REAL (rprec), DIMENSION(3), INTENT(in) :: dxcart
624  REAL (rprec), INTENT(in) :: length
625  REAL (rprec), INTENT(in) :: dx
626  END FUNCTION
627  END INTERFACE
628 
629 ! local variables
630  REAL (rprec), DIMENSION(3) :: xcart
631  REAL (rprec), DIMENSION(3) :: xcart1
632  REAL (rprec), DIMENSION(3) :: dxcart
633  REAL (rprec) :: length
634  REAL (rprec) :: lengthp
635  REAL (rprec) :: temp_length
636  INTEGER :: i
637  INTEGER :: j
638  INTEGER :: nsteps
639  REAL (rprec) :: int_length
640  REAL (rprec) :: integratej
641 
642 ! Start of executable code
643 
644 ! Determine the number of integration steps to take by dividing the path length
645 ! by the step length and rounding to the nearest integer.
646  dxcart = vertex2%position - vertex1%position
647  length = sqrt(dot_product(dxcart, dxcart))
648  dxcart(:) = dxcart(:)/length
649 
650  nsteps = int(length/this%length)
651  int_length = length/real(nsteps, rprec)
652 
653 ! Integrate the length using Gauss-Legendre quadrature
654  integrate_hp_gleg = 0.0
655  xcart = vertex1%position
656  lengthp = 0.0
657  DO j = 1, nsteps
658  integratej = 0.0
659  DO i = 1, SIZE(this%weights)
660  xcart1 = xcart + this%absc(i)*dxcart*int_length
661  temp_length = lengthp + this%absc(i)*int_length
662  integratej = integratej + &
663  & this%weights(i)*integration_function(context, xcart1, &
664  & dxcart, temp_length, &
665  & 1.0_rprec)
666  END DO
667 
668  integrate_hp_gleg = integrate_hp_gleg + integratej*int_length
669  lengthp = lengthp + int_length
670  xcart = xcart + dxcart*int_length
671  END DO
672 
673  END FUNCTION
674 
675 !-------------------------------------------------------------------------------
689 !-------------------------------------------------------------------------------
690  FUNCTION search(context, vertex1, vertex2, search_function, xcart)
691 
692  IMPLICIT NONE
693 
694 ! Declare Arguments
695  LOGICAL :: search
696  CHARACTER(len=1), INTENT(in) :: context(:)
697  TYPE (vertex), INTENT(in) :: vertex1
698  TYPE (vertex), INTENT(in) :: vertex2
699  REAL (rprec), DIMENSION(3), INTENT(out) :: xcart
700  INTERFACE
701  FUNCTION search_function(context, xcart1, xcart2)
702  USE stel_kinds
703  LOGICAL :: search_function
704  CHARACTER (len=1), INTENT(in) :: context(:)
705  REAL (rprec), DIMENSION(3), INTENT(in) :: xcart1
706  REAL (rprec), DIMENSION(3), INTENT(in) :: xcart2
707  END FUNCTION
708  END INTERFACE
709 
710 ! local variables
711  REAL (rprec), DIMENSION(3) :: dxcart
712  REAL (rprec) :: dx
713  INTEGER :: i
714  INTEGER :: nsteps
715 
716 ! local parameters
717  REAL (rprec), PARAMETER :: dxcourse = 0.01
718  REAL (rprec), PARAMETER :: dxfine = 1.0e-20
719 
720 ! Start of executable code
721 ! Determine the number of integration steps to take by dividing the path length
722 ! by the step length and rounding to the nearest integer.
723  dxcart = vertex2%position - vertex1%position
724  nsteps = int(sqrt(dot_product(dxcart, dxcart))/dxcourse)
725 
726 ! Choose the actual step size.
727  dxcart = dxcart/nsteps
728 
729  search = .false.
730  xcart = vertex1%position
731 
732 ! Linearly search the line until an interval containing the point is detected.
733  DO i = 1, nsteps
734  search = search_function(context, xcart, xcart + dxcart)
735  IF (search) THEN
736 
737 ! Found an interval. Bisect the interval until the length is machine precision.
738  DO WHILE (sqrt(dot_product(dxcart, dxcart)) .gt. dxfine)
739  dxcart = dxcart/2.0
740  IF (.not.search_function(context, xcart, &
741  & xcart + dxcart)) THEN
742  xcart = xcart + dxcart
743  END IF
744  END DO
745 
746  xcart = xcart + dxcart/2.0
747  RETURN
748 
749  END IF
750 
751  xcart = xcart + dxcart
752  END DO
753 
754  END FUNCTION
755 
756 !-------------------------------------------------------------------------------
768 !-------------------------------------------------------------------------------
769  SUBROUTINE path_get_gaussqad_weights(a, b, abscissas, weights)
770 
771  IMPLICIT NONE
772 
773 ! Declare Arguments
774  REAL (rprec), INTENT(in) :: a
775  REAL (rprec), INTENT(in) :: b
776  REAL(rprec), DIMENSION(:), INTENT(out) :: abscissas
777  REAL(rprec), DIMENSION(:), INTENT(out) :: weights
778 
779 ! local variables
780  INTEGER :: i
781  INTEGER :: j
782  INTEGER :: m
783  REAL (rprec) :: p1
784  REAL (rprec) :: p2
785  REAL (rprec) :: p3
786  REAL (rprec) :: pp
787  REAL (rprec) :: xl
788  REAL (rprec) :: xm
789  REAL (rprec) :: current
790  REAL (rprec) :: last
791 
792 ! local parameters
793  REAL (rprec), PARAMETER :: eps = 1.0e-14
794 
795 ! Start of executable code
796  m = (SIZE(weights) + 1)/2
797 
798 ! Change basis from [-1,1] to [a.b]
799  xm = 0.5*(b + a)
800  xl = 0.5*(b - a)
801 
802  DO i = 1, m
803 ! current guess for i-th root of P_n(x)
804  current = cos(pi*(real(i, rprec) - 0.25_rprec)/ &
805  & (real(SIZE(weights), rprec) + 0.5_rprec))
806 
807 ! Newtons methods to find the roots of P_n(x)
808  DO
809 
810 ! Calculate P_n(x)
811  p1 = 1.0
812  p2 = 0.0
813  DO j = 1, SIZE(weights)
814  p3 = p2
815  p2 = p1
816  p1 = (real(2.0*j - 1, rprec)*current*p2 - &
817  & real(j - 1, rprec)*p3)/real(j, rprec)
818  END DO
819 
820 ! P'_n = n(xP_n - P_(n-1)/(x^2-1)
821  pp = real(SIZE(weights), rprec)*(current*p1 - p2) &
822  & / (current**2 - 1._rprec)
823 
824 ! Previous guess for i-th root
825  last = current
826  current = last - p1/pp
827  IF (abs(current - last) .le. eps) EXIT
828  END DO
829 
830  abscissas(i) = xm - xl*current
831  abscissas(SIZE(weights) + 1 - i) = xm + xl*current
832 
833 ! wi = 2/[(1 - x^2)P'n(x)^2]
834  weights(i) = 2.0*xl/((1.0 - current*current)*pp*pp)
835  weights(SIZE(weights) + 1 - i) = weights(i)
836  END DO
837 
838  END SUBROUTINE
839 
840 !*******************************************************************************
841 ! UNIT TESTS
842 !*******************************************************************************
843 !-------------------------------------------------------------------------------
849 !-------------------------------------------------------------------------------
850  FUNCTION path_test()
851 
852  IMPLICIT NONE
853 
854 ! Declare Arguments
855  LOGICAL :: path_test
856 
857 ! local variables
858  REAL (rprec) :: result
859  TYPE (vertex), POINTER :: test_path => null()
860  TYPE (path_int_class), POINTER :: int_params
861  CHARACTER(len=1), ALLOCATABLE :: context(:)
862  REAL(rprec), DIMENSION(3) :: absc
863  REAL(rprec), DIMENSION(3) :: wgts
864  REAL(rprec) :: work
865 
866 ! Start of executable code
867 ! Test to make sure the vertices begin in an unallocated state.
868  path_test = check(.false., ASSOCIATED(test_path), 1, 'ASSOCIATED')
869  IF (.not.path_test) RETURN
870 
871 ! Test to make sure first vertex is created.
872  CALL path_append_vertex(test_path, &
873  & (/ 1.0_rprec, 2.0_rprec, 3.0_rprec /))
874  path_test = check(.true., ASSOCIATED(test_path), 1, &
875  & 'path_append_vertex')
876  IF (.not.path_test) RETURN
877  path_test = check(1.0_rprec, test_path%position(1), 2, &
878  & 'path_append_vertex')
879  IF (.not.path_test) RETURN
880  path_test = check(2.0_rprec, test_path%position(2), 3, &
881  & 'path_append_vertex')
882  IF (.not.path_test) RETURN
883  path_test = check(3.0_rprec, test_path%position(3), 4, &
884  & 'path_append_vertex')
885  IF (.not.path_test) RETURN
886 
887 ! Test to make sure second vertex is appended.
888  CALL path_append_vertex(test_path, &
889  & (/ 4.0_rprec, 5.0_rprec, 6.0_rprec /))
890  path_test = check(.true., ASSOCIATED(test_path%next), 5, &
891  & 'path_append_vertex')
892  IF (.not.path_test) RETURN
893  path_test = check(1.0_rprec, test_path%position(1), 6, &
894  & 'path_append_vertex')
895  IF (.not.path_test) RETURN
896  path_test = check(2.0_rprec, test_path%position(2), 7, &
897  & 'path_append_vertex')
898  IF (.not.path_test) RETURN
899  path_test = check(3.0_rprec, test_path%position(3), 8, &
900  & 'path_append_vertex')
901  IF (.not.path_test) RETURN
902  path_test = check(4.0_rprec, test_path%next%position(1), 9, &
903  & 'path_append_vertex')
904  IF (.not.path_test) RETURN
905  path_test = check(5.0_rprec, test_path%next%position(2), 10, &
906  & 'path_append_vertex')
907  IF (.not.path_test) RETURN
908  path_test = check(6.0_rprec, test_path%next%position(3), 11, &
909  & 'path_append_vertex')
910  IF (.not.path_test) RETURN
911 
912 ! Test to make sure third vertex is appended.
913  CALL path_append_vertex(test_path, &
914  & (/ 7.0_rprec, 8.0_rprec, 9.0_rprec /))
915  path_test = check(.true., ASSOCIATED(test_path%next%next), &
916  & 12, 'path_append_vertex')
917  IF (.not.path_test) RETURN
918  path_test = check(7.0_rprec, test_path%next%next%position(1), 12, &
919  & 'path_append_vertex')
920  IF (.not.path_test) RETURN
921  path_test = check(8.0_rprec, test_path%next%next%position(2), 13, &
922  & 'path_append_vertex')
923  IF (.not.path_test) RETURN
924  path_test = check(9.0_rprec, test_path%next%next%position(3), 14, &
925  & 'path_append_vertex')
926  IF (.not.path_test) RETURN
927 
928 ! Test to make sure path object is destroyed.
929  CALL path_destruct(test_path)
930  path_test = check(.false., ASSOCIATED(test_path), 1, &
931  & 'path_destruct')
932  IF (.not.path_test) RETURN
933 
934 ! Test to make sure a path_in_class object was allocated.
935  int_params => path_construct('add', 0, 0.0_rprec)
936  path_test = check(.true., ASSOCIATED(int_params), 1, &
937  & 'path_construct_int')
938  IF (.not.path_test) RETURN
939 
940 ! Test path integration.
941  CALL path_append_vertex(test_path, &
942  & (/ 2.0_rprec, 0.0_rprec, 0.0_rprec /))
943  CALL path_append_vertex(test_path, &
944  & (/ 0.0_rprec, 0.0_rprec, 0.0_rprec /))
945  result = path_integrate(int_params, test_path, test_function, &
946  & context)
947  path_test = check(2.0_rprec, result, 1, 'path_integrate')
948  CALL path_destruct(int_params)
949 
950 ! Test gll integrator weights and abscissas
951  CALL path_get_gaussqad_weights(-1.0_rprec, 1.0_rprec, absc, wgts)
952  work = sqrt(0.6_rprec)
953  path_test = check(-1.0_rprec*work, absc(1), 1, 'get_gaussqad_a')
954  path_test = check( 0.0_rprec, absc(2), 2, 'get_gaussqad_a')
955  path_test = check( work, absc(3), 3, 'get_gaussqad_a')
956  path_test = check(5._rprec/9._rprec, wgts(1), 4, 'get_gaussqad_w')
957  path_test = check(8._rprec/9._rprec, wgts(2), 5, 'get_gaussqad_w')
958  path_test = check(5._rprec/9._rprec, wgts(3), 6, 'get_gaussqad_w')
959  IF (.not.path_test) RETURN
960 
961 ! Test gleg integrator.
962  int_params => path_construct('gleg', 100, 0.0_rprec)
963  result = path_integrate(int_params, test_path, test_function, &
964  & context)
965  path_test = check(2.0_rprec, result, 1, 'path_integrate_gleg')
966  IF (.not.path_test) RETURN
967  CALL path_destruct(int_params)
968 
969 ! Test hp_gleg integrator.
970  int_params => path_construct('hp_gleg', 100, 0.01_rprec)
971  result = path_integrate(int_params, test_path, test_function, &
972  & context)
973  path_test = check(2.0_rprec, result, 1, &
974  & 'path_integrate_hp_gleg')
975  IF (.not.path_test) RETURN
976  CALL path_destruct(test_path)
977  CALL path_destruct(int_params)
978 
979  END FUNCTION
980 
981 !*******************************************************************************
982 ! CHECK FUNCTIONS
983 !*******************************************************************************
984 !-------------------------------------------------------------------------------
995 !-------------------------------------------------------------------------------
996  FUNCTION check_log(expected, received, testNum, name)
997 
998  IMPLICIT NONE
999 
1000 ! Declare Arguments
1001  LOGICAL :: check_log
1002  LOGICAL, INTENT(in) :: expected
1003  LOGICAL, INTENT(in) :: received
1004  INTEGER, INTENT(in) :: testnum
1005  CHARACTER (LEN=*), INTENT(in) :: name
1006 
1007 ! Start of executable code
1008  check_log = expected .eqv. received
1009  IF (.not.check_log) THEN
1010  WRITE(*,*) "integration_path.f: ", name, " test", testnum, &
1011  & "failed."
1012  WRITE(*,*) "Expected", expected, "Received", received
1013  END IF
1014 
1015  END FUNCTION
1016 
1017 !-------------------------------------------------------------------------------
1028 !-------------------------------------------------------------------------------
1029  FUNCTION check_real(expected, received, testNum, name)
1030 
1031  IMPLICIT NONE
1032 
1033 ! Declare Arguments
1034  LOGICAL :: check_real
1035  REAL(rprec), INTENT(in) :: expected
1036  REAL(rprec), INTENT(in) :: received
1037  INTEGER, INTENT(in) :: testnum
1038  CHARACTER (LEN=*), INTENT(in) :: name
1039 
1040 ! local variables
1041  REAL(rprec) :: eps = 1.0e-8
1042 
1043 ! Start of executable code
1044  check_real = (abs(expected - received) .lt. eps)
1045  IF (.not.check_real) THEN
1046  WRITE(*,*) "integration_path.f: ", name, " test", testnum, &
1047  & "failed."
1048  WRITE(*,*) "Expected", expected, "Received", received
1049  END IF
1050 
1051  END FUNCTION
1052 
1053 !-------------------------------------------------------------------------------
1064 !-------------------------------------------------------------------------------
1065  FUNCTION check_int(expected, received, testNum, name)
1066 
1067  IMPLICIT NONE
1068 
1069 ! Declare Arguments
1070  LOGICAL :: check_int
1071  INTEGER, INTENT(in) :: expected
1072  INTEGER, INTENT(in) :: received
1073  INTEGER, INTENT(in) :: testnum
1074  CHARACTER (LEN=*), INTENT(in) :: name
1075 
1076 ! Start of executable code
1077  check_int = expected .eq. received
1078  IF (.not.check_int) THEN
1079  WRITE(*,*) "integration_path.f: ", name, " test", testnum, &
1080  & "failed."
1081  WRITE(*,*) "Expected", expected, "Received", received
1082  END IF
1083 
1084  END FUNCTION
1085 
1086 !-------------------------------------------------------------------------------
1098 !-------------------------------------------------------------------------------
1099  FUNCTION test_function(context, xcart, dxcart, length, dx)
1100 
1101  IMPLICIT NONE
1102 
1103 ! Declare Arguments
1104  REAL (rprec) :: test_function
1105  CHARACTER (len=1), INTENT(in) :: context(:)
1106  REAL (rprec), DIMENSION(3), INTENT(in) :: xcart
1107  REAL (rprec), DIMENSION(3), INTENT(in) :: dxcart
1108  REAL (rprec), INTENT(in) :: length
1109  REAL (rprec), INTENT(in) :: dx
1110 
1111 ! Start of executable code
1112  test_function = dx
1113 
1114  END FUNCTION
1115 
1116  END MODULE
integration_path::path_destruct_vertex
recursive subroutine path_destruct_vertex(this)
Deconstruct a vertex object.
Definition: integration_path.f:222
integration_path::path_construct_vertex
type(vertex) function, pointer path_construct_vertex(position)
Construct a single vertex.
Definition: integration_path.f:167
integration_path::path_get_gaussqad_weights
subroutine path_get_gaussqad_weights(a, b, abscissas, weights)
@breif Calculate the weights and abscissas for Gauss-Legendre integration.
Definition: integration_path.f:770
integration_path::path_add_type
integer, parameter path_add_type
Integrate the length in addition.
Definition: integration_path.f:25
integration_path::path_construct
Construction interface using either path_construct_int or path_construct_vertex.
Definition: integration_path.f:75
integration_path::integrate_hp_gleg
real(rprec) function, private integrate_hp_gleg(this, context, vertex1, vertex2, integration_function)
Integrate between to points using hp Guass-Legendre quadrature method.
Definition: integration_path.f:607
integration_path::integrate_add
real(rprec) function, private integrate_add(this, context, vertex1, vertex2, integration_function)
Line integrate between to points.
Definition: integration_path.f:460
integration_path::path_default_dx
real(rprec), parameter path_default_dx
Default step size of the integration.
Definition: integration_path.f:32
integration_path::check_int
logical function, private check_int(expected, received, testNum, name)
Check an integer value.
Definition: integration_path.f:1066
vertex
A vertex.
Definition: vertex.hpp:21
integration_path::path_gleg_type
integer, parameter path_gleg_type
Integrate the path using Gauss Legendre Quadrature.
Definition: integration_path.f:27
mpi_inc
Umbrella module avoid multiple inlcudes of the mpif.h header.
Definition: mpi_inc.f:11
integration_path::path_construct_int
type(path_int_class) function, pointer path_construct_int(method, npoints, length)
Construct a single path_int_class.
Definition: integration_path.f:116
integration_path
Module is part of the LIBSTELL. This modules contains code to define and integrate along an arbitray ...
Definition: integration_path.f:12
integration_path::integrate
real(rprec) function, private integrate(this, context, vertex1, vertex2, integration_function)
Line integrate between to points.
Definition: integration_path.f:396
integration_path::path_hp_glep_type
integer, parameter path_hp_glep_type
Integrate the path using a hp approch.
Definition: integration_path.f:29
integration_path::check_log
logical function, private check_log(expected, received, testNum, name)
Check a logical value.
Definition: integration_path.f:997
integration_path::path_integrate
recursive real(rprec) function path_integrate(this, path, integration_function, context)
Integrate along the path.
Definition: integration_path.f:293
integration_path::path_append_vertex
recursive subroutine path_append_vertex(this, position)
Append a vertex to a path.
Definition: integration_path.f:253
integration_path::path_search
recursive real(rprec) function, dimension(3) path_search(path, search_function, context, found)
Search along the path.
Definition: integration_path.f:344
integration_path::path_test
logical function path_test()
Path unit test function.
Definition: integration_path.f:851
integration_path::test_function
real(rprec) function, private test_function(context, xcart, dxcart, length, dx)
Call back function to test the integration.
Definition: integration_path.f:1100
integration_path::check_real
logical function, private check_real(expected, received, testNum, name)
Check a real value.
Definition: integration_path.f:1030
integration_path::integrate_gleg
real(rprec) function, private integrate_gleg(this, context, vertex1, vertex2, integration_function)
Integrate between to points using Guass-Legendre quadrature.
Definition: integration_path.f:531
integration_path::check
Interface for checking the results of the unit tests.
Definition: integration_path.f:92
integration_path::path_destruct
Destruct interface using either path_destruct_int or path_destruct_vertex.
Definition: integration_path.f:84
integration_path::search
logical function, private search(context, vertex1, vertex2, search_function, xcart)
Search line between to points.
Definition: integration_path.f:691
integration_path::path_destruct_int
subroutine path_destruct_int(this)
Deconstruct a path_int_class object.
Definition: integration_path.f:192
integration_path::path_none_type
integer, parameter path_none_type
No qaudrature type.
Definition: integration_path.f:23
integration_path::path_int_class
Class containing the parameters of the integration method to use.
Definition: integration_path.f:42
Generated on Thu Mar 5 2020 15:49:23 for V3FIT by   doxygen 1.8.17