ode_SL_ms_mod.f90

Go to the documentation of this file.

module ode_sl_ms_mod

  !$ use OMP_LIB
  use choral_constants
  use choral_variables
  use real_type
  use algebra_set
  use basic_tools
  use r1d_mod, only: xpay, copy, axpby
  use ode_def
  use ode_problem_mod
  use ode_solution_mod
  use krylov_mod

  implicit none
  private

  public :: solve_ode_sl_ms
  public :: create_ode_sl_ms_sol
  public :: set_solver_ode_sl_ms
  public :: memsize_ode_sl_ms
  public :: check_ode_method_sl_ms


  real(RP), parameter :: c4_3    =   4._rp /  3._rp
  real(RP), parameter :: c2_3    =   2._rp /  3._rp
  real(RP), parameter :: c_3     =   1._rp /  3._rp
  real(RP), parameter :: c2_11   =   2._rp / 11._rp
  real(RP), parameter :: c6_11   =   6._rp / 11._rp
  real(RP), parameter :: c9_11   =   9._rp / 11._rp
  real(RP), parameter :: c18_11  =  18._rp / 11._rp
  real(RP), parameter :: c3_25   =   3._rp / 25._rp
  real(RP), parameter :: c12_25  =  12._rp / 25._rp
  real(RP), parameter :: c16_25  =  16._rp / 25._rp
  real(RP), parameter :: c36_25  =  36._rp / 25._rp
  real(RP), parameter :: c48_25  =  48._rp / 25._rp
  real(RP), parameter :: c72_25  =  72._rp / 25._rp
  real(RP), parameter :: c12_137 =  12._rp /137._rp
  real(RP), parameter :: c60_137 =  60._rp /137._rp
  real(RP), parameter :: c75_137 =  75._rp /137._rp
  real(RP), parameter :: c200_137= 200._rp /137._rp
  real(RP), parameter :: c300_137= 300._rp /137._rp
  real(RP), parameter :: c600_137= 600._rp /137._rp

contains

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

    b = .false.

    select case(method)
    case(ode_fbe, ode_fe, ode_cnab2, ode_mcnab2, &
         & ode_bdfsbdf2, ode_bdfsbdf3, ode_bdfsbdf4, ode_bdfsbdf5)
   
       b = .true.       
    end select

  end function check_ode_method_sl_ms


  subroutine memsize_ode_sl_ms(n_V, n_FY, method)
    integer, intent(out) :: n_V, n_FY
    integer, intent(in)  :: method

    select case(method)
    
    case(ode_fbe, ode_fe)
       n_v = 1 ; n_fy = 1

    case(ode_cnab2)
       n_v = 1 ; n_fy = 2

    case(ode_mcnab2, ode_bdfsbdf2)
       n_v = 2 ; n_fy = 2

    case(ode_bdfsbdf3)
       n_v = 3 ; n_fy = 3

    case(ode_bdfsbdf4)
       n_v = 4 ; n_fy = 4

    case(ode_bdfsbdf5)
       n_v = 5 ; n_fy = 5

    case default
       call quit( "ode_SL_ms_mod: memSize_ode_SL_ms: error")
       
    end select

  end subroutine memsize_ode_sl_ms
  

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

    select case(method)
    case(ode_fbe)
       slv => sl_ms_fbe
    case(ode_fe)
       slv => sl_ms_fe
    case(ode_bdfsbdf2)
       slv => sl_ms_bdfsbdf2 
    case(ode_bdfsbdf3)
       slv =>  sl_ms_bdfsbdf3
    case(ode_bdfsbdf4)
       slv => sl_ms_bdfsbdf4
    case(ode_bdfsbdf5)
       slv =>  sl_ms_bdfsbdf5
    case(ode_cnab2)
       slv => sl_ms_cnab2
    case(ode_mcnab2)
       slv => sl_ms_mcnab2

    case default
       call quit( "ode_SL_ms_mod:&
    & set_solver_ode_SL_ms: uncorrect method")
    end select
    
  end subroutine set_solver_ode_sl_ms


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

    integer :: n_V, n_FY
    logical :: bool

    bool = check_ode_method_sl_ms(method)
    if (.NOT.bool) call quit(&
         & "ode_SL_ms_mod: create_ode_SL_ms_sol: uncorrect method")

    if (pb%N/=1) call quit(&
         & "ode_SL_ms_mod: create_ode_SL_ms_sol: wrong problem size")

    call memsize_ode_sl_ms(n_v, n_fy, method)
    sol = ode_solution(pb, nv=n_v, nfy=n_fy, ny=1)

  end subroutine create_ode_sl_ms_sol


  subroutine solve_ode_sl_ms(sol, pb, t0, T, dt, method, &
       & out, check_overflow, 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
    logical            , intent(in)       :: check_overflow
    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
    real(RP)     :: tn, CS
    real(DP)     :: t_SL, t_out, t_reac, t_tot, cpu
    logical      :: ierr, exit_comp
    integer      :: ns, ii, jj, n0, n0_F, n0_V
    real(RP), dimension(0) :: ay0

    if (choral_verb>1) write(*,*) &
         &"ode_SL_ms_mod   : solve_ode_SL_ms"
    if (choral_verb>2) write(*,*) &
         &"  SemiLin multistep 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 elementary solver
    !!
    call set_solver_ode_sl_ms(slv, method)

    !! initialise the solution indexes
    !!
    call ode_solution_init_indexes(sol)

    !! timer initialisation
    t_sl   = 0.0_dp
    t_out  = 0.0_dp
    t_reac = .0_dp
    t_tot  = clock()

    !! ode resolution
    !!
    ns = int( (t-t0) / dt) 
    exit_comp = .false. 
    sol%ierr  = 0

    do jj=1, ns
       tn = t0 + re(jj-1)*dt

       ! output
       cpu = clock()
       call out(tn, sol, exit_comp)
       if (exit_comp) then
          if (choral_verb>2) write(*,*) &
               &"ode_SL_ms_mod   : solve&
               & time     = ", tn, ': EXIT_COMPp'
          return
       end if
       t_out = t_out + clock() - cpu

       !! perform one time-step
       !!
       cpu = clock()
       call slv(sol, ierr, dt, pb, kinv_1)
       t_sl = t_sl + clock() - cpu

       !! check linear system resolution
       !!
       if (ierr) then
          sol%ierr = 1
          if (choral_verb>1) write(*,*) &
               &"ode_SL_ms_mod   : solve,&
               &    time = ", tn, ': linear system solver error'
          return
       end if

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

         call circperm(sol%Y_i)
         call circperm(sol%F_i)
         call circperm(sol%V_i)
         
         n0   = sol%Y_i(1)
         n0_f = sol%F_i(1)
         n0_v = sol%V_i(1)

         ! 1. recast V into Y (note that pb%N = 1 here)
         !
         ! 2. compute reaction terms
         !
         if (pb%Na==0) then
            !$OMP PARALLEL PRIVATE(ay0)
            !$OMP DO
            do ii=1, sol%dof
               y(1, n0, ii) =  v(ii, n0_v)

               call pb%AB(ay0, by(:, n0_f,ii), &
                    & x(:,ii), tn+dt, y(:,n0,ii), 1, 0 )
            end do
            !$OMP END DO
            !$OMP END PARALLEL
         else
            !$OMP PARALLEL 
            !$OMP DO
            do ii=1, sol%dof
               y(1, n0, ii) =  v(ii, n0_v)

               call pb%AB(  ay(:,n0_f,ii), by(:,n0_f,ii), &
                    & x(:,ii), tn+dt, y(:,n0,ii), 1, 1 )

               by(1,n0_f, ii) = by(1,n0_f, ii) &
                    &         + ay(1,n0_f, ii)*v(ii, n0_v)
            end do
            !$OMP END DO
            !$OMP END PARALLEL
         end if

         !! Check overflow
         !!
         if (check_overflow) then
            ierr = overflow(sol%Y(:,n0,:))
            if (ierr) then
               sol%ierr = 10
               if (choral_verb>2) write(*,*) &
                    &"ode_NL_ms_mod   : solve,&
                    &    time = ", tn, ': OVERFLOW Y'
               return
            end if

            ierr = overflow(sol%BY(:,n0_f,:))
            if (ierr) then
               sol%ierr = 10
               if (choral_verb>2) write(*,*) &
                    &"ode_NL_ms_mod   : solve,&
                    &    time = ", tn, ': OVERFLOW BY'
               return
            end if

            if (pb%Na>0) then
               ierr = overflow(sol%AY(:,n0_f,:))
               if (ierr) then
                  sol%ierr = 10
                  if (choral_verb>2) write(*,*) &
                       &"ode_NL_ms_mod   : solve,&
                       &    time = ", tn, ': OVERFLOW AY'
                  return
               end if
            end if
         end if

       end associate
       t_reac = t_reac + clock() - cpu

    end do

    if (choral_verb>1) then
       write(*,*)"ode_SL_NL_mod   : solve, end     = timing"
       write(*,*)"  Semilin. eq.    integration    =", &
            & real(t_sl, sp)
       write(*,*)"  Reaction term evaluation       =", &
            & real(t_reac, sp)
       write(*,*)"  Output                         =", &
            & real(t_out, sp)
       t_tot  = clock() - t_tot
       write(*,*)"  Total CPU                      =", &
            & real(t_tot, sp)

    end if

  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_sl_ms

  
  subroutine sl_ms_fbe(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

    integer :: n0_V, nl_V, n0_F
    
    n0_v = sol%V_i(1)
    nl_v = sol%V_i(sol%nV)
    n0_f = sol%F_i(1)

    associate( by=>sol%BY, v=>sol%V, aux=>sol%aux, rhs=>sol%rhs)

      call xpay(aux, v(:, n0_v), dt, by(pb%N, n0_f, :))
      call pb%M(rhs, aux)
      
      !! initial guess
      if (n0_v /= nl_v) then
         call copy(v(:,nl_v), v(:,n0_v))
      end if
      
      !! solve the linear system
      call kinv(v(:,nl_v), ierr, rhs)
      
    end associate

  end subroutine sl_ms_fbe


  subroutine sl_ms_fe(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

    integer :: n0_V, nl_V, n0_F
    
    n0_v = sol%V_i(1)
    nl_v = sol%V_i(sol%nV)
    n0_f = sol%F_i(1)

    associate( by=>sol%BY, v=>sol%V, aux=>sol%aux, rhs=>sol%rhs)

      call xpay(aux, v(:, n0_v), dt, by(pb%N, n0_f, :))
      call pb%M(rhs, aux)
      call pb%S(aux, v(:, n0_v))
      call xpay(rhs, -dt, aux)

      
      !! initial guess
      if (n0_v /= nl_v) then
         call copy(v(:,nl_v), v(:,n0_v))
      end if

      !! solve the linear system
      call kinv(v(:,nl_v), ierr, rhs)
      
    end associate

  end subroutine sl_ms_fe

  
  subroutine sl_ms_cnab2(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

    integer  :: ii, N, n0_V, nl_V, n0_F, n1_F
    real(RP) :: h3_2, h1_2

    h3_2 = dt * 1.5_rp
    h1_2 = dt * 0.5_rp

    n    = pb%N
    n0_v = sol%V_i(1)
    nl_v = sol%V_i(sol%nV)
    n0_f = sol%F_i(1)
    n1_f = sol%F_i(2)

    associate( by=>sol%BY, v=>sol%V, aux=>sol%aux, rhs=>sol%rhs)

      !$OMP PARALLEL 
      !$OMP DO
      do ii=1, sol%dof

         aux(ii) = v(ii,n0_v) &
              &  + h3_2 * by(n,n0_f,ii) - h1_2 * by(n,n1_f,ii)

      end do
      !$OMP END DO
      !$OMP END PARALLEL

      call pb%M(rhs, aux)
      call pb%S(aux, v(:,n0_v))
      call xpay(rhs, -h1_2, aux)

      !! initial guess
      if (sol%NV>1) call copy( v(:,nl_v), v(:,n0_v))
      
      !! solve the linear system
      call kinv(v(:,nl_v), ierr, rhs)

    end associate

  end subroutine sl_ms_cnab2
    

  subroutine sl_ms_mcnab2(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

    integer  :: ii, N, n0_V, n1_V, nl_V, n0_F, n1_F
    real(RP) :: h3_2, h1_2

    h3_2 = dt * 1.5_rp
    h1_2 = dt * 0.5_rp

    n    = pb%N
    n0_v = sol%V_i(1)
    n1_v = sol%V_i(2)
    nl_v = sol%V_i(sol%nV)
    n0_f = sol%F_i(1)
    n1_f = sol%F_i(2)

    associate( by=>sol%BY, v=>sol%V, aux=>sol%aux, rhs=>sol%rhs)

      !$OMP PARALLEL 
      !$OMP DO
      do ii=1, sol%dof

         aux(ii) = v(ii,n0_v) &
              &  + h3_2 * by(n,n0_f,ii) - h1_2 * by(n,n1_f,ii)

      end do
      !$OMP END DO
      !$OMP END PARALLEL

      call pb%M(rhs, aux)
      call pb%S(aux, v(:,n0_v))
      call xpay(rhs, -re(3,8)*dt, aux)
      call pb%S(aux, v(:,n1_v))
      call xpay(rhs, -re(1,16)*dt, aux)

      !! initial guess
      call copy( v(:,nl_v), v(:,n0_v))
      
      !! solve the linear system
      call kinv(v(:,nl_v), ierr, rhs)

    end associate

  end subroutine sl_ms_mcnab2

    
    ! call copy(aux, V(:,n0_V))
    ! call add_aux_AY_p_B(1, dt*1.5_RP)
    ! call add_aux_AY_p_B(2,-dt/2._RP)

    ! call M(rhs, aux)
    ! call S(aux, V(:,n0_V))
    ! call xpay(rhs, -re(3,8)*dt, aux)
    ! call S(aux, V(:,n1_V))
    ! call xpay(rhs, -re(1,16)*dt, aux)
    
    ! !! initial guess
    ! if (size(V)>1) call copy(V(NV), V(:,n0_V))

    ! !! solve the linear system
    ! call KInv(V(NV), ierr, rhs)



  subroutine sl_ms_bdfsbdf2(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

    integer  :: ii, N, n0_V, n1_V, nl_V, n0_F, n1_F
    real(RP) :: h4_3, h2_3

    h4_3 = dt * c4_3
    h2_3 = dt * c2_3

    n    = pb%N
    n0_v = sol%V_i(1)
    n1_v = sol%V_i(2)
    nl_v = sol%V_i(sol%nV)
    n0_f = sol%F_i(1)
    n1_f = sol%F_i(2)

    associate( by=>sol%BY, v=>sol%V, aux=>sol%aux, rhs=>sol%rhs)

      !$OMP PARALLEL 
      !$OMP DO
      do ii=1, sol%dof

         aux(ii) = c4_3 * v(ii,n0_v) - c_3 * v(ii,n1_v)   &
              &  + h4_3 * by(n,n0_f,ii) - h2_3 * by(n,n1_f,ii)

      end do
      !$OMP END DO
      !$OMP END PARALLEL

      call pb%M(rhs, aux)

      !! initial guess
      call copy( v(:,nl_v), v(:,n0_v))
      
      !! solve the linear system
      call kinv(v(:,nl_v), ierr, rhs)

    end associate

  end subroutine sl_ms_bdfsbdf2


  subroutine sl_ms_bdfsbdf3(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

    integer  :: ii, N, n0_V, n1_V, n2_V, nl_V, n0_F, n1_F, n2_F
    real(RP) :: a, h18_11, h6_11

    h18_11 = dt * c18_11
    h6_11  = dt * c6_11

    n    = pb%N
    n0_v = sol%V_i(1)
    n1_v = sol%V_i(2)
    n2_v = sol%V_i(3)
    nl_v = sol%V_i(sol%nV)
    n0_f = sol%F_i(1)
    n1_f = sol%F_i(2)
    n2_f = sol%F_i(3)

    associate( by=>sol%BY, v=>sol%V, aux=>sol%aux, rhs=>sol%rhs)

      !$OMP PARALLEL PRIVATE(a) 
      !$OMP DO
      do ii=1, sol%dof

         a =       c18_11 * v(ii,n0_v) - c9_11  * v(ii,n1_v)   &
              &  + c2_11  * v(ii,n2_v)

         a = a  + h18_11 * by(n,n0_f,ii) - h18_11 * by(n,n1_f,ii) &
              & + h6_11  * by(n,n2_f,ii)

          aux(ii) = a

      end do
      !$OMP END DO
      !$OMP END PARALLEL

      call pb%M(rhs, aux)

      !! initial guess
      call copy( v(:,nl_v), v(:,n0_v))
      
      !! solve the linear system
      call kinv(v(:,nl_v), ierr, rhs)

    end associate

  end subroutine sl_ms_bdfsbdf3



  subroutine sl_ms_bdfsbdf4(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

    integer  :: ii, N, n0_V, n1_V, n2_V, nl_V, n0_F, n1_F, n2_F
    integer  :: n3_V, n3_F
    real(RP) :: a, h12_25, h48_25, h72_25

    h12_25 = dt * c12_25
    h48_25 = dt * c48_25
    h72_25 = dt * c72_25

    n    = pb%N
    n0_v = sol%V_i(1)
    n1_v = sol%V_i(2)
    n2_v = sol%V_i(3)
    n3_v = sol%V_i(4)
    nl_v = sol%V_i(sol%nV)
    n0_f = sol%F_i(1)
    n1_f = sol%F_i(2)
    n2_f = sol%F_i(3)
    n3_f = sol%F_i(4)

    associate( by=>sol%BY, v=>sol%V, aux=>sol%aux, rhs=>sol%rhs)

      !$OMP PARALLEL PRIVATE(a) 
      !$OMP DO
      do ii=1, sol%dof

         a =       c48_25 * v(ii,n0_v) - c36_25  * v(ii,n1_v)   &
              &  + c16_25 * v(ii,n2_v) - c3_25   * v(ii,n3_v)   

         a = a  + h48_25 * by(n,n0_f,ii) - h72_25 * by(n,n1_f,ii) &
              & + h48_25 * by(n,n2_f,ii) - h12_25 * by(n,n3_f,ii)

          aux(ii) = a

      end do
      !$OMP END DO
      !$OMP END PARALLEL

      call pb%M(rhs, aux)

      !! initial guess
      call copy( v(:,nl_v), v(:,n0_v))
      
      !! solve the linear system
      call kinv(v(:,nl_v), ierr, rhs)

    end associate

  end subroutine sl_ms_bdfsbdf4
    

  subroutine sl_ms_bdfsbdf5(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

    integer  :: ii, N, n0_V, n1_V, n2_V, nl_V, n0_F, n1_F, n2_F
    integer  :: n3_V, n3_F, n4_V, n4_F
    real(RP) :: a, h60_137, h300_137, h600_137

    h60_137  = dt * c60_137
    h300_137 = dt * c300_137
    h600_137 = dt * c600_137

    n    = pb%N
    n0_v = sol%V_i(1)
    n1_v = sol%V_i(2)
    n2_v = sol%V_i(3)
    n3_v = sol%V_i(4)
    n4_v = sol%V_i(5)
    nl_v = sol%V_i(sol%nV)
    n0_f = sol%F_i(1)
    n1_f = sol%F_i(2)
    n2_f = sol%F_i(3)
    n3_f = sol%F_i(4)
    n4_f = sol%F_i(5)

    associate( by=>sol%BY, v=>sol%V, aux=>sol%aux, rhs=>sol%rhs)

      !$OMP PARALLEL PRIVATE(a) 
      !$OMP DO
      do ii=1, sol%dof

         a =       c300_137 * v(ii,n0_v) - c300_137  * v(ii,n1_v) &
              &  + c200_137 * v(ii,n2_v) - c75_137   * v(ii,n3_v) &
              &  + c12_137  * v(ii,n4_v)

         a = a  + h300_137* by(n,n0_f,ii) - h600_137* by(n,n1_f,ii) &
              & + h600_137* by(n,n2_f,ii) - h300_137* by(n,n3_f,ii) &
              & + h60_137 * by(n,n4_f,ii)

          aux(ii) = a

      end do
      !$OMP END DO
      !$OMP END PARALLEL

      call pb%M(rhs, aux)

      !! initial guess
      call copy( v(:,nl_v), v(:,n0_v))
      
      !! solve the linear system
      call kinv(v(:,nl_v), ierr, rhs)

    end associate

  end subroutine sl_ms_bdfsbdf5
    

end module ode_sl_ms_mod