ode_Lin_0d.f90

Go to the documentation of this file.

program ode_lin_0d

  use choral_constants
  use choral

  implicit none

  !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !!
  !!      VARAIBLE DEFINITION : BEGIN
  !!
  !!  

  !! verbosity level
  integer, parameter :: verb = 2

  !!         TIME DISCRETISATION
  !!
  !! pb      = definition of the ode problem
  !! slv     = definition of the ode solver
  !! sol     = data structure for the ode solution
  !!
  type(ode_solution) :: sol
  type(ode_problem)  :: pb 
  type(ode_solver)   :: slv
  !! t0      = initial time
  !! T       = final time
  !! dt      = time step
  !! L_meth   = method for the non-ilinear system
  !!
  real(RP), parameter :: t0       = 0.0_rp
  real(RP), parameter :: T        = 1.5_rp
  integer , parameter :: L_meth   = ode_cn
  real(RP), parameter :: dt       = 0.1_rp
  !!
  !! OUTPUT DEF.
  !!   ode_out = definition of the output
  type(ode_output) :: ode_out


  !!
  !!   VARAIBLE DEFINITION : END
  !!
  !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!


  call choral_init(verb=verb)
  write(*,*) "ode_Lin_0d      : start"


  !!
  !! output settings
  !!
  write(*,*) ""
  write(*,*) "==================================  OUTPUT SETTINGS"
  ode_out = ode_output(t0, t, n=1)
  call set(ode_out, verb=5)
  call set(ode_out, vxn_period = 0.1_rp, vxn_plot=.true.)
  if (verb>0) call print(ode_out)


  !! !!!!!!!!!!!!!!!!!!!!! ODE PROBLEM DEF.
  !!
  write(*,*) ""
  write(*,*) "==================================  TIME DISCRETISATION"
  !!
  pb = ode_problem(ode_pb_lin, dim=0, m=id, s=id_m) 
  if (verb>0) call print(pb)

  !! !!!!!!!!!!!!!!!!!!!!! ODE SOLVER DEF.
  !!
  slv = ode_solver(pb, ode_slv_1s, l_meth=l_meth, verb=verb)
  if (verb>0) call print(slv)


  !! !!!!!!!!!!!!!!!!!!!!! ODE SOLUTION DEF.
  !!
  sol = ode_solution(slv, pb)
  if (verb>0) call print(sol)

  !!!!!!!!!!!!!!!!!!!!!!!  SOLVE
  !!
  write(*,*) ""
  write(*,*) "==================================  NUMERICAL RESOLUTION"
  write(*,*) "==================================  1- USER-DEFINED OUTPUT"
  !!
  !! initial condition
  call initialcond(sol, pb, slv, t0, (/1.0_rp/))
  !!
  !! set the output to the routine "user_defined_output" below
  call set_ode_solver_output(slv, user_defined_output)
  !!
  !! numerical resolution
  call solve(sol, slv, pb, t0, t, dt)!, output=ode_out)

  write(*,*) ""
  write(*,*) "==================================  NUMERICAL RESOLUTION"
  write(*,*) "==================================  2- PRE-DEFINED OUTPUT&
       & USING 'ode-output' type"
  !!
  !! initial condition
  call initialcond(sol, pb, slv, t0, (/1.0_rp/))
  !!
  !! finalise the output definition
  call assemble(ode_out, dt)
  !!
  !! numerical resolution
  call solve(sol, slv, pb, t0, t, dt, output=ode_out)
  !!
  !! Numerical error at time T
  write(*,*) " Error = |y_h(T) - y(T)| / |y(T)| = ", &
       & abs(sol%V(1,1) - exp(t)) / exp(t)

  write(*,*) "ode_Lin_0d      : end"

contains

  subroutine id(yy,xx)
    real(RP), dimension(:), intent(out) :: yy
    real(RP), dimension(:), intent(in)  :: xx
    
    yy = xx
    
  end subroutine id

  subroutine id_m(yy,xx)
    real(RP), dimension(:), intent(out) :: yy
    real(RP), dimension(:), intent(in)  :: xx
    
    yy = -xx
    
  end subroutine id_m

  subroutine user_defined_output(tn, s, stop)
    real(RP)          , intent(in)    :: tn
    type(ode_solution), intent(in)    :: s
    logical           , intent(inout) :: stop

    if (tn<=real_tol) then 
       print*, 'Solving with a user_defined_output:'
       print*, 'Resolution method = ', name_ode_method(l_meth)
       print*, ' Time t          | y_h(t)        | errro |y_h(t) - y(t)|/y(t)'
    end if
    print*, real(tn, SP), real(s%v, SP), real(abs(s%v-exp(tn))/exp(tn), sp)

    stop = .false.
  end subroutine user_defined_output

end program ode_lin_0d