timedisc_rkfehlberg.f90

Go to the documentation of this file.

!#############################################################################
!#                                                                           #
!# fosite - 3D hydrodynamical simulation program                             #
!# module: timedisc_rkfehlberg .f90                                          #
!#                                                                           #
!# Copyright (C) 2011,2014                                                   #
!# Björn Sperling   <sperling@astrophysik.uni-kiel.de>                       #
!# Tobias Illenseer <tillense@astrophysik.uni-kiel.de>                       #
!#                                                                           #
!# This program is free software; you can redistribute it and/or modify      #
!# it under the terms of the GNU General Public License as published by      #
!# the Free Software Foundation; either version 2 of the License, or (at     #
!# your option) any later version.                                           #
!#                                                                           #
!# This program is distributed in the hope that it will be useful, but       #
!# WITHOUT ANY WARRANTY; without even the implied warranty of                #
!# MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE, GOOD TITLE or        #
!# NON INFRINGEMENT.  See the GNU General Public License for more            #
!# details.                                                                  #
!#                                                                           #
!# You should have received a copy of the GNU General Public License         #
!# along with this program; if not, write to the Free Software               #
!# Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.                 #
!#                                                                           #
!#############################################################################

!----------------------------------------------------------------------------!
!----------------------------------------------------------------------------!
MODULE timedisc_rkfehlberg_mod
  USE timedisc_base_mod
  USE mesh_base_mod
  USE fluxes_base_mod
  USE boundary_base_mod
  USE physics_base_mod
  USE timedisc_modeuler_mod
  USE marray_base_mod
  USE marray_compound_mod
  USE sources_base_mod
  USE common_dict
#ifdef PARALLEL
#ifdef HAVE_MPI_MOD
  USE mpi
#endif
#endif
  IMPLICIT NONE
#ifdef PARALLEL
#ifdef HAVE_MPIF_H
  include 'mpif.h'
#endif
#endif
  !--------------------------------------------------------------------------!
  PRIVATE
  TYPE coeff_type
    CLASS(marray_compound), POINTER :: p
  END TYPE coeff_type
  TYPE, EXTENDS (timedisc_modeuler) :: timedisc_rkfehlberg
     TYPE(coeff_type), DIMENSION(:), ALLOCATABLE :: coeff
     REAL, DIMENSION(:), POINTER        :: b_low,b_high,c
     REAL, DIMENSION(:,:), POINTER      :: a
  CONTAINS
    PROCEDURE :: inittimedisc
    PROCEDURE :: inittimedisc_rkfehlberg
    PROCEDURE :: finalize
    PROCEDURE :: solveode
    PROCEDURE :: computecvar_rkfehlberg
    generic   :: computecvar => computecvar_rkfehlberg
    PROCEDURE :: setbutchertableau
    PROCEDURE :: showbutchertableau
  END TYPE timedisc_rkfehlberg
  CHARACTER(LEN=32), PARAMETER :: odesolver_name = "Runge-Kutta Fehlberg"

  !--------------------------------------------------------------------------!
  PUBLIC :: &
       ! types
       coeff_type, &
       timedisc_rkfehlberg
  !--------------------------------------------------------------------------!

CONTAINS

  SUBROUTINE inittimedisc(this,Mesh,Physics,config,IO,ttype,tname)
    USE physics_eulerisotherm_mod, ONLY : physics_eulerisotherm
    USE physics_euler_mod, ONLY : physics_euler
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(timedisc_rkfehlberg), INTENT(INOUT) :: this
    CLASS(mesh_base),     INTENT(INOUT) :: Mesh
    CLASS(physics_base),  INTENT(IN)    :: Physics
    TYPE(Dict_TYP),       POINTER       :: config,IO
    INTEGER,              INTENT(IN)    :: ttype
    CHARACTER(LEN=32),    INTENT(IN)    :: tname
    !------------------------------------------------------------------------!
    INTEGER :: err,k,ShowBT
    !------------------------------------------------------------------------!
    CALL this%timedisc_modeuler%InitTimedisc(mesh,physics,config,io,ttype,tname)

    ALLOCATE(this%b_high(this%m),this%b_low(this%m),&
             this%c(this%m),this%a(this%m,this%m), &
             this%coeff(this%m), &
             stat = err)
    IF (err.NE.0) THEN
       CALL this%Error("timedisc_rkfehlberg::InitTimedisc", "Unable to allocate memory.")
    END IF

    ! init Butcher tableau
    CALL this%SetButcherTableau()
    CALL getattr(config, "ShowButcherTableau", showbt, 0)
    IF (showbt.EQ.1) CALL this%ShowButcherTableau()

    ! init coefficient compounds
    this%coeff(1)%p => this%rhs
    DO k=2,this%m
      CALL physics%new_statevector(this%coeff(k)%p,conservative)
      this%coeff(k)%p%data1d(:) = 0.0
    END DO

    IF ((this%tol_rel.LT.0.0).OR.minval(this%tol_abs(:)).LT.0.0) &
         CALL this%Error("timedisc_rkfehlberg::InitTimedisc", &
         "error tolerance levels must be greater than 0")

    IF (this%tol_rel.GT.1.0) THEN
         CALL this%Warning("timedisc_rkfehlberg::InitTimedisc", &
            "adaptive step size control disabled (tol_rel>1)",0)
    ELSE IF(this%tol_rel.GE.0.01 .AND. this%order .GE. 5) THEN
         CALL this%Warning("timedisc_rkfehlberg::InitTimedisc", &
             "You chose a relatively high tol_rel (in comparison to order)",0)
    END IF

  END SUBROUTINE

  SUBROUTINE inittimedisc_rkfehlberg(this,Mesh,Physics,config,IO)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(Timedisc_rkfehlberg), INTENT(INOUT) :: this
    CLASS(Mesh_base),           INTENT(INOUT) :: Mesh
    CLASS(Physics_base),        INTENT(IN)    :: Physics
    TYPE(Dict_TYP), POINTER                   :: config,IO
    !------------------------------------------------------------------------!
    ! set default order
    CALL getattr(config, "order", this%order, 5)

    ! set number of coefficients
    SELECT CASE(this%GetOrder())
    CASE(3)
      this%m = 3
    CASE(5)
      this%m = 6
    CASE DEFAULT
       CALL this%Error("InitTimedisc_rkfehlberg","time order must be 3 or 5")
    END SELECT

    CALL this%InitTimedisc(mesh,physics,config,io,rk_fehlberg,odesolver_name)
  END SUBROUTINE inittimedisc_rkfehlberg


  SUBROUTINE solveode(this,Mesh,Physics,Sources,Fluxes,time,dt,err)
  IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(Timedisc_rkfehlberg), INTENT(INOUT) :: this
    CLASS(Mesh_base),           INTENT(IN)    :: Mesh
    CLASS(Physics_base),        INTENT(INOUT) :: Physics
    CLASS(Fluxes_base),         INTENT(INOUT) :: Fluxes
    CLASS(sources_base),        POINTER       :: Sources
    REAL,                       INTENT(IN)    :: time
    REAL,                       INTENT(INOUT) :: dt,err
    !------------------------------------------------------------------------!
    INTEGER            :: i,j,k,l,m
    REAL               :: t
    !------------------------------------------------------------------------!
    t = time
    ! ATTENTION: rhs and therefore coeff(1)%p should already be up to date.
    ! In the very first step this is done in FirstStepFosite, any later step
    ! must update coeff(1)%p at the end of the whole time step update. This is
    ! done in timedisc_base::AcceptTimestep or timedisc_base::RejectTimestep.
    ! Don't forget: this%rhs => this%coeff(1)%p.
!NEC$ SHORTLOOP
    DO m=2,this%m
       ! time step update of cell mean values
       CALL this%ComputeCVar_rkfehlberg(mesh,physics,fluxes,dt,m,this%coeff,this%cvar,this%ctmp)
       ! compute right-hand-side
       ! coeff_m is k_m/dt from Butcher tableau
       CALL this%ComputeRHS(mesh,physics,sources,fluxes,t+this%c(m)*dt,dt,&
                           this%ptmp,this%ctmp,check_nothing,this%coeff(m)%p)
    END DO

    !reset ctmp
    this%ctmp = this%cvar
!NEC$ SHORTLOOP
    DO m=1,this%m
        ! compute two solutions with different numerical orders
        ! y_n+1 = y_n + SUM(b_i*k_i) = y_n + SUM(b_i*dt*coeff_i)
        this%ctmp%data1d(:) = this%ctmp%data1d(:) &
                           - this%b_low(m)*dt*this%coeff(m)%p%data1d(:)
        this%cvar%data1d(:) = this%cvar%data1d(:) &
                           - this%b_high(m)*dt*this%coeff(m)%p%data1d(:)
    END DO

    ! at the boundary the this%rhs contains the boundary fluxes
    ! (see subroutine ComputeRHS_modeuler )
!NEC$ SHORTLOOP
    DO m=1,this%m
!NEC$ SHORTLOOP
      DO l=1,physics%VNUM
      IF(mesh%INUM.GT.1) THEN
        ! western and eastern
!NEC$ IVDEP
        DO k=mesh%KMIN,mesh%KMAX
!NEC$ IVDEP
          DO j=mesh%JMIN,mesh%JMAX
             fluxes%bxflux(j,k,1,l) = fluxes%bxflux(j,k,1,l) &
                               - dt*this%b_high(m)*this%coeff(m)%p%data4d(mesh%IMIN-mesh%IP1,j,k,l)
             fluxes%bxflux(j,k,2,l) = fluxes%bxflux(j,k,2,l) &
                               - dt*this%b_high(m)*this%coeff(m)%p%data4d(mesh%IMAX+mesh%IP1,j,k,l)
          END DO
        END DO
      END IF
        ! southern and northern
      IF(mesh%JNUM.GT.1) THEN
!NEC$ IVDEP
        DO k=mesh%KMIN,mesh%KMAX
!NEC$ IVDEP
          DO i=mesh%IMIN,mesh%IMAX
            fluxes%byflux(k,i,1,l) = fluxes%byflux(k,i,1,l) &
                              - dt*this%b_high(m)*this%coeff(m)%p%data4d(i,mesh%JMIN-mesh%JP1,k,l)
            fluxes%byflux(k,i,2,l) = fluxes%byflux(k,i,2,l) &
                              - dt*this%b_high(m)*this%coeff(m)%p%data4d(i,mesh%JMAX+mesh%JP1,k,l)
          END DO
        END DO
      END IF
        ! bottom and top
      IF(mesh%KNUM.GT.1) THEN
!NEC$ IVDEP
        DO j=mesh%JMIN,mesh%JMAX
!NEC$ IVDEP
          DO i=mesh%IMIN,mesh%IMAX
            fluxes%bzflux(i,j,1,l) = fluxes%bzflux(i,j,1,l) &
                              - dt*this%b_high(m)*this%coeff(m)%p%data4d(i,j,mesh%KMIN-mesh%KP1,l)
            fluxes%bzflux(i,j,2,l) = fluxes%bzflux(i,j,2,l) &
                              - dt*this%b_high(m)*this%coeff(m)%p%data4d(i,j,mesh%KMAX+mesh%KP1,l)
          END DO
        END DO
      END IF
      END DO
    END DO

    ! compute an error estimate based on the two independent numerical
    ! solutions err and dt
    IF (this%tol_rel.LE.1.0) THEN
      err = this%ComputeError(mesh,physics,this%cvar,this%ctmp)
      dt = this%AdjustTimestep(err,dt)
    END IF

  END SUBROUTINE solveode

  SUBROUTINE computecvar_rkfehlberg(this,Mesh,Physics,Fluxes,dt,m,coeff,cvar,cnew)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(Timedisc_rkfehlberg),INTENT(INOUT) :: this
    CLASS(Mesh_base)          ,INTENT(IN)    :: Mesh
    CLASS(Physics_base)       ,INTENT(INOUT) :: Physics
    CLASS(Fluxes_base)        ,INTENT(INOUT) :: Fluxes
    REAL                      ,INTENT(IN)    :: dt
    INTEGER                   ,INTENT(IN)    :: m
    TYPE(coeff_type), DIMENSION(m), INTENT(IN) :: coeff
    CLASS(marray_compound)    ,INTENT(INOUT) :: cvar,cnew
    !------------------------------------------------------------------------!
    INTEGER :: mm
    !------------------------------------------------------------------------!
    ! compute y + SUM(a_mi*k_i) = y + SUM(a_mi*dt*rhs_i)  : i in [1,m-1]
    ! see Runge-Kutta method (Butcher tableau)
    IF (m.EQ.1) THEN
      cnew = cvar
    ELSE
      cnew%data1d(:) = cvar%data1d(:) - this%a(m,1)*dt*coeff(1)%p%data1d(:)
!NEC$ SHORTLOOP
      DO mm=2,m-1
        cnew%data1d(:) = cnew%data1d(:) - this%a(m,mm)*dt*coeff(mm)%p%data1d(:)
      END DO
    END IF
  END SUBROUTINE computecvar_rkfehlberg

  SUBROUTINE setbutchertableau(this)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(timedisc_rkfehlberg)   :: this
    !------------------------------------------------------------------------!
    SELECT CASE(this%GetOrder())
    CASE(3) ! RK-Fehlberg rkf23
      this%b_high = (/ 1.0/6.0, 2.0/3.0, 1.0/6.0 /)
      this%b_low  = (/ 0.0, 1.0, 0.0 /)
      this%c  = (/ 0.0, 0.5, 1.0 /)
      this%a  = transpose(reshape((/ &
                0.0, 0.0, 0.0, &
                0.5, 0.0, 0.0, &
              -1.0, 2.0, 0.0 /),(/this%m,this%m/)))
    CASE(5) ! RK-Fehlberg rkf45
      this%b_high = (/ 16.0/135.0, 0.0, 6656.0/12825.0, &
                      28561.0/56430.0, -9.0/50.0, 2.0/55.0 /)
      this%b_low  = (/ 25.0/216.0, 0.0, 1408.0/2565.0, 2197.0/4104.0, -0.2, 0.0 /)
      this%c  = (/ 0.0, 0.25, 3.0/8.0, 12.0/13.0, 1.0, 0.5 /)
      this%a  = transpose(reshape((/ &
                0.0, 0.0, 0.0, 0.0, 0.0, 0.0, &
                0.25, 0.0, 0.0, 0.0, 0.0, 0.0, &
                3.0/32.0, 9.0/32.0, 0.0, 0.0, 0.0, 0.0, &
                1932.0/2197.0, -7200.0/2197.0, 7296.0/2197.0, 0.0, 0.0, 0.0, &
                439.0/216.0, -8.0, 3680.0/513.0, -845.0/4104.0, 0.0, 0.0, &
                -8.0/27.0, 2.0, -3544.0/2565.0, 1859.0/4104.0, -11.0/40.0, 0.0/),&
                (/this%m,this%m/)))
    CASE DEFAULT
       CALL this%Error("timedisc_rkfehlberg::SetButcherTableau","only order 3 or 5 supported")
    END SELECT
  END SUBROUTINE setbutchertableau


  SUBROUTINE showbutchertableau(this)
  IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(Timedisc_rkfehlberg) :: this
    !------------------------------------------------------------------------!
    INTEGER            :: i
    CHARACTER(LEN=64)  :: sformat
    CHARACTER(LEN=128)  :: buffer
    !------------------------------------------------------------------------!
    INTENT(INOUT)      :: this
    !------------------------------------------------------------------------!
    CALL this%Info(repeat("+",(this%m+1)*12+4))
    CALL this%Info("Butcher tableau")
    DO i=1,this%m
      IF (i == 1) THEN
        WRITE (sformat,'(A)') '(ES12.3,A)'
      ELSE
        WRITE (sformat,'(A,I1,A)') '(ES12.3,A,',i-1,'ES12.3)'
      END IF
      WRITE (buffer,fmt=sformat) this%c(i), " | " , this%a(i,1:i-1)
      CALL this%Info(buffer)
    END DO
    CALL this%Info(repeat("-",(this%m+1)*12+4))
    WRITE (sformat,'(A,I1,A)') '(A,',this%m,'ES12.3)'
    WRITE (buffer,fmt=sformat) " high order  | ", this%b_high(:)
    CALL this%Info(buffer)
    WRITE (buffer,fmt=sformat) "  low order  | ", this%b_low(:)
    CALL this%Info(buffer)
    CALL this%Info(repeat("+",(this%m+1)*12+4))
    WRITE (buffer,fmt=sformat) "c_i-SUM(a_ij)| ", this%c(:)-sum(this%a(:,:),dim=2)
    CALL this%Info(buffer)
  END SUBROUTINE showbutchertableau


  SUBROUTINE finalize(this)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(Timedisc_rkfehlberg)   :: this
    !------------------------------------------------------------------------!
    INTEGER :: k
    !------------------------------------------------------------------------!
    DEALLOCATE(this%b_high,this%b_low,this%c,this%a)
    DO k=2,this%m
      CALL this%coeff(k)%p%Destroy()
      DEALLOCATE(this%coeff(k)%p)
    END DO
    DEALLOCATE(this%coeff)
    CALL this%timedisc_modeuler%Finalize()
  END SUBROUTINE finalize

END MODULE timedisc_rkfehlberg_mod