advance_FOeqn_vars

public subroutine advance_FOeqn_vars(tt, a, q_cache, m_cache, params, X_X, X_Y, X_Z, V_X, V_Y, V_Z, B_X, B_Y, B_Z, E_X, E_Y, E_Z, P, F, g, flagCon, flagCol, PSIp)

Calls radiation_force_p in korc_ppusher.

Arguments

Type	Intent	Attributes		Name
integer(kind=ip),	intent(in)		::	tt	Time step used in the leapfrog step ( $\Delta t$ ).
real(kind=rp),	intent(in)		::	a	This variable is used to simplify notation in the code, and is given by $a=q\Delta t/m$ ,
real(kind=rp),	intent(in)		::	q_cache	Time step used in the leapfrog step ( $\Delta t$ ).
real(kind=rp),	intent(in)		::	m_cache	Time step used in the leapfrog step ( $\Delta t$ ).
type(KORC_PARAMS),	intent(in)		::	params	Core KORC simulation parameters.
real(kind=rp),	intent(inout),	DIMENSION(params%pchunk)	::	X_X
real(kind=rp),	intent(inout),	DIMENSION(params%pchunk)	::	X_Y
real(kind=rp),	intent(inout),	DIMENSION(params%pchunk)	::	X_Z
real(kind=rp),	intent(inout),	DIMENSION(params%pchunk)	::	V_X
real(kind=rp),	intent(inout),	DIMENSION(params%pchunk)	::	V_Y
real(kind=rp),	intent(inout),	DIMENSION(params%pchunk)	::	V_Z
real(kind=rp),	intent(in),	DIMENSION(params%pchunk)	::	B_X
real(kind=rp),	intent(in),	DIMENSION(params%pchunk)	::	B_Y
real(kind=rp),	intent(in),	DIMENSION(params%pchunk)	::	B_Z
real(kind=rp),	intent(in),	DIMENSION(params%pchunk)	::	E_X
real(kind=rp),	intent(in),	DIMENSION(params%pchunk)	::	E_Y
real(kind=rp),	intent(in),	DIMENSION(params%pchunk)	::	E_Z
type(PROFILES),	intent(in)		::	P
type(FIELDS),	intent(in)		::	F
real(kind=rp),	intent(inout),	DIMENSION(params%pchunk)	::	g
integer(kind=is),	intent(inout),	DIMENSION(params%pchunk)	::	flagCon
integer(kind=is),	intent(inout),	DIMENSION(params%pchunk)	::	flagCol
real(kind=rp),	intent(in),	DIMENSION(params%pchunk)	::	PSIp

Source Code

  subroutine advance_FOeqn_vars(tt,a,q_cache,m_cache,params,X_X,X_Y,X_Z, &
       V_X,V_Y,V_Z,B_X,B_Y,B_Z,E_X,E_Y,E_Z,P,F,g,flagCon,flagCol,PSIp)
    TYPE(PROFILES), INTENT(IN)                                 :: P
    TYPE(FIELDS), INTENT(IN)      :: F
    TYPE(KORC_PARAMS), INTENT(IN)                              :: params
    !! Core KORC simulation parameters.

    INTEGER(ip), INTENT(IN)                                       :: tt
    !! Time step used in the leapfrog step (\(\Delta t\)).
    REAL(rp)                                      :: dt
    !! Time step used in the leapfrog step (\(\Delta t\)).
    REAL(rp), INTENT(IN)                       :: m_cache,q_cache
    !! Time step used in the leapfrog step (\(\Delta t\)).

    REAL(rp),DIMENSION(params%pchunk)                                  :: Bmag



    REAL(rp),INTENT(in)                                       :: a
    !! This variable is used to simplify notation in the code, and
    !! is given by \(a=q\Delta t/m\),
    REAL(rp),DIMENSION(params%pchunk)                                    :: sigma
    !! This variable is \(\sigma = \gamma'^2 - \tau^2\) in the above equations.
    REAL(rp),DIMENSION(params%pchunk)                               :: us
    !! This variable is \(u^{*} = p^{*}/m\) where \( p^{*} =
    !! \mathbf{p}'\cdot \mathbf{\tau}/mc\).
    !! Variable 'u^*' in Vay, J.-L. PoP (200params%pchunk).
    REAL(rp),DIMENSION(params%pchunk),INTENT(INOUT)                 :: g
    REAL(rp),DIMENSION(params%pchunk) :: gp,g0
    !! Relativistic factor \(\gamma\).
    REAL(rp),DIMENSION(params%pchunk)                                 :: s
    !! This variable is \(s = 1/(1+t^2)\) in the equations above.
    !! Variable 's' in Vay, J.-L. PoP (2008).
    REAL(rp),DIMENSION(params%pchunk)                            :: U_hs_X,U_hs_Y,U_hs_Z
    !! Is \(\mathbf{u}=\mathbf{p}/m\) at half-time step (\(i+1/2\)) in
    !! the absence of radiation losses or collisions. \(\mathbf{u}^{i+1/2} =
    !! \mathbf{u}^i + \frac{q\Delta t}{2m}\left( \mathbf{E}^{i+1/2} +
    !! \mathbf{v}^i\times \mathbf{B}^{i+1/2} \right)\).
    REAL(rp),DIMENSION(params%pchunk)                           :: tau_X,tau_Y,tau_Z
    !! This variable is \(\mathbf{\tau} = (q\Delta t/2)\mathbf{B}^{i+1/2}\).
    REAL(rp),DIMENSION(params%pchunk)                            :: up_X,up_Y,up_Z
    !! This variable is \(\mathbf{u}'= \mathbf{p}'/m\), where \(\mathbf{p}'
    !! = \mathbf{p}^i + q\Delta t \left( \mathbf{E}^{i+1/2} +
    !! \frac{\mathbf{v}^i}{2} \times \mathbf{B}^{i+1/2} \right)\).
    REAL(rp),DIMENSION(params%pchunk)                                     :: t_X,t_Y,t_Z
    !! This variable is \(\mathbf{t} = {\mathbf \tau}/\gamma^{i+1}\).
    REAL(rp),DIMENSION(params%pchunk),INTENT(INOUT)                     :: X_X,X_Y,X_Z
    REAL(rp),DIMENSION(params%pchunk),INTENT(INOUT)                      :: V_X,V_Y,V_Z
    REAL(rp),DIMENSION(params%pchunk),INTENT(IN)                      :: B_X,B_Y,B_Z
    REAL(rp),DIMENSION(params%pchunk),INTENT(IN)          :: E_X,E_Y,E_Z,PSIp
    REAL(rp),DIMENSION(params%pchunk)                     :: U_L_X,U_L_Y,U_L_Z
    REAL(rp),DIMENSION(params%pchunk)                     :: U_X,U_Y,U_Z
    REAL(rp),DIMENSION(params%pchunk)                     :: U_RC_X,U_RC_Y,U_RC_Z
    REAL(rp),DIMENSION(params%pchunk)                     :: U_os_X,U_os_Y,U_os_Z
    !! This variable is \(\mathbf{u}^{i+1}= \mathbf{p}^{i+1}/m\).
    REAL(rp),DIMENSION(params%pchunk)                          :: cross_X,cross_Y,cross_Z

    REAL(rp), DIMENSION(params%pchunk)                       :: Frad_X,Frad_Y,Frad_Z
    !! Synchrotron radiation reaction force of each particle.

    REAL(rp),DIMENSION(params%pchunk) :: ne,Te,Zeff,Y_R,Y_PHI,Y_Z

    INTEGER                                      :: cc,pchunk
    !! Chunk iterator.

    INTEGER(is),DIMENSION(params%pchunk),intent(inout)             :: flagCon,flagCol

    dt=params%dt

    pchunk=params%pchunk

    !write(6,*) 'X',X_X,X_Y,X_Z
    !write(6,*) 'V',V_X,V_Y,V_Z
    !write(6,*) 'B',B_X,B_Y,B_Z
    
    !$OMP SIMD
    !    !$OMP& aligned(g0,g,U_X,U_Y,U_Z,V_X,V_Y,V_Z,Bmag,B_X,B_Y,B_Z, &
    !    !$OMP& U_L_X,U_L_Y,U_L_Z,U_RC_X,U_RC_Y,U_RC_Z, &
    !    !$OMP& cross_X,cross_Y,cross_Z,U_hs_X,U_hs_Y,U_hs_Z,E_X,E_Y,E_Z, &
    !    !$OMP& tau_X,tau_Y,tau_Z,up_X,up_Y,up_Z,gp,sigma,us,t_X,t_Y,t_Z,s, &
    !    !$OMP& U_os_X,U_os_Y,U_os_Z,Frad_X,Frad_Y,Frad_Z)
    do cc=1_idef,pchunk

       g0(cc)=g(cc)

       U_X(cc) = g(cc)*V_X(cc)
       U_Y(cc) = g(cc)*V_Y(cc)
       U_Z(cc) = g(cc)*V_Z(cc)


       ! Magnitude of magnetic field
       Bmag(cc) = SQRT(B_X(cc)*B_X(cc)+B_Y(cc)*B_Y(cc)+B_Z(cc)*B_Z(cc))

       U_L_X(cc)=U_X(cc)
       U_L_Y(cc)=U_Y(cc)
       U_L_Z(cc)=U_Z(cc)

       U_RC_X(cc)=U_X(cc)
       U_RC_Y(cc)=U_Y(cc)
       U_RC_Z(cc)=U_Z(cc)

       ! LEAP-FROG SCHEME FOR LORENTZ FORCE !

       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)

       !write(6,*) 'vcrossB',cross_X,cross_Y,cross_Z

       U_hs_X(cc) = U_L_X(cc) + 0.5_rp*a*(E_X(cc) +cross_X(cc))
       U_hs_Y(cc) = U_L_Y(cc) + 0.5_rp*a*(E_Y(cc) +cross_Y(cc))
       U_hs_Z(cc) = U_L_Z(cc) + 0.5_rp*a*(E_Z(cc) +cross_Z(cc))


       !write(6,*) 'half step',0.5_rp*a*cross_X(cc),0.5_rp*a*cross_Y(cc),0.5_rp*a*cross_Z(cc)

       tau_X(cc) = 0.5_rp*a*B_X(cc)
       tau_Y(cc) = 0.5_rp*a*B_Y(cc)
       tau_Z(cc) = 0.5_rp*a*B_Z(cc)



       up_X(cc) = U_hs_X(cc) + 0.5_rp*a*E_X(cc)
       up_Y(cc) = U_hs_Y(cc) + 0.5_rp*a*E_Y(cc)
       up_Z(cc) = U_hs_Z(cc) + 0.5_rp*a*E_Z(cc)

       gp(cc) = SQRT( 1.0_rp + up_X(cc)*up_X(cc)+up_Y(cc)*up_Y(cc)+ &
            up_Z(cc)*up_Z(cc) )

       sigma(cc) = gp(cc)*gp(cc) - (tau_X(cc)*tau_X(cc)+ &
            tau_Y(cc)*tau_Y(cc)+tau_Z(cc)*tau_Z(cc))

       us(cc) = up_X(cc)*tau_X(cc)+up_Y(cc)*tau_Y(cc)+ &
            up_Z(cc)*tau_Z(cc)

       ! variable 'u^*' in Vay, J.-L. PoP (2008)
       g(cc) = SQRT( 0.5_rp*(sigma(cc) + SQRT(sigma(cc)*sigma(cc) + &
            4.0_rp*(tau_X(cc)*tau_X(cc)+tau_Y(cc)*tau_Y(cc)+ &
            tau_Z(cc)*tau_Z(cc) + us(cc)*us(cc)))) )

       t_X(cc) = tau_X(cc)/g(cc)
       t_Y(cc) = tau_Y(cc)/g(cc)
       t_Z(cc) = tau_Z(cc)/g(cc)


       s(cc) = 1.0_rp/(1.0_rp + t_X(cc)*t_X(cc)+t_Y(cc)*t_Y(cc)+ &
            t_Z(cc)*t_Z(cc))
       ! variable 's' in Vay, J.-L. PoP (2008)

       cross_X(cc)=up_Y(cc)*t_Z(cc)-up_Z(cc)*t_Y(cc)
       cross_Y(cc)=up_Z(cc)*t_X(cc)-up_X(cc)*t_Z(cc)
       cross_Z(cc)=up_X(cc)*t_Y(cc)-up_Y(cc)*t_X(cc)

       U_L_X(cc) = s(cc)*(up_X(cc) + (up_X(cc)*t_X(cc)+ &
            up_Y(cc)*t_Y(cc)+up_Z(cc)*t_Z(cc))*t_X(cc) + cross_X(cc))
       U_L_Y(cc) = s(cc)*(up_Y(cc) + (up_X(cc)*t_X(cc)+ &
            up_Y(cc)*t_Y(cc)+up_Z(cc)*t_Z(cc))*t_Y(cc) + cross_Y(cc))
       U_L_Z(cc) = s(cc)*(up_Z(cc) + (up_X(cc)*t_X(cc)+ &
            up_Y(cc)*t_Y(cc)+up_Z(cc)*t_Z(cc))*t_Z(cc) + cross_Z(cc))
       ! LEAP-FROG SCHEME FOR LORENTZ FORCE !      

       U_os_X(cc) = 0.5_rp*(U_L_X(cc) + U_X(cc))
       U_os_Y(cc) = 0.5_rp*(U_L_Y(cc) + U_Y(cc))
       U_os_Z(cc) = 0.5_rp*(U_L_Z(cc) + U_Z(cc))
       ! Splitting operator for including radiation

       if (params%radiation) then
          !! Calls [[radiation_force_p]] in [[korc_ppusher]].
          call radiation_force_p(pchunk,q_cache,m_cache,U_os_X,U_os_Y,U_os_Z, &
               E_X,E_Y,E_Z,B_Z,B_Y,B_Z,Frad_X,Frad_Y,Frad_Z)
          U_RC_X(cc) = U_RC_X(cc) + a*Frad_X(cc)/q_cache
          U_RC_Y(cc) = U_RC_Y(cc) + a*Frad_Y(cc)/q_cache
          U_RC_Z(cc) = U_RC_Z(cc) + a*Frad_Z(cc)/q_cache
       end if
       ! Splitting operator for including radiation

       U_X(cc) = U_L_X(cc) + U_RC_X(cc) - U_X(cc)
       U_Y(cc) = U_L_Y(cc) + U_RC_Y(cc) - U_Y(cc)
       U_Z(cc) = U_L_Z(cc) + U_RC_Z(cc) - U_Z(cc)

    end do
    !$OMP END SIMD


    if (params%collisions) then

       call include_CoulombCollisions_FO_p(tt,params,X_X,X_Y,X_Z, &
            U_X,U_Y,U_Z,B_X,B_Y,B_Z,m_cache,P,F,flagCon,flagCol,PSIp)

    end if

    if (params%radiation.or.params%collisions) then

       !$OMP SIMD
       !       !$OMP& aligned(g,U_X,U_Y,U_Z)
       do cc=1_idef,pchunk
          g(cc)=sqrt(1._rp+U_X(cc)*U_X(cc)+U_Y(cc)*U_Y(cc)+U_Z(cc)*U_Z(cc))
       end do
       !$OMP END SIMD

    end if

    !$OMP SIMD
    !    !$OMP& aligned(g,g0,V_X,V_Y,V_Z,U_X,U_Y,U_Z,X_X,X_Y,X_Z,flagCon,flagCol)
    do cc=1_idef,pchunk

       if ((flagCon(cc).eq.0_is).or.(flagCol(cc).eq.0_is)) then
          g(cc)=g0(cc)
       else
          V_X(cc) = U_X(cc)/g(cc)
          V_Y(cc) = U_Y(cc)/g(cc)
          V_Z(cc) = U_Z(cc)/g(cc)
       end if

       X_X(cc) = X_X(cc) + dt*V_X(cc)*REAL(flagCon(cc))*REAL(flagCol(cc))
       X_Y(cc) = X_Y(cc) + dt*V_Y(cc)*REAL(flagCon(cc))*REAL(flagCol(cc))
       X_Z(cc) = X_Z(cc) + dt*V_Z(cc)*REAL(flagCon(cc))*REAL(flagCol(cc))
    end do
    !$OMP END SIMD

  end subroutine advance_FOeqn_vars

advance_FOeqn_vars Subroutine

Contents

Source Code

public subroutine advance_FOeqn_vars(tt, a, q_cache, m_cache, params, X_X, X_Y, X_Z, V_X, V_Y, V_Z, B_X, B_Y, B_Z, E_X, E_Y, E_Z, P, F, g, flagCon, flagCol, PSIp)

Arguments

Calls

Called by

Contents

Source Code

Source Code