ode_def.f90

Go to the documentation of this file.

module ode_def

  use choral_constants
  use choral_variables
  use real_type
  use basic_tools

  implicit none
  private

  public :: s_prefactor
  public :: name_ode_method
  public :: order_ode_method
  public :: overflow

  !! ABSTRACT INTERFACES
  public :: ode_reaction
  public :: linsystem_solver

  !       %----------------------------------------%
  !       |                                        |
  !       |        ABSTRACT INTERFACES             |
  !       |                                        |
  !       %----------------------------------------%
  abstract interface

     subroutine ode_reaction(a, b, x, t, y, N, Na)
       import :: rp
       real(RP), dimension(Na), intent(out) :: a
       real(RP), dimension(N) , intent(out) :: b
       real(RP), dimension(3) , intent(in)  :: x
       real(RP)               , intent(in)  :: t
       real(RP), dimension(N) , intent(in)  :: y
       integer                , intent(in)  :: N, Na
     end subroutine ode_reaction

     subroutine linsystem_solver(x, ierr, y)
       import :: rp
       real(RP), dimension(:), intent(inout) :: x
       logical               , intent(out)   :: ierr
       real(RP), dimension(:), intent(in)    :: y
     end subroutine linsystem_solver

  end interface


contains


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

    select case(method)

    case(ode_be)
       name="Bwd-Euler"

    case(ode_fe)
       name="Fwd-Euler"

    case(ode_cn)
       name="CN"

    case(ode_sdirk4)
       name="SDIRK4"

    case(ode_rk2)
       name="RK2"

    case(ode_erk1)
       name="Exp. Euler"

    case(ode_erk2_a)
       name="ERK2 (A)"

    case(ode_erk2_b)
       name="ERK2 (B)"

    case(ode_modif_erk2_b)
       name="Modif. ERK2 (B)"

    case(ode_rk4)
       name="RK4"

    case(ode_fbe)
       name="Fwd-Bwd-Euler"

    case(ode_bdfsbdf2)
       name="BDF2-SBDF2"

    case(ode_rl2)
       name="RL2"

    case(ode_eab2)
       name="Exp-Adams-2"

    case(ode_cnab2)
       name="CN-AB2"

    case(ode_mcnab2)
       name="Mod-CN-AB2"

    case(ode_bdfsbdf3)
       name="BDF3-SBDF3"

    case(ode_rl3)
       name="RL3"

    case(ode_eab3)
       name="Exp-Adams-3"

    case(ode_bdfsbdf4)
       name="BDF4-SBDF4"

    case(ode_rl4)
       name="RL4"

    case(ode_eab4)
       name="Exp-Adams-4"

    case(ode_bdfsbdf5)
       name="BDF5-SBDF5"

    case(ode_dc_2)
       name="DC-2"

    case(ode_dc_3)
       name="DC-3"

    case default
       name="invalid"

    end select

  end function name_ode_method


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

    !! invalid method
    if ( (method<=0).OR.(method>ode_tot_nb) ) then
       o = -1
       return
    end if

    !! default
    o = 2

    select case(method)

    case(ode_be,ode_fe,ode_fbe, ode_erk1)
       o = 1

    case(ode_rl3,ode_eab3,ode_bdfsbdf3)
       o = 3

    case(ode_sdirk4,ode_rk4,ode_rl4,ode_eab4,ode_bdfsbdf4)
       o = 4

    case(ode_bdfsbdf5)
       o = 5

    case(ode_dc_2)
       o = 2

    case(ode_dc_3)
       o = 3

    end select

  end function order_ode_method



  function s_prefactor(method, dt) result(cs)
    integer , intent(in) :: method
    real(RP), intent(in) :: dt
    real(RP)             :: Cs

    select case(method)
    case(ode_fe)
       cs = 0.0_rp

    case(ode_fbe, ode_be, ode_dc_2, ode_dc_3)
       cs = dt

    case(ode_bdfsbdf2)
       cs = dt * 2._rp / 3._rp

    case(ode_cnab2)
       cs = dt / 2._rp

    case(ode_mcnab2)
       cs = dt * 9._rp / 16._rp

    case(ode_bdfsbdf3)
       cs = dt * 6._rp / 11._rp

    case(ode_bdfsbdf4)
       cs = dt * 12._rp / 25._rp

    case(ode_bdfsbdf5)
       cs = dt * 60._rp / 137._rp

    case(ode_cn)
       cs = dt / 2._rp

    case(ode_sdirk4)
       cs = 1._rp/sqrt(3._rp) * cos(pi/18._rp) 
       cs = cs + 0.5_rp
       cs = dt * cs

    case default
       cs = 0.0_rp
       call quit( "ode_def: S_prefactor: unknown method" )

    end select

  end function s_prefactor


  function overflow(yy) result(bool)
    logical                              :: bool
    real(RP), dimension(:,:), intent(in) :: yy

    real(RP):: nrm

    bool = .false.

    nrm   = maxval( abs(yy) )

    if ((nrm>1.e4_rp) .OR. (isnan(nrm)) ) then
       bool = .true.

       if (choral_verb>0) then
          write(*,*)"ode_def         : OVERFLOW" 
          if (choral_verb>2) call print(yy)
       end if

    end if

  end function overflow

end module ode_def