derivs_gca Subroutine

private subroutine derivs_gca(t_s, y, dydt)

Uses

  • proc~~derivs_gca~~UsesGraph proc~derivs_gca derivs_gca module~mod_global_parameters mod_global_parameters proc~derivs_gca->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~~CallsGraph proc~derivs_gca derivs_gca proc~cross cross proc~derivs_gca->proc~cross proc~get_bfield get_bfield proc~derivs_gca->proc~get_bfield proc~get_efield get_efield proc~derivs_gca->proc~get_efield proc~get_vec get_vec proc~derivs_gca->proc~get_vec proc~mpistop mpistop proc~derivs_gca->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~~CalledByGraph proc~derivs_gca derivs_gca proc~gca_get_particle_dt gca_get_particle_dt proc~gca_get_particle_dt->proc~derivs_gca proc~gca_integrate_particles gca_integrate_particles proc~gca_integrate_particles->proc~derivs_gca proc~gca_integrate_particles->proc~gca_get_particle_dt

Source Code

  subroutine derivs_gca(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

    ! Here the normal interpolation is done for the terms in the GCA equations of motion

    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
      call mpistop("ERROR IN DERIVS_GCA: NaNs IN X! 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)

    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) vE(1:ndir)   = vE(1:ndir) / absb

    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