Time integration

 

Example: we want to integrate du_p/dt = \lambda u_p on the particles using the 4th order Runge-Kutta of the ODE suite. We'll assume that we have np particles with positions xp and values up already allocated and initialized.

Local variables:

REAL(mk), DIMENSION(:,:),  POINTER :: dup     ! du/dt on particles
REAL(mk), DIMENSION(:,:),  POINTER :: bfr     ! storage space for the stages
REAL(mk), DIMENSION(4)             :: time    ! time things
INTEGER                            :: istage  ! stage counter
INTEGER                            :: nstages ! number of stages
INTEGER                            :: bfrsz   ! size of the buffer "bfr"
INTEGER                            :: scheme  ! which scheme to use
INTEGER                            :: odeid   ! handle on the solver
LOGICAL                            :: adapt   ! use adaptive time step
INTEGER                            :: lda     ! leading dimension of our mode
INTEGER, EXTERNAL                  :: MyRHS   ! your implementation of the RHS

 

Initialize the ODE suite:

!-------------------------------------------------------------
!  Initialize the Ode solver
!-------------------------------------------------------------
CALL ppm_ode_init (info)

 

Create Mode:

scheme = PPM_PARAM_ODE_SCHEME_RK4  ! we want the 4th order RK scheme
odeid  = -1 ! let the PPM choose an ID for us
adapt  = .FALSE.! don't need adaptive time stepping
lda    = 2
!-------------------------------------------------------------
!  Create the mode
!-------------------------------------------------------------
CALL ppm_ode_create_ode(odeid, bfrsz, nstages, scheme, scheme, adapt, info)

 

Allocate space for the stages:

ALLOCATE(bfr(bfrsz*lda,np))

 

Set the time:

dt      = 0.1
time(1) = 0.0  ! set the start time
time(2) = 1.0  ! set the end time
time(3) = 0.0  ! set the current time
time(4) = dt   ! set the time step size

 

Start the ODE solver:

CALL ppm_ode_start(info)

 

Start the time integration:

DO WHILE(.NOT.ppm_ode_alldone(info))
 
DO istage=1,nstages
 
CALL ppm_ode_step(odeid, xp, up, dup, lda, np, &
&            bfr, istage, time, MyRHS, info=info)
 
!-- say particles move, then we need to map after each stage
maptype = ppm_param_map_partial
CALL ppm_map_part(xp,3,np,mpart,topo_id,maptype,info)
maptype = ppm_param_map_push
CALL ppm_map_part(up,lda,np,mpart,topo_id,maptype,info)
CALL ppm_map_part(dup,lda,np,mpart,topo_id,maptype,info)
!-- now have the ode suite map the stages
CALL ppm_ode_map_push(odeid,bfr,lda,np,mpart,info)
!-- send
maptype = ppm_param_map_send
CALL ppm_map_part(dup,lda,np,mpart,topo_id,maptype,info)
!-- pop in the reverse order
CALL ppm_ode_map_pop(odeid,bfr,lda,np,mpart,info)
maptype = ppm_param_map_pop
CALL ppm_map_part(dup,lda,np,mpart,topo_id,maptype,info)
CALL ppm_map_part(up,lda,np,mpart,topo_id,maptype,info)
CALL ppm_map_part(xp,3,np,mpart,topo_id,maptype,info)
 
END DO
 
END DO
 
CALL ppm_ode_finalize(info)

 

A possible implementation for the right-hand side function:

FUNCTION MyRHS(vxp, vup, vdup, vdime, vnp, rpack, ipack, lpack, info)
USE myGlobalData
 
!-- Arguments
INTEGER, INTENT(in)               :: vdime, vnp
REAL(mk),DIMENSION(:,:),POINTER     :: vxp, vup, vdup
REAL(MK),DIMENSION(:,:), POINTER,OPTIONAL  :: rpack
INTEGER, DIMENSION(:,:), POINTER,OPTIONAL  :: ipack
LOGICAL, DIMENSION(:,:), POINTER,OPTIONAL  :: lpack
INTEGER, INTENT(inout)            :: info
 
INTEGER :: MyRHS
 
!-- Local variables
INTEGER :: p
 
!-- Compute the right-hand side
!   assuming the parameter REAL(mk), DIMENSION(2) :: lambda
!   is specified in the module myGlobalData
DO p=1,vnp
dup(1,p) = -lambda(1) * up(1,p)
dup(2,p) = -lambda(2) * up(2,p)
END DO
 
!-- bogus return value
MyRHS = 123456
RETURN
END FUNCTION MyRHS