advect_integrate_particles Subroutine

private subroutine advect_integrate_particles(end_time)

Uses

  • proc~~advect_integrate_particles~~UsesGraph proc~advect_integrate_particles advect_integrate_particles module~mod_global_parameters mod_global_parameters proc~advect_integrate_particles->module~mod_global_parameters module~mod_odeint mod_odeint proc~advect_integrate_particles->module~mod_odeint module~mod_usr_methods mod_usr_methods proc~advect_integrate_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_comm_lib mod_comm_lib module~mod_odeint->module~mod_comm_lib module~mod_variables->module~mod_basic_types

Arguments

Type IntentOptional Attributes Name
double precision, intent(in) :: end_time

Calls

proc~~advect_integrate_particles~~CallsGraph proc~advect_integrate_particles advect_integrate_particles proc~advect_get_particle_dt advect_get_particle_dt proc~advect_integrate_particles->proc~advect_get_particle_dt proc~advect_update_payload advect_update_payload proc~advect_integrate_particles->proc~advect_update_payload proc~derivs_advect derivs_advect proc~advect_integrate_particles->proc~derivs_advect proc~get_vec_advect get_vec_advect proc~advect_integrate_particles->proc~get_vec_advect proc~odeint odeint proc~advect_integrate_particles->proc~odeint proc~advect_get_particle_dt->proc~derivs_advect proc~limit_dt_endtime limit_dt_endtime proc~advect_get_particle_dt->proc~limit_dt_endtime proc~interpolate_var interpolate_var proc~advect_update_payload->proc~interpolate_var proc~derivs_advect->proc~get_vec_advect proc~get_vec_advect->proc~interpolate_var proc~mpistop mpistop proc~interpolate_var->proc~mpistop mpi_abort mpi_abort proc~mpistop->mpi_abort

Source Code

  subroutine advect_integrate_particles(end_time)
    ! this solves dx/dt=v for particles
    use mod_odeint
    use mod_global_parameters
    use mod_usr_methods, only: usr_create_particles, usr_update_payload
    double precision, intent(in) :: end_time

    double precision, dimension(1:ndir) :: v, x
    double precision                 :: defpayload(ndefpayload)
    double precision                 :: usrpayload(nusrpayload)
    double precision                 :: tloc, tlocnew, dt_p, h1, tk
    double precision, dimension(1:ndir) :: yk, k1, k2, k3, k4
    double precision,parameter       :: eps=1.0d-6, hmin=1.0d-8
    integer                          :: ipart, iipart, igrid
    integer                          :: nok, nbad, ierror

    do iipart=1,nparticles_active_on_mype
      ipart                   = particles_active_on_mype(iipart);
      igrid                   = particle(ipart)%igrid
      igrid_working           = igrid
      dt_p                    = advect_get_particle_dt(particle(ipart),&
          end_time)
      particle(ipart)%self%dt = dt_p

      tloc                    = particle(ipart)%self%time
      x(1:ndir)               = particle(ipart)%self%x(1:ndir)
      tlocnew                 = tloc+dt_p

      ! Position update
      select case (integrator)
      case (Euler) 
        ! Simple forward Euler 
        call derivs_advect(tloc,x,v)
        particle(ipart)%self%x(1:ndir) = particle(ipart)%self%x(1:ndir) + dt_p &
           * v(1:ndir)

      case (RK4)
        ! Runge-Kutta order 4
        tk = tloc
        yk = x
        call derivs_advect(tk,yk,k1)
        tk = tloc + dt_p/2.d0
        yk = x + dt_p*k1/2.d0
        call derivs_advect(tk,yk,k2)
        tk = tloc + dt_p/2.d0
        yk = x + dt_p*k2/2.d0
        call derivs_advect(tk,yk,k3)
        tk = tloc + dt_p
        yk = x + dt_p*k3
        call derivs_advect(tloc,yk,k4)
        x = x + dt_p/6.d0*(k1 + 2.d0*k2 + 2.d0*k3 + k4)
      
      case (ARK4)
        ! Adaptive stepwidth RK4:
        h1 = dt_p/2.0d0
        call odeint(x,ndir,tloc,tlocnew,eps,h1,hmin,nok,nbad,derivs_advect,&
           rkqs,ierror)
  
        if (ierror /= 0) then
          print *, "odeint returned error code", ierror
          print *, "1 means hmin too small, 2 means MAXSTP exceeded"
          print *, "Having a problem with particle", iipart
        end if
      end select
      particle(ipart)%self%x(1:ndir) = x(1:ndir)

      ! Velocity update
      call get_vec_advect(igrid,x,tlocnew,v,vp(1),vp(ndir))
      particle(ipart)%self%u(1:ndir) = v(1:ndir)

      ! Time update
      particle(ipart)%self%time = tlocnew

      ! Update payload
      if(ndefpayload>0) then
        call advect_update_payload(igrid,x,v,0.d0,0.d0,defpayload,ndefpayload,&
           tlocnew)
        particle(ipart)%payload(1:ndefpayload) = defpayload
      end if
      if (associated(usr_update_payload)) then
        call usr_update_payload(igrid,x,v,0.d0,0.d0,usrpayload,nusrpayload,&
           tlocnew)
        particle(ipart)%payload(ndefpayload+1:npayload) = usrpayload
      end if

    end do

  end subroutine advect_integrate_particles