GCEoM1_p – KORC

private subroutine GCEoM1_p(pchunk, tt, P, F, params, RHS_R, RHS_PHI, RHS_Z, RHS_PLL, RHS_MU, B_R, B_PHI, B_Z, E_R, E_PHI, E_Z, curlb_R, curlb_PHI, curlb_Z, gradB_R, gradB_PHI, gradB_Z, V_PLL, V_MU, Y_R, Y_PHI, Y_Z, q_cache, m_cache, PSIp, ne, flag_cache)

Uses

- omp_lib
Help
Graph Key

Nodes of different colours represent the following:

Solid arrows point from a submodule to the (sub)module which it is descended from. Dashed arrows point from a module or program unit to modules which it uses.

Arguments

Type	Intent	Attributes		Name
integer,	intent(in)		::	pchunk
integer(kind=ip),	intent(in)		::	tt
type(PROFILES),	intent(in)		::	P
type(FIELDS),	intent(in)		::	F
type(KORC_PARAMS),	intent(inout)		::	params	Core KORC simulation parameters.
real(kind=rp),	intent(out),	DIMENSION(pchunk)	::	RHS_R
real(kind=rp),	intent(out),	DIMENSION(pchunk)	::	RHS_PHI
real(kind=rp),	intent(out),	DIMENSION(pchunk)	::	RHS_Z
real(kind=rp),	intent(out),	DIMENSION(pchunk)	::	RHS_PLL
real(kind=rp),	intent(out),	DIMENSION(pchunk)	::	RHS_MU
real(kind=rp),	intent(in),	DIMENSION(pchunk)	::	B_R
real(kind=rp),	intent(in),	DIMENSION(pchunk)	::	B_PHI
real(kind=rp),	intent(in),	DIMENSION(pchunk)	::	B_Z
real(kind=rp),	intent(in),	DIMENSION(pchunk)	::	E_R
real(kind=rp),	intent(in),	DIMENSION(pchunk)	::	E_PHI
real(kind=rp),	intent(in),	DIMENSION(pchunk)	::	E_Z
real(kind=rp),	intent(in),	DIMENSION(pchunk)	::	curlb_R
real(kind=rp),	intent(in),	DIMENSION(pchunk)	::	curlb_PHI
real(kind=rp),	intent(in),	DIMENSION(pchunk)	::	curlb_Z
real(kind=rp),	intent(in),	DIMENSION(pchunk)	::	gradB_R
real(kind=rp),	intent(in),	DIMENSION(pchunk)	::	gradB_PHI
real(kind=rp),	intent(in),	DIMENSION(pchunk)	::	gradB_Z
real(kind=rp),	intent(in),	DIMENSION(pchunk)	::	V_PLL
real(kind=rp),	intent(in),	DIMENSION(pchunk)	::	V_MU
real(kind=rp),	intent(in),	DIMENSION(pchunk)	::	Y_R
real(kind=rp),	intent(in),	DIMENSION(pchunk)	::	Y_PHI
real(kind=rp),	intent(in),	DIMENSION(pchunk)	::	Y_Z
real(kind=rp),	intent(in)		::	q_cache
real(kind=rp),	intent(in)		::	m_cache
real(kind=rp),	intent(in),	DIMENSION(pchunk)	::	PSIp
real(kind=rp),	intent(out),	DIMENSION(pchunk)	::	ne
integer(kind=is),	intent(out),	DIMENSION(pchunk)	::	flag_cache

Source Code

  subroutine GCEoM1_p(pchunk,tt,P,F,params,RHS_R,RHS_PHI,RHS_Z,RHS_PLL,RHS_MU, &
       B_R,B_PHI,B_Z,E_R,E_PHI,E_Z,curlb_R,curlb_PHI,curlb_Z, &
       gradB_R,gradB_PHI,gradB_Z,V_PLL,V_MU,Y_R,Y_PHI,Y_Z,q_cache,m_cache,PSIp,ne,flag_cache)

    USE omp_lib
    IMPLICIT NONE
    
    TYPE(KORC_PARAMS), INTENT(INOUT)                           :: params
    !! Core KORC simulation parameters.
    TYPE(FIELDS), INTENT(IN)      :: F
    TYPE(PROFILES), INTENT(IN)                                 :: P
    INTEGER, INTENT(IN) :: pchunk
    REAL(rp),DIMENSION(pchunk)  :: Bmag,bhat_R,bhat_PHI,bhat_Z,Bst_R,Bst_PHI
    REAL(rp),DIMENSION(pchunk)  :: BstdotE,BstdotgradB,EcrossB_R,EcrossB_PHI,bdotBst
    REAL(rp),DIMENSION(pchunk)  :: bcrossgradB_R,bcrossgradB_PHI,bcrossgradB_Z,gamgc
    REAL(rp),DIMENSION(pchunk)  :: EcrossB_Z,Bst_Z
    REAL(rp),DIMENSION(pchunk)  :: pm,xi,tau_R
    REAL(rp),DIMENSION(pchunk)  :: SR_PLL,SR_MU,BREM_PLL,BREM_MU,BREM_P
    REAL(rp),DIMENSION(pchunk),INTENT(in) :: gradB_R,gradB_PHI,gradB_Z,curlb_R
    REAL(rp),DIMENSION(pchunk),INTENT(in) :: curlb_Z,B_R,B_PHI,B_Z,E_R,E_PHI,E_Z
    REAL(rp),DIMENSION(pchunk),INTENT(OUT) :: RHS_R,RHS_PHI,RHS_Z
    REAL(rp),DIMENSION(pchunk),INTENT(OUT) :: RHS_PLL,RHS_MU
    REAL(rp),DIMENSION(pchunk),INTENT(IN) :: V_PLL,V_MU,Y_R,Y_PHI,Y_Z,curlb_PHI
    REAL(rp),DIMENSION(pchunk),INTENT(IN) :: PSIp
    REAL(rp),INTENT(in) :: q_cache,m_cache
    INTEGER  :: cc
    INTEGER(ip),INTENT(IN)  :: tt
    INTEGER(is),DIMENSION(pchunk),INTENT(OUT)  :: flag_cache
    REAL(rp)  :: time,re_cache,alpha_cache
    REAL(rp), DIMENSION(pchunk) 			:: Zeff,Te
    REAL(rp), DIMENSION(pchunk),INTENT(OUT) 		:: ne
    INTEGER             :: thread_num


    thread_num = OMP_GET_THREAD_NUM()

    !$OMP SIMD
    !    !$OMP& aligned(gradB_R,gradB_PHI,gradB_Z,curlb_R,curlb_Z, &
    !    !$OMP& B_R,B_PHI,B_Z,E_R,E_PHI,E_Z,RHS_R,RHS_PHI,RHS_Z,RHS_PLL,RHS_MU, &
    !    !$OMP& V_PLL,V_MU,Y_R,curlb_PHI,tau_R)
    do cc=1_idef,pchunk

       ne(cc)=-1._rp
       Te(cc)=-1._rp
       Zeff(cc)=-1._rp
       
       Bmag(cc) = SQRT(B_R(cc)*B_R(cc)+B_PHI(cc)*B_PHI(cc)+B_Z(cc)*B_Z(cc))
       
       bhat_R(cc) = B_R(cc)/Bmag(cc)
       bhat_PHI(cc) = B_PHI(cc)/Bmag(cc)
       bhat_Z(cc) = B_Z(cc)/Bmag(cc)

       Bst_R(cc)=q_cache*B_R(cc)+V_PLL(cc)*curlb_R(cc)
       Bst_PHI(cc)=q_cache*B_PHI(cc)+V_PLL(cc)*curlb_PHI(cc)
       Bst_Z(cc)=q_cache*B_Z(cc)+V_PLL(cc)*curlb_Z(cc)

      ! write(output_unit_write,*) 'bmag',Bmag(cc),'bhat',bhat_R(cc),bhat_PHI(cc),bhat_Z(cc),'Bst',Bst_R(cc),Bst_PHI(cc),Bst_Z(cc)
       
       bdotBst(cc)=bhat_R(cc)*Bst_R(cc)+bhat_PHI(cc)*Bst_PHI(cc)+ &
            bhat_Z(cc)*Bst_Z(cc)
       BstdotE(cc)=Bst_R(cc)*E_R(cc)+Bst_PHI(cc)*E_PHI(cc)+Bst_Z(cc)*E_Z(cc)   
       BstdotgradB(cc)=Bst_R(cc)*gradB_R(cc)+Bst_PHI(cc)*gradB_PHI(cc)+ &
            Bst_Z(cc)*gradB_Z(cc)

       !write(output_unit_write,*) 'bdotBst',bdotBst(cc),BstdotE(cc),BstdotgradB(cc)
       
       Ecrossb_R(cc)=E_PHI(cc)*bhat_Z(cc)-E_Z(cc)*bhat_PHI(cc)
       Ecrossb_PHI(cc)=E_Z(cc)*bhat_R(cc)-E_R(cc)*bhat_Z(cc)
       Ecrossb_Z(cc)=E_R(cc)*bhat_PHI(cc)-E_PHI(cc)*bhat_R(cc)

       !write(output_unit_write,*) 'Ecrossb',Ecrossb_R(cc),Ecrossb_PHI(cc),Ecrossb_Z(cc)
       
       bcrossgradB_R(cc)=bhat_PHI(cc)*gradB_Z(cc)-bhat_Z(cc)*gradB_PHI(cc)
       bcrossgradB_PHI(cc)=bhat_Z(cc)*gradB_R(cc)-bhat_R(cc)*gradB_Z(cc)
       bcrossgradB_Z(cc)=bhat_R(cc)*gradB_PHI(cc)-bhat_PHI(cc)*gradB_R(cc)

      ! write(output_unit_write,*) 'bcrossgradB',bcrossgradB_R(cc),bcrossgradB_PHI(cc),bcrossgradB_Z(cc)
       
       gamgc(cc)=sqrt(1+V_PLL(cc)*V_PLL(cc)+2*V_MU(cc)*Bmag(cc))
       
       pm(cc)=sqrt(gamgc(cc)**2-1)
       xi(cc)=V_PLL(cc)/pm(cc)

       RHS_R(cc)=(q_cache*Ecrossb_R(cc)+(m_cache*V_MU(cc)* &
            bcrossgradB_R(cc)+V_PLL(cc)*Bst_R(cc))/(m_cache*gamgc(cc)))/ &
            bdotBst(cc)
       RHS_PHI(cc)=(q_cache*Ecrossb_PHI(cc)+(m_cache*V_MU(cc)* &
            bcrossgradB_PHI(cc)+V_PLL(cc)*Bst_PHI(cc))/(m_cache*gamgc(cc)))/ &
            (Y_R(cc)*bdotBst(cc))
       RHS_Z(cc)=(q_cache*Ecrossb_Z(cc)+(m_cache*V_MU(cc)* &
            bcrossgradB_Z(cc)+V_PLL(cc)*Bst_Z(cc))/(m_cache*gamgc(cc)))/ &
            bdotBst(cc)
       RHS_PLL(cc)=(q_cache*BstdotE(cc)-V_MU(cc)*BstdotgradB(cc)/gamgc(cc))/ &
            bdotBst(cc)
       RHS_MU(cc)=0._rp

    end do
    !$OMP END SIMD

    !write(output_unit_write,*) 'bmag',Bmag(1),'bhat',bhat_R(1),bhat_PHI(1),bhat_Z(1),'Bst',Bst_R(1),Bst_PHI(1),Bst_Z(1)
    !write(output_unit_write,*) 'bdotBst',bdotBst(1),BstdotE(1),BstdotgradB(1)
    !write(output_unit_write,*) 'Ecrossb',Ecrossb_R(1),Ecrossb_PHI(1),Ecrossb_Z(1)
    !write(output_unit_write,*) 'bcrossgradB',bcrossgradB_R(1),bcrossgradB_PHI(1),bcrossgradB_Z(1)

    if (params%radiation.and.(params%GC_rad_model.eq.'SDE')) then

       !       write(output_unit_write,*) 'RHS_PLL',RHS_PLL(1)

       re_cache=C_RE/params%cpp%length
       alpha_cache=C_a

       if (.not.params%LargeCollisions) then
          time=params%init_time+(params%it-1+tt)*params%dt
       else
          time=params%init_time+params%it*params%dt+ &
               tt*params%coll_per_dump_dt
       end if

       if (params%profile_model(1:10).eq.'ANALYTICAL') then
          call analytical_profiles_p(pchunk,time,params,Y_R,Y_Z,P,F,ne,Te,Zeff,PSIp)
       else
#ifdef PSPLINE
          call interp_FOcollision_p(pchunk,Y_R,Y_PHI,Y_Z,ne,Te,Zeff,flag_cache)
#endif
       endif

       !write(6,*) 'Y: ',Y_R,Y_PHI,Y_Z
       !write(6,*) 'Profs: ',ne,Te,Zeff
       
       
       !$OMP SIMD
       !       !$OMP& aligned(tau_R,Bmag,RHS_PLL,V_PLL,xi,gamgc,RHS_MU,V_MU)
       do cc=1_idef,pchunk

          tau_R(cc)=6*C_PI*E0/(Bmag(cc)*Bmag(cc))

          SR_PLL(cc)=V_PLL(cc)*(1._rp-xi(cc)*xi(cc))/tau_R(cc)* &
               (1._rp/gamgc(cc)-gamgc(cc))
          SR_MU(cc)=-2._rp*V_MU(cc)/tau_R(cc)* &
               (gamgc(cc)*(1-xi(cc)*xi(cc))+xi(cc)*xi(cc)/gamgc(cc))

          !Normalizations done here
          BREM_P(cc)=-4._rp*re_cache**2*ne(cc)* &
               Zeff(cc)*(Zeff(cc)+1._rp)*alpha_cache* &
               (gamgc(cc)-1._rp)*(log(2._rp*gamgc(cc))-1._rp/3._rp)
          BREM_PLL(cc)=xi(cc)*BREM_P(cc)
          BREM_MU(cc)=(1._rp-xi(cc)*xi(cc))*V_PLL(cc)/ &
               (Bmag(cc)*xi(cc))*BREM_P(cc)

          RHS_PLL(cc)=RHS_PLL(cc)+SR_PLL(cc)+BREM_PLL(cc)
          RHS_MU(cc)=SR_MU(cc)+BREM_MU(cc)

       end do
       !$OMP END SIMD

    end if   

    
#if DBG_CHECK    
    !$OMP SIMD
    do cc=1_idef,pchunk 
       if(isnan(gamgc(cc)).and.flag_cache(cc)==1) then           
          write(6,*) 'gamgc is NaN'
          write(6,*) 'V_PLL',V_PLL(cc)
          write(6,*) 'V_MU',V_MU(cc)
          write(6,*) 'Bmag',Bmag(cc)
          
          stop 'gamgc is a NaN'
       endif
       if(isnan(RHS_R(cc)).and.flag_cache(cc)==1) then  

          write(6,*) thread_num,'flags',flag_cache(cc)
          write(6,*) thread_num,'Y',Y_R(cc)*params%cpp%length,Y_PHI(cc),Y_Z(cc)*params%cpp%length
          write(6,*) thread_num,'B',B_R(cc),B_PHI(cc),B_Z(cc)
          write(6,*) thread_num,'E',E_R(cc),E_PHI(cc),E_Z(cc)
          write(6,*) thread_num,'gradB',gradB_R(cc),gradB_PHI(cc),gradB_Z(cc)
          write(6,*) thread_num,'curlb',curlb_R(cc),curlb_PHI(cc),curlb_Z(cc)
          write(6,*) thread_num,'V',V_PLL(cc),V_MU(cc)
          write(6,*) thread_num,'Exb',Ecrossb_R(cc)
          write(6,*) thread_num,'bxgradB',bcrossgradB_R(cc)
          write(6,*) thread_num,'Bst',Bst_R(cc)
          write(6,*) thread_num,'bdotBst',bdotBst(cc)
          write(6,*) thread_num,'gamma',gamgc(cc)
          
          stop 'RHS_R1 is a NaN'
       endif
       if(isnan(RHS_PHI(cc)).and.flag_cache(cc)==1) then
          stop 'RHS_PHI1 is a NaN'
       end if
       if(isnan(RHS_Z(cc)).and. flag_cache(cc)==1) then
          stop 'RHS_Z1 is a NaN'
       endif
       if(isnan(RHS_PLL(cc)).and. flag_cache(cc)==1) then

          write(6,*) params%mpi_params%rank,thread_num,'flags',flag_cache(cc)
          write(6,*) params%mpi_params%rank,thread_num,'Y',Y_R(cc)*params%cpp%length,Y_PHI(cc),Y_Z(cc)*params%cpp%length
          write(6,*) params%mpi_params%rank,thread_num,'V',V_PLL(cc),V_MU(cc)
          write(6,*) params%mpi_params%rank,thread_num,'gamma',gamgc(cc)
          write(6,*) params%mpi_params%rank,thread_num,'xi',xi(cc)
          stop 'RHS_PLL1 is a NaN'
       endif
       if(isnan(RHS_MU(cc)).and. flag_cache(cc)==1) then
          stop 'RHS_MU1 is a NaN'
       endif
    end do
    !$OMP END SIMD

    if (params%radiation.and.(params%GC_rad_model.eq.'SDE')) then
       !$OMP SIMD
       do cc=1_idef,pchunk
          if(isnan(ne(cc)).and. flag_cache(cc)==1) then
             stop 'ne is a NaN'
          endif
          if(isnan(Zeff(cc)).and. flag_cache(cc)==1) then
             stop 'Zeff is a NaN'
          endif
          if(isnan(BREM_P(cc)).and. flag_cache(cc)==1) then
             stop 'BREM_P is a NaN'
          endif
          if(isnan(BREM_PLL(cc)).and. flag_cache(cc)==1) then
             stop 'BREM_PLL is a NaN'
          endif
          if(isnan(BREM_MU(cc)).and. flag_cache(cc)==1) then
             stop 'BREM_MU is a NaN'
          endif
          if(isnan(SR_PLL(cc)).and. flag_cache(cc)==1) then
             stop 'SR_PLL is a NaN'
          endif
          if(isnan(SR_MU(cc)).and. flag_cache(cc)==1) then
             stop 'SR_MU is a NaN'
          endif
       end do
       !$OMP END SIMD
    end if
#endif
    
    !    write(output_unit_write,*) 'RHS_R: ',RHS_R(1)
    !    write(output_unit_write,*) 'RHS_PHI: ',RHS_PHI(1)
    !    write(output_unit_write,*) 'RHS_Z: ',RHS_Z(1)
    !    write(output_unit_write,*) 'RHS_PLL: ',RHS_PLL(1)
    !    write(output_unit_write,*) 'RHS_MU: ',RHS_MU(1)

  end subroutine GCEoM1_p

GCEoM1_p Subroutine

Contents

Source Code

private subroutine GCEoM1_p(pchunk, tt, P, F, params, RHS_R, RHS_PHI, RHS_Z, RHS_PLL, RHS_MU, B_R, B_PHI, B_Z, E_R, E_PHI, E_Z, curlb_R, curlb_PHI, curlb_Z, gradB_R, gradB_PHI, gradB_Z, V_PLL, V_MU, Y_R, Y_PHI, Y_Z, q_cache, m_cache, PSIp, ne, flag_cache)

Uses

Arguments

Calls

Called by

Contents

Source Code

Source Code