Computes the auxiliary fields and that are used in the RHS of the evolution equations for the GC orbit model.
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| type(FIELDS), | intent(inout) | :: | F | An instance of the KORC derived type FIELDS. |
subroutine initialize_GC_fields_3D(F)
!! Computes the auxiliary fields \(\nabla|{\bf B}|\) and
!! \(\nabla\times\hat{b}\) that are used in the RHS of the
!! evolution equations for the GC orbit model.
TYPE(FIELDS), INTENT(INOUT) :: F
!! An instance of the KORC derived type FIELDS.
INTEGER :: ii,jj
!! Iterator across F%dim
REAL(rp), DIMENSION(:,:,:),ALLOCATABLE :: Bmag
!! Magnetic field magnitude
REAL(rp), DIMENSION(:,:,:,:),ALLOCATABLE :: bhat
!! Magnetic field unit vector
Bmag=SQRT(F%B_3D%R**2+F%B_3D%PHI**2+F%B_3D%Z**2)
ALLOCATE(bhat(F%dims(1),F%dims(2),F%dims(3),3))
bhat(:,:,:,1)=F%B_3D%R/Bmag
bhat(:,:,:,2)=F%B_3D%PHI/Bmag
bhat(:,:,:,3)=F%B_3D%Z/Bmag
! Single-sided difference for axiliary fields at edge nodes
! Differential over R on first index, differential over Z
! on second index.
F%gradB_3D%R=0._rp
F%gradB_3D%PHI=0._rp
F%gradB_3D%Z=0._rp
F%curlb_3D%R=0._rp
F%curlb_3D%PHI=0._rp
F%curlb_3D%Z=0._rp
! gradB
! edge nodes at minimum R,Z
F%gradB_3D%R(1,:,:)=F%gradB_3D%R(1,:,:)+ &
(Bmag(2,:,:)-Bmag(1,:,:))/(F%X%R(2)-F%X%R(1))
do ii=1_idef,F%dims(1)
F%gradB_3D%PHI(ii,1,:)=F%gradB_3D%PHI(ii,1,:)+ &
(Bmag(ii,2,:)-Bmag(ii,F%dims(2),:))/ &
(F%X%R(ii)*(F%X%PHI(2)-F%X%PHI(F%dims(2))))
end do
F%gradB_3D%Z(:,:,1)=F%gradB_3D%Z(:,:,1)+ &
(Bmag(:,:,2)-Bmag(:,:,1))/(F%X%Z(2)-F%X%Z(1))
! edge nodes at maximum R,Z
F%gradB_3D%R(F%dims(1),:,:)=F%gradB_3D%R(F%dims(1),:,:)+ &
(Bmag(F%dims(1),:,:)-Bmag(F%dims(1)-1,:,:))/ &
(F%X%R(F%dims(1))-F%X%R(F%dims(1)-1))
do ii=1_idef,F%dims(1)
F%gradB_3D%PHI(ii,F%dims(2),:)=F%gradB_3D%PHI(ii,F%dims(2),:)+ &
(Bmag(ii,1,:)-Bmag(ii,F%dims(2)-1,:))/ &
(F%X%R(ii)*(F%X%PHI(1)-F%X%PHI(F%dims(2)-1)))
end do
F%gradB_3D%Z(:,:,F%dims(3))=F%gradB_3D%Z(:,:,F%dims(3))+ &
(Bmag(:,:,F%dims(3))-Bmag(:,:,F%dims(3)-1))/ &
(F%X%Z(F%dims(3))-F%X%Z(F%dims(3)-1))
! curlb
! edge nodes at minimum R,PHI,Z
! R component has differential over PHI and Z
do ii=1_idef,F%dims(1)
F%curlb_3D%R(ii,1,:)=F%curlb_3D%R(ii,1,:)+ &
(bhat(ii,2,:,3)-bhat(ii,F%dims(2),:,3))/ &
(F%X%R(ii)*(F%X%PHI(2)-F%X%PHI(F%dims(2))))
end do
F%curlb_3D%R(:,:,1)=F%curlb_3D%R(:,:,1)-&
(bhat(:,:,2,2)-bhat(:,:,1,2))/(F%X%Z(2)-F%X%Z(1))
! PHI component has differentials over R and Z
F%curlb_3D%PHI(1,:,:)=F%curlb_3D%PHI(1,:,:)- &
(bhat(2,:,:,3)-bhat(1,:,:,3))/ &
(F%X%R(2)-F%X%R(1))
F%curlb_3D%PHI(:,:,1)=F%curlb_3D%PHI(:,:,1)+ &
((bhat(:,:,2,1)-bhat(:,:,1,1))/(F%X%Z(2)-F%X%Z(1)))
! Z component has differentials over R and PHI
F%curlb_3D%Z(1,:,:)=F%curlb_3D%Z(1,:,:)+ &
((bhat(2,:,:,2)*F%X%R(2)- &
bhat(1,:,:,2)*F%X%R(1))/(F%X%R(2)-F%X%R(1)))/F%X%R(1)
do ii=1_idef,F%dims(1)
F%curlb_3D%Z(ii,1,:)=F%curlb_3D%Z(ii,1,:)-&
(bhat(ii,2,:,1)-bhat(ii,F%dims(2),:,1))/ &
(F%X%R(ii)*(F%X%PHI(2)-F%X%PHI(F%dims(2))))
end do
! edge nodes at maximum R,PHI,Z
! R component has differential over PHI and Z
do ii=1_idef,F%dims(1)
F%curlb_3D%R(ii,F%dims(2),:)=F%curlb_3D%R(ii,F%dims(2),:)+ &
(bhat(ii,1,:,3)-bhat(ii,F%dims(2)-1,:,3))/ &
(F%X%R(ii)*(F%X%PHI(1)-F%X%PHI(F%dims(2)-1)))
end do
F%curlb_3D%R(:,:,F%dims(3))=F%curlb_3D%R(:,:,F%dims(3)) &
-(bhat(:,:,F%dims(3),2)-bhat(:,:,F%dims(3)-1,2))/ &
(F%X%Z(F%dims(3))-F%X%Z(F%dims(3)-1))
! PHI component has differentials over R and Z
F%curlb_3D%PHI(F%dims(1),:,:)=F%curlb_3D%PHI(F%dims(1),:,:)- &
(bhat(F%dims(1),:,:,3)-bhat(F%dims(1)-1,:,:,3))/ &
(F%X%R(F%dims(1))-F%X%R(F%dims(1)-1))
F%curlb_3D%PHI(:,:,F%dims(3))=F%curlb_3D%PHI(:,:,F%dims(3))+ &
((bhat(:,:,F%dims(3),1)-bhat(:,:,F%dims(3)-1,1))/ &
(F%X%Z(F%dims(3))-F%X%Z(F%dims(3)-1)))
! Z component has differentials over R and PHI
F%curlb_3D%Z(F%dims(1),:,:)=F%curlb_3D%Z(F%dims(1),:,:)+ &
((bhat(F%dims(1),:,:,2)*F%X%R(F%dims(1))- &
bhat(F%dims(1)-1,:,:,2)*F%X%R(F%dims(1)-1))/(F%X%R(F%dims(1))- &
F%X%R(F%dims(1)-1)))/F%X%R(F%dims(1))
do ii=1_idef,F%dims(1)
F%curlb_3D%Z(ii,F%dims(2),:)=F%curlb_3D%Z(ii,F%dims(2),:)-&
(bhat(ii,1,:,1)-bhat(ii,F%dims(2)-1,:,1))/ &
(F%X%R(ii)*(F%X%PHI(1)-F%X%PHI(F%dims(2)-1)))
end do
do ii=2_idef,F%dims(1)-1
! central difference over R for interior nodes
F%gradB_3D%R(ii,:,:)=F%gradB_3D%R(ii,:,:)+ &
(Bmag(ii+1,:,:)-Bmag(ii-1,:,:))/ &
(F%X%R(ii+1)-F%X%R(ii-1))
F%curlb_3D%Z(ii,:,:)=F%curlb_3D%Z(ii,:,:)+ &
((bhat(ii+1,:,:,2)*F%X%R(ii+1)- &
bhat(ii-1,:,:,2)*F%X%R(ii-1))/(F%X%R(ii+1)-F%X%R(ii-1)))/ &
F%X%R(ii)
F%curlb_3D%PHI(ii,:,:)=F%curlb_3D%PHI(ii,:,:)- &
(bhat(ii+1,:,:,3)-bhat(ii-1,:,:,3))/ &
(F%X%R(ii+1)-F%X%R(ii-1))
end do
do ii=2_idef,F%dims(2)-1
! central difference over PHI for interior nodes
do jj=1_idef,F%dims(1)
F%gradB_3D%PHI(jj,ii,:)=F%gradB_3D%PHI(jj,ii,:)+ &
(Bmag(jj,ii+1,:)-Bmag(jj,ii-1,:))/ &
(F%X%R(jj)*(F%X%PHI(ii+1)-F%X%PHI(ii-1)))
F%curlb_3D%Z(jj,ii,:)=F%curlb_3D%Z(jj,ii,:)-&
(bhat(jj,ii+1,:,1)-bhat(jj,ii-1,:,1))/ &
(F%X%R(jj)*(F%X%PHI(ii+1)-F%X%PHI(ii-1)))
F%curlb_3D%R(jj,ii,:)=F%curlb_3D%R(jj,ii,:)+ &
(bhat(jj,ii+1,:,3)-bhat(jj,ii-1,:,3))/ &
(F%X%R(jj)*(F%X%PHI(ii+1)-F%X%PHI(ii-1)))
end do
end do
do ii=2_idef,F%dims(3)-1
! central difference over Z for interior nodes
F%gradB_3D%Z(:,:,ii)=F%gradB_3D%Z(:,:,ii)+ &
(Bmag(:,:,ii+1)-Bmag(:,:,ii-1))/ &
(F%X%Z(ii+1)-F%X%Z(ii-1))
F%curlb_3D%R(:,:,ii)=F%curlb_3D%R(:,:,ii)- &
(bhat(:,:,ii+1,2)-bhat(:,:,ii-1,2))/ &
(F%X%Z(ii+1)-F%X%Z(ii-1))
F%curlb_3D%PHI(:,:,ii)=F%curlb_3D%PHI(:,:,ii)+ &
((bhat(:,:,ii+1,1)-bhat(:,:,ii-1,1))/(F%X%Z(ii+1)-F%X%Z(ii-1)))
end do
DEALLOCATE(Bmag)
DEALLOCATE(bhat)
end subroutine initialize_GC_fields_3D