-
Notifications
You must be signed in to change notification settings - Fork 1
/
bnormal.f90
205 lines (197 loc) · 11.4 KB
/
bnormal.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
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
!title (bnormal) ! Calculate total bnormal and its derivatives.
!latex \briefly{Calculate the total bnormal of all coils on plasma surface and the derivatives with respect to coil geometry and currents, including the first and second dirivatives.
!latex Calling \emph{bnormal(0), bnormal(1), bnormal(2)} calculates the $0-order$, $1^{st}-order$ and $2^{nd}-order$ derivatives respectively.}
!latex \calledby{\link{costfun}}
!latex \calls{\link{bfield}}
!latex \tableofcontents
!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
subroutine bnormal( ideriv )
!------------------------------------------------------------------------------------------------------
! DATE: 04/02/2017;
! Calculate the Bn surface integral and its derivatives;
! ideriv = 0 -> only calculate the Bn surface integral;
! ideriv = 1 -> calculate the Bn surface integral and its first derivatives;
! ideriv = 2 -> calculate the Bn surface integral and its first & second derivatives;
!------------------------------------------------------------------------------------------------------
use globals, only: dp, zero, half, one, pi2, sqrtmachprec, bsconstant, ncpu, myid, ounit, &
coil, DoF, surf, Ncoils, Nteta, Nzeta, discretefactor, plasma, &
bnorm, t1B, bn, Ndof, Cdof, weight_bharm, case_bnormal, &
weight_bnorm, ibnorm, mbnorm, ibharm, mbharm, LM_fvec, LM_fjac, &
bharm, t1H, Bmnc, Bmns, wBmn, tBmnc, tBmns, Bmnim, Bmnin, NBmn, MPI_COMM_FOCUS
use bnorm_mod
use bharm_mod
use mpi
implicit none
INTEGER, INTENT(in) :: ideriv
!--------------------------------------------------------------------------------------------
INTEGER :: astat, ierr, request
INTEGER :: icoil, iteta, jzeta, idof, ND, NumGrid, isurf
!--------------------------initialize and allocate arrays-------------------------------------
isurf = plasma
NumGrid = Nteta*Nzeta
! reset to zero;
bnorm = zero
surf(isurf)%Bx = zero; surf(isurf)%By = zero; surf(isurf)%Bz = zero; surf(isurf)%Bn = zero
dBx = zero; dBy = zero; dBz = zero; Bm = zero
bn = zero
call MPI_BARRIER( MPI_COMM_FOCUS, ierr )
!-------------------------------calculate Bn--------------------------------------------------
if( ideriv >= 0 ) then
do jzeta = 0, Nzeta - 1
do iteta = 0, Nteta - 1
if( myid.ne.modulo(jzeta*Nteta+iteta,ncpu) ) cycle ! parallelization loop;
do icoil = 1, Ncoils
call bfield0(icoil, surf(isurf)%xx(iteta, jzeta), surf(isurf)%yy(iteta, jzeta), &
& surf(isurf)%zz(iteta, jzeta), dBx(0,0), dBy(0,0), dBz(0,0))
surf(isurf)%Bx(iteta, jzeta) = surf(isurf)%Bx(iteta, jzeta) + dBx( 0, 0)
surf(isurf)%By(iteta, jzeta) = surf(isurf)%By(iteta, jzeta) + dBy( 0, 0)
surf(isurf)%Bz(iteta, jzeta) = surf(isurf)%Bz(iteta, jzeta) + dBz( 0, 0)
enddo ! end do icoil
surf(isurf)%Bn(iteta, jzeta) = surf(isurf)%Bx(iteta, jzeta)*surf(isurf)%nx(iteta, jzeta) &
& + surf(isurf)%By(iteta, jzeta)*surf(isurf)%ny(iteta, jzeta) &
& + surf(isurf)%Bz(iteta, jzeta)*surf(isurf)%nz(iteta, jzeta) &
& - surf(isurf)%pb(iteta, jzeta)
select case (case_bnormal)
case (0) ! no normalization over |B|;
bnorm = bnorm + surf(isurf)%Bn(iteta, jzeta) * surf(isurf)%Bn(iteta, jzeta) * surf(isurf)%ds(iteta, jzeta)
case (1) ! normalized over |B|;
Bm(iteta, jzeta) = surf(isurf)%Bx(iteta, jzeta)*surf(isurf)%Bx(iteta, jzeta) &
& + surf(isurf)%By(iteta, jzeta)*surf(isurf)%By(iteta, jzeta) &
& + surf(isurf)%Bz(iteta, jzeta)*surf(isurf)%Bz(iteta, jzeta)
bnorm = bnorm + surf(isurf)%Bn(iteta, jzeta) * surf(isurf)%Bn(iteta, jzeta) &
& / Bm(iteta, jzeta) * surf(isurf)%ds(iteta, jzeta)
case default
FATAL( bnorm, .true., case_bnormal can only be 0/1 )
end select
enddo ! end do iteta
enddo ! end do jzeta
! gather data
call MPI_BARRIER( MPI_COMM_FOCUS, ierr )
call MPI_ALLREDUCE( MPI_IN_PLACE, surf(isurf)%Bx, NumGrid, MPI_DOUBLE_PRECISION, MPI_SUM, MPI_COMM_FOCUS, ierr )
call MPI_ALLREDUCE( MPI_IN_PLACE, surf(isurf)%By, NumGrid, MPI_DOUBLE_PRECISION, MPI_SUM, MPI_COMM_FOCUS, ierr )
call MPI_ALLREDUCE( MPI_IN_PLACE, surf(isurf)%Bz, NumGrid, MPI_DOUBLE_PRECISION, MPI_SUM, MPI_COMM_FOCUS, ierr )
call MPI_ALLREDUCE( MPI_IN_PLACE, surf(isurf)%Bn, NumGrid, MPI_DOUBLE_PRECISION, MPI_SUM, MPI_COMM_FOCUS, ierr )
call MPI_ALLREDUCE( MPI_IN_PLACE, bnorm, 1 , MPI_DOUBLE_PRECISION, MPI_SUM, MPI_COMM_FOCUS, ierr )
call MPI_BARRIER( MPI_COMM_FOCUS, ierr )
bnorm = bnorm * half * discretefactor
bn = surf(isurf)%Bn + surf(isurf)%pb ! bn is B.n from coils
! bn = surf(isurf)%Bx * surf(isurf)%nx + surf(isurf)%By * surf(isurf)%ny + surf(isurf)%Bz * surf(isurf)%nz
!! if (case_bnormal == 0) bnorm = bnorm * bsconstant * bsconstant ! take bsconst back
! collect |B|
if (case_bnormal == 1) then
call MPI_ALLREDUCE( MPI_IN_PLACE, Bm, NumGrid, MPI_DOUBLE_PRECISION, MPI_SUM, MPI_COMM_FOCUS, ierr )
!! bm = bm * bsconstant * bsconstant
endif
! LM required discrete cost functions
if (mbnorm > 0) then
select case (case_bnormal)
case (0) ! no normalization over |B|;
LM_fvec(ibnorm+1:ibnorm+mbnorm) = weight_bnorm &
& * reshape(surf(isurf)%bn(0:Nteta-1, 0:Nzeta-1) , (/Nteta*Nzeta/))
case (1) ! normalized over |B|;
LM_fvec(ibnorm+1:ibnorm+mbnorm) = weight_bnorm &
& * reshape(surf(isurf)%bn(0:Nteta-1, 0:Nzeta-1)/sqrt(bm(0:Nteta-1, 0:Nzeta-1)), (/Nteta*Nzeta/))
case default
FATAL( bnorm, .true., case_bnormal can only be 0/1 )
end select
endif
! Bn harmonics related
if (weight_bharm > sqrtmachprec) then
call twodft( bn, Bmns, Bmnc, Bmnim, Bmnin, NBmn ) ! Bn from coils
bharm = half * sum( wBmn * ((Bmnc - tBmnc)**2 + (Bmns - tBmns)**2) )
if (mbharm > 0) then
LM_fvec(ibharm+1:ibharm+mbharm/2) = weight_bharm * wBmn * (Bmnc - tBmnc)
LM_fvec(ibharm+mbharm/2+1:ibharm+mbharm) = weight_bharm * wBmn * (Bmns - tBmns)
endif
endif
endif
!-------------------------------calculate Bn/x------------------------------------------------
if ( ideriv >= 1 ) then
! reset data
t1B = zero ; d1B = zero
dBn = zero ; dBm = zero
do jzeta = 0, Nzeta - 1
do iteta = 0, Nteta - 1
if( myid.ne.modulo(jzeta*Nteta+iteta,ncpu) ) cycle ! parallelization loop;
idof = 0
do icoil = 1, Ncoils
ND = DoF(icoil)%ND
! derivatives w.r.t currents
if ( coil(icoil)%Ic /= 0 ) then
call bfield0(icoil, surf(isurf)%xx(iteta, jzeta), surf(isurf)%yy(iteta, jzeta), &
& surf(isurf)%zz(iteta, jzeta), dBx(0,0), dBy(0,0), dBz(0,0))
if (coil(icoil)%type == 3) dBz(0,0) = zero ! Bz doesn't change in type=3
dBn(idof+1) = ( dBx(0,0)*surf(isurf)%nx(iteta,jzeta) &
& + dBy(0,0)*surf(isurf)%ny(iteta,jzeta) &
& + dBz(0,0)*surf(isurf)%nz(iteta,jzeta) ) / coil(icoil)%I
if (case_bnormal == 1) then ! normalized over |B|;
dBm(idof+1) = ( dBx(0,0)*surf(isurf)%Bx(iteta,jzeta) &
& + dBy(0,0)*surf(isurf)%By(iteta,jzeta) &
& + dBz(0,0)*surf(isurf)%Bz(iteta,jzeta) ) / coil(icoil)%I
endif
idof = idof +1
endif
! derivatives w.r.t geometries
if ( coil(icoil)%Lc /= 0 ) then
call bfield1(icoil, surf(isurf)%xx(iteta, jzeta), surf(isurf)%yy(iteta, jzeta), &
& surf(isurf)%zz(iteta, jzeta), dBx(1:ND,0), dBy(1:ND,0), dBz(1:ND,0), ND)
dBn(idof+1:idof+ND) = ( dBx(1:ND,0)*surf(isurf)%nx(iteta,jzeta) &
& + dBy(1:ND,0)*surf(isurf)%ny(iteta,jzeta) &
& + dBz(1:ND,0)*surf(isurf)%nz(iteta,jzeta) )
if (case_bnormal == 1) then ! normalized over |B|;
dBm(idof+1:idof+ND) = ( dBx(1:ND,0)*surf(isurf)%Bx(iteta,jzeta) &
& + dBy(1:ND,0)*surf(isurf)%By(iteta,jzeta) &
& + dBz(1:ND,0)*surf(isurf)%Bz(iteta,jzeta) )
endif
idof = idof + ND
endif
enddo !end icoil;
FATAL( bnormal , idof .ne. Ndof, counting error in packing )
select case (case_bnormal)
case (0) ! no normalization over |B|;
t1B(1:Ndof) = t1B(1:Ndof) + surf(isurf)%bn(iteta,jzeta) * surf(isurf)%ds(iteta,jzeta) * dBn(1:Ndof)
d1B(1:Ndof, iteta, jzeta) = d1B(1:Ndof, iteta, jzeta) + dBn(1:Ndof)
case (1) ! normalized over |B|;
t1B(1:Ndof) = t1B(1:Ndof) + ( surf(isurf)%Bn(iteta,jzeta) * dBn(1:Ndof) &
& / bm(iteta, jzeta) &
& - surf(isurf)%Bn(iteta,jzeta) * surf(isurf)%Bn(iteta,jzeta) &
& / (bm(iteta, jzeta)*bm(iteta, jzeta)) &
& * dBm(1:Ndof) ) * surf(isurf)%ds(iteta,jzeta)
d1B(1:Ndof, iteta, jzeta) = d1B(1:Ndof, iteta, jzeta) + dBn(1:Ndof) &
& / sqrt(bm(iteta, jzeta)) &
& - surf(isurf)%Bn(iteta,jzeta) * dBm(1:Ndof) &
& / (bm(iteta, jzeta) * sqrt(bm(iteta, jzeta)))
case default
FATAL( bnorm, .true., case_bnormal can only be 0/1 )
end select
enddo !end iteta;
enddo !end jzeta;
! gather data
call MPI_BARRIER( MPI_COMM_FOCUS, ierr )
call MPI_ALLREDUCE( MPI_IN_PLACE, t1B, Ndof , MPI_DOUBLE_PRECISION, MPI_SUM, MPI_COMM_FOCUS, ierr )
call MPI_ALLREDUCE( MPI_IN_PLACE, d1B, Ndof*NumGrid, MPI_DOUBLE_PRECISION, MPI_SUM, MPI_COMM_FOCUS, ierr )
t1B = t1B * discretefactor
! LM discrete derivatives
if (mbnorm > 0) then
do idof = 1, Ndof
LM_fjac(ibnorm+1:ibnorm+mbnorm, idof) = weight_bnorm &
& * reshape(d1B(idof, 0:Nteta-1, 0:Nzeta-1), (/Nteta*Nzeta/))
enddo
endif
! derivatives for Bn harmonics
if (weight_bharm > sqrtmachprec) then
dBc = zero ; dBs = zero
do idof = 1, Ndof
call twodft( d1B(idof, 0:Nteta-1, 0:Nzeta-1), dBs, dBc, Bmnim, Bmnin, NBmn )
t1H(idof) = sum( wBmn * ( (Bmnc - tBmnc)*dBc + (Bmns - tBmns)*dBs ) )
if (mbharm > 0) then
LM_fjac(ibharm+1 :ibharm+mbharm/2, idof) = weight_bharm * wBmn * dBc
LM_fjac(ibharm+mbharm/2+1:ibharm+mbharm , idof) = weight_bharm * wBmn * dBs
endif
enddo
endif
endif
!--------------------------------------------------------------------------------------------
call MPI_BARRIER( MPI_COMM_FOCUS, ierr )
return
end subroutine bnormal