ode_Lin_1s_mod.f90

Go to the documentation of this file.

module ode_lin_1s_mod

  use choral_constants
  use choral_variables
  use real_type
  use basic_tools
  use r1d_mod, only: xpay, scale, axpby
  use ode_def
  use ode_problem_mod
  use ode_solution_mod
  use krylov_mod

  implicit none
  private

  public :: solve_ode_lin_1s
  public :: create_ode_lin_1s_sol
  public :: set_solver_ode_lin_1s
  public :: memsize_ode_lin_1s
  public :: check_ode_method_lin_1s

contains

  function check_ode_method_lin_1s(method) result(b)
    logical :: b
    integer, intent(in) :: method

    b = .false.

    select case(method)   
    case(ode_be, ode_cn, ode_sdirk4)
       b = .true.
       
    end select

  end function check_ode_method_lin_1s


  subroutine memsize_ode_lin_1s(n_V, method)
    integer, intent(out) :: n_V
    integer, intent(in)  :: method 

    select case(method)
    
    case(ode_be, ode_cn)
       n_v  =  1

    case(ode_sdirk4)
       n_v  =  4

    case default
       call quit("ode_Lin_1s_mod: memSize_ode_Lin_1s;&
            & uncorrect method")
       
    end select

  end subroutine memsize_ode_lin_1s

  subroutine set_solver_ode_lin_1s(slv, method) 
    procedure(ode_lin_solver), pointer    :: slv
    integer                  , intent(in) :: method

    select case(method)
    case(ode_be)
       slv => lin_1s_be
    case(ode_cn)
       slv => lin_1s_cn
    case(ode_sdirk4)
       slv => lin_1s_sdirk4
       
    case default
       call quit("ode_Lin_1s_mod: set_solver_ode_Lin_1s;&
            & uncorrect method")
    end select
    
  end subroutine set_solver_ode_lin_1s


  subroutine create_ode_lin_1s_sol(sol, pb, method)
    type(ode_solution), intent(inout) :: sol
    type(ode_problem) , intent(in)    :: pb
    integer           , intent(in)    :: method

    integer :: n_V
    logical :: bool

    bool = check_ode_method_lin_1s(method)
    if (.NOT.bool) call quit(&
         & "ode_Lin_1s_mod: create_ode_Lin_1s_sol: uncorrect method")
    call memsize_ode_lin_1s(n_v, method)

    sol = ode_solution(pb, nv=n_v)

  end subroutine create_ode_lin_1s_sol
  

  subroutine solve_ode_lin_1s(sol, pb, t0, T, dt, method, &
       & out, Kinv, kry)
    type(ode_solution) , intent(inout)    :: sol
    type(ode_problem)  , intent(in)       :: pb
    real(RP)           , intent(in)       :: t0, T, dt
    integer            , intent(in)       :: method
    procedure(ode_output_proc)            :: out
    procedure(linsystem_solver), optional :: Kinv
    type(krylov)               , optional :: kry

    procedure(ode_lin_solver)  , pointer :: slv=>null()
    procedure(linsystem_solver), pointer :: KInv_1
    type(krylov) :: kry_1
    logical  :: ierr, exit_comp
    integer  :: ns, jj
    real(RP) :: tn, CS


    if (choral_verb>1) write(*,*) &
         & "ode_Lin_1s_mod  : solve_ode_Lin_1s"
    if (choral_verb>2) write(*,*) &
         & "  Linear  onestep solver         = ",&
         & name_ode_method(method)
    if (choral_verb>2) write(*,*) &
         & "  K_inv provided                 =",&
         & present(kinv)
    if (.NOT.present(kinv)) then
       if (choral_verb>2) write(*,*) &
            &"  Krylov settings provided       = ",&
            & present(kry)
    end if

    !! Set CS to define K = M + CS*S
    !!
    cs = s_prefactor(method, dt)

    !! set the linear system solver for the system K*x = rhs
    !!
    if (present(kinv))  then
       kinv_1 => kinv
    else
       kinv_1 => kinv_default
       if (present(kry)) kry_1 = kry
    end if

    !! set the ode solver depending on the method
    !!
    call set_solver_ode_lin_1s(slv, method)

    !! ode resolution
    !!
    ns = int( (t-t0) / dt) 
    sol%ierr = 0 
    exit_comp = .false. 
    do jj=1, ns
       tn = t0 + re(jj-1)*dt

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

       call slv(sol, ierr, dt, pb, kinv_1)

       if (ierr) then
          sol%ierr = 1
          return
       end if
    end do

  contains

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

      call solve(x, kry_1, b, k_default)
      bool = kry_1%ierr
      
    end subroutine kinv_default

    !! matrix-vector product x --> K*x, K = M + CS*S
    !!
    subroutine k_default(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_default

  end subroutine solve_ode_lin_1s


  subroutine lin_1s_be(sol, ierr, dt, pb, KInv)
    type(ode_solution) , intent(inout) :: sol
    logical            , intent(out)   :: ierr  
    real(RP)           , intent(in)    :: dt
    type(ode_problem)  , intent(in)    :: pb
    procedure(linsystem_solver)        :: KInv

    call pb%M(sol%rhs, sol%V(:,1))
    call kinv(sol%V(:,1), ierr, sol%rhs)

  end subroutine lin_1s_be


  subroutine lin_1s_cn(sol, ierr, dt, pb, KInv)
    type(ode_solution) , intent(inout) :: sol
    logical            , intent(out)   :: ierr  
    real(RP)           , intent(in)    :: dt
    type(ode_problem)  , intent(in)    :: pb
    procedure(linsystem_solver)        :: KInv

    call pb%M(sol%rhs, sol%V(:,1))
    call pb%S(sol%aux, sol%V(:,1))
    call xpay(sol%rhs, -dt/2._rp, sol%aux)

    call kinv(sol%V(:,1), ierr, sol%rhs)

  end subroutine lin_1s_cn


  subroutine lin_1s_sdirk4(sol, ierr, dt, pb, KInv)
    type(ode_solution) , intent(inout) :: sol
    logical            , intent(out)   :: ierr  
    real(RP)           , intent(in)    :: dt
    type(ode_problem)  , intent(in)    :: pb
    procedure(linsystem_solver)        :: KInv

    real(RP) :: gamma, a

    gamma = 1._rp/sqrt(3._rp) * cos(pi/18._rp) 
    gamma = gamma + 0.5_rp

    ! sol%V(:,2) := k1
    call scale(sol%aux, -1._rp, sol%V(:,1))
    call pb%S(sol%rhs, sol%aux)
    call kinv(sol%V(:,2), ierr, sol%rhs)
    if (ierr) return
    
    ! sol%V(:,3) := k2
    a = dt * (0.5_rp - gamma)
    call axpby(sol%aux, -1._rp, sol%V(:,1), -a, sol%V(:,2))
    call pb%S(sol%rhs, sol%aux)
    call kinv(sol%V(:,3), ierr, sol%rhs)
    if (ierr) return

    ! sol%V(:,4) := k3
    a = dt * 2._rp * gamma
    call axpby(sol%aux, -1._rp, sol%V(:,1), -a, sol%V(:,2))
    a = dt * (1._rp - 4._rp*gamma)
    call xpay(sol%aux, -a, sol%V(:,3))
    call pb%S(sol%rhs, sol%aux)
    call kinv(sol%V(:,4), ierr, sol%rhs)
    if (ierr) return

    a = 6._rp * ( 2._rp*gamma - 1._rp )**2
    a = dt / a
    call xpay(sol%V(:,1), a, sol%V(:,2) )
    call xpay(sol%V(:,1), a, sol%V(:,4) )

    a = dt - 2._rp*a
    call xpay(sol%V(:,1), a, sol%V(:,3) )
    
  end subroutine lin_1s_sdirk4

  
end module ode_lin_1s_mod