gca_create_particles Subroutine

public subroutine gca_create_particles()

Uses

  • proc~~gca_create_particles~~UsesGraph proc~gca_create_particles gca_create_particles module~mod_global_parameters mod_global_parameters proc~gca_create_particles->module~mod_global_parameters module~mod_usr_methods mod_usr_methods proc~gca_create_particles->module~mod_usr_methods 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

None

Calls

proc~~gca_create_particles~~CallsGraph proc~gca_create_particles gca_create_particles mpi_allreduce mpi_allreduce proc~gca_create_particles->mpi_allreduce mpi_bcast mpi_bcast proc~gca_create_particles->mpi_bcast proc~find_particle_ipe find_particle_ipe proc~gca_create_particles->proc~find_particle_ipe proc~gca_update_payload gca_update_payload proc~gca_create_particles->proc~gca_update_payload proc~get_lfac_from_velocity get_lfac_from_velocity proc~gca_create_particles->proc~get_lfac_from_velocity proc~get_vec get_vec proc~gca_create_particles->proc~get_vec proc~push_particle_into_particles_on_mype push_particle_into_particles_on_mype proc~gca_create_particles->proc~push_particle_into_particles_on_mype proc~unif_01 rng_t%unif_01 proc~gca_create_particles->proc~unif_01 proc~get_igslice get_igslice proc~find_particle_ipe->proc~get_igslice proc~particle_in_domain particle_in_domain proc~find_particle_ipe->proc~particle_in_domain proc~gca_update_payload->proc~get_vec proc~cross cross proc~gca_update_payload->proc~cross proc~get_bfield get_bfield proc~gca_update_payload->proc~get_bfield proc~get_efield get_efield proc~gca_update_payload->proc~get_efield proc~interpolate_var interpolate_var proc~get_vec->proc~interpolate_var proc~next rng_t%next proc~unif_01->proc~next proc~mpistop mpistop proc~cross->proc~mpistop proc~get_bfield->proc~get_vec proc~get_bfield->proc~interpolate_var proc~get_efield->proc~get_vec proc~get_efield->proc~get_bfield proc~get_efield->proc~interpolate_var proc~get_igslice->proc~mpistop proc~interpolate_var->proc~mpistop proc~rotl rotl proc~next->proc~rotl mpi_abort mpi_abort proc~mpistop->mpi_abort

Called by

proc~~gca_create_particles~~CalledByGraph proc~gca_create_particles gca_create_particles proc~particles_create particles_create proc~particles_create->proc~gca_create_particles proc~main main proc~main->proc~particles_create program~agile agile program~agile->proc~main

Source Code

  subroutine gca_create_particles()
    ! initialise the particles
    use mod_global_parameters
    use mod_usr_methods, only: usr_create_particles, usr_update_payload,&
        usr_check_particle

    double precision :: b(ndir), u(ndir), magmom
    double precision :: bnorm, lfac, vnorm, vperp, vpar
    integer          :: igrid, ipe_particle
    integer          :: n, idir, nparticles_local
    double precision :: x(3, num_particles)
    double precision :: v(3, num_particles)
    double precision :: q(num_particles)
    double precision :: m(num_particles)
    double precision :: rrd(num_particles,ndir)
    double precision :: defpayload(ndefpayload)
    double precision :: usrpayload(nusrpayload)
    logical          :: follow(num_particles), check

    if (mype==0) then
      if (.not. associated(usr_create_particles)) then
        ! Randomly distributed
        do idir=1,ndir
          do n = 1, num_particles
            rrd(n,idir) = rng%unif_01()
          end do
        end do
        do n=1, num_particles
          x(1,n) = xprobmin1 + rrd(n+1,1) * (xprobmax1 - xprobmin1)
          x(2,n) = xprobmin2 + rrd(n+1,2) * (xprobmax2 - xprobmin2)
          x(3,n) = xprobmin3 + rrd(n+1,3) * (xprobmax3 - xprobmin3)
        end do
      else
        call usr_create_particles(num_particles, x, v, q, m, follow)
      end if
    end if

    call MPI_BCAST(x,3*num_particles,MPI_DOUBLE_PRECISION,0,icomm,ierrmpi)
    call MPI_BCAST(v,3*num_particles,MPI_DOUBLE_PRECISION,0,icomm,ierrmpi)
    call MPI_BCAST(q,num_particles,MPI_DOUBLE_PRECISION,0,icomm,ierrmpi)
    call MPI_BCAST(m,num_particles,MPI_DOUBLE_PRECISION,0,icomm,ierrmpi)
    call MPI_BCAST(follow,num_particles,MPI_LOGICAL,0,icomm,ierrmpi)

    nparticles_local = 0

    ! first find ipe and igrid responsible for particle
    do n = 1, num_particles
      call find_particle_ipe(x(:, n),igrid,ipe_particle)

      particle(n)%igrid = igrid
      particle(n)%ipe   = ipe_particle

      if(ipe_particle == mype) then
        check = .true.

        ! Check for user-defined modifications or rejection conditions
        if (associated(usr_check_particle)) call usr_check_particle(igrid, x(:,&
           n), v(:,n), q(n), m(n), follow(n), check)
        if (check) then
          call push_particle_into_particles_on_mype(n)
        else
          cycle
        end if

        nparticles_local = nparticles_local + 1

        call get_lfac_from_velocity(v(:, n), lfac)

        allocate(particle(n)%self)
        particle(n)%self%x      = x(:, n)
        particle(n)%self%q      = q(n)
        particle(n)%self%m      = m(n)
        particle(n)%self%follow = follow(n)
        particle(n)%self%index  = n
        particle(n)%self%time   = global_time
        particle(n)%self%dt     = 0.0d0

        call get_vec(bp, igrid, x(:, n), particle(n)%self%time, b)

        bnorm = norm2(b(:))
        vnorm = norm2(v(:, n))
        vpar  = sum(v(:, n) * b/bnorm)
        vperp = sqrt(vnorm**2 - vpar**2)

        ! The momentum vector u(1:3) is filled with the following components

        ! parallel momentum component (gamma v||)
        particle(n)%self%u(1) = lfac * vpar

        ! Mr: the conserved magnetic moment
        magmom = m(n) * (vperp * lfac)**2 / (2.0d0 * bnorm)
        particle(n)%self%u(2) = magmom

        ! Lorentz factor
        particle(n)%self%u(3) = lfac

        ! initialise payloads for GCA module
        allocate(particle(n)%payload(npayload))
        call gca_update_payload(igrid,x(:,n),particle(n)%self%u(:),q(n),m(n),&
           defpayload,ndefpayload,0.d0)
        particle(n)%payload(1:ndefpayload) = defpayload
        if (associated(usr_update_payload)) then
          call usr_update_payload(igrid,x(:,n),particle(n)%self%u(:),q(n),m(n),&
             usrpayload,nusrpayload,0.d0)
          particle(n)%payload(ndefpayload+1:npayload) = usrpayload
        end if
      end if
    end do

    call MPI_ALLREDUCE(nparticles_local,nparticles,1,MPI_INTEGER,MPI_SUM,icomm,&
       ierrmpi)

  end subroutine gca_create_particles