Subroutine that initializes the analytical or externally calculated electric and magnetic fields.
In this subroutine we load the parameters of the electric and magnetic fields from the namelists 'analytical_fields_params' and 'externalPlasmaModel' in the input file. Sign convention in analytical fields corresponds to DIII-D fields with and . Sign convention in analytical fields corresponds to DIII-D fields with and .
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| type(KORC_PARAMS), | intent(inout) | :: | params | Core KORC simulation parameters. |
||
| type(FIELDS), | intent(out) | :: | F | An instance of the KORC derived type FIELDS. |
subroutine initialize_fields(params,F)
!! @note Subroutine that initializes the analytical or externally
!! calculated electric and magnetic fields. @endnote
!! In this subroutine we load the parameters of the electric and
!! magnetic fields from the namelists 'analytical_fields_params' and
!! 'externalPlasmaModel' in the input file.
TYPE(KORC_PARAMS), INTENT(INOUT) :: params
!! Core KORC simulation parameters.
TYPE(FIELDS), INTENT(OUT) :: F
!! An instance of the KORC derived type FIELDS.
REAL(rp) :: Bo
!! Magnetic field at magnetic axis for an 'ANALYTICAL' magnetic field,
!! or the magnitude of the magnetic field for a 'UNFIROM' plasma.
REAL(rp) :: minor_radius
!! Plasma edge \(r_{edge}\) as measured from the magnetic axis.
REAL(rp) :: major_radius
!! Radial position of the magnetic axis \(R_0\)
REAL(rp) :: qa
!! Safety factor at the plasma edge.
REAL(rp) :: qo
!! Safety factor at the magnetic axis \(q_0\).
CHARACTER(MAX_STRING_LENGTH) :: current_direction
!! String with information about the direction of the plasma current,
!! 'PARALLEL' or 'ANTI-PARALLEL' to the toroidal magnetic field.
CHARACTER(MAX_STRING_LENGTH) :: E_model
REAL(rp) :: Eo,E_dyn,E_pulse,E_width
!! Electric field at the magnetic axis.
LOGICAL :: Efield
!! Logical variable that specifies if the electric field is
!! going to be used on in a given simulation.
LOGICAL :: dBfield
LOGICAL :: Bfield
!! Logical variable that specifies if the magnetic field is
!! going to be used on in a given simulation.
LOGICAL :: Bflux
LOGICAL :: Bflux3D
LOGICAL :: Dim2x1t
LOGICAL :: E_2x1t,ReInterp_2x1t
!! Logical variable that specifies if the poloidal magnetic
!! flux is going to be used on in a given simulation.
LOGICAL :: axisymmetric_fields
!! Logical variable that specifies if the plasma is axisymmetric.
INTEGER :: ii
!! Iterators for creating mesh for GC model with analytic fields
INTEGER :: kk
!! Iterators for creating mesh for GC model with analytic fields
INTEGER :: nR
!! Number of mesh points in R for grid in GC model of analytical field
INTEGER :: nZ,nPHI
!! Number of mesh points in Z for grid in GC model of analytical field
real(rp) :: rm
!! Minor radius at each position in the grid for
!! GC model of analytical field
real(rp) :: qr
!! Safety factor at each position in the grid for
!! GC model of analytical field
real(rp) :: theta
!! Poloidal angle at each position in the grid for
!! GC model of analytical field
logical :: test
integer :: res_double
real(rp) :: RMAX,RMIN,ZMAX,ZMIN
integer :: dim_1D,ind0_2x1t
real(rp) :: dt_E_SC,Ip_exp,PSIp_lim
real(rp) :: t0_2x1t
NAMELIST /analytical_fields_params/ Bo,minor_radius,major_radius,&
qa,qo,Eo,current_direction,nR,nZ,nPHI,dim_1D,dt_E_SC,Ip_exp, &
E_dyn,E_pulse,E_width
NAMELIST /externalPlasmaModel/ Efield, Bfield, Bflux,Bflux3D,dBfield, &
axisymmetric_fields, Eo,E_dyn,E_pulse,E_width,res_double, &
dim_1D,dt_E_SC,Ip_exp,PSIp_lim,Dim2x1t,t0_2x1t,E_2x1t,ReInterp_2x1t, &
ind0_2x1t
if (params%mpi_params%rank .EQ. 0) then
write(6,'(/,"* * * * * * * * INITIALIZING FIELDS * * * * * * * *")')
end if
! SELECT CASE (TRIM(params%field_model))
if (params%field_model(1:10).eq.'ANALYTICAL') then
! CASE('ANALYTICAL')
! Load the parameters of the analytical magnetic field
open(unit=default_unit_open,file=TRIM(params%path_to_inputs), &
status='OLD',form='formatted')
read(default_unit_open,nml=analytical_fields_params)
close(default_unit_open)
F%AB%Bo = Bo
F%AB%a = minor_radius
F%AB%Ro = major_radius
F%Ro = major_radius
F%Zo = 0.0_rp
F%AB%qa = qa
F%AB%qo = qo
F%AB%lambda = F%AB%a/SQRT(qa/qo - 1.0_rp)
F%AB%Bpo = F%AB%lambda*F%AB%Bo/(F%AB%qo*F%AB%Ro)
F%AB%current_direction = TRIM(current_direction)
SELECT CASE (TRIM(F%AB%current_direction))
CASE('PARALLEL')
F%AB%Bp_sign = 1.0_rp
CASE('ANTI-PARALLEL')
F%AB%Bp_sign = -1.0_rp
CASE DEFAULT
END SELECT
F%Eo = Eo
F%Bo = F%AB%Bo
F%E_dyn = E_dyn
F%E_pulse = E_pulse
F%E_width = E_width
F%res_double=res_double
if (params%mpi_params%rank .EQ. 0) then
write(6,'("ANALYTIC")')
write(6,'("Magnetic field: ",E17.10)') F%Bo
write(6,'("Electric field: ",E17.10)') F%Eo
end if
if (params%field_eval.eq.'interp') then
F%dims(1) = nR
F%dims(2) = nPHI
F%dims(3) = nZ
if (params%field_model(12:14).eq.'PSI') then
F%axisymmetric_fields = .TRUE.
F%Bfield=.TRUE.
F%Efield=.TRUE.
F%Bflux=.TRUE.
call ALLOCATE_2D_FIELDS_ARRAYS(params,F,F%Bfield,F%Bflux, &
.false.,F%Efield)
do ii=1_idef,F%dims(1)
F%X%R(ii)=(F%Ro-F%AB%a)+(ii-1)*2*F%AB%a/(F%dims(1)-1)
end do
do ii=1_idef,F%dims(3)
F%X%Z(ii)=(F%Zo-F%AB%a)+(ii-1)*2*F%AB%a/(F%dims(3)-1)
end do
! write(6,*) F%X%R
do ii=1_idef,F%dims(1)
do kk=1_idef,F%dims(3)
rm=sqrt((F%X%R(ii)-F%Ro)**2+(F%X%Z(kk)-F%Zo)**2)
qr=F%AB%qo*(1+(rm/F%AB%lambda)**2)
theta=atan2(F%X%Z(kk)-F%Zo,F%X%R(ii)-F%Ro)
F%B_2D%R(ii,kk)=(rm/F%X%R(ii))* &
(F%AB%Bo/qr)*sin(theta)
F%B_2D%PHI(ii,kk)=-(F%Ro/F%X%R(ii))*F%AB%Bo
F%B_2D%Z(ii,kk)=-(rm/F%X%R(ii))* &
(F%AB%Bo/qr)*cos(theta)
F%E_2D%R(ii,kk)=0.0_rp
F%E_2D%PHI(ii,kk)=-(F%Ro/F%X%R(ii))*F%Eo
F%E_2D%Z(ii,kk)=0.0_rp
F%PSIp(ii,kk)=F%X%R(ii)*F%AB%lambda**2*F%Bo/ &
(2*F%AB%qo*(F%Ro+rm*cos(theta)))* &
log(1+(rm/F%AB%lambda)**2)
!! Sign convention in analytical fields corresponds to
!! DIII-D fields with \(B_\phi<0\) and \(B_\theta<0\).
F%FLAG2D=1.
end do
end do
F%FLAG2D(1:2,:)=0.
F%FLAG2D(F%dims(1)-1:F%dims(1),:)=0.
F%FLAG2D(:,1:2)=0.
F%FLAG2D(:,F%dims(3)-1:F%dims(3))=0.
else if (params%field_model(12:13).eq.'3D') then
F%axisymmetric_fields = .FALSE.
F%Bfield=.TRUE.
F%Efield=.TRUE.
call ALLOCATE_3D_FIELDS_ARRAYS(params,F,F%Bfield,F%Efield,.false.)
do ii=1_idef,F%dims(1)
F%X%R(ii)=(F%Ro-F%AB%a)+(ii-1)*2*F%AB%a/(F%dims(1)-1)
end do
do ii=1_idef,F%dims(2)
F%X%PHI(ii)=0._rp+(ii-1)*2*C_PI/(F%dims(1)-1)
end do
do ii=1_idef,F%dims(3)
F%X%Z(ii)=(F%Zo-F%AB%a)+(ii-1)*2*F%AB%a/(F%dims(3)-1)
end do
!write(6,*) size(F%B_3D%R)
do ii=1_idef,F%dims(1)
do kk=1_idef,F%dims(3)
!write(6,*) ii,kk
rm=sqrt((F%X%R(ii)-F%Ro)**2+(F%X%Z(kk)-F%Zo)**2)
qr=F%AB%qo*(1+(rm/F%AB%lambda)**2)
theta=atan2(F%X%Z(kk)-F%Zo,F%X%R(ii)-F%Ro)
F%B_3D%R(ii,:,kk)=(rm/F%X%R(ii))* &
(F%AB%Bo/qr)*sin(theta)
F%B_3D%PHI(ii,:,kk)=-(F%Ro/F%X%R(ii))*F%AB%Bo
!write(6,*) F%B_3D%PHI(ii,1,kk)
F%B_3D%Z(ii,:,kk)=-(rm/F%X%R(ii))* &
(F%AB%Bo/qr)*cos(theta)
F%E_3D%R(ii,:,kk)=0.0_rp
F%E_3D%PHI(ii,:,kk)=-(F%Ro/F%X%R(ii))*F%Eo
F%E_3D%Z(ii,:,kk)=0.0_rp
!! Sign convention in analytical fields corresponds to
!! DIII-D fields with \(B_\phi<0\) and \(B_\theta<0\).
F%FLAG3D=1.
end do
end do
F%FLAG3D(1:2,:,:)=0.
F%FLAG3D(F%dims(1)-1:F%dims(1),:,:)=0.
F%FLAG3D(:,:,1:2)=0.
F%FLAG3D(:,:,F%dims(3)-1:F%dims(3))=0.
end if
if (params%orbit_model(3:5).eq.'pre') then
if (params%mpi_params%rank .EQ. 0) then
write(6,'("Initializing GC fields from analytic EM fields")')
end if
if (params%field_model(12:13).eq.'2D') then
call initialize_GC_fields(F)
else if (params%field_model(12:13).eq.'3D') then
call initialize_GC_fields_3D(F)
end if
end if
!F%Bfield= .FALSE.
!F%axisymmetric_fields = .TRUE.
!F%Bflux=.TRUE.
!F%Efield=.FALSE.
end if
if (params%SC_E) then
F%dim_1D=dim_1D
F%dt_E_SC=dt_E_SC
F%Ip_exp=Ip_exp
ALLOCATE(F%E_SC_1D%PHI(F%dim_1D))
ALLOCATE(F%A1_SC_1D%PHI(F%dim_1D))
ALLOCATE(F%A2_SC_1D%PHI(F%dim_1D))
ALLOCATE(F%A3_SC_1D%PHI(F%dim_1D))
ALLOCATE(F%J1_SC_1D%PHI(F%dim_1D))
ALLOCATE(F%J2_SC_1D%PHI(F%dim_1D))
ALLOCATE(F%J3_SC_1D%PHI(F%dim_1D))
ALLOCATE(F%r_1D(F%dim_1D))
F%E_SC_1D%PHI=0._rp
F%A1_SC_1D%PHI=0._rp
F%A2_SC_1D%PHI=0._rp
F%A3_SC_1D%PHI=0._rp
F%J1_SC_1D%PHI=0._rp
F%J2_SC_1D%PHI=0._rp
F%J3_SC_1D%PHI=0._rp
F%r_1D=0._rp
do ii=1_idef,F%dim_1D
F%r_1D(ii)=(ii-1)*F%AB%a/(F%dim_1D-1)
end do
end if
! CASE('EXTERNAL')
else if (params%field_model(1:8).eq.'EXTERNAL') then
! Load the magnetic field from an external HDF5 file
open(unit=default_unit_open,file=TRIM(params%path_to_inputs), &
status='OLD',form='formatted')
read(default_unit_open,nml=externalPlasmaModel)
close(default_unit_open)
F%Bfield = Bfield
F%dBfield = dBfield
F%Bflux = Bflux
F%Bflux3D = Bflux3D
F%Efield = Efield
F%axisymmetric_fields = axisymmetric_fields
F%Dim2x1t = Dim2x1t
F%ReInterp_2x1t = ReInterp_2x1t
F%t0_2x1t = t0_2x1t
F%ind0_2x1t = ind0_2x1t
if (params%proceed) then
call load_prev_iter(params)
F%ind_2x1t=params%prev_iter_2x1t
else
F%ind_2x1t=F%ind0_2x1t
end if
F%E_2x1t = E_2x1t
F%E_dyn = E_dyn
F%E_pulse = E_pulse
F%E_width = E_width
F%PSIp_lim=PSIp_lim
F%res_double=res_double
! write(6,'("E_dyn: ",E17.10)') E_dyn
! write(6,'("E_pulse: ",E17.10)') E_pulse
! write(6,'("E_width: ",E17.10)') E_width
call load_dim_data_from_hdf5(params,F)
!sets F%dims for 2D or 3D data
!write(6,*) F%dims
call which_fields_in_file(params,F%Bfield_in_file,F%Efield_in_file, &
F%Bflux_in_file,F%dBfield_in_file)
if (F%Bflux.AND..NOT.F%Bflux_in_file) then
write(6,'("ERROR: Magnetic flux to be used but no data in file!")')
call KORC_ABORT()
end if
if (F%Bfield.AND..NOT.F%Bfield_in_file) then
write(6,'("ERROR: Magnetic field to be used but no data in file!")')
call KORC_ABORT()
end if
if (F%dBfield.AND..NOT.F%dBfield_in_file) then
write(6,'("ERROR: differential Magnetic field to be used &
but no data in file!")')
! call KORC_ABORT()
end if
if (F%Efield.AND..NOT.F%Efield_in_file) then
if (params%mpi_params%rank.EQ.0_idef) then
write(6,'(/,"* * * * * * * * * * FIELDS * * * * * * * * * *")')
write(6,'("MESSAGE: Analytical electric field will be used.")')
write(6,'("* * * * * * * * * * * * ** * * * * * * * * * * *",/)')
end if
end if
if (F%axisymmetric_fields) then
if (F%Dim2x1t) then
call ALLOCATE_2D_FIELDS_ARRAYS(params,F,F%Bfield, &
F%Bflux,F%dBfield,F%Efield.AND.F%Efield_in_file)
call ALLOCATE_3D_FIELDS_ARRAYS(params,F,F%Bfield, &
F%Efield,F%dBfield)
else
F%Efield_in_file=.TRUE.
call ALLOCATE_2D_FIELDS_ARRAYS(params,F,F%Bfield, &
F%Bflux,F%dBfield,F%Efield.AND.F%Efield_in_file)
F%Efield_in_file=.FALSE.
end if
else
call ALLOCATE_3D_FIELDS_ARRAYS(params,F,F%Bfield,F%Efield,F%dBfield)
end if
!allocates 2D or 3D data arrays (fields and spatial)
call load_field_data_from_hdf5(params,F)
! write(6,*) F%PSIp
!write(6,*) F%E_3D%PHI(:,F%ind0_2x1t,:)
! end if
if (F%Bflux) then
F%PSIP_min=minval(F%PSIp)
else if(F%Bflux3D) then
F%PSIP_min=minval(F%PSIp3D(:,1,:))
end if
if ((.not.F%Efield_in_file).and.(.not.F%Dim2x1t)) then
F%Eo = Eo
if (F%axisymmetric_fields) then
F%E_2D%R=0._rp
do ii=1_idef,F%dims(1)
F%E_2D%PHI(ii,:)=F%Eo*F%Ro/F%X%R(ii)
end do
F%E_2D%Z=0._rp
else
F%E_3D%R=0._rp
do ii=1_idef,F%dims(1)
F%E_3D%PHI(ii,:,:)=F%Eo*F%Ro/F%X%R(ii)
end do
F%E_3D%Z=0._rp
end if
end if
if(F%dBfield.and..not.F%dBfield_in_file) then
if (F%axisymmetric_fields) then
F%dBdR_2D%R=0._rp
F%dBdR_2D%PHI=0._rp
F%dBdR_2D%Z=0._rp
F%dBdPHI_2D%R=0._rp
F%dBdPHI_2D%PHI=0._rp
F%dBdPHI_2D%Z=0._rp
F%dBdZ_2D%R=0._rp
F%dBdZ_2D%PHI=0._rp
F%dBdZ_2D%Z=0._rp
else
F%dBdR_3D%R=0._rp
F%dBdR_3D%PHI=0._rp
F%dBdR_3D%Z=0._rp
F%dBdPHI_3D%R=0._rp
F%dBdPHI_3D%PHI=0._rp
F%dBdPHI_3D%Z=0._rp
F%dBdZ_3D%R=0._rp
F%dBdZ_3D%PHI=0._rp
F%dBdZ_3D%Z=0._rp
end if
end if
if (params%mpi_params%rank .EQ. 0) then
write(6,'("EXTERNAL")')
write(6,'("Magnetic field: ",E17.10)') F%Bo
write(6,'("Electric field: ",E17.10)') F%Eo
end if
if (params%SC_E) then
F%dim_1D=dim_1D
F%dt_E_SC=dt_E_SC
F%Ip_exp=Ip_exp
call allocate_1D_FS_arrays(params,F)
call load_1D_FS_from_hdf5(params,F)
! write(6,*) F%PSIP_1D
end if
! test=.true.
! if (F%Bflux.and.(.not.test)) then
! call initialize_fields_interpolant(params,F)
! F%Bfield=.TRUE.
! F%Efield=.TRUE.
! F%Efield_in_file=.TRUE.
! RMIN=F%X%R(1)
! RMAX=F%X%R(F%dims(1))
! ZMIN=F%X%Z(1)
! ZMAX=F%X%Z(F%dims(3))
! do ii=1_idef,res_double
! F%dims(1)=2*F%dims(1)-1
! F%dims(3)=2*F%dims(3)-1
! end do
! if (res_double>0) then
! DEALLOCATE(F%X%R)
! DEALLOCATE(F%X%Z)
! DEALLOCATE(F%PSIp)
! end if
! call ALLOCATE_2D_FIELDS_ARRAYS(params,F,F%Bfield, &
! F%Bflux,F%Efield.AND.F%Efield_in_file)
! do ii=1_idef,F%dims(1)
! F%X%R(ii)=RMIN+REAL(ii-1)/REAL(F%dims(1)-1)*(RMAX-RMIN)
! end do
! do ii=1_idef,F%dims(3)
! F%X%Z(ii)=ZMIN+REAL(ii-1)/REAL(F%dims(3)-1)*(ZMAX-ZMIN)
! end do
! call calculate_initial_magnetic_field(F)
! F%E_2D%R=0._rp
! do ii=1_idef,F%dims(1)
! F%E_2D%PHI(ii,:)=F%Eo*F%Ro/F%X%R(ii)
! end do
! F%E_2D%Z=0._rp
! end if
! if (F%Bflux.and.test) then
! F%Bfield=.TRUE.
! F%Efield=.TRUE.
! F%Efield_in_file=.TRUE.
! call ALLOCATE_2D_FIELDS_ARRAYS(params,F,F%Bfield, &
! F%Bflux,F%Efield.AND.F%Efield_in_file)
! B
! edge nodes at minimum R,Z
! F%B_2D%Z(1,:)=-(F%PSIp(2,:)-F%PSIp(1,:))/(F%X%R(2)-F%X%R(1))/F%X%R(1)
! F%B_2D%R(:,1)=(F%PSIp(:,2)-F%PSIp(:,1))/(F%X%Z(2)-F%X%Z(1))/F%X%R(:)
! edge nodes at maximum R,Z
! F%B_2D%Z(F%dims(1),:)=-(F%PSIp(F%dims(1),:)-F%PSIp(F%dims(1)-1,:))/ &
! (F%X%R(F%dims(1))-F%X%R(F%dims(1)-1))/F%X%R(F%dims(1))
! F%B_2D%R(:,F%dims(3))=(F%PSIp(:,F%dims(3))-F%PSIp(:,F%dims(3)-1))/ &
! (F%X%Z(F%dims(3))-F%X%Z(F%dims(3)-1))/F%X%R(:)
! do ii=2_idef,F%dims(1)-1
! central difference over R for interior nodes for BZ
! F%B_2D%Z(ii,:)=-(F%PSIp(ii+1,:)-F%PSIp(ii-1,:))/ &
! (F%X%R(ii+1)-F%X%R(ii-1))/F%X%R(ii)
! end do
! do ii=2_idef,F%dims(3)-1
! central difference over Z for interior nodes for BR
! F%B_2D%R(:,ii)=(F%PSIp(:,ii+1)-F%PSIp(:,ii-1))/ &
! (F%X%Z(ii+1)-F%X%Z(ii-1))/F%X%R(:)
! end do
! do ii=1_idef,F%dims(1)
! F%B_2D%PHI(ii,:)=-F%Bo*F%Ro/F%X%R(ii)
! end do
! F%E_2D%R=0._rp
! do ii=1_idef,F%dims(1)
! F%E_2D%PHI(ii,:)=F%Eo*F%Ro/F%X%R(ii)
! end do
! F%E_2D%Z=0._rp
! F%Bfield=.FALSE.
if (params%mpi_params%rank.EQ.0) then
if (F%axisymmetric_fields) then
if (F%Bflux) then
write(6,'("PSIp(r=0)",E17.10)') F%PSIp(F%dims(1)/2,F%dims(3)/2)
write(6,'("BPHI(r=0)",E17.10)') F%Bo
write(6,'("EPHI(r=0)",E17.10)') F%Eo
else if (F%Bflux3D) then
write(6,'("PSIp(r=0)",E17.10)') F%PSIp3D(F%dims(1)/2,1,F%dims(3)/2)
else
write(6,'("BR(r=0)",E17.10)') F%B_2D%R(F%dims(1)/2,F%dims(3)/2)
write(6,'("BPHI(r=0)",E17.10)') &
F%B_2D%PHI(F%dims(1)/2,F%dims(3)/2)
write(6,'("BZ(r=0)",E17.10)') F%B_2D%Z(F%dims(1)/2,F%dims(3)/2)
write(6,'("EPHI(r=0)",E17.10)') &
F%E_2D%PHI(F%dims(1)/2,F%dims(3)/2)
end if
end if
end if
if (params%orbit_model(3:5).EQ.'pre') then
if (params%mpi_params%rank.eq.0) then
write(6,'("Initializing GC fields from external EM fields")')
end if
if (params%field_model(10:12).eq.'2DB') then
if (F%axisymmetric_fields) then
call initialize_GC_fields(F)
else
call initialize_GC_fields_3D(F)
end if
end if
end if
! write(6,'("gradBR",E17.10)') F%gradB_2D%R(F%dims(1)/2,F%dims(3)/2)
! write(6,'("gradBPHI",E17.10)') F%gradB_2D%PHI(F%dims(1)/2,F%dims(3)/2)
! write(6,'("gradBZ",E17.10)') F%gradB_2D%Z(F%dims(1)/2,F%dims(3)/2)
end if
if (params%mpi_params%rank.eq.0) then
write(6,'("* * * * * * * * * * * * * * * * * * * * * * * * *",/)')
end if
end subroutine initialize_fields