ode_solver_mod.f90

Go to the documentation of this file.

module ode_solver_mod

  use choral_constants
  use choral_variables
  use real_type
  use basic_tools
  use krylov_mod
  use ode_def
  use ode_problem_mod
  use ode_solution_mod
  use ode_nl_1s_mod
  use ode_nl_ms_mod
  use ode_lin_1s_mod
  use ode_sl_ms_mod
  use ode_sl_nl_ms_mod
  use ode_sl_nl_dc_mod
  use ode_opsplt_mod
  use ode_output_mod

  implicit none
  private

  !       %----------------------------------------%
  !       |                                        |
  !       |          PUBLIC DATA                   |
  !       |                                        |
  !       %----------------------------------------%

  !! DERIVED TYPE
  public :: ode_solver

  !! GENERIC SUBROUTINES
  public :: clear
  public :: ode_solution
  public :: valid
  public :: print
  public :: solve
  public :: initialcond

  !! NON GENERIC SUBROUTINES
  public :: name_ode_solver_type
  public :: check_ode_method
  public :: set_ode_solver_output
  
  !       %----------------------------------------%
  !       |                                        |
  !       |          DERIVED TYPE                  |
  !       |                                        |
  !       %----------------------------------------%

  type :: ode_solver

     integer :: verb = 0

     integer :: type = -1

     integer :: pb = -1

     type(ode_opsplt) :: os 

     integer :: l_meth = -1

     integer :: sl_meth = -1

     integer :: dc_meth = -1

     integer :: nl_meth = -1

     type(krylov) :: kry        

     logical :: check_overflow = .false.

     procedure(ode_output_proc), nopass, pointer :: output => void_ode_output

   contains

     final :: ode_solver_clear

  end type ode_solver


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

  interface clear
     module procedure ode_solver_clear
  end interface clear

  interface ode_solver
     module procedure ode_solver_create
  end interface ode_solver

  interface valid
     module procedure ode_solver_valid
  end interface valid

  interface print
     module procedure ode_solver_print
  end interface print

  interface ode_solution
     module procedure ode_solution_create
  end interface ode_solution

  interface solve
     module procedure ode_solver_solve
  end interface solve

  interface initialcond
     module procedure homogeneous_initialcond
  end interface initialcond

contains


  subroutine ode_solver_clear(slv)
    type(ode_solver), intent(inout) :: slv

    slv%verb    = 0
    slv%type    = -1
    slv%pb      = -1
    slv%L_meth  = -1
    slv%SL_meth = -1
    slv%DC_meth = -1
    slv%NL_meth = -1
    slv%check_overflow = .false.
    slv%output => void_ode_output
    call clear(slv%os)
    call clear(slv%kry)

  end subroutine ode_solver_clear


  function ode_solver_create(pb, type, &
       &               os_meth, SL_meth, NL_meth, L_meth, &
       &               DC_meth, &
       &               check_overflow, verb) result(slv)
    type(ode_solver)              :: slv
    type(ode_problem), intent(in)  :: pb
    integer          , intent(in)  :: type
    integer, optional, intent(in)  :: os_meth, SL_meth, DC_meth
    integer, optional, intent(in)  :: NL_meth, L_meth, verb
    logical, optional, intent(in)  :: check_overflow

    logical :: bool

    if (choral_verb>0) write(*,*) &
         &"ode_solver_mod  : ode_solver_create"

    call clear(slv)

    !! default krylov
    !!
    slv%kry = krylov(type=kry_cg, tol=1e-6_rp, itmax=100, verb=0)

    !! Asociated ODE problem
    !!
    bool = valid(pb)  
    if (.NOT.bool) call quit(&
         & "ode_solver_mod: ode_solver_create: 'pb' not valid")
    slv%pb = pb%type

    !! Type
    !!
    bool = (type>0) .AND. (type <= ode_slv_tot_nb)
    if (.NOT.bool) call quit(&
         &"ode_solver_mod: ode_solver_create:&
         &  wrong argument 'type'")
    slv%type = type

    if (present(nl_meth)) slv%NL_meth =  nl_meth
    if (present(sl_meth)) slv%SL_meth =  sl_meth
    if (present(dc_meth)) slv%DC_meth =  dc_meth
    if (present(l_meth )) slv%L_meth  =  l_meth
    if (present(verb))    slv%verb    = verb

    !! Opêrator splitting
    if (present(os_meth)) then
       slv%os = ode_opsplt(os_meth, slv%L_meth, slv%NL_meth, &
            & check_overflow)  
    end if

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

    if (.NOT.valid(slv)) call quit(&
         & "ode_solver_mod: ode_solver_create: not valid")

  end function ode_solver_create


  function ode_solver_valid(slv) result(bool)
    type(ode_solver) , intent(in) :: slv

    logical :: bool

    !! check type
    bool = (slv%type>0) .AND. (slv%type <= ode_slv_tot_nb)
    if (.NOT.bool) return

    !! check operator splitting
    if (slv%type==ode_slv_os) bool = valid(slv%os)
    if (.NOT.bool) return

    !! check ODE problem
    bool = (slv%pb>0) .AND. (slv%pb<=ode_pb_tot_nb)
    if (.NOT.bool) return

    !! check ODE solver depending on the ODE type
    select case(slv%pb)

    case(ode_pb_lin)

       select case(slv%type)
       case(ode_slv_1s)
          bool = check_ode_method_lin_1s(slv%L_meth)
       end select

    case(ode_pb_nl)

       select case(slv%type)
       case(ode_slv_1s)

          bool = check_ode_method_nl_1s(slv%NL_meth)

       case(ode_slv_ms)

          bool = check_ode_method_nl_ms(slv%NL_meth)

       end select

    case(ode_pb_sl)

       select case(slv%type)
       case(ode_slv_ms)
          
          bool = check_ode_method_sl_ms(slv%SL_meth)

       end select
       
    case(ode_pb_sl_nl)

       select case(slv%type)

       case(ode_slv_ms)
          bool = check_ode_method_sl_ms(slv%SL_meth)
          bool = bool .AND. check_ode_method_nl_ms(slv%NL_meth)

       case(ode_slv_dc)
          bool = check_ode_method_sl_nl_dc(slv%DC_meth)
       end select

    end select

  end function ode_solver_valid


  subroutine ode_solver_print(slv) 
    type(ode_solver) , intent(in) :: slv

    write(*,*)"ode_solver_mod  : print"

    write(*,*)"  Associated ODE problem         = ",&
         & trim(name_ode_problem(slv%pb))

    write(*,*)"  ODE solver type                = ",&
         & trim(name_ode_solver_type(slv%type))

    if ( valid(slv)) then
       write(*,*)"  Status                         = valid"
    else
       write(*,*)"  Status                         = invalid"
    end if

    select case(slv%type)
    case(ode_slv_os)
       call print(slv%os)

    case(ode_slv_1s)
       if (slv%L_meth /= -1 ) then
          write(*,*)"  one-step linear solver         = ",&
               & trim(name_ode_method(slv%L_meth))
       end if
       if (slv%NL_meth /= -1 ) then
          write(*,*)"  one-step non-linear solver     = ",&
               & trim(name_ode_method(slv%NL_meth))
       end if

    case(ode_slv_ms)
       if (slv%NL_meth /= -1 ) then
          write(*,*)"  multistep non-linear solver    = ",&
               & trim(name_ode_method(slv%NL_meth))
       end if
       if (slv%SL_meth /= -1 ) then
          write(*,*)"  multistep semilinear solver    = ",&
               & trim(name_ode_method(slv%SL_meth))
       end if

    case(ode_slv_dc)
       write(*,*)"  deferred correction ode solver = ",&
            & trim(name_ode_method(slv%DC_meth))

    end select
    write(*,*)"  check overfow                  =", slv%check_overflow
    write(*,*)"  verbosity                      =", slv%verb

  end subroutine ode_solver_print



  function name_ode_solver_type(type) result(name)
    integer, intent(in) :: type
    character(len=20)   :: name

    select case(type)
    case(ode_slv_1s)
       name="one-step"

    case(ode_slv_ms)
       name="multistep"

    case(ode_slv_dc)
       name="deferred corrections"

    case(ode_slv_os)
       name="op splitting"

    case default
       name = "invalid"

    end select

  end function name_ode_solver_type

  function check_ode_method(method, pb_type, slv_type) result(b)
    logical             :: b
    integer, intent(in) :: method, pb_type, slv_type

    b = .false.
    select case(pb_type)

    case(ode_pb_lin)

       select case(slv_type)
       case(ode_slv_1s)
          b = check_ode_method_lin_1s(method)
       
       end select

    case(ode_pb_nl)

       select case(slv_type)
       case(ode_slv_ms)
          b = check_ode_method_nl_ms(method)

       case(ode_slv_1s)
          b = check_ode_method_nl_1s(method)

       end select

    case(ode_pb_sl)

       select case(slv_type)
       case(ode_slv_ms)
          b = check_ode_method_sl_ms(method)

       end select

    case(ode_pb_sl_nl)

       select case(slv_type)
       case(ode_slv_dc)
          b = check_ode_method_sl_nl_dc(method)

       end select

    end select

  end function check_ode_method



  function ode_solution_create(slv, pb) result(sol)
    type(ode_solution)                :: sol
    type(ode_solver)  , intent(in)    :: slv
    type(ode_problem) , intent(in)    :: pb

    logical :: bool

    if (choral_verb>0) write(*,*) &
         &"ode_solver_mod  : ode_solution_create"

    bool = valid(pb)  
    if (.NOT.bool) call quit(&
         & "ode_solver_mod: ode_solution_create: 'pb' not valid")

    bool = valid(slv)  
    if (.NOT.bool) call quit(&
         & "ode_solver_mod: ode_solution_create: 'slv' not valid")

    !! define ODE solution
    select case(pb%type)

    case(ode_pb_lin)

       select case(slv%type)
       case(ode_slv_1s)
          call create_ode_lin_1s_sol(sol, pb, slv%L_meth)

       end select

    case(ode_pb_nl)

       select case(slv%type)
       case(ode_slv_1s)
          call create_ode_nl_1s_sol(sol, pb, slv%NL_meth)

       case(ode_slv_ms)
          call create_ode_nl_ms_sol(sol, pb, slv%NL_meth)

       end select

    case(ode_pb_sl)

       select case(slv%type)
       case(ode_slv_os)
          call create_ode_opsplt_sol(sol, pb, slv%os)

       case(ode_slv_ms)
          call create_ode_sl_ms_sol(sol, pb, slv%SL_meth)

       end select

    case(ode_pb_sl_nl)

       select case(slv%type)
       case(ode_slv_os)
          call create_ode_opsplt_sol(sol, pb, slv%os)

       case(ode_slv_ms)
          call create_ode_sl_nl_ms_sol(sol, pb, &
               & slv%SL_meth, slv%NL_meth)

       case(ode_slv_dc)
          call create_ode_sl_nl_dc_sol(sol, pb, slv%DC_meth)

       end select

    end select

  end function ode_solution_create
  

  subroutine homogeneous_initialcond(sol, pb, slv, t0, y0)
    type(ode_solution)    , intent(inout):: sol
    type(ode_problem )    , intent(in)   :: pb
    type(ode_solver)      , intent(in)   :: slv
    real(RP)              , intent(in)   :: t0
    real(RP), dimension(:), intent(in)   :: y0

    integer  :: jj, ii, N, Na

    if (choral_verb>2) write(*,*) &
         &"ode_solver_mod  :&
         & initialCond    = homogeneous"

    if (size(y0)/=sol%N) call quit( &
         & "ode_solver_mod: homogeneous_initialCond:&
         & wrong size for Y0" )

    n  = sol%N
    na = sol%Na

    associate( y=>sol%Y, by=>sol%BY, ay=>sol%AY, &
         & x=>pb%X, v=>sol%V)

      !! initialise Y
      do ii=1, sol%dof
         do jj=1, sol%nY
            y(:,jj,ii) = y0(:) 
         end do
      end do

      !! initialise V
      do jj=1, sol%nV
         v(:,jj) = y0(n)
      end do

      !! initialise AY and BY
      if (sol%NFY==0) return

      do ii=1, sol%dof
         call pb%AB(ay(:,1,ii), by(:,1,ii), x(:,ii), &
              & t0, y0, n, na)
      end do

      do jj=2, sol%nFY
         ay(:,jj,:) = ay(:,1,:)
         by(:,jj,:) = by(:,1,:)
      end do
      
      if (na == n) then
         if ( (pb%type==ode_pb_sl).OR.(pb%type==ode_pb_sl_nl) ) then
            if (slv%type==ode_slv_ms) then
               
               do ii=1, sol%dof
                  do jj=1, sol%nFY
                     
                     by(n,jj,ii) = by(n,jj,ii) &
                          &         + ay(n,jj,ii)*y0(n)
                  end do
               end do
               
            end if
         end if
      end if

    end associate

  end subroutine homogeneous_initialcond


  subroutine ode_solver_solve(sol, slv, pb, t0, T, dt, KInv, output)
    type(ode_solution) , intent(inout) :: sol
    type(ode_problem)  , intent(in)    :: pb
    type(ode_solver)   , intent(inout) :: slv
    real(RP)           , intent(in)    :: dt, t0, T

    procedure(linsystem_solver), optional :: KInv       
    type(ode_output)           , optional :: output

    procedure(ode_output_proc), pointer :: out => null()
    logical  :: bool
    real(RP) :: time
    real(DP) :: cpu

    if (slv%verb>0) write(*,*) &
         &"ode_solver_mod  : solve    dt    =", dt

    bool = valid(pb)  
    if (.NOT.bool) call quit(&
         &"ode_solver_mod: ode_solver_solve: arg 'pb' not valid")

    bool = valid(slv)  
    if (.NOT.bool) call quit(&
         &"ode_solver_mod: ode_solver_solve: arg 'slv' not valid")

    !! output defined with 'out'
    !!
    if (present(output)) then

       !! check 
       if(    .NOT.output%assembled         &
            & .OR. (output%nbDof /= pb%dof) &
            & .OR. (output%N     /= pb%N  ) ) then
          call quit(&
               & "ode_solver_mod: ode_solver_solve: arg 'out' not valid")
       end if

       out => local_output

    else
       out => slv%output

    end if



    cpu = clock()
    sol%ierr = 0

    select case(pb%type)

    case(ode_pb_lin)
       select case(slv%type)
       case(ode_slv_1s)

          call solve_ode_lin_1s(sol, pb, t0, t, dt, &
               & slv%L_meth, out, kinv, slv%kry)

       end select

    case(ode_pb_nl)
       select case(slv%type)
       case(ode_slv_1s)

          call solve_ode_nl_1s(sol, pb, t0, t, dt, &
               &               slv%NL_meth, out, &
               &               slv%check_overflow)

       case(ode_slv_ms)

          call solve_ode_nl_ms(sol, pb, t0, t, dt, &
               &               slv%NL_meth, out, &
               &               slv%check_overflow)

       end select

    case(ode_pb_sl)
       select case(slv%type)
       case(ode_slv_ms)

          call solve_ode_sl_ms(sol, pb, t0, t, dt, &
               &               slv%SL_meth, out, &
               &               slv%check_overflow, kinv, slv%kry)

       case(ode_slv_os)
          call solve_ode_opsplt(sol, t0, t, dt, &
               &                slv%os, pb, slv%kry, out)
       end select

    case(ode_pb_sl_nl)

       select case(slv%type)
       case(ode_slv_ms)

          call solve_ode_sl_nl_ms(sol, pb, t0, t, dt, &
               & slv%SL_meth, slv%NL_meth, out, &
               & slv%check_overflow, kinv, slv%kry)

       case(ode_slv_os)
          call solve_ode_opsplt(sol, t0, t, dt, &
               &slv%os, pb, slv%kry, out)

       case(ode_slv_dc)

          call solve_ode_sl_nl_dc(sol, pb, t0, t, dt, &
               & slv%DC_meth, out, &
               & slv%check_overflow, kinv, slv%kry)
       end select

    end select

    if (slv%verb>1) then
       cpu = clock() - cpu
       write(*,*)"ode_solver_mod  : solve end, CPU =",&
         & real(cpu, sp)
    end if

    if (sol%ierr/=0) then
       if (slv%verb>0) write(*,*) &
            &"ode_solver_mod  :&
            & solve          = ERROR DETECTED, IERR = ", sol%ierr

    else
       !! external post treatment 
       !! special output at computation end
       time = t
       bool = .true.
       call out(time, sol, bool)
    end if

  contains

    subroutine local_output(tn, s, stop)
      real(RP)          , intent(in)    :: tn
      type(ode_solution), intent(in)    :: s
      logical           , intent(inout) :: stop
      
      call output_proc(output, tn, s, stop)
      
    end subroutine local_output
    
  end subroutine ode_solver_solve



  subroutine set_ode_solver_output(slv, output)
    type(ode_solver), intent(inout) :: slv
    procedure(ode_output_proc)      :: output

    slv%output => output

  end subroutine set_ode_solver_output
  
end module ode_solver_mod