ode_NL_1s_mod.f90

Go to the documentation of this file.

module ode_nl_1s_mod

  use choral_constants
  use choral_variables
  use basic_tools
  use real_type
  use ode_def
  use ode_problem_mod
  use ode_solution_mod
  !$ use OMP_LIB

  implicit none
  private

  public :: solve_ode_nl_1s
  public :: create_ode_nl_1s_sol
  public :: set_solver_ode_nl_1s
  public :: memsize_ode_nl_1s
  public :: check_ode_method_nl_1s

  
contains   


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

    b = .false.

    select case(method)   
    case(ode_fe, ode_fbe, ode_rk2, ode_rk4, &
         & ode_erk1, ode_erk2_a, ode_erk2_b, &
         & ode_modif_erk2_b )
       b = .true.
       
    end select

  end function check_ode_method_nl_1s


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

    select case(method)
    
    case(ode_fe, ode_fbe, ode_rk2, ode_rk4, &
         & ode_erk1, ode_erk2_a, ode_erk2_b, &
         & ode_modif_erk2_b )
       n_y  =  1; n_fy = 0  

    case default
       call quit("ode_NL_1s_mod: memSize_ode_NL_1s: uncorrect method")
       
    end select

  end subroutine memsize_ode_nl_1s


  subroutine set_solver_ode_nl_1s(slv, method) 
    integer                    , intent(in) :: method
    procedure(ode_nl_1s_solver), pointer    :: slv

    select case(method)
    case(ode_fbe)
       slv => nl_1s_fbe
    case(ode_fe)
       slv => nl_1s_fe
    case(ode_rk2)
       slv => nl_1s_rk2
    case(ode_rk4)
       slv => nl_1s_rk4
    case(ode_erk1)
       slv => nl_1s_erk1
    case(ode_erk2_a)
       slv => nl_1s_erk2_a
    case(ode_erk2_b)
       slv => nl_1s_erk2_b
    case(ode_modif_erk2_b)
       slv => nl_1s_modif_erk2_a
       
    case default
       call quit("ode_NL_1s_mod: set_solver_ode_NL_1s;&
            & uncorrect method")
    end select
      
  end subroutine set_solver_ode_nl_1s


  subroutine create_ode_nl_1s_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
    logical :: bool

    bool = check_ode_method_nl_1s(method)
    if (.NOT.bool) call quit(&
         & "ode_NL_1s_mod: create_ode_NL_1s_sol: uncorrect method")

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

  end subroutine create_ode_nl_1s_sol



  subroutine solve_ode_nl_1s(sol, pb, t0, T, dt, method, &
       & 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)    :: method
    procedure(ode_output_proc)                 :: out
    logical           , intent(in)    :: check_overflow

    procedure(ode_nl_1s_solver), pointer :: slv=>null()
    real(RP) :: tn
    logical  :: ierr
    integer  :: ns, jj
    logical  :: exit_comp

    if (choral_verb>1) write(*,*) &
         &"ode_NL_1s_mod   : solve_ode_NL_1s"
    if (choral_verb>2) write(*,*) &
         & "  NonLin  onestep solver         = ",&
         & name_ode_method(method)

    call set_solver_ode_nl_1s(slv, method)

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

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

       call slv(sol, dt, tn, pb)

       if (check_overflow) then
          ierr = overflow(sol%Y(:,1,:))
          if (ierr) then
             sol%ierr = 10
             if (choral_verb > 1) write(*,*) &
                  & "ode_NL_1s_mod   : solve,&
                  &    time = ", tn, ': OVERFLOW'
             return
          end if
       end if

    end do

  end subroutine solve_ode_nl_1s


  
  subroutine nl_1s_fbe(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)  :: BY
    real(RP), dimension(pb%Na) :: AY
    integer :: ii, Na, N
    
    na = pb%Na
    n  = pb%N

    !$OMP PARALLEL PRIVATE(AY, BY) 
    !$OMP DO
    do ii=1, sol%dof
       
       call pb%AB(ay, by,&
            & pb%X(:,ii), t, sol%Y(:,1,ii), n, na)
       
       sol%Y(:,1,ii) = sol%Y(:,1,ii) + dt*by
       
       sol%Y(1:na,1,ii) = sol%Y(1:na,1,ii) / &
            & (1._rp - dt * ay )                 
    end do
    !$OMP END DO
    !$OMP END PARALLEL

  end subroutine nl_1s_fbe


  subroutine nl_1s_rk2(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)  :: K, Y
    real(RP), dimension(pb%Na) :: AY
    real(RP) :: tph2, h2
    integer  :: ii, Na, N
    
    h2 = dt / 2.0_rp

    tph2 = t + h2

    na = pb%Na
    n  = pb%N

    !$OMP PARALLEL PRIVATE(K, Y, AY) 
    !$OMP DO
    do ii=1, sol%dof
       
       y = sol%Y(:,1,ii)
       call pb%AB(ay, k, pb%X(:,ii), t, y, n, na)
       k(1:na) = k(1:na) + ay*y(1:na)
       
       y  = sol%Y(:,1,ii) + h2*k
       call pb%AB(ay, k, pb%X(:,ii), tph2, y, n, na)
       k(1:na) = k(1:na) + ay*y(1:na)
       
       sol%Y(:,1,ii) = sol%Y(:,1,ii) + dt * k
       
    end do
    !$OMP END DO
    !$OMP END PARALLEL

  end subroutine nl_1s_rk2




  elemental subroutine phi_1(z)
    real(RP), intent(inout) :: z
    
    if (abs(z)> real_tol) then
       z = (exp(z) - 1._rp)/z       
    else
       z = 1._rp       
    end if

  end subroutine phi_1


  elemental subroutine phi_2(z)
    real(RP), intent(inout) :: z

    real(RP) :: z2

    z2 = z**2

    if (abs(z2)>real_tol) then
       z = (exp(z) - z - 1._rp)/z2

    else
       z = 0.5_rp

    end if

  end subroutine phi_2


  subroutine nl_1s_fe(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)  :: y, b
    real(RP), dimension(pb%Na) :: a
    integer  :: ii, Na, N
    
    na = pb%Na
    n  = pb%N

    !$OMP PARALLEL PRIVATE(y, b, a) 
    !$OMP DO
    do ii=1, sol%dof
       
       y = sol%Y(:,1,ii)
       call pb%AB(a, b, pb%X(:,ii), t, y, n, na)
       b(1:na) = b(1:na) + a*y(1:na)
       sol%Y(:,1,ii) = y + dt*b

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

  end subroutine nl_1s_fe

  subroutine nl_1s_erk1(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)  :: y, b
    real(RP), dimension(pb%Na) :: a
    integer  :: ii, Na, N
    
    na = pb%Na
    n  = pb%N

    !$OMP PARALLEL PRIVATE(y, b, a) 
    !$OMP DO
    do ii=1, sol%dof
       
       y = sol%Y(:,1,ii)
       call pb%AB(a, b,&
            & pb%X(:,ii), t, y, n, na)

       b(1:na) = b(1:na) + a*y(1:na)

       a = a*dt
       call phi_1(a)
       b(1:na) = a * b(1:na)
       y = y + dt*b

       sol%Y(:,1,ii) = y

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

  end subroutine nl_1s_erk1


  subroutine nl_1s_erk2_a(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)  :: y, y2, G1, G2
    real(RP), dimension(pb%Na) :: A, phi
    real(RP) :: t2, h2
    integer  :: ii, Na, N
    
    real(RP), parameter :: C2 = 0.5_rp !! ERK_2 parameter c_2

    real(RP), parameter :: S2 = 1._rp / c2

    h2 = dt * c2
    t2 = t + h2       
    na = pb%Na
    n  = pb%N

    !$OMP PARALLEL PRIVATE(y, y2, G1, G2, A, phi) 
    !$OMP DO
    do ii=1, sol%dof
       
       !! y  <-- yn
       !! A  <-- A1 := a(tn, yn)
       !! G1 <-- B1 := b(tn, yn)
       !!
       y = sol%Y(:,1,ii)
       call pb%AB(a, g1, pb%X(:,ii), t, y, n, na)

       !! G1 <-- Ayn + B1
       !!
       g1(1:na) = g1(1:na) + a*y(1:na)

       !! y2 <-- yn + h/2 phi_1(h/2 A) G1
       !!
       phi      = a * h2
       call phi_1(phi)
       y2       = g1
       y2(1:na) = phi * y2(1:na)
       y2       = y + h2*y2

       !! phi <-- A2 := a(t2, y2)
       !! G2  <-- B2 := b(t2, y2)
       !!
       call pb%AB(phi, g2, pb%X(:,ii), t2, y2, n, na)

       !! G2  <-- G^n_2(t2, y2) = B2 + (A2 - A) y2
       !!
       phi = phi - a
       g2(1:na) = g2(1:na) + phi * y2(1:na)

       !! G2  <-- G^n_2(t2, y2) +  A yn 
       !!
       g2(1:na) = g2(1:na) + a * y(1:na)

       !! phi <-- phi_1(hA)
       !! A   <-- phi_2(hA)
       !!
       a   = a * dt
       phi = a 
       call phi_1(phi)
       call phi_2(a)
       
       !! A   <-- 1/c_2 phi_2(hA)
       !! phi <-- phi_1(hA) - 1/c_2 phi_2(hA) 
       !!
       a   = s2  * a
       phi = phi - a

       !! compute y_{n+1}
       !!
       g1(1:na)  = phi * g1(1:na)
       g1(na+1:) = 0._rp

       g2(1:na)  = a * g2(1:na)

       sol%Y(:,1,ii) = y + dt * (g1 + g2)
     
    end do
    !$OMP END DO
    !$OMP END PARALLEL

  end subroutine nl_1s_erk2_a


  subroutine nl_1s_erk2_b(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)  :: y, y2, G1, G2
    real(RP), dimension(pb%Na) :: A, phi
    real(RP) :: t2, h2
    integer  :: ii, Na, N

    real(RP), parameter :: C2 = 0.5_rp !! ERK_2 parameter c_2

    real(RP), parameter :: S2 = 1._rp / (2._rp * c2)
    real(RP), parameter :: S1 = 1._rp - s2
    
    h2 = dt * c2
    t2 = t + h2
    na = pb%Na
    n  = pb%N

    !$OMP PARALLEL PRIVATE(y, y2, G1, G2, A, phi) 
    !$OMP DO
    do ii=1, sol%dof
       
       
       !! y  <-- yn
       !! A  <-- A1 := a(tn, yn)
       !! G1 <-- B1 := b(tn, yn)
       !!
       y = sol%Y(:,1,ii)
       call pb%AB(a, g1, pb%X(:,ii), t, y, n, na)

       !! G1 <-- Ayn + B1
       !!
       g1(1:na) = g1(1:na) + a*y(1:na)

       !! y2 <-- yn + c_2 h phi_1(c_2 h A) G1
       !!
       phi      = a * h2
       call phi_1(phi)
       y2       = g1
       y2(1:na) = phi * y2(1:na)
       y2       = y + h2*y2

       !! phi <-- A2 := a(t2, y2)
       !! G2  <-- B2 := b(t2, y2)
       !!
       call pb%AB(phi, g2, pb%X(:,ii), t2, y2, n, na)

       !! G2  <-- G^n_2(t2, y2) = B2 + (A2 - A) y2
       !!
       phi = phi - a
       g2(1:na) = g2(1:na) + phi * y2(1:na)

       !! G2  <-- G^n_2(t2, y2) +  A yn 
       !!
       g2(1:na) = g2(1:na) + a * y(1:na)

       !! compute y_{n+1}
       !!
       phi      = a * dt
       call phi_1(phi)

       g1(1:na) = phi * g1(1:na)
       g2(1:na) = phi * g2(1:na)
       g1       = s1*g1 + s2*g2

       sol%Y(:,1,ii) = sol%Y(:,1,ii) + dt*g1
       
    end do
    !$OMP END DO
    !$OMP END PARALLEL

  end subroutine nl_1s_erk2_b




  subroutine nl_1s_modif_erk2_b(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)  :: y, b, y2
    real(RP), dimension(pb%Na) :: a
    real(RP) :: t2, h2
    integer  :: ii, Na, N
    
    h2 = dt * 0.5_rp
    t2 = t + h2
    na = pb%Na
    n  = pb%N

    !$OMP PARALLEL PRIVATE(y, b, y2, a) 
    !$OMP DO
    do ii=1, sol%dof
       
       !! compute a_n, b_n
       !!
       y = sol%Y(:,1,ii)
       call pb%AB(a, b, pb%X(:,ii), t, y, n, na)

       !! compute y2 = y_{n+1/2} with ERK1, half a time step
       !!
       b(1:na) = b(1:na) + a*y(1:na)
       a       = a*h2
       call phi_1(a)
       b(1:na) = a * b(1:na)
       y2      = y + h2*b

       !! compute a_{n+1/2}, b_{n+1/2}
       !!
       call pb%AB(a, b, pb%X(:,ii), t2, y2, n, na)

       !! compute y_{n+1}
       !!
       b(1:na) = b(1:na) + a*y(1:na)
       a       = a*dt
       call phi_1(a)
       b(1:na) = a * b(1:na)

       sol%Y(:,1,ii) = sol%Y(:,1,ii) + dt*b
       
    end do
    !$OMP END DO
    !$OMP END PARALLEL

  end subroutine nl_1s_modif_erk2_b



  subroutine nl_1s_modif_erk2_a(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)  :: y1, G1, y2, G2
    real(RP), dimension(pb%Na) :: a1, a2, phi1, phi2
    real(RP) :: t2, h2
    integer  :: ii, Na, N
    
    h2 = dt * 0.5_rp
    t2 = t + h2
    na = pb%Na
    n  = pb%N

    !$OMP PARALLEL PRIVATE(y1, G1, y2, G2, a1, a2, phi1, phi2) 
    !$OMP DO
    do ii=1, sol%dof
       
       !! compute a1, G1
       !!
       y1 = sol%Y(:,1,ii)
       call pb%AB(a1, g1, pb%X(:,ii), t, y1, n, na)

       !! compute y2 
       !!
       y2       = g1
       y2(1:na) = y2(1:na) + a1*y1(1:na)
       phi1     = a1*h2
       call phi_1(phi1)
       y2(1:na) = phi1 * y2(1:na)
       y2       = y1 + h2*y2

       !! compute a2, G2
       !!
       call pb%AB(a2, g2, pb%X(:,ii), t2, y2, n, na)

       !! compute y_{n+1}
       !!
       g1(1:na)  = g1(1:na) + a1*y1(1:na)
       g1(na+1:) = 0._rp
       phi1      = a1*dt
       phi2      = phi1
       call phi_1(phi1)
       call phi_2(phi2)
       phi1      = phi1 - 2._rp*phi2
       g1(1:na)  = g1(1:na) * phi1

       g2(1:na)  = g2(1:na) + a2*y1(1:na)
       phi2      = a2*dt
       call phi_2(phi2)
       phi2      = phi2 * 2._rp
       g2(1:na)  = g2(1:na) * phi2

       g1            = g1 + g2
       sol%Y(:,1,ii) = y1 + dt*g1 
       
    end do
    !$OMP END DO
    !$OMP END PARALLEL

  end subroutine nl_1s_modif_erk2_a




  subroutine nl_1s_rk4(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)  :: K1, K2, K3, K4, Y
    real(RP), dimension(pb%Na) :: AY
    real(RP) :: h6, h3, h2, tph, tph2
    integer  :: ii, Na, N
    
    h2 = dt / 2.0_rp
    h6 = dt / 6.0_rp
    h3 = dt / 3.0_rp
    
    tph  = t  + dt
    tph2 = t + h2

    na = pb%Na
    n  = pb%N

    !$OMP PARALLEL PRIVATE(K1, K2, K3, K4, Y, AY)
    !$OMP DO
    do ii=1, sol%dof
       
       y = sol%Y(:,1,ii)
       call pb%AB(ay, k1, pb%X(:,ii), t, y, n, na)
       k1(1:na) = k1(1:na) + ay*y(1:na)
       
       y  = sol%Y(:,1,ii) + h2*k1
       call pb%AB(ay, k2, pb%X(:,ii), tph2, y, n, na)
       k2(1:na) = k2(1:na) + ay*y(1:na)

       y  = sol%Y(:,1,ii) + h2*k2
       call pb%AB(ay, k3, pb%X(:,ii), tph2, y, n, na)
       k3(1:na) = k3(1:na) + ay*y(1:na)

       y  = sol%Y(:,1,ii) + dt*k3
       call pb%AB(ay, k4, pb%X(:,ii), tph, y, n, na)
       k4(1:na) = k4(1:na) + ay*y(1:na)

       y = k1*h6 + k2*h3 + k3*h3 + k4*h6
       sol%Y(:,1,ii) = sol%Y(:,1,ii) + y
       
    end do
    !$OMP END DO
    !$OMP END PARALLEL

  end subroutine nl_1s_rk4

  
end module ode_nl_1s_mod