@brief Subroutine that calculates the axisymmetric magnetic field to the particles' position using the poloidal magnetic flux. @details When the poloidal magnetic flux is used in a KORC simulation, the magnetic field components are calculated as it follows:
where and are the radial position of the magnetic axis and the magnetic field as measured at the magnetic axis, respectively. First, the derivatives of the poloidal magnetic flux are calculated at the particles' position using the PSPLINE interpolant of the poloidal magnetic flux. Then, we calculate the cylindrical components of the magnetic field, and finally we calculate its Cartesian components that will be used in the particle pusher.
@param[in] Y Particles' position in cylindrical coordinates, Y(1,:) = , Y(2,:) = , and Y(3,:) = . @param[in] F An instance of KORC's derived type FIELDS containing all the information about the fields used in the simulation. See korc_types.f90 and korc_fields.f90. @param[in,out] B Cartesian components of interpolated magnetic field components. B(1,:)=, B(2,:)=, and B(3,:)=. @param flag Flag that indicates whether particles are followed in the simulation (flag=1), or not (flag=0). @param A Variable containing the partial derivatives of the poloidal magnetic flux and the cylindrical components of the magnetic field (its value changes through the subroutine). @param pp Particle iterator. @param ss Species iterator.
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| type(KORC_PARAMS), | intent(in) | :: | params | |||
| real(kind=rp), | intent(in), | DIMENSION(:,:), ALLOCATABLE | :: | Y | ||
| type(FIELDS), | intent(in) | :: | F | |||
| real(kind=rp), | intent(inout), | DIMENSION(:,:), ALLOCATABLE | :: | B | ||
| real(kind=rp), | intent(inout), | DIMENSION(:,:), ALLOCATABLE | :: | E | ||
| real(kind=rp), | intent(inout), | DIMENSION(:), ALLOCATABLE | :: | PSI_P | ||
| integer(kind=is), | intent(inout), | DIMENSION(:), ALLOCATABLE | :: | flag |
subroutine calculate_magnetic_field(params,Y,F,B,E,PSI_P,flag)
TYPE(KORC_PARAMS), INTENT(IN) :: params
REAL(rp), DIMENSION(:,:), ALLOCATABLE, INTENT(IN) :: Y
TYPE(FIELDS), INTENT(IN) :: F
REAL(rp), DIMENSION(:,:), ALLOCATABLE, INTENT(INOUT) :: B
REAL(rp), DIMENSION(:,:), ALLOCATABLE, INTENT(INOUT) :: E
REAL(rp), DIMENSION(:), ALLOCATABLE, INTENT(INOUT) :: PSI_P
INTEGER(is), DIMENSION(:), ALLOCATABLE, INTENT(INOUT) :: flag
REAL(rp), DIMENSION(:,:), ALLOCATABLE :: A
INTEGER :: pp
INTEGER :: ss
REAL(rp) :: psip_conv
psip_conv=F%psip_conv
if (size(Y,1).eq.1) then
ss = size(Y,1)
else
if (Y(2,1).eq.0) then
ss=1_idef
else
ss = size(Y,1)
end if
endif
ALLOCATE(A(ss,3))
A=0._rp
if(F%Dim2x1t.and.(.not.F%ReInterp_2x1t)) then
!$OMP PARALLEL DO FIRSTPRIVATE(ss) PRIVATE(pp,ezerr) &
!$OMP& SHARED(F,Y,A,B,flag,bfield_2X1T,PSI_P)
do pp=1_idef,ss
! write(output_unit_write,'("pp: ",I16)') pp
! write(output_unit_write,'("Y_R: ",E17.10)') Y(:,1)
! write(output_unit_write,'("Y_PHI: ",E17.10)') Y(:,2)
! write(output_unit_write,'("Y_Z: ",E17.10)') Y(:,3)
call EZspline_interp(bfield_2X1T%A, Y(pp,1), F%t0_2x1t, Y(pp,3), &
PSI_P(pp), ezerr)
call EZspline_error(ezerr)
! write(output_unit_write,'("PSI_P: ",E17.10)') PSI_P(1)
! FR = (dA/dZ)/R
call EZspline_derivative(bfield_2X1T%A, 0, 0, 1, Y(pp,1), F%t0_2x1t, &
Y(pp,3), A(pp,1), ezerr)
! call EZspline_error(ezerr)
if (ezerr .NE. 0) then ! We flag the particle as lost
flag(pp) = 0_is
else
! write(output_unit_write,'("R*B_R: ",E17.10)') A(pp,1)
!if(params%SC_E) A(pp,1)=A(pp,1)/(2*C_PI)
A(pp,1) = psip_conv*A(pp,1)/Y(pp,1)
! FPHI = Fo*Ro/R
A(pp,2) = -F%Bo*F%Ro/Y(pp,1)
! FR = -(dA/dR)/R
call EZspline_derivative(bfield_2X1T%A, 1, 0, 0, Y(pp,1), &
F%t0_2x1t, Y(pp,3), A(pp,3), ezerr)
call EZspline_error(ezerr)
! write(output_unit_write,'("R*B_Z: ",E17.10)') A(pp,3)
!if(params%SC_E) A(pp,3)=A(pp,3)/(2*C_PI)
call EZspline_derivative(bfield_2X1T%A, 0, 1, 0, Y(pp,1), &
F%t0_2x1t, Y(pp,3), E(pp,2), ezerr)
E(pp,2) = E(pp,2)/(2*C_PI*Y(pp,1))
A(pp,3) = -psip_conv*A(pp,3)/Y(pp,1)
if (.not.params%GC_coords) then
B(pp,1) = A(pp,1)*COS(Y(pp,2)) - A(pp,2)*SIN(Y(pp,2))
B(pp,2) = A(pp,1)*SIN(Y(pp,2)) + A(pp,2)*COS(Y(pp,2))
B(pp,3) = A(pp,3)
E(pp,1) = -E(pp,2)*sin(Y(pp,2))
E(pp,2) = E(pp,2)*cos(Y(pp,2))
else
B(pp,1) = A(pp,1)
B(pp,2) = A(pp,2)
B(pp,3) = A(pp,3)
end if
end if
end do
!$OMP END PARALLEL DO
else
!$OMP PARALLEL DO FIRSTPRIVATE(ss) PRIVATE(pp,ezerr) &
!$OMP& SHARED(F,Y,A,B,flag,bfield_2d,PSI_P)
do pp=1_idef,ss
! write(output_unit_write,'("pp: ",I16)') pp
! write(output_unit_write,'("Y_R: ",E17.10)') Y(:,1)
! write(output_unit_write,'("Y_PHI: ",E17.10)') Y(:,2)
! write(output_unit_write,'("Y_Z: ",E17.10)') Y(:,3)
call EZspline_interp(bfield_2d%A, Y(pp,1), Y(pp,3), &
PSI_P(pp), ezerr)
call EZspline_error(ezerr)
! write(output_unit_write,'("PSI_P: ",E17.10)') PSI_P(1)
! FR = (dA/dZ)/R
call EZspline_derivative(bfield_2d%A, 0, 1, Y(pp,1), Y(pp,3), &
A(pp,1), ezerr)
! call EZspline_error(ezerr)
if (ezerr .NE. 0) then ! We flag the particle as lost
flag(pp) = 0_is
else
! write(output_unit_write,'("R*B_R: ",E17.10)') A(pp,1)
!if(params%SC_E) A(pp,1)=A(pp,1)/(2*C_PI)
A(pp,1) = psip_conv*A(pp,1)/Y(pp,1)
! FPHI = Fo*Ro/R
A(pp,2) = -F%Bo*F%Ro/Y(pp,1)
! FR = -(dA/dR)/R
call EZspline_derivative(bfield_2d%A, 1, 0, Y(pp,1), Y(pp,3), &
A(pp,3), ezerr)
call EZspline_error(ezerr)
! write(output_unit_write,'("R*B_Z: ",E17.10)') A(pp,3)
!if(params%SC_E) A(pp,3)=A(pp,3)/(2*C_PI)
A(pp,3) = -psip_conv*A(pp,3)/Y(pp,1)
if (.not.params%GC_coords) then
B(pp,1) = A(pp,1)*COS(Y(pp,2)) - A(pp,2)*SIN(Y(pp,2))
B(pp,2) = A(pp,1)*SIN(Y(pp,2)) + A(pp,2)*COS(Y(pp,2))
B(pp,3) = A(pp,3)
else
B(pp,1) = A(pp,1)
B(pp,2) = A(pp,2)
B(pp,3) = A(pp,3)
end if
end if
end do
!$OMP END PARALLEL DO
end if
! write(output_unit_write,'("calculate_fields")')
! write(output_unit_write,'("B_R: ",E17.10)') A(:,1)
! write(output_unit_write,'("B_PHI: ",E17.10)') A(:,2)
! write(output_unit_write,'("B_Z: ",E17.10)') A(:,3)
! write(output_unit_write,'("B_X: ",E17.10)') B(:,1)
! write(output_unit_write,'("B_Y: ",E17.10)') B(:,2)
! write(output_unit_write,'("B_Z: ",E17.10)') B(:,3)
DEALLOCATE(A)
end subroutine calculate_magnetic_field