-
Notifications
You must be signed in to change notification settings - Fork 0
/
demo_abscab.f90
140 lines (109 loc) · 4.19 KB
/
demo_abscab.f90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
module mod_demo_abscab
use abscab
implicit none
real(wp) :: radius, rCorr, omega
contains
subroutine vertexSupplierStd(i, pointData)
integer, intent(in) :: i
real(wp), intent(out), dimension(3) :: pointData
real(wp) :: phi
phi = omega * real(i, kind=wp)
pointData(1) = radius * cos(phi)
pointData(2) = radius * sin(phi)
pointData(3) = 0.0_wp
end subroutine ! vertexSupplierStd
subroutine vertexSupplierMcG(i, pointData)
integer, intent(in) :: i
real(wp), intent(out), dimension(3) :: pointData
real(wp) :: phi
phi = omega * real(i, kind=wp)
pointData(1) = rCorr * cos(phi)
pointData(2) = rCorr * sin(phi)
pointData(3) = 0.0_wp
end subroutine ! vertexSupplierMcG
subroutine demo_McGreivy()
use mod_testutil, only: errorMetric
!$ifdef _OPENMP
use omp_lib
!$endif
real(wp), parameter :: current = 17.0_wp ! A
real(wp), dimension(3) :: center, normal
real(wp), dimension(3,1) :: evalPos, magneticField
logical :: useCompensatedSummation
integer :: i, numCases, numPhi, numProcessors
real(wp) :: bZRef, bZStd, bZMcG, dPhi
integer, dimension(*), parameter :: allNumPhi = (/ &
10, 30, 100, 300, 1000, 3000, &
10000, 30000, 100000, 300000, &
1000000, 3000000, &
10000000, 30000000, &
100000000, 300000000, 1000000000 /)
real(wp), dimension(:), allocatable :: &
allBzStdErr, allBzMcGErr
real(wp), dimension(:,:), allocatable :: &
resultTable
useCompensatedSummation = .true.
!$ifdef _OPENMP
numProcessors = omp_get_max_threads()
print *, "using OpenMP"
!$else
! numProcessors = 1
!$endif
print *, "numProcessors = ", numProcessors
radius = 1.23_wp ! m
center = (/ 0.0_wp, 0.0_wp, 0.0_wp /)
normal = (/ 0.0_wp, 0.0_wp, 1.0_wp /)
evalPos = reshape((/ 10.0_wp, 5.0_wp, 0.0_wp /), shape(evalPos))
call magneticFieldCircularFilament(center, normal, radius, current, &
1, evalPos, magneticField)
bZRef = magneticField(3, 1)
print *, "ref B_z = ", bZRef
! mimic circular wire loop as:
! a) Polygon with points on the circule to be mimiced
! b) Polygon with points slightly offset radially outward (McGreivy correction)
! --> a) should have 2nd-order convergence;
! b) should have 4th-order convergence wrt. number of Polygon points
numCases = size(allNumPhi)
print *, "number of cases: ", numCases
allocate(allBzStdErr(numCases), &
allBzMcGErr(numCases), &
resultTable(3, numCases))
do i = 1, numCases
numPhi = allNumPhi(i)
print *, "case ",i,"/",numCases," numPhi=",numPhi
omega = 2.0_wp * PI / real(numPhi - 1, kind=wp)
call magneticFieldPolygonFilamentVertexSupplier( &
numPhi, vertexSupplierStd, current, &
1, evalPos, magneticField, &
numProcessors, useCompensatedSummation)
bZStd = magneticField(3,1)
allBzStdErr(i) = errorMetric(bZRef, bZStd)
print *, " on-poly B_z: ", bZStd, " (err ", allBzStdErr(i), ")"
!> McGreivy radius correction
!> spacing between points
dPhi = 2.0 * PI / real(numPhi - 1, kind=wp)
!> TODO: understand derivation of alpha for special case of closed circle
!> |dr/ds| = 2*pi
!> --> alpha = 1/R * (dr)^2 / 12
!> == 4 pi^2 / (12 R)
rCorr = radius * (1.0_wp + dPhi * dPhi/ 12.0_wp)
call magneticFieldPolygonFilamentVertexSupplier( &
numPhi, vertexSupplierMcG, current, &
1, evalPos, magneticField, &
numProcessors, useCompensatedSummation)
bZMcG = magneticField(3,1)
allBzMcGErr(i) = errorMetric(bZRef, bZMcG)
print *, "McGreivy B_z: ", bZMcG, " (err ", allBzMcGErr(i), ")"
resultTable(1, i) = real(numPhi, kind=wp)
resultTable(2, i) = allBzStdErr(i)
resultTable(3, i) = allBzMcGErr(i)
end do ! i = 1, numCases
! TODO: dump resultTable to output file
deallocate(allBzStdErr, allBzMcGErr, resultTable)
end subroutine ! demo_McGreivy
end module ! mod_demo_abscab
program demo_abscab
use mod_demo_abscab
implicit none
call demo_McGreivy()
end program ! demo_abscab