ode_solution_mod.f90

Go to the documentation of this file.

module ode_solution_mod

  use real_type
  use basic_tools
  use ode_def
  use ode_problem_mod

  implicit none
  private

  public :: ode_solution
  public :: clear, print

  public :: ode_lin_solver
  public :: ode_nl_1s_solver
  public :: ode_nl_ms_solver
  public :: ode_solution_init_indexes
  public :: ode_output_proc, void_ode_output

  !       %----------------------------------------%
  !       |                                        |
  !       |          DERIVED TYPE                  |
  !       |                                        |
  !       %----------------------------------------%

  type :: ode_solution

     integer :: type = -1

     integer :: n    = -1

     integer :: na   = -1

     integer :: dof  = -1

     real(RP), dimension(:,:,:), allocatable :: y

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

     real(RP), dimension(:,:,:), allocatable :: by

     real(RP), dimension(:,:)  , allocatable :: v 

     real(RP), dimension(:), allocatable :: rhs

     real(RP), dimension(:), allocatable :: aux

     integer :: nv=0, ny=0, nfy=0

     integer :: ierr = 0

     integer, dimension(:), allocatable :: y_i

     integer, dimension(:), allocatable :: f_i

     integer, dimension(:), allocatable :: v_i
     
   contains

     final :: ode_solution_clear

  end type ode_solution


  !       %----------------------------------------%
  !       |                                        |
  !       |       GENERIc SUBROUTINES              |
  !       |                                        |
  !       %----------------------------------------%
  interface clear
     module procedure ode_solution_clear
  end interface clear

  interface ode_solution
     module procedure ode_solution_create
  end interface ode_solution

  interface print
     module procedure ode_solution_print
  end interface print


  !       %----------------------------------------%
  !       |                                        |
  !       |        ABSTRACT INTERFACES             |
  !       |                                        |
  !       %----------------------------------------%
  abstract interface

     subroutine ode_output_proc(tn, s, stop)
       import :: rp, ode_solution
       real(RP)          , intent(in)    :: tn
       type(ode_solution), intent(in)    :: s
       logical           , intent(inout) :: stop
     end subroutine ode_output_proc

     subroutine ode_lin_solver(sol, ierr, dt, pb, KInv)
       import :: rp, ode_problem, ode_solution, linsystem_solver 
       type(ode_solution) , intent(inout) :: sol
       logical            , intent(out)   :: ierr  
       real(RP)           , intent(in)    :: dt
       type(ode_problem)  , intent(in)    :: pb
       procedure(linsystem_solver)        :: KInv
     end subroutine ode_lin_solver

     subroutine ode_nl_1s_solver(sol, dt, t, pb)
       import :: rp, ode_problem, ode_solution 
       type(ode_solution), intent(inout) :: sol
       real(RP)          , intent(in)    :: dt, t
       type(ode_problem) , intent(in)    :: pb
     end subroutine ode_nl_1s_solver

     subroutine ode_nl_ms_solver(sol, dt, N, Na)
       import :: rp, ode_solution 
       type(ode_solution), intent(inout) :: sol
       real(RP)          , intent(in)    :: dt
       integer           , intent(in)    :: N, Na  
     end subroutine ode_nl_ms_solver

  end interface

contains

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

    stop = .false.
  end subroutine void_ode_output


  subroutine ode_solution_clear(sol)
    type(ode_solution), intent(inout) :: sol

    call freemem( sol%Y )
    call freemem( sol%AY )
    call freemem( sol%BY )
    call freemem( sol%V )
    call freemem( sol%aux )
    call freemem( sol%rhs )

    call freemem( sol%Y_i )
    call freemem( sol%F_i )
    call freemem( sol%V_i )

    sol%N   = -1
    sol%Na  = -1
    sol%dof = -1

    sol%nV  = 0
    sol%nY  = 0
    sol%nFY = 0

    sol%ierr = 0

  end subroutine ode_solution_clear


  function ode_solution_create(pb, nV, NY, NFY) result(sol)
    type(ode_solution)              :: sol
    type(ode_problem)  , intent(in) :: pb
    integer, intent(in), optional   :: nV, NY, NFY

    call clear(sol)

    sol%type = pb%type
    sol%N    = pb%N
    sol%Na   = pb%Na
    sol%dof  = pb%dof

    if (present(nv)) then
       if ( nv < 0) call quit(&
            & "ode_solution_mod: ode_solution_create:&
            &  parameter 'nV' must be >= 0" )

       sol%nV = nv

       call allocmem(sol%V  , pb%dof, nv)
       call allocmem(sol%aux, pb%dof)
       call allocmem(sol%rhs, pb%dof)
       call allocmem(sol%V_i, nv)

    end if

    if (present(ny)) then

       if ( ny < 0) call quit(&
            & "ode_solution_mod: ode_solution_create::&
            &  parameter 'nY' must be >= 0" )

       sol%nY = ny
       call allocmem(sol%Y, pb%N, ny, pb%dof)
       call allocmem(sol%Y_i, ny)

    end if

    if (present(nfy)) then

       if ( nfy < 0) call quit(&
            & "ode_solution_mod: ode_solution_create::&
            &  parameter 'nFY' must be >= 0" )

       sol%nFY = nfy
       call allocmem(sol%BY, pb%N , nfy, pb%dof)
       call allocmem(sol%AY, pb%Na, nfy, pb%dof)
       call allocmem(sol%F_i, nfy)

    end if

    call ode_solution_init_indexes(sol)

  end function ode_solution_create


  subroutine ode_solution_init_indexes(sol)
    type(ode_solution) , intent(inout) :: sol

    integer :: ii

    if (sol%nY > 0) then
       do ii=1, sol%nY
          sol%Y_i(ii) = ii
       end do
    end if

    if (sol%nFY > 0) then
       do ii=1, sol%nFY
          sol%F_i(ii) = ii
       end do
    end if

    if (sol%nV > 0) then
       do ii=1, sol%nV
          sol%V_i(ii) = ii
       end do
    end if


  end subroutine ode_solution_init_indexes

  subroutine ode_solution_print(sol)
    type(ode_solution) , intent(in) :: sol

    write(*,*)"ode_solution_mod: ode_solution_print"

    if (allocated(sol%Y)) then
       write(*,*)"  Array Y ,              shape   =", &
            & shape(sol%Y)
    end if

    if (allocated(sol%AY)) then
       write(*,*)"  Array AY,              shape   =", &
            & shape(sol%AY)
    end if

    if (allocated(sol%BY)) then
       write(*,*)"  Array BY,              shape   =", &
            & shape(sol%BY)
    end if


    if (allocated(sol%V)) then
       write(*,*)"  Array V ,              shape   =", &
            & "            ", shape(sol%V)
    end if

  end subroutine ode_solution_print

end module ode_solution_mod