ode_SL_NL_DC_mod.f90

Go to the documentation of this file.

module ode_sl_nl_dc_mod

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

  implicit none
  private

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

  public :: solve_ode_sl_nl_dc
  public :: create_ode_sl_nl_dc_sol
  public :: check_ode_method_sl_nl_dc

contains


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

    b = .false.

    select case(method)   
    case(ode_dc_2, ode_dc_3)
       b = .true.
    end select

  end function check_ode_method_sl_nl_dc
  

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

    select case(method)
    
    case(ode_dc_2)
       n_y  =  5; n_fy = 2; n_v = 1

    case(ode_dc_3)
       n_y  =  9; n_fy = 4; n_v = 1  

    case default
       call quit("ode_SL_NL_DC_mod: memSize_ode_SL_NL_DC: uncorrect method")
       
    end select

  end subroutine memsize_ode_sl_nl_dc

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

    integer :: n_Y, n_FY, n_V

    if (pb%Na/=0) call quit("ode_SL_NL_DC_mod: create_ode_SL_NL_DC:&
         & Na muste be ezero")

    call memsize_ode_sl_nl_dc(n_y, n_fy, n_v, method)
    sol = ode_solution(pb, ny=n_y, nfy=n_fy, nv=n_v)

  end subroutine create_ode_sl_nl_dc_sol



  subroutine init_nl_dc(sol, pb, meth, t0, dt, KInv)
    type(ode_solution), intent(inout) :: sol
    type(ode_problem) , intent(in)    :: pb
    integer           , intent(in)    :: meth
    real(RP)          , intent(in)    :: t0, dt
    procedure(linsystem_solver)       :: Kinv

    select case(meth)
    case(ode_dc_2)
       call init_dc_2(sol, pb, t0, dt, kinv)

    case(ode_dc_3)
       call init_dc_3(sol, pb, t0, dt, kinv)

    case default
       call quit("ode_NL_DC_mod: init_NL_DC: wrong method")

    end select

  end subroutine init_nl_dc


  subroutine init_dc_2(sol, pb, t0, dt, KInv)
    type(ode_solution), intent(inout) :: sol
    type(ode_problem) , intent(in)    :: pb
    real(RP)          , intent(in)    :: t0, dt
    procedure(linsystem_solver)       :: Kinv

    real(RP), dimension(0) :: a

    integer :: ii, N
    logical :: ierr

    n  = pb%N

    associate( y=>sol%Y, f=>sol%BY, x=>pb%X, &
         & v=>sol%V, aux=>sol%aux, rhs=>sol%rhs)

      !! Y(:,2,:) = Y_0^0 = Y^0
      y(:,2,:) = y(:,1,:)      
      
      !! Y(:,4,:) = Y_1^0 = 0
      y(:,4,:) = 0.0_rp

      !! F(:,1,:) = F_0^0
      do ii=1, sol%dof
         call pb%AB(a, f(:,1,ii), x(:,ii), t0, y(:,1,ii), n, 0 )
      end do

      !! V = V^0 (initial guess to compute V^1)
      v(:,1) = y(n,1,:)

      !! Y(:,3,:) = Y_0^{1}
      y(:,3,:) = y(:,1,:) + dt*f(:,1,:)

      !! V = V_0^{n+2}
      aux = y(n,3,:)
      call pb%M(rhs, aux)
      call kinv(v(:,1), ierr, rhs)

      !! update  Y(N,yi(2),ii) = V_0^{n+2}
      y(n,3,:) = v(:,1)

      !! restore V(:,1) = V^0
      v(:,1) = y(n,1,:)
      
    end associate

  end subroutine init_dc_2


  subroutine init_dc_3(sol, pb, t0, dt, KInv)
    type(ode_solution), intent(inout) :: sol
    type(ode_problem) , intent(in)    :: pb
    real(RP)          , intent(in)    :: t0, dt
    procedure(linsystem_solver)       :: Kinv

    real(RP), dimension(0) :: a

    integer :: ii, N
    logical :: ierr

    n  = pb%N

    associate( y=>sol%Y, f=>sol%BY, x=>pb%X, &
         & v=>sol%V, aux=>sol%aux, rhs=>sol%rhs)

      !! Y(:,2,:) = Y_0^0 = Y^0
      !! Y(:,5,:) = Y_1^0 = 0
      !! Y(:,7,:) = Y_2^0 = 0
      y(:,2,:) = y(:,1,:)            
      y(:,5,:) = 0.0_rp
      y(:,7,:) = 0.0_rp

      !! F(:,1,:) = F_0^0
      do ii=1, sol%dof
         call pb%AB(a, f(:,1,ii), x(:,ii), t0, y(:,1,ii), n, 0 )
      end do

      !! Y(:,3,:) = Y_0^{1}
      y(:,3,:) = y(:,1,:) + dt*f(:,1,:)
      aux = y(n,3,:)
      call pb%M(rhs, aux)
      v(:,1)   = y(n,1,:)         ! initial guess
      call kinv(v(:,1), ierr, rhs)
      y(n,3,:) = v(:,1)

      !! F(:,2,:) = F_0^1
      do ii=1, sol%dof
         call pb%AB(a, f(:,2,ii), x(:,ii), t0+dt, y(:,3,ii), n, 0 )
      end do

      !! Y(:,4,:) = Y_0^{2}
      y(:,4,:) = y(:,3,:) + dt*f(:,2,:)
      aux = y(n,4,:)
      call pb%M(rhs, aux)
      v(:,1)   = y(n,3,:)         ! initial guess
      call kinv(v(:,1), ierr, rhs)
      y(n,4,:) = v(:,1)

      !! F(:,4,:) = F_1^1 = F_0^1
      f(:,4,:) = f(:,2,:)

      !! Y(:,6,:) = Y_1^{1}
      y(:,6,:) = (y(:,4,:) - 2._rp*y(:,3,:) + y(:,2,:))*(-0.5_rp/dt) &
           &   +  f(:,4,:) - f(:,1,:) 
      aux = y(n,6,:)
      call pb%M(rhs, aux)
      v(:,1)   = 0._rp             ! initial guess
      call kinv(v(:,1), ierr, rhs)
      y(n,6,:) = v(:,1)

      !! restore V(:,1) = V^0
      v(:,1) = y(n,1,:)
      
    end associate

  end subroutine init_dc_3


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

    procedure(linsystem_solver), pointer :: KInv_1
    type(krylov) :: kry_1
    real(RP)     :: tn, CS
    integer      :: nS, jj, ii
    logical      :: exit_comp, bool

    if (choral_verb>1) write(*,*) &
         &"ode_SL_NL_DC_mod: solve"
    if (choral_verb>2) write(*,*) &
         &"  Method                         = ",&
         & name_ode_method(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

    bool = check_ode_method_sl_nl_dc(meth)
    if (.NOT.bool) call quit("ode_SL_NL_DC_mod: solve_ode_SL_NL_DC:&
         & wrong method")

    if (pb%Na/=0) call quit("ode_SL_NL_DC_mod: solve_ode_SL_NL_DC:&
         & Na muste be ezero")

    !! Set CS to define K = M + CS*S
    !!
    cs = s_prefactor(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

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

    !! Initialise
    !!
    call init_nl_dc(sol, pb, meth, t0, dt, kinv_1)

    !! ode resolution
    !!
    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_DC_mod: solve&
               & time     = ", tn, ': EXIT_COMPp'
          return
       end if

       select case(meth)
       case(ode_dc_2)
          call sl_nl_dc_2(sol, dt, tn, pb, kinv_1)

       case(ode_dc_3)
          call sl_nl_dc_3(sol, dt, tn, pb, kinv_1)

       end select

       !! exit in case of linear system inversion error
       if (sol%ierr/=0) then
          return
       end if

       !! check overflow
       !! WARNING : for the monodoain model
       if (check_overflow) then
          
          bool = ( minval(sol%V) < -200.0_rp ) .OR. &
               & ( maxval(sol%V) >  200.0_rp ) 
          if (bool) then
             sol%ierr = 10
             return
          end if
          
       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_sl_nl_dc


  subroutine sl_nl_dc_2(sol, dt, t, pb, KInv)
    type(ode_solution), intent(inout) :: sol
    real(RP)          , intent(in)    :: dt, t
    type(ode_problem) , intent(in)    :: pb    
    procedure(linsystem_solver)       :: Kinv

    integer                   :: ii, N
    real(RP), dimension(pb%N) :: F1, yy
    real(RP), dimension(0)    :: a
    logical                   :: ierr
    real(RP)                  :: dt_inv, t1

    dt_inv = 1.0_rp / dt
    t1     = t + dt

    n  = pb%N

    associate( y=>sol%Y, f=>sol%BY, x=>pb%X, &
         & v=>sol%V, aux=>sol%aux, rhs=>sol%rhs, &
         & yi=>sol%Y_i, fi=>sol%F_i)

      !$OMP PARALLEL PRIVATE(a)
      !$OMP DO
      do ii=1, sol%dof
         
         !! Y(:,5,ii) = d Y_0^{n+1}
         y(:,5,ii) = (y(:,yi(3),ii)-y(:,yi(2),ii)) * dt_inv

         !! V = V_0^{n+1} (initial guess to compute V_0^{n+2})
         v(ii,1) = y(n,yi(3),ii)

         !! F(:,fi(2),ii) = F_0^{n+1}
         call pb%AB(a,f(:,fi(2),ii),x(:,ii),t1,y(:,yi(3),ii),n,0)

         !! Y(:,yi(2),ii) = Y_0^{n+2}
         y(:,yi(2),ii) = y(:,yi(3),ii) + dt*f(:,fi(2),ii)

         !! aux = [ Y_0^{n+2} ]_N
         aux(ii) = y(n,yi(2),ii)
         
      end do
      !$OMP END DO
      !$OMP END PARALLEL

      !! V = V_0^{n+2}
      call pb%M(rhs, aux)
      call kinv(v(:,1), ierr, rhs)
      if (ierr) then
         sol%ierr = 1
         return
      end if


      !$OMP PARALLEL PRIVATE(a, yy, F1)
      !$OMP DO
      do ii=1, sol%dof
         
         !! update  Y(N,yi(2),ii) = V_0^{n+2}
         y(n,yi(2),ii) = v(ii,1)

         !! V = V_1^{n} (initial guess to compute V_1^{n+1})
         v(ii,1) = y(n,4,ii)

         !! F1 = F_1^{n+1}
         yy = y(:,yi(3),ii) + dt*y(:,4,ii)
         call pb%AB(a, f1, x(:,ii), t1, yy, n, 0 )

         !! yy = - d^2 Y_0^{n+2} * dt / 2
         yy = y(:,5,ii) - ( y(:,yi(2),ii) - y(:,yi(3),ii) ) * dt_inv
         yy = yy * 0.5_rp

         !! Y(:,4,ii) = Y_1^{n+1}
         y(:,4,ii) = y(:,4,ii) + yy + f1 - f(:,fi(1),ii)

         !! aux = [ Y_1^{n+1} ]_N
         aux(ii) = y(n,4,ii)

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

      !! V = V_1^{n+1}
      call pb%M(rhs, aux)
      call kinv(v(:,1), ierr, rhs)
      if (ierr) then
         sol%ierr = 1
         return
      end if

      !$OMP PARALLEL 
      !$OMP DO
      do ii=1, sol%dof
         
         !! update  Y(N,4,ii) = V_1^{n+1}
         y(n,4,ii) = v(ii,1)

         !! Compute Y_{n+1}
         y(:,1,ii) = y(:,yi(3),ii) + dt * y(:,4,ii)

         !! Update V = V_{n+1} = [ Y_{n+1} ]_N
         v(ii,1) = y(n,1,ii)

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

      !! index permutation for sol%Y_i and sol%FY_i
      call perm2(yi(2:3))
      call perm2(fi)

    end associate

  end subroutine sl_nl_dc_2





  subroutine sl_nl_dc_3(sol, dt, t, pb, KInv)
    type(ode_solution), intent(inout) :: sol
    real(RP)          , intent(in)    :: dt, t
    type(ode_problem) , intent(in)    :: pb    
    procedure(linsystem_solver)       :: Kinv

    integer                   :: ii, N
    real(RP), dimension(pb%N) :: yy, F2
    real(RP), dimension(0)    :: a
    logical                   :: ierr
    real(RP)                  :: t1, t2 
    real(RP)                  :: dt_inv, dt2_inv_6, dt_inv_2

    dt_inv    = 1.0_rp / dt
    dt_inv_2  = 1.0_rp / (dt    * 2._rp)
    dt2_inv_6 = 1.0_rp / (dt**2 * 6._rp)

    t1 = t  + dt
    t2 = t1 + dt

    n  = pb%N

    associate( y=>sol%Y, f=>sol%BY, x=>pb%X, &
         & v=>sol%V, aux=>sol%aux, rhs=>sol%rhs, &
         & yi=>sol%Y_i, fi=>sol%F_i)

      !$OMP PARALLEL PRIVATE(a)
      !$OMP DO
      do ii=1, sol%dof
         
         !! Y(:,8,ii) = d^2 Y_0^{n+2} * dt**2
         y(:,8,ii) = y(:,yi(4),ii) - 2._rp*y(:,yi(3),ii) &
              &    + y(:,yi(2),ii)

         !! F(:,fi(1),ii) = F_0^{n+2}
         call pb%AB(a,f(:,fi(1),ii),x(:,ii),t2,y(:,yi(4),ii),n,0)

         !! Y(:,yi(2),ii) = Y_0^{n+3}
         y(:,yi(2),ii) = y(:,yi(4),ii) + dt*f(:,fi(1),ii)

         !! aux = [ Y_0^{n+3} ]_N
         aux(ii) = y(n,yi(2),ii)
         
         !! V = V_0^{n+2} (initial guess to compute V_0^{n+3})
         v(ii,1) = y(n,yi(4),ii)

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

      !! V = V_0^{n+3}
      call pb%M(rhs, aux)
      call kinv(v(:,1), ierr, rhs)
      if (ierr) then
         sol%ierr = 1
         return
      end if

      !$OMP PARALLEL PRIVATE(a, yy)
      !$OMP DO
      do ii=1, sol%dof
         
         !! update  Y(N,yi(2),ii) = V_0^{n+3}
         y(n,yi(2),ii) = v(ii,1)

         !! Y(:,9,ii) = Y_1^{n}
         y(:,9,ii) = y(:,yi(5),ii) 

         !! F(:,fi(3),ii) = F_1^{n+2}
         yy = y(:,yi(4),ii) + dt*y(:,yi(6),ii)
         call pb%AB(a,f(:,fi(3),ii),x(:,ii),t2,yy,n,0)
         
         !! yy = d^2 Y_0^{n+3} * dt**2
         yy = y(:,yi(2),ii) - 2._rp*y(:,yi(4),ii) &
              &    +  y(:,yi(3),ii) 

         !! Y(:,8,ii) = d^3 Y_0^{n+3} * dt * 1/6
         y(:,8,ii) = ( yy - y(:,8,ii) ) * dt2_inv_6

         !! yy = -d^2 Y_0^{n+3} * dt/2 * 1/2
         yy = - yy * dt_inv_2

         !! Y(:,yi(5),ii) = Y_1^{n+2}
         y(:,yi(5),ii) = y(:,yi(6),ii) + yy &
              & + f(:,fi(3),ii) - f(:,fi(2),ii)

         !! V = V_1^{n+1} (initial guess to compute V_1^{n+2})
         v(ii,1) = y(n,yi(6),ii)

         !! aux = [ Y_1^{n+2} ]_N
         aux(ii) = y(n,yi(5),ii)

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

      !! V = V_1^{n+2}
      call pb%M(rhs, aux)
      call kinv(v(:,1), ierr, rhs)
      if (ierr) then
         sol%ierr = 1
         return
      end if

      !$OMP PARALLEL PRIVATE(a, yy, F2)
      !$OMP DO
      do ii=1, sol%dof
         
         !! update  Y(N,yi(5),ii) = V_1^{n+2}
         y(n,yi(5),ii) = v(ii,1)

         !! V = V_2^{n} (initial guess to compute V_2^{n+1})
         v(ii,1) = y(n,7,ii)

         !! F2 = F_2^{n+1}
         yy = y(:,yi(3),ii) + dt*y(:,yi(6),ii) + dt**2*y(:,7,ii)
         call pb%AB(a,f2,x(:,ii),t1,yy,n,0)

         !! yy = - d^2 Y_1{n+2}  * dt * 1/2
         yy = -y(:,9,ii) + 2._rp*y(:,yi(6),ii) - y(:,yi(5),ii) 
         yy = yy * dt_inv_2

         !! Y(:,7,:) = Y_2^{n+1}         
         yy = yy + ( f2 - f(:,fi(4),ii) ) * dt_inv 
         y(:,7,ii) = y(:,7,ii) + y(:,8,ii) + yy   

         !! aux = [ Y_2^{n+1} ]_N
         aux(ii) = y(n,7,ii)

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

      !! V = V_2^{n+1}
      call pb%M(rhs, aux)
      call kinv(v(:,1), ierr, rhs)
      if (ierr) then
         sol%ierr = 1
         return
      end if

      !$OMP PARALLEL 
      !$OMP DO
      do ii=1, sol%dof
         
         !! update  Y(N,7,ii) = V_2^{n+1}
         y(n,7,ii) = v(ii,1)

         !! Compute Y_{n+1}
         y(:,1,ii) = y(:,yi(3),ii) + dt    * y(:,yi(6),ii) &
              &                    + dt**2 * y(:,7,ii)    

         !! Update V = V_{n+1} = [ Y_{n+1} ]_N
         v(ii,1) = y(n,1,ii)

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


      !! index permutation for sol%Y_i and sol%FY_i
      call perm3(yi(2:4))
      call perm2(yi(5:6))
      call perm2(fi(1:2))
      call perm2(fi(3:4))

    end associate

  end subroutine sl_nl_dc_3


  subroutine perm2(T)
    integer, dimension(2), intent(inout) :: T

    integer :: t1

    t1 = t(1)
    t(1) = t(2)
    t(2) = t1

  end subroutine perm2


  subroutine perm3(T)
    integer, dimension(3), intent(inout) :: T

    integer :: t1

    t1 = t(1)
    t(1) = t(2)
    t(2) = t(3)
    t(3) = t1

  end subroutine perm3

end module ode_sl_nl_dc_mod