ode_NL_DC_mod.f90

Go to the documentation of this file.

module ode_nl_dc_mod

  use choral_variables
  use choral_constants
  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 krylov_mod
  !$ use OMP_LIB

  implicit none
  private

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

  public :: solve_ode_nl_dc
  public :: create_ode_nl_dc_sol
  public :: check_ode_method_nl_dc

contains


  function check_ode_method_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_nl_dc

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

    select case(method)
    
    case(ode_dc_2)
       n_y  =  3; n_fy = 1  

    case(ode_dc_3)
       n_y  =  6; n_fy = 3  

    case default
       call quit("ode_NL_DC_mod: memSize_ode_NL_DC: uncorrect method")
       
    end select

  end subroutine memsize_ode_nl_dc


  subroutine create_ode_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

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

    call memsize_ode_nl_dc(n_y, n_fy, method)
    sol = ode_solution(pb, ny=n_y, nfy=n_fy)

  end subroutine create_ode_nl_dc_sol
  

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

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

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

    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)
    type(ode_solution), intent(inout) :: sol
    type(ode_problem) , intent(in)    :: pb
    real(RP)          , intent(in)    :: t0

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

    integer :: ii, N

    n  = pb%N

    associate( y=>sol%Y, by=>sol%BY, x=>pb%X)

      !$OMP PARALLEL PRIVATE(ay)
      !$OMP DO
      do ii=1, sol%dof
         y(:,2,ii) = y(:,1,ii)      
         y(:,3,ii) = 0.0_rp
         call pb%AB(  ay, by(:,1,ii), &
              & x(:,ii), t0, y(:,1,ii), n, 0 )
      end do
      !$OMP END DO
      !$OMP END PARALLEL
      
    end associate
  end subroutine init_dc_2


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

    integer                   :: ii, N
    real(RP), dimension(0)    :: ay

    n  = pb%N

    associate( y=>sol%Y, fy=>sol%BY, x=>pb%X)

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

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

         ! F_0^0
         call pb%AB(  ay, fy(:,1,ii), &
              & x(:,ii), t0, y(:,2,ii), n, 0 )

         ! Y_0^1
         y(:,5,ii) = y(:,2,ii) + dt*fy(:,1,ii)

         ! F_0^1
         call pb%AB(  ay, fy(:,2,ii), &
              & x(:,ii), t0+dt, y(:,5,ii), n, 0 )

         ! F_1^1
         fy(:,3,ii) = fy(:,2,ii)

         ! Y_1^1
         y(:,6,ii) = 0.5_rp*(fy(:,2,ii)-fy(:,1,ii))

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

    end associate
  end subroutine init_dc_3

  subroutine solve_ode_nl_dc(sol, pb, t0, T, dt, &
       & meth, out, check_overflow)
    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

    real(RP)     :: tn
    integer      :: nS, jj
    logical      :: exit_comp, b

    if (choral_verb>1) write(*,*) &
         &"ode_NL_DC_mod: solve"
    if (choral_verb>2) write(*,*) &
         &"  Method                         = ",&
         & name_ode_method(meth)

    b = check_ode_method_nl_dc(meth)
    if (.NOT.b) call quit("ode_NL_DC_mod: solve_ode_NL_DC:&
         & wrong method")

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

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

    !! 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
       call out(tn, sol, exit_comp)
       if (exit_comp) then
          if (choral_verb>2) write(*,*) &
               &"ode_NL_DC_mod: solve&
               & time     = ", tn, ': EXIT_COMPp'
          return
       end if

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

       case(ode_dc_3)
          call nl_dc_3(sol, dt, tn, pb)

       end select
       
    end do

  end subroutine solve_ode_nl_dc



  subroutine nl_dc_2(sol, dt, t, pb)
    type(ode_solution), intent(inout) :: sol
    real(RP)          , intent(in)    :: dt, t
    type(ode_problem) , intent(in)    :: pb    
    
    integer                   :: ii, N
    real(RP), dimension(pb%N) :: y0, F0, y1, F1
    real(RP), dimension(0)    :: ay

    n  = pb%N

    associate( y=>sol%Y, fy=>sol%BY, x=>pb%X)

      !$OMP PARALLEL PRIVATE(y0, F0, y1, F1, ay) 
      !$OMP DO
      do ii=1, sol%dof

         !! y0 = y_0^{n+1} := y_0^{n} + dt * F(t_n, y_0^{n})
         y0 = y(:,2,ii) + dt*fy(:,1,ii)
         
         !! F0 = F_0^{n+1} := F(y0^{n+1}, t_{n+1})
         call pb%AB(ay, f0, x(:,ii), t+dt, y0, n, 0 )

         !! F1 = F_1^{n+1} := F(t_{n+1}, y_0^{n+1} + dt*y_1^{n})
         y1 = y0 + dt*y(:,3,ii)
         call pb%AB(ay, f1, x(:,ii), t+dt, y1, n, 0 )

         !! y1 y_1^{n+1} 
         y1 = y(:,3,ii) + f1 - 0.5_rp*( f0 + fy(:,1,ii) )

         !! end, updates
         y(:,1,ii)  = y0 + dt*y1
         y(:,2,ii)  = y0
         y(:,3,ii)  = y1
         fy(:,1,ii) = f0

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

    end associate

  end subroutine nl_dc_2



  subroutine nl_dc_3(sol, dt, t, pb)
    type(ode_solution), intent(inout) :: sol
    real(RP)          , intent(in)    :: dt, t
    type(ode_problem) , intent(in)    :: pb    
    
    real(RP), dimension(pb%N) :: y0, y1, F0, F1
    real(RP), dimension(pb%N) :: F2, y2
    real(RP), dimension(0)    :: ay
    real(RP), parameter       :: C1_6 = 1.0_rp/6.0_rp
    integer                   :: ii, N
    real(RP)                  :: t1, t2, dt_inv

    n  = pb%N
    t1 = t+dt
    t2 = t+dt+dt
    dt_inv = 1.0_rp / dt

    associate( y=>sol%Y, fy=>sol%BY, x=>pb%X)

      !$OMP PARALLEL PRIVATE(y0, F0, y1, F1, ay, F2, y2) 
      !$OMP DO
      do ii=1, sol%dof

         ! !!!!!!!!!!!!!!!!!!!!!!!
         ! START DC-2 ar stage n+2

         ! y0 := Y_0^{n+2}
         y0 = y(:,5,ii) + dt*fy(:,2,ii)

         ! F0 := F_0^{n+2}
         call pb%AB(ay, f0, x(:,ii), t2, y0, n, 0 )

         ! F1 = F_1^{n+2}
         y1 = y0 + dt*y(:,6,ii)
         call pb%AB(ay, f1, x(:,ii), t2, y1, n, 0 )

         ! y1 = Y_1^{n+2}
         y1 = y(:,6,ii) + f1 - (f0 + fy(:,2,ii))*0.5_rp

         ! END DC-2 
         ! !!!!!!!!!!!!!!!!!!!!!!!

         ! F2 = F_2^{n+1}
         y2 = y(:,5,ii) + y(:,6,ii)*dt + y(:,4,ii)*dt**2
         call pb%AB(ay, f2, x(:,ii), t1, y2, n, 0 )

         ! Computation of y2 = y_2^{n+1}
         !
         ! y2 = (1/6) * d^3Y_0^{n+3} * dt**2
         y2 = f0 - 2.0_rp*fy(:,2,ii) + fy(:,1,ii)
         y2 = y2 * c1_6
         !
         ! y2 = [ (1/6) * d^3Y_0^{n+3} - (1/2) * d^2Y_1^{n+2} ]
         !      * dt**2
         y2 = y2 - ( y1 - 2.0_rp*y(:,6,ii) + y(:,3,ii) ) * 0.5_rp
         !
         ! y2 = [ (1/6) * d^3Y_0^{n+3} - (1/2) * d^2Y_1^{n+2} 
         !        + F_2^{n+1} - F_1^{n+1} ] * dt**2
         y2 =  y2 + f2 - fy(:,3,ii) 
         !
         ! y2 = y_2^{n+1}
         y2 = y(:,4,ii) + y2 * dt_inv

         !! end, updates
         y(:,1,ii) = y(:,5,ii) + dt*y(:,6,ii) + dt**2*y2
         y(:,2,ii) = y(:,5,ii)
         y(:,3,ii) = y(:,6,ii)
         y(:,4,ii) = y2
         y(:,5,ii) = y0
         y(:,6,ii) = y1

         fy(:,1,ii) = fy(:,2,ii)
         fy(:,2,ii) = f0          
         fy(:,3,ii) = f1

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

    end associate

  end subroutine nl_dc_3

end module ode_nl_dc_mod