Added ode buffers for multi-stage schemes

src/ode/integrator_interfaces.f

 
 
-subroutine integrator_create_(this,fields,rhsfunc,rhs_fields_discr,info)
-  import ppm_t_integrator_, ppm_v_main_abstr, ppm_p_rhsfunc, ppm_v_field_discr_pair
-  class(ppm_t_integrator_)               :: this
-  class(ppm_v_main_abstr), pointer       :: fields
-  procedure(ppm_p_rhsfunc)               :: rhsfunc
-  class(ppm_v_field_discr_pair), pointer :: rhs_fields_discr
-  integer,                 intent(  out) :: info
+subroutine integrator_create_(this,fields,rhsfunc,rhs_fields_discr,info,options)
+  import ppm_v_main_abstr, ppm_v_field_discr_pair, ppm_t_options
+  import ppm_t_integrator_, ppm_p_rhsfunc
+  class(ppm_t_integrator_)                    :: this
+  class(ppm_v_main_abstr), pointer            :: fields
+  procedure(ppm_p_rhsfunc)                    :: rhsfunc
+  class(ppm_v_field_discr_pair), pointer      :: rhs_fields_discr
+  integer,                      intent(  out) :: info
+  class(ppm_t_options),target,optional,intent(in   ) :: options
 end subroutine
 
 subroutine integrator_destroy_(this,info)
@@ -15,10 +17,11 @@ subroutine integrator_destroy_(this,info)
   integer,             intent(  out) :: info
 end subroutine
 
-subroutine integrator_step_(this,t,dt,info)
+subroutine integrator_step_(this,t,dt,istage,info)
   import ppm_t_integrator_, ppm_kind_double
   class(ppm_t_integrator_)                 :: this
   real(ppm_kind_double)            ,intent(inout)    :: t
   real(ppm_kind_double)            ,intent(in   )    :: dt
+  integer,                          intent(in   )    :: istage
   integer,                          intent(  out)    :: info
 end subroutine

src/ode/ode_interfaces.f

 
-subroutine ode_create_(this,scheme,fields,rhsfunc,rhs_fields_discr,info,kickoff_scheme)
-  import ppm_t_ode_,ppm_v_main_abstr, ppm_p_rhsfunc, ppm_v_field_discr_pair
+subroutine ode_create_(this,scheme,fields,rhsfunc,rhs_fields_discr,info,options,kickoff_scheme)
+  import ppm_v_main_abstr, ppm_v_field_discr_pair, ppm_t_options
+  import ppm_t_ode_, ppm_p_rhsfunc
   class(ppm_t_ode_)                      :: this
   integer,          intent(in   )        :: scheme
   class(ppm_v_main_abstr), pointer       :: fields
   procedure(ppm_p_rhsfunc)               :: rhsfunc
   class(ppm_v_field_discr_pair), pointer :: rhs_fields_discr
   integer,          intent(  out)        :: info
+  class(ppm_t_options),target,optional,intent(in   ) :: options
   integer,optional, intent(in   )        :: kickoff_scheme
 end subroutine
 
@@ -16,11 +18,23 @@ subroutine ode_destroy_(this,info)
   integer,          intent(  out)        :: info
 end subroutine
 
-subroutine ode_step_(this,t,dt,info)
+subroutine ode_step_(this,t,dt,istage,info)
   import ppm_t_ode_, ppm_kind_double
   class(ppm_t_ode_)                      :: this
   real(ppm_kind_double), intent(inout)   :: t
   real(ppm_kind_double), intent(in   )   :: dt
+  integer,               intent(in   )   :: istage
   integer,               intent(  out)   :: info
 end subroutine
 
+subroutine ode_map_push_(this,info)
+  import ppm_t_ode_
+  class(ppm_t_ode_)                      :: this
+  integer,               intent(  out)   :: info
+end subroutine
+
+subroutine ode_map_pop_(this,info)
+  import ppm_t_ode_
+  class(ppm_t_ode_)                      :: this
+  integer,               intent(  out)   :: info
+end subroutine

src/ode/test/ppm_module_ode_typedef.fun

@@ -185,7 +185,7 @@ integer                                        :: nterms
         Assert_Equal(info,0)
         t = 0.0_mk
         dt = 0.1_mk
-        call ode%step(t,dt,info)
+        call ode%step(t,dt,1,info)
         Assert_Equal(info,0)
         
         call Part1%get_field(Field1,wp_1r,info)
@@ -274,7 +274,7 @@ integer                                        :: nterms
         Assert_Equal(info,0)
         t = 0.0_mk
         dt = 0.1_mk
-        call ode%step(t,dt,info)
+        call ode%step(t,dt,1,info)
         Assert_Equal(info,0)
        
         call Part1%get_xp(moved_xp,info)
@@ -400,7 +400,7 @@ integer                                        :: nterms
         Assert_Equal(info,0)
         t = 0.0_mk
         dt = 0.1_mk
-        call ode%step(t,dt,info)
+        call ode%step(t,dt,1,info)
         Assert_Equal(info,0)
         
         IF (ndim.EQ.2) THEN

src/ppm_module_data_ode.f

@@ -50,10 +50,14 @@
         INTEGER, DIMENSION(7)       :: ppm_ode_scheme_memsize
         INTEGER, DIMENSION(7)       :: ppm_ode_scheme_nstages
         INTEGER, DIMENSION(7)       :: ppm_ode_scheme_kickoff
-        DATA ppm_ode_scheme_order /1,2,2,4,2,3,1/
+        DATA ppm_ode_scheme_order   /1,2,2,4,2,3,1/
         DATA ppm_ode_scheme_memsize /1,2,2,4,2,1,0/
         DATA ppm_ode_scheme_nstages /1,2,2,4,2,3,999999/
         DATA ppm_ode_scheme_kickoff /1,2,2,4,2,6,7/
+        ! suggest DATA ppm_ode_scheme_order    /1,2,2,4,1/
+        ! suggest: DATA ppm_ode_scheme_memsize /0,1,1,4,0/
+        ! suggest: DATA ppm_ode_scheme_nstages /1,2,2,4,999999/
+        ! suggest DATA ppm_ode_scheme_kickoff  /1,2,2,4,7/
         
         
         !-----------------------------------------------------------------------

src/ppm_module_numerics_interfaces.f

@@ -62,12 +62,12 @@ type,abstract :: ppm_t_integrator_
     integer                     :: scheme_kickoff
     !!! suitable kick off scheme
     
-    procedure(ppm_p_rhsfunc), pointer, nopass   :: rhsfunc => null()
-    class(ppm_v_main_abstr), pointer       :: fields => null()
-    class(ppm_v_discr_kind), pointer       :: discretizations => null()
-    class(ppm_v_field_discr_pair), pointer :: rhs_fields_discr => null()
-    class(ppm_v_field), pointer            :: changes => null()
-    
+    procedure(ppm_p_rhsfunc),         pointer, nopass :: rhsfunc          => null()
+    class(ppm_v_main_abstr),          pointer         :: fields           => null()
+    class(ppm_v_discr_kind),          pointer         :: discretizations  => null()
+    class(ppm_v_field_discr_pair),    pointer         :: rhs_fields_discr => null()
+    class(ppm_v_field),               pointer         :: changes          => null()
+    class(ppm_v_field), dimension(:), pointer         :: buffers          => null() 
     contains
     procedure(integrator_create_),       deferred :: create
     procedure(integrator_destroy_),      deferred :: destroy
@@ -103,6 +103,8 @@ type,abstract :: ppm_t_ode_
     procedure(ode_create_),        deferred :: create
     procedure(ode_destroy_),       deferred :: destroy
     procedure(ode_step_),          deferred :: step
+    procedure(ode_map_push_),      deferred :: map_push
+    procedure(ode_map_push_),      deferred :: map_pop
 end type ppm_t_ode_
 
 !----------------------------------------------------------------------

src/ppm_module_ode_typedef.f

@@ -23,6 +23,8 @@ type,extends(ppm_t_ode_) :: ppm_t_ode
   procedure :: create => ode_create
   procedure :: destroy => ode_destroy
   procedure :: step => ode_step
+  procedure :: map_push => ode_map_push
+  procedure :: map_pop => ode_map_pop
 end type ppm_t_ode
         
 
@@ -31,7 +33,7 @@ end type ppm_t_ode
 !----------------------------------------------------------------------
 contains
 
-subroutine ode_create(this,scheme,fields,rhsfunc,rhs_fields_discr,info,kickoff_scheme)
+subroutine ode_create(this,scheme,fields,rhsfunc,rhs_fields_discr,info,options,kickoff_scheme)
   use ppm_module_integrator_typedef
   implicit none
 
@@ -41,19 +43,26 @@ subroutine ode_create(this,scheme,fields,rhsfunc,rhs_fields_discr,info,kickoff_s
   procedure(ppm_p_rhsfunc)               :: rhsfunc
   class(ppm_v_field_discr_pair), pointer :: rhs_fields_discr
   integer,          intent(  out)        :: info
+  class(ppm_t_options),target,optional,intent(in   ) :: options
   integer,optional, intent(in   )        :: kickoff_scheme
   !----------------------------------------------------------------------
   !  Variables
   !----------------------------------------------------------------------
   integer                          :: kickoff
 
-  start_subroutine("ode_add_integrator")
+  start_subroutine("ode_create")
 
   select case(scheme)
   case(ppm_param_ode_scheme_eulerf)
     ! allocate changes array
     allocate(ppm_t_eulerf::this%integrator,STAT=info)
-    call this%integrator%create(fields,rhsfunc,rhs_fields_discr,info)
+    call this%integrator%create(fields,rhsfunc,rhs_fields_discr,info,options)
+    or_fail("Creating eulerf failed")
+  case(ppm_param_ode_scheme_sts)
+    ! allocate changes array
+    allocate(ppm_t_sts::this%integrator,STAT=info)
+    call this%integrator%create(fields,rhsfunc,rhs_fields_discr,info,options)
+    or_fail("Creating STS failed")
   case default
     ppm_fail("Integrator not implemented")
   end select
@@ -72,8 +81,13 @@ subroutine ode_create(this,scheme,fields,rhsfunc,rhs_fields_discr,info,kickoff_s
     case(ppm_param_ode_scheme_eulerf)
       ! allocate changes array
       allocate(ppm_t_eulerf::this%kickoff,STAT=info)
-      call this%kickoff%create(fields,rhsfunc,rhs_fields_discr,info)
+      call this%kickoff%create(fields,rhsfunc,rhs_fields_discr,info,options)
      or_fail("Creating eulerf failed")
+    case(ppm_param_ode_scheme_sts)
+      ! allocate changes array
+      allocate(ppm_t_sts::this%kickoff,STAT=info)
+      call this%kickoff%create(fields,rhsfunc,rhs_fields_discr,info,options)
+     or_fail("Creating STS failed")
     case default
       ppm_fail("Integrator not implemented")
     end select
@@ -93,28 +107,47 @@ subroutine ode_destroy(this,info)
   end_subroutine()
 end subroutine ode_destroy
 
-subroutine ode_step(this,t,dt,info)
+subroutine ode_step(this,t,dt,istage,info)
   implicit none
   class(ppm_t_ode)                                   :: this
   real(ppm_kind_double),            intent(inout)    :: t
   real(ppm_kind_double),            intent(in   )    :: dt
+  integer,                          intent(in   )    :: istage
   integer,                          intent(  out)    :: info
   start_subroutine("ode_step")
 
 
 
   if (this%state.EQ.ode_state_init) then
-    call this%integrator%step(t,dt,info)
+    call this%integrator%step(t,dt,istage,info)
     this%state = ode_state_running
   else if (this%state.EQ.ode_state_kickoff) then
-    call this%kickoff%step(t,dt,info)
+    call this%kickoff%step(t,dt,istage,info)
     this%state = ode_state_running
   else if(this%state.EQ.ode_state_running) then
-    call this%integrator%step(t,dt,info)
+    call this%integrator%step(t,dt,istage,info)
   end if
   
   end_subroutine()
 
 end subroutine ode_step
 
+subroutine ode_map_push(this,info)
+  IMPLICIT NONE
+  class(ppm_t_ode)      :: this
+  integer,               intent(  out)   :: info
+  start_subroutine("ode_map_push")
+  
+  end_subroutine()
+end subroutine ode_map_push
+
+subroutine ode_map_pop(this,info)
+  IMPLICIT NONE
+  class(ppm_t_ode)      :: this
+  integer,               intent(  out)   :: info
+  start_subroutine("ode_map_pop")
+  
+  end_subroutine()
+end subroutine ode_map_pop
+
 end module ppm_module_ode_typedef

src/ppm_numerics.h

@@ -61,8 +61,8 @@
       INTEGER, PARAMETER :: ppm_param_ode_scheme_tvdrk2   = 2
       INTEGER, PARAMETER :: ppm_param_ode_scheme_midrk2   = 3
       INTEGER, PARAMETER :: ppm_param_ode_scheme_rk4      = 4
-      INTEGER, PARAMETER :: ppm_param_ode_scheme_trapez   = 5
-      INTEGER, PARAMETER :: ppm_param_ode_scheme_tvdrk3   = 6
+      !INTEGER, PARAMETER :: ppm_param_ode_scheme_trapez   = 5
+      !INTEGER, PARAMETER :: ppm_param_ode_scheme_tvdrk3   = 6
       INTEGER, PARAMETER :: ppm_param_ode_scheme_sts      = 7