ode_SL_NL_ms_mod.f90

Go to the documentation of this file.

module ode_sl_nl_ms_mod

  use choral_variables
  use real_type
  use basic_tools
  use algebra_set
  use r1d_mod, only: xpay
  use ode_def
  use ode_problem_mod
  use ode_solution_mod
  use ode_nl_ms_mod
  use krylov_mod
  use ode_sl_ms_mod
  !$ use OMP_LIB

  implicit none
  private

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

  public :: solve_ode_sl_nl_ms
  public :: create_ode_sl_nl_ms_sol

contains

  subroutine create_ode_sl_nl_ms_sol(sol, pb, SL_meth, NL_meth)
    type(ode_solution), intent(inout) :: sol
    type(ode_problem) , intent(in)    :: pb
    integer           , intent(in)    :: SL_meth, NL_meth
    
    integer :: nV, nY, nFY, nFY2

    call memsize_ode_nl_ms(ny, nfy , nl_meth)
    call memsize_ode_sl_ms(nv, nfy2, sl_meth)
    nfy = max(nfy, nfy2)
    sol = ode_solution(pb, nv=nv, ny=ny, nfy=nfy)

  end subroutine create_ode_sl_nl_ms_sol
  

  subroutine solve_ode_sl_nl_ms(sol, pb, t0, T, dt, &
       & SL_meth, NL_meth, 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)       :: SL_meth, NL_meth
    procedure(ode_output_proc)                     :: out
    logical            , intent(in)       :: check_overflow
    procedure(linsystem_solver), optional :: Kinv
    type(krylov)               , optional :: kry

    procedure(ode_nl_ms_solver), pointer :: slv_NL
    procedure(ode_lin_solver)  , pointer :: slv_SL
    procedure(linsystem_solver), pointer :: KInv_1
    type(krylov) :: kry_1
    real(RP)     :: tn, CS
    real(DP)     :: t_NL, t_SL, t_out, t_reac, cpu, t_tot
    integer      :: nS, ii, jj, N, Na, P, n0, n0_F, n0_V
    logical      :: exit_comp, ierr

    if (choral_verb>1) write(*,*) &
         &"ode_SL_NL_ms_mod: solve"
    if (choral_verb>2) write(*,*) &
         &"  SemiLin multistep solver       = ",&
         & name_ode_method(sl_meth)
    if (choral_verb>2) write(*,*) &
         &"  NonLin  multistep solver       = ",&
         & name_ode_method(nl_meth)
    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(sl_meth, 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 solvers 
    !!
    call set_solver_ode_sl_ms(slv_sl, sl_meth)
    call set_solver_ode_nl_ms(slv_nl, nl_meth)

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

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

    !! ode resolution
    !!
    n  = pb%N
    na = pb%Na
    if (na < n) then
       p = na
    else
       p = n-1
    end if
    ns = int( (t-t0) / dt) 
    exit_comp = .false. 
    sol%ierr  = 0
    !!
    !! Time loop
    !!
    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_NL_ms_mod: solve&
               & time     = ", tn, ': EXIT_COMPp'
          return
       end if
       t_out = t_out + clock() - cpu

       !! solve the non linear ODE system
       cpu = clock()
       call slv_nl(sol, dt, n, p)
       t_nl = t_nl +  clock() - cpu

       ! solve the semilinear equation
       cpu = clock()
       call slv_sl(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>2) write(*,*) &
               &"ode_SL_NL_ms_mod: solve,&
               &    time = ", tn, ': lin pb inv 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
         !
         ! 2. compute reaction terms
         !
         if (na<n) then
            !$OMP PARALLEL 
            !$OMP DO
            do ii=1, sol%dof
               y(n, n0, ii) =  v(ii, n0_v)

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

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

               by(n, n0_f, ii) = by(n, n0_f, ii) &
                    &         + ay(n, 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

            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 associate
       t_reac = t_reac + clock() - cpu

    end do

    if (choral_verb>1) then
       write(*,*)"ode_SL_NL_ms_mod: solve, end     = timing"
       write(*,*)"  Non Lin. system integration    =", &
            & real(t_nl, sp)
       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_nl_ms


end module ode_sl_nl_ms_mod