ode_opSplt_mod.f90

Go to the documentation of this file.

module ode_opsplt_mod

  use choral_constants
  use choral_variables
  use real_type
  use basic_tools
  use r1d_mod, only: copy
  use krylov_mod
  use r1d_mod, only: xpay
  use ode_def
  use ode_problem_mod
  use ode_solution_mod
  use ode_nl_1s_mod
  use ode_lin_1s_mod

  implicit none
  private

  !       %----------------------------------------%
  !       |                                        |
  !       |          PUBLIC DATA                   |
  !       |                                        |
  !       %----------------------------------------%
  !! DERIVED TYPE
  public :: ode_opsplt

  !! GENERIC SUBROUTINES
  public :: clear
  public :: valid
  public :: print

  !! NON GENERIC SUBROUTINES
  public :: solve_ode_opsplt
  public :: name_ode_method_opsplt
  public :: order_ode_method_opsplt
  public :: create_ode_opsplt_sol



  !       %----------------------------------------%
  !       |                                        |
  !       |          DERIVED TYPE                  |
  !       |                                        |
  !       %----------------------------------------%
  type :: ode_opsplt

     integer :: method = -1

     integer :: ns = -1

     integer :: l_meth = -1

     integer :: nl_meth = -1

     real(RP), dimension(:), allocatable :: a_nl 

     real(RP), dimension(:), allocatable :: a_l 

     logical , dimension(:), allocatable :: nl_first 
     
     procedure(ode_lin_solver  ), nopass, pointer :: slv_l =>null()

     procedure(ode_nl_1s_solver), nopass, pointer :: slv_nl=>null()

     logical :: check_overflow = .false.

   contains

     final :: ode_opsplt_clear

  end type ode_opsplt


  !       %----------------------------------------%
  !       |                                        |
  !       |       GENERIc SUBROUTINES              |
  !       |                                        |
  !       %----------------------------------------%

  interface clear
     module procedure ode_opsplt_clear
  end interface clear

  interface ode_opsplt
     module procedure ode_opsplt_create
  end interface ode_opsplt

  interface valid
     module procedure ode_opsplt_valid
  end interface valid

  interface print
     module procedure ode_opsplt_print
  end interface print


contains

  subroutine ode_opsplt_clear(os)
    type(ode_opsplt), intent(inout):: os

     os%method  = -1
     os%ns      = -1
     os%L_meth  = -1
     os%NL_meth = -1
     os%check_overflow = .false.

     call freemem(os%a_NL)
     call freemem(os%a_L)
     call freemem(os%NL_first)

     os%slv_L =>null()
     os%slv_NL=>null()


  end subroutine ode_opsplt_clear


  function ode_opsplt_create(method, L_meth, NL_meth, &
       &                    check_overflow) result(os)
    type(ode_opsplt)              :: os
    integer          , intent(in) :: method, L_meth, NL_meth
    logical, optional, intent(in) :: check_overflow

    logical :: bool

    if (choral_verb>0) write(*,*) &
         &"ode_opSplt_mod  : ode_opSplt_create"
    call clear(os)
    
    os%method = method
    os%L_meth  = l_meth
    os%NL_meth = nl_meth 

    bool = valid(os)
    if (.NOT.bool) call quit("ode_opSplt_mod: ode_opSplt_set")

    call set_solver_ode_lin_1s(os%slv_L, l_meth)
    call set_solver_ode_nl_1s(os%slv_NL, nl_meth)
    call def_splitting(os)

    !! check overflow when solving
    if (present(check_overflow)) then
       os%check_overflow = check_overflow
    end if

    if (.NOT.valid(os)) call quit(&
         & "ode_opSplt_mod: ode_opSplt_create: not valid")

  end function ode_opsplt_create


  subroutine def_splitting(os)
    type(ode_opsplt), intent(inout) :: os
    
    real(RP) :: theta
    
    !! number of stages
    !!
    select case(os%method)
    case(ode_os_lie)
       os%ns = 1
    case(ode_os_strang)
       os%ns = 2
    case(ode_os_ruth, ode_os_aks3)
       os%ns = 3
    case(ode_os_yoshida)        
       os%ns = 4
    case default
       call quit( 'rdEq_mod: def_splitting: 1' )
    end select
    
    call allocmem(os%a_NL, os%ns)
    call allocmem(os%a_L , os%ns)
    call allocmem(os%NL_first, os%ns)
    
    !! Definition of a_L(:), a_NL(:) and NL_first(:)
    !!
    os%NL_first = .true.       ! default = reaction first
    select case(os%method)
    case(ode_os_lie)
       os%a_NL=(/1._rp/)
       os%a_L=(/1._rp/)
       
    case(ode_os_strang)
       os%a_NL=(/.5_rp, .5_rp/)
       os%a_L=(/.5_rp, .5_rp/)
       os%NL_first(2) = .false.
       
    case(ode_os_ruth)
       os%a_L=(/ re(2,3) ,-re(2,3), 1._rp/)
       os%a_NL=(/ re(7,24), re(3,4),-re(1,24)/)
       os%NL_first = .true.
       
    case(ode_os_aks3)
       os%a_L(1)= 0.91966152301739985_rp
       os%a_L(2)=-0.18799161879915978_rp
       os%a_L(3)= 0.26833009578175992_rp
       
       os%a_NL(1)=os%a_L(3)
       os%a_NL(2)=os%a_L(2)
       os%a_NL(3)=os%a_L(1)
       
       os%NL_first = .true.
       
    case(ode_os_yoshida)
       theta = 1._rp/(2._rp - 2._rp ** ( 1._rp / 3._rp ) )
       
       os%a_NL=(/ theta/ 2._rp, (1._rp-theta)/2._rp, &
            &            (1._rp-theta)/2._rp, theta/2._rp /)
       
       os%a_L=(/ theta, 1._rp-2._rp*theta, theta, 0._rp/)
       
    end select
    
  end subroutine def_splitting


  !! check ode_opSplt
  !!
  function ode_opsplt_valid(os) result(bool)
    type(ode_opsplt), intent(in):: os

    logical :: bool

    bool = (os%method>0) .AND. (os%method <= ode_os_tot_nb)
    bool = bool .AND. check_ode_method_lin_1s(os%L_meth)
    bool = bool .AND. check_ode_method_nl_1s(os%NL_meth)

  end function ode_opsplt_valid


  !! print ode_opSplt
  !!
  subroutine ode_opsplt_print(os)
    type(ode_opsplt), intent(in):: os

    write(*,*)"ode_opSplt_mod  : print"

    write(*,*)"  operator spltiting Method      = ",&
         & trim( name_ode_method_opsplt(os%method) )

    write(*,*)"  one step Non-linear ODE solver = ",&
         & trim(name_ode_method(os%NL_meth))
    
    write(*,*)"  one step linear ODE solver     = ",&
         & trim(name_ode_method(os%L_meth))

    write(*,*)"  check overfow                  =", os%check_overflow

  end subroutine ode_opsplt_print


  function name_ode_method_opsplt(method) result(name)
    integer, intent(in) :: method
    character(len=15)   :: name

    select case(method)
    case(ode_os_lie)
       name="Lie OS" 

    case(ode_os_strang)
       name="Strang OS" 

    case(ode_os_ruth)
       name="Ruth OS" 

    case(ode_os_aks3)
       name="AKS3 OS" 

    case(ode_os_yoshida)
       name="Yoshida OS" 

    case default
       name="invalid"

    end select

  end function name_ode_method_opsplt


  function order_ode_method_opsplt(method) result(o)
    integer, intent(in) :: method
    integer             :: o

    select case(method)
    case(ode_os_lie)
       o = 1

    case(ode_os_strang)
       o = 2

    case(ode_os_ruth, ode_os_aks3)
       o = 3

    case(ode_os_yoshida)
       o = 4

    case default
       o= -1

    end select

  end function order_ode_method_opsplt



  subroutine create_ode_opsplt_sol(sol, pb, os)
    type(ode_solution), intent(inout) :: sol
    type(ode_problem) , intent(in)    :: pb
    type(ode_opsplt)  , intent(in)    :: os
    
    integer :: n_V, n_Y, n_FY
    logical :: bool

    bool = valid(os)
    if (.NOT.bool) call quit(&
         & "ode_opSplt_mod: create_ode_opSplt_sol: 1")

    call memsize_ode_nl_1s(n_y, n_fy, os%NL_meth)
    call memsize_ode_lin_1s(n_v, os%L_meth)

    sol = ode_solution(pb, nv=n_v, ny=n_y, nfy=n_fy)

  end subroutine create_ode_opsplt_sol
  

  subroutine solve_ode_opsplt(sol, t0, T, dt, os, pb, kry, out)
    type(ode_solution), intent(inout) :: sol
    real(RP)          , intent(in)    :: dt, t0, T
    type(ode_opsplt)  , intent(in)    :: os
    type(ode_problem) , intent(in)    :: pb
    type(krylov)      , intent(inout) :: kry
    procedure(ode_output_proc)                 :: out

    real(RP) :: CS, tn, ts, h
    integer  :: ii, jj, nStep
    logical  :: exit_comp, ierr

    if (choral_verb>1) write(*,*) &
         &"ode_opSplt_mod  : solve"
    if (choral_verb>2) write(*,*) &
         &"  Operator splitting method      = ",&
         & name_ode_method_opsplt(os%method)
    if (choral_verb>2) write(*,*) &
         &"  NonLin  onestep solver         = ",&
         & name_ode_method(os%NL_meth)
    if (choral_verb>2) write(*,*) &
         &"  Linear  onestep solver         = ",&
         & name_ode_method(os%L_meth)

    sol%ierr  = 0
    nstep     = int( (t-t0) / dt) 
    exit_comp = .false. 

    !! Time loop
    !!
    time_loop: do ii=1, nstep
       tn = t0 + re(ii-1)*dt

       ! output
       call out(tn, sol, exit_comp)
       if (exit_comp) then
          if (choral_verb>2) &
               &write(*,*)"ode_opSplt_mod   :&
               & solve  = ", tn, ': EXIT_COMPp'
          return
       end if

       !! Loop on operator splitting stages
       !!
       ts = tn
       stage_loop: do jj=1, os%ns

          if (os%NL_first(jj)) then    
             !! Stage starting with the non-linear solver

             !! non-linear equation
             h = dt * os%a_NL(jj)
             call os%slv_NL(sol, h, ts, pb)
             call copy(sol%V(:,1), sol%Y(sol%N, 1, :))
             if (os%check_overflow) then
                ierr = overflow(sol%Y(:,1,:))
                if (ierr) then
                   sol%ierr = 10
                   exit time_loop
                end if
             end if

             !! linear equation
             h  = dt * os%a_L(jj)
             cs = s_prefactor(os%L_meth, h)
             call os%slv_L(sol, ierr, h, pb, kinv)
             call copy(sol%Y(sol%N,1,:), sol%V(:,1))
             if (ierr) then
                sol%ierr = 1
                exit time_loop
             end if

          else        
             !! Stage starting with the linear equationN
       
             !! linear equation
             h  = dt * os%a_L(jj)
             cs = s_prefactor(os%L_meth, h)
             call os%slv_L(sol, ierr, h, pb, kinv)
             call copy(sol%Y(sol%N,1,:), sol%V(:,1))

             if (ierr) then
                sol%ierr = 1
                exit time_loop
             end if

             ! non-linear equation       
             h = dt * os%a_NL(jj)
             call os%slv_NL(sol, h, ts, pb)
             call copy(sol%V(:,1), sol%Y(sol%N,1,:))

             if (os%check_overflow) then
                ierr = overflow(sol%Y(:,1,:))
                if (ierr) then
                   sol%ierr = 10
                   exit time_loop
                end if
             end if

          end if

          ts = ts + os%a_NL(jj)*dt
       end do stage_loop

    end do time_loop

    if (sol%ierr>0) then

       if (choral_verb>1) then
          select case(sol%ierr)
          case(10)
             write(*,*) "ode_opSplt_mod  : solve,&
                  & time =", real(tn, SP), ", stage =", int(jj,1), &
                  & ': OVERFLOW Y'

          case(1) 
             write(*,*) "ode_opSplt_mod  : solve,&
                  & time =", real(tn, SP), ", stage =", int(jj,1), &
                  & ': linear pb inversion error'

          case default
             write(*,*) "ode_opSplt_mod  : solve,&
                  & time =", real(tn, SP), ", stage =", int(jj,1), &
                  & ': ERROR'

          end select
       end if
    end if

  contains

    subroutine kinv(x, bool, b)    
      real(RP), dimension(:), intent(inout) :: x
      logical               , intent(out)   :: bool
      real(RP), dimension(:), intent(in)    :: b

      call solve(x, kry, b, k)
      bool = kry%ierr
      
    end subroutine kinv

    !! matrix-vector product x --> K*x, K = M + CS*S
    !!
    subroutine k(y, x)
      real(RP), dimension(:), intent(out) :: y
      real(RP), dimension(:), intent(in)  :: x

      call pb%M(y , x)
      call pb%S(sol%aux, x)
      call xpay(y, cs, sol%aux)

    end subroutine k

  end subroutine solve_ode_opsplt


end module ode_opsplt_mod