derivs_gca_rk Subroutine

private subroutine derivs_gca_rk(t_s, y, dydt)

Uses

  • proc~~derivs_gca_rk~~UsesGraph proc~derivs_gca_rk derivs_gca_rk module~mod_global_parameters mod_global_parameters proc~derivs_gca_rk->module~mod_global_parameters module~mod_basic_types mod_basic_types module~mod_global_parameters->module~mod_basic_types module~mod_connectivity mod_connectivity module~mod_global_parameters->module~mod_connectivity module~mod_constants mod_constants module~mod_global_parameters->module~mod_constants module~mod_physicaldata mod_physicaldata module~mod_global_parameters->module~mod_physicaldata module~mod_variables mod_variables module~mod_global_parameters->module~mod_variables mpi mpi module~mod_global_parameters->mpi module~mod_variables->module~mod_basic_types

Arguments

Type IntentOptional Attributes Name
double precision :: t_s
double precision :: y(ndir+2)
double precision :: dydt(ndir+2)

Calls

proc~~derivs_gca_rk~~CallsGraph proc~derivs_gca_rk derivs_gca_rk proc~cross cross proc~derivs_gca_rk->proc~cross proc~get_bfield get_bfield proc~derivs_gca_rk->proc~get_bfield proc~get_efield get_efield proc~derivs_gca_rk->proc~get_efield proc~get_vec get_vec proc~derivs_gca_rk->proc~get_vec proc~mpistop mpistop proc~derivs_gca_rk->proc~mpistop proc~cross->proc~mpistop proc~get_bfield->proc~get_vec proc~interpolate_var interpolate_var proc~get_bfield->proc~interpolate_var proc~get_efield->proc~get_bfield proc~get_efield->proc~get_vec proc~get_efield->proc~interpolate_var proc~get_vec->proc~interpolate_var mpi_abort mpi_abort proc~mpistop->mpi_abort proc~interpolate_var->proc~mpistop

Called by

proc~~derivs_gca_rk~~CalledByGraph proc~derivs_gca_rk derivs_gca_rk proc~gca_get_particle_dt gca_get_particle_dt proc~gca_get_particle_dt->proc~derivs_gca_rk proc~gca_integrate_particles gca_integrate_particles proc~gca_integrate_particles->proc~derivs_gca_rk proc~gca_integrate_particles->proc~gca_get_particle_dt

Source Code

  subroutine derivs_gca_rk(t_s,y,dydt)
    use mod_global_parameters

    double precision                :: t_s, y(ndir+2)
    double precision                :: dydt(ndir+2)

    double precision,dimension(ndir):: vE, b, e, x, bhat, bdotgradb,&
        vEdotgradb, gradkappaB, vfluid, current
    double precision,dimension(ndir):: bdotgradvE, vEdotgradvE, u, utmp1,&
        utmp2, utmp3
    double precision                :: upar, Mr, gamma, absb, q, m, epar,&
        kappa
    integer                         :: ic1,ic2,ic3

    ! Here the terms in the guiding centre equations of motion are interpolated for
    ! the RK integration. The interpolation is also done in the ghost cells such
    ! that the RK integration does not give an error

    q = particle(ipart_working)%self%q
    m = particle(ipart_working)%self%m

    x(1:ndir) = y(1:ndir)
    upar      = y(ndir+1) ! gamma v||
    Mr        = y(ndir+2)
    !gamma     = y(ndir+3)

    if (any(x .ne. x)) then
      write(*,*) "ERROR IN DERIVS_GCA_RK: NaNs IN X OR Y!"
      write(*,*) "x",x
      write(*,*) "y",y(ndir+1:ndir+2)
      call mpistop("ABORTING...")
    end if

    call get_bfield(igrid_working, x, t_s, b)
    call get_efield(igrid_working, x, t_s, e)
!    call get_vec(bp, igrid_working,x,t_s,b)
!    call get_vec(vp, igrid_working,x,t_s,vfluid)
!    call get_vec(jp, igrid_working,x,t_s,current)
!    e(1) = -vfluid(2)*b(3)+vfluid(3)*b(2) + particles_eta*current(1)
!    e(2) = vfluid(1)*b(3)-vfluid(3)*b(1) + particles_eta*current(2)
!    e(3) = -vfluid(1)*b(2)+vfluid(2)*b(1) + particles_eta*current(3)
    call get_vec(b_dot_grad_b, igrid_working,x,t_s,bdotgradb)
    call get_vec(vE_dot_grad_b, igrid_working,x,t_s,vEdotgradb)
    call get_vec(grad_kappa_B, igrid_working,x,t_s,gradkappaB)
    call get_vec(b_dot_grad_vE, igrid_working,x,t_s,bdotgradvE)
    call get_vec(vE_dot_grad_vE, igrid_working,x,t_s,vEdotgradvE)

    if (any(b .ne. b) .or. any(e .ne. e) .or. any(bdotgradb .ne. bdotgradb) &
       .or. any(vEdotgradb .ne. vEdotgradb) .or. any(gradkappaB .ne. &
       gradkappaB) .or. any(bdotgradvE .ne. bdotgradvE) .or. any(vEdotgradvE &
       .ne. vEdotgradvE)) then
      write(*,*) "ERROR IN DERIVS_GCA_RK: NaNs IN FIELD QUANTITIES!"
      write(*,*) "b",b
      write(*,*) "e",e
      write(*,*) "bdotgradb",bdotgradb
      write(*,*) "vEdotgradb",vEdotgradb
      write(*,*) "gradkappaB",gradkappaB
      write(*,*) "bdotgradvE",bdotgradvE
      write(*,*) "vEdotgradvE",vEdotgradvE
      call mpistop("ABORTING...")
    end if

    absb  = sqrt(sum(b(:)**2))
    if (absb .gt. 0.d0) then
      bhat(1:ndir) = b(1:ndir) / absb
    else
      bhat = 0.d0
    end if
    epar = sum(e(:)*bhat(:))

    call cross(e,bhat,vE)
    if (absb .gt. 0.d0) then
      vE(1:ndir) = vE(1:ndir) / absb
    else
      vE(1:ndir) = 0.d0
    end if

    if (relativistic) then
      kappa = 1.d0/sqrt(1.0d0 - sum(vE(:)**2)/c_norm**2)
      gamma = sqrt(1.0d0+upar**2/c_norm**2+2.0d0*Mr*absb/m/c_norm**2)*kappa
    else
      kappa = 1.d0
      gamma = 1.d0
    end if

    if (absb .gt. 0.d0) then
      utmp1(1:ndir) = bhat(1:ndir)/(absb/kappa**2)
    else
      utmp1 = 0.d0
    end if
    utmp2(1:ndir) = Mr/(gamma*q)*gradkappaB(1:ndir) + m/q* &
       (upar**2/gamma*bdotgradb(1:ndir) + upar*vEdotgradb(1:ndir) + &
       upar*bdotgradvE(1:ndir) + gamma*vEdotgradvE(1:ndir))
    if (relativistic) then
      utmp2(1:ndir) = utmp2(1:ndir) + upar*epar/(gamma)*vE(1:ndir)
    end if

    call cross(utmp1,utmp2,utmp3)
    u(1:ndir) = vE(1:ndir) + utmp3(1:ndir)

    ! done assembling the terms, now write rhs:
    dydt(1:ndir) = ( u(1:ndir) + upar/gamma * bhat(1:ndir) )
    dydt(ndir+1) = q/m*epar - Mr/(m*gamma) * sum(bhat(:)*gradkappaB(:)) + &
       sum(vE(:)*(upar*bdotgradb(:)+gamma*vEdotgradb(:)))
    dydt(ndir+2) = 0.0d0 ! magnetic moment is conserved

  end subroutine derivs_gca_rk