FO_init – KORC

public subroutine FO_init(params, F, spp, output, step)

Arguments

Type	Intent	Attributes		Name
type(KORC_PARAMS),	intent(in)		::	params	Core KORC simulation parameters.
type(FIELDS),	intent(in)		::	F	An instance of the KORC derived type FIELDS.
type(SPECIES),	intent(inout),	DIMENSION(:), ALLOCATABLE	::	spp	An instance of the derived type SPECIES containing all the parameters and simulation variables of the different species in the simulation.
logical,	intent(in)		::	output
logical,	intent(in)		::	step
Calls

Help
Called by

Help
Source Code

FO_init
Source Code

  subroutine FO_init(params,F,spp,output,step)
    TYPE(KORC_PARAMS), INTENT(IN)                              :: params
    !! Core KORC simulation parameters.
    TYPE(FIELDS), INTENT(IN)                                   :: F
    !! An instance of the KORC derived type FIELDS.
    TYPE(SPECIES), DIMENSION(:), ALLOCATABLE, INTENT(INOUT)    :: spp
    !! An instance of the derived type SPECIES containing all the parameters
    !! and simulation variables of the different species in the simulation.

    REAL(rp)                                      :: dt
    !! Time step used in the leapfrog step (\(\Delta t\)).

    REAL(rp), DIMENSION(params%pchunk)               :: Bmag
    REAL(rp), DIMENSION(params%pchunk)               :: b_unit_X,b_unit_Y,b_unit_Z
    REAL(rp), DIMENSION(params%pchunk)               :: v,vpar,vperp
    REAL(rp), DIMENSION(params%pchunk)               :: tmp
    REAL(rp), DIMENSION(params%pchunk)               :: g
    REAL(rp), DIMENSION(params%pchunk)               :: cross_X,cross_Y,cross_Z
    REAL(rp), DIMENSION(params%pchunk)               :: vec_X,vec_Y,vec_Z
    INTEGER                                      :: ii
    !! Species iterator.
    INTEGER                                      :: pp
    !! Particles iterator.
    INTEGER                                      :: cc,pchunk
    !! Chunk iterator.

    LOGICAL,intent(in) :: output
    LOGICAL,intent(in) :: step   

    REAL(rp),DIMENSION(params%pchunk) :: X_X,X_Y,X_Z
    REAL(rp),DIMENSION(params%pchunk) :: V_X,V_Y,V_Z
    REAL(rp),DIMENSION(params%pchunk) :: Y_R,Y_PHI,Y_Z
    REAL(rp),DIMENSION(params%pchunk) :: B_X,B_Y,B_Z
    REAL(rp),DIMENSION(params%pchunk) :: E_X,E_Y,E_Z
    REAL(rp),DIMENSION(params%pchunk) :: PSIp
    REAL(rp) :: m_cache,q_cache,B0,EF0,lam,R0,q0,ar
    TYPE(C_PTR),DIMENSION(params%pchunk) :: hint
    INTEGER(is) ,DIMENSION(params%pchunk) :: flagCon,flagCol

    pchunk=params%pchunk

    B0=F%Bo
    EF0=F%Eo
    lam=F%AB%lambda
    R0=F%AB%Ro
    q0=F%AB%qo
    ar=F%AB%a

    do ii = 1_idef,params%num_species
       
       m_cache=spp(ii)%m
       q_cache=spp(ii)%q
       
       if(output) then

          !$OMP PARALLEL DO default(none) &
          !$OMP firstprivate(E0,m_cache,q_cache,B0,EF0,lam,R0,q0,ar,pchunk) &
          !$OMP& shared(params,ii,spp,F) &
          !$OMP& PRIVATE(pp,cc,X_X,X_Y,X_Z,B_X,B_Y,B_Z,V_X,V_Y,V_Z, &
          !$OMP& E_X,E_Y,E_Z,Y_R,Y_PHI,Y_Z,flagCon,flagCol,PSIp,hint,Bmag, &
          !$OMP& b_unit_X,b_unit_Y,b_unit_Z,v,vpar,vperp,tmp, &
          !$OMP& cross_X,cross_Y,cross_Z,vec_X,vec_Y,vec_Z,g)
          do pp=1_idef,spp(ii)%ppp,pchunk

             !$OMP SIMD
             do cc=1_idef,pchunk
                X_X(cc)=spp(ii)%vars%X(pp-1+cc,1)
                X_Y(cc)=spp(ii)%vars%X(pp-1+cc,2)
                X_Z(cc)=spp(ii)%vars%X(pp-1+cc,3)

                V_X(cc)=spp(ii)%vars%V(pp-1+cc,1)
                V_Y(cc)=spp(ii)%vars%V(pp-1+cc,2)
                V_Z(cc)=spp(ii)%vars%V(pp-1+cc,3)

                Y_R(cc)=0._rp
                Y_PHI(cc)=0._rp
                Y_Z(cc)=0._rp
                
                B_X(cc)=0._rp
                B_Y(cc)=0._rp
                B_Z(cc)=0._rp

                E_X(cc)=0._rp
                E_Y(cc)=0._rp
                E_Z(cc)=0._rp

                PSIp(cc)=100._rp
                
                flagCon(cc)=spp(ii)%vars%flagCon(pp-1+cc)
                flagCol(cc)=spp(ii)%vars%flagCol(pp-1+cc)
             end do
             !$OMP END SIMD

#ifdef FIO
             if (TRIM(params%field_model).eq.'M3D_C1'.or. &
                  (TRIM(params%field_model).eq.'NIMROD')) then
                !$OMP SIMD
                do cc=1_idef,pchunk
                   hint(cc)=spp(ii)%vars%hint(pp-1+cc)
                end do
                !$OMP END SIMD
             end if
#endif

             call cart_to_cyl_p(pchunk,X_X,X_Y,X_Z,Y_R,Y_PHI,Y_Z)             

             if (params%field_model(1:3).eq.'ANA') then
                call analytical_fields_p(params,pchunk,F, &
                     X_X,X_Y,X_Z,B_X,B_Y,B_Z,E_X,E_Y,E_Z,flagCon)
             else if (params%field_model(1:3).eq.'UNI') then
                call uniform_fields_p(pchunk,F, &
                    B_X,B_Y,B_Z,E_X,E_Y,E_Z)
#ifdef PSPLINE
             else if (F%axisymmetric_fields.and. &
                  (params%orbit_model(3:3).eq.'B')) then
                call interp_FOfields_p(pchunk,F,Y_R,Y_PHI,Y_Z,B_X,B_Y,B_Z, &
                     E_X,E_Y,E_Z,PSIp,flagCon)
             else if ((.not.F%axisymmetric_fields).and. &
                  (params%orbit_model(3:3).eq.'B')) then
                call interp_FO3Dfields_p(pchunk,F,Y_R,Y_PHI,Y_Z, &
                     B_X,B_Y,B_Z,E_X,E_Y,E_Z,flagCon)
             else if (params%orbit_model(3:5).eq.'psi') then
                call interp_FOfields1_p(pchunk,F,Y_R,Y_PHI,Y_Z,B_X,B_Y,B_Z, &
                     E_X,E_Y,E_Z,PSIp,flagCon)
             else if (params%field_model(10:13).eq.'MARS') then
                call interp_FOfields_mars_p(params,pchunk,F,Y_R,Y_PHI,Y_Z, &
                     B_X,B_Y,B_Z,PSIp,flagCon)
             else if (params%field_model(10:14).eq.'AORSA') then
                call interp_FOfields_aorsa_p(0._rp,params,pchunk,F, &
                     Y_R,Y_PHI,Y_Z,B_X,B_Y,B_Z,E_X,E_Y,E_Z,PSIp,flagCon)
#endif
#ifdef FIO
             else if (TRIM(params%field_model).eq.'M3D_C1'.or. &
                  TRIM(params%field_model).eq.'NIMROD') then
                call get_fio_FOmagnetic_fields_p(params,F,Y_R,Y_PHI,Y_Z, &
                     B_X,B_Y,B_Z,flagCon,hint)
                if (F%FIO_E .ge. 0) then
                   call get_fio_FOelectric_fields_p(params,F, &
                        Y_R,Y_PHI,Y_Z,E_X,E_Y,E_Z,flagCon,hint)
                end if
                if (F%FIO_A .ge. 0) then
                   call get_fio_vector_potential_p(params,F,Y_R,Y_PHI,Y_Z, &
                        PSIp,flagCon,hint)
                end if
#endif                
             end if

             !write(6,'("Y_R: ",E17.10)') Y_R*params%cpp%length
             !write(6,'("Y_PHI: ",E17.10)') Y_PHI
             !write(6,'("Y_Z: ",E17.10)') Y_Z*params%cpp%length
             !write(6,*) 'r_sam',sqrt((Y_R-spp(ii)%Ro)**2+Y_Z**2)*params%cpp%length
             
             !write(6,'("B_X: ",E17.10)') B_X*params%cpp%Bo
             !write(6,'("B_Y: ",E17.10)') B_Y*params%cpp%Bo
             !write(6,'("B_Z: ",E17.10)') B_Z*params%cpp%Bo
             
             !$OMP SIMD
             do cc=1_idef,pchunk
                spp(ii)%vars%B(pp-1+cc,1) = B_X(cc)
                spp(ii)%vars%B(pp-1+cc,2) = B_Y(cc)
                spp(ii)%vars%B(pp-1+cc,3) = B_Z(cc)

                spp(ii)%vars%E(pp-1+cc,1) = E_X(cc)
                spp(ii)%vars%E(pp-1+cc,2) = E_Y(cc)
                spp(ii)%vars%E(pp-1+cc,3) = E_Z(cc)

                spp(ii)%vars%PSI_P(pp-1+cc) = PSIp(cc)
             end do
             !$OMP END SIMD

#ifdef FIO
             if (TRIM(params%field_model).eq.'M3D_C1'.or. &
                  TRIM(params%field_model).eq.'NIMROD') then
                !$OMP SIMD
                do cc=1_idef,pchunk
                   spp(ii)%vars%hint(pp-1+cc) = hint(cc)
                end do
                !$OMP END SIMD
             end if
#endif

             !$OMP SIMD
             !          !$OMP& aligned(Bmag,B_X,B_Y,B_Z, &
             !          !$OMP& b_unit_X,b_unit_Y,b_unit_Z,v,V_X,V_Y,V_Z,vpar, &
             !          !$OMP& vperp,tmp,cross_X,cross_Y,cross_Z, &
             !          !$OMP& vec_X,vec_Y,vec_Z,E_X,E_Y,E_Z)
             do cc=1_idef,pchunk
                !Derived output data
                Bmag(cc) = SQRT(B_X(cc)*B_X(cc)+B_Y(cc)*B_Y(cc)+B_Z(cc)*B_Z(cc))

                ! Parallel unit vector
                b_unit_X(cc) = B_X(cc)/Bmag(cc)
                b_unit_Y(cc) = B_Y(cc)/Bmag(cc)
                b_unit_Z(cc) = B_Z(cc)/Bmag(cc)

                !write(6,*) 'X',X_X,X_Y,X_Z
                !write(6,*) 'b_unit',b_unit_X,b_unit_Y,b_unit_Z
                
                v(cc) = SQRT(V_X(cc)*V_X(cc)+V_Y(cc)*V_Y(cc)+V_Z(cc)*V_Z(cc))
                if (v(cc).GT.korc_zero) then
                   ! Parallel and perpendicular components of velocity
                   vpar(cc) = (V_X(cc)*b_unit_X(cc)+V_Y(cc)*b_unit_Y(cc)+ &
                        V_Z(cc)*b_unit_Z(cc))

                   vperp(cc) =  v(cc)**2 - vpar(cc)**2
                   if ( vperp(cc) .GE. korc_zero ) then
                      vperp(cc) = SQRT( vperp(cc) )
                   else
                      vperp(cc) = 0.0_rp
                   end if

                   !write(6,*) 'v,vpar,vperp',v(cc),vpar(cc),vperp(cc)
                   
                   ! Pitch angle
                   spp(ii)%vars%eta(pp-1+cc) = 180.0_rp* &
                        MODULO(ATAN2(vperp(cc),vpar(cc)),2.0_rp*C_PI)/C_PI

                   ! Magnetic moment
                   spp(ii)%vars%mu(pp-1+cc) = 0.5_rp*m_cache* &
                        g(cc)**2*vperp(cc)**2/Bmag(cc)
                   ! See Northrop's book (The adiabatic motion of charged
                   ! particles)

                   ! Radiated power
                   tmp(cc) = q_cache**4/(6.0_rp*C_PI*E0*m_cache**2)

                   cross_X(cc) = V_Y(cc)*B_Z(cc)-V_Z(cc)*B_Y(cc)
                   cross_Y(cc) = V_Z(cc)*B_X(cc)-V_X(cc)*B_Z(cc)
                   cross_Z(cc) = V_X(cc)*B_Y(cc)-V_Y(cc)*B_X(cc)

                   vec_X(cc) = E_X(cc) + cross_X(cc)
                   vec_Y(cc) = E_Y(cc) + cross_Y(cc)
                   vec_Z(cc) = E_Z(cc) + cross_Z(cc)

                   spp(ii)%vars%Prad(pp-1+cc) = tmp(cc)* &
                        ( E_X(cc)*E_X(cc)+E_Y(cc)*E_Y(cc)+E_Z(cc)*E_Z(cc) + &
                        cross_X(cc)*E_X(cc)+cross_Y(cc)*E_Y(cc)+ &
                        cross_Z(cc)*E_Z(cc) + g(cc)**2* &
                        ((E_X(cc)*V_X(cc)+E_Y(cc)*V_Y(cc)+E_Z(cc)*V_Z(cc))**2 &
                        - vec_X(cc)*vec_X(cc)-vec_Y(cc)*vec_Y(cc)- &
                        vec_Z(cc)*vec_Z(cc)) )

                   ! Input power due to electric field
                   spp(ii)%vars%Pin(pp-1+cc) = q_cache*(E_X(cc)*V_X(cc)+ &
                        E_Y(cc)*V_Y(cc)+E_Z(cc)*V_Z(cc))
                else
                   spp(ii)%vars%eta(pp-1+cc) = 0.0_rp
                   spp(ii)%vars%mu(pp-1+cc) = 0.0_rp
                   spp(ii)%vars%Prad(pp-1+cc) = 0.0_rp
                   spp(ii)%vars%Pin(pp-1+cc) = 0.0_rp
                end if

             end do
             !$OMP END SIMD
             
          end do
          !$OMP END PARALLEL DO

       end if !(if output)

       if(step.and.(.not.params%FokPlan)) then
          dt=0.5_rp*params%dt

          !$OMP PARALLEL DO FIRSTPRIVATE(dt) PRIVATE(pp,cc) &
          !$OMP& SHARED(ii,spp,params)
          do pp=1_idef,spp(ii)%ppp,pchunk

             !$OMP SIMD
             do cc=1_idef,pchunk
                spp(ii)%vars%X(pp-1+cc,1) = spp(ii)%vars%X(pp-1+cc,1) + &
                     dt*spp(ii)%vars%V(pp-1+cc,1)
                spp(ii)%vars%X(pp-1+cc,2) = spp(ii)%vars%X(pp-1+cc,2) + &
                     dt*spp(ii)%vars%V(pp-1+cc,2)
                spp(ii)%vars%X(pp-1+cc,3) = spp(ii)%vars%X(pp-1+cc,3) + &
                     dt*spp(ii)%vars%V(pp-1+cc,3)
             end do
             !$OMP END SIMD

          end do
          !$OMP END PARALLEL DO

       end if !(if step)

    end do ! over species

  end subroutine FO_init
FO_init Subroutine

Contents

Source Code

public subroutine FO_init(params, F, spp, output, step)

Arguments

Calls

Called by

Contents

Source Code

Source Code
