Time integration |
Example: we want to integrate 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
|