sources_rotframe.f90

Go to the documentation of this file.

!#############################################################################
!#                                                                           #
!# fosite - 3D hydrodynamical simulation program                             #
!# module: sources_rotframe.f90                                              #
!#                                                                           #
!# Copyright (C) 2010-2018                                                   #
!# Tobias Illenseer <tillense@astrophysik.uni-kiel.de>                       #
!# Jannes Klee      <jklee@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 sources_rotframe_mod
  USE sources_c_accel_mod
  USE sources_base_mod
  USE physics_base_mod
  USE fluxes_base_mod
  USE mesh_base_mod
  USE marray_base_mod
  USE marray_compound_mod
  USE common_dict
  !--------------------------------------------------------------------------!
  PRIVATE
  CHARACTER(LEN=32), PARAMETER :: source_name = "inertial forces"
  TYPE, EXTENDS(sources_c_accel) :: sources_rotframe
    TYPE(marray_base) :: cent
    TYPE(marray_base) :: centproj
    TYPE(marray_base) :: cos1
  CONTAINS
    PROCEDURE :: initsources_rotframe
    PROCEDURE :: infosources
    PROCEDURE :: externalsources_single
    PROCEDURE :: convert2rotatingframe
    PROCEDURE :: finalize
  END TYPE

  PUBLIC :: &
       ! classes
       sources_rotframe
  !--------------------------------------------------------------------------!

CONTAINS

  SUBROUTINE initsources_rotframe(this,Mesh,Physics,Fluxes,config,IO)
    USE physics_euler_mod, ONLY: physics_euler
    USE physics_eulerisotherm_mod, ONLY: physics_eulerisotherm
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(sources_rotframe)          :: this
    CLASS(mesh_base),     INTENT(IN) :: Mesh
    CLASS(physics_base),  INTENT(IN) :: Physics
    CLASS(fluxes_base),   INTENT(IN) :: Fluxes
    TYPE(Dict_TYP),          POINTER :: config, IO
    !------------------------------------------------------------------------!
    INTEGER           :: stype
    !------------------------------------------------------------------------!
    CALL getattr(config, "stype", stype)
    CALL this%InitLogging(stype,source_name)
    CALL this%InitSources(mesh,fluxes,physics,config,io)

    SELECT TYPE(physics)
    TYPE IS(physics_euler)
      ! do nothing
    TYPE IS(physics_eulerisotherm)
      ! do nothing
    CLASS DEFAULT
       CALL this%Error("sources_rotframe::InitSources","physics not supported")
    END SELECT

    IF (mesh%FARGO.GT.0) &
      CALL this%Error("sources_rotframe::InitSources", &
         "rotating frame with FARGO seems broken, don't use it together")

    IF (physics%VDIM.NE.2) &
      CALL this%Error("sources_rotframe::InitSources", &
         "Only 2D simulations working with rotating frame at the moment." )

    this%accel = marray_base(physics%VDIM)
    this%accel%data1d(:) = 0.
    this%cent = marray_base(3)
    this%cent%data1d(:) = 0.
    this%centproj = marray_base(3)
    this%centproj%data1d(:) = 0.
    this%cos1 = marray_base(3)
    this%cos1%data1d(:) = 0.

    CALL getattr(config, "gparam", this%gparam, 1.0)

    ! define position vectors
    ! for bianglespherical no shift of axis possible
    ! (for a planet reasonable)
!    IF ((Mesh%Geometry%GetType()).EQ.BIANGLESPHERICAL) THEN
!
!      this%centproj%data4d(:,:,:,1) = this%gparam*SIN(Mesh%bcenter(:,:,:,1))*&
!                              COS(Mesh%bcenter(:,:,:,1))
!      ! for better performance
!      this%cos1%data4d(:,:,:,1)  = COS(Mesh%bcenter(:,:,:,1))
!      this%cos1%data4d(:,:,:,2)  = COS(Mesh%bcenter(:,:,:,2))
!    ELSE
      IF (abs(mesh%rotcent(1)).LE.tiny(mesh%rotcent(1)) &
       .AND.abs(mesh%rotcent(2)).LE.tiny(mesh%rotcent(2))) THEN
        ! no shift of point mass: set position vector to Mesh defaults
        this%cent%data4d(:,:,:,:) = mesh%posvec%bcenter(:,:,:,:)
      ELSE
        ! shifted center of rotation:
        ! compute curvilinear components of shift vector
        this%cent%data4d(:,:,:,1) = mesh%rotcent(1)
        this%cent%data4d(:,:,:,2) = mesh%rotcent(2)
        this%cent%data4d(:,:,:,3) = mesh%rotcent(3)
        CALL mesh%Geometry%Convert2Curvilinear(mesh%bcenter,this%cent%data4d,this%cent%data4d)
        ! subtract the result from the position vector:
        ! this gives you the curvilinear components of all vectors pointing
        ! from the center of rotation to the bary center of any cell on the mesh
        this%cent%data4d(:,:,:,:) = mesh%posvec%bcenter(:,:,:,:) - this%cent%data4d(:,:,:,:)
      END IF
!    END IF
  END SUBROUTINE initsources_rotframe


  SUBROUTINE infosources(this,Mesh)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(sources_rotframe), INTENT(IN) :: this
    CLASS(mesh_base),        INTENT(IN) :: Mesh
    !------------------------------------------------------------------------!
    CHARACTER(LEN=32) :: omega_str
    !------------------------------------------------------------------------!
    WRITE (omega_str,'(ES8.2)') mesh%OMEGA
    CALL this%Info("            angular velocity:  " // trim(omega_str))
  END SUBROUTINE infosources


  SUBROUTINE externalsources_single(this,Mesh,Physics,Fluxes,Sources,time,dt,pvar,cvar,sterm)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(sources_rotframe), INTENT(INOUT) :: this
    CLASS(mesh_base),        INTENT(IN)    :: Mesh
    CLASS(physics_base),     INTENT(INOUT) :: Physics
    CLASS(fluxes_base),      INTENT(IN)    :: Fluxes
    CLASS(sources_base),     INTENT(INOUT) :: Sources
    REAL,                    INTENT(IN)    :: time, dt
    CLASS(marray_compound),  INTENT(INOUT) :: pvar,cvar,sterm
    !------------------------------------------------------------------------!
    INTEGER       :: i,j,k
    !------------------------------------------------------------------------!
    ! Two cases due to different angles between angular velocity to mesh
    ! 1. only a projected part plays role for bianglespherical geometry
!    IF ((Mesh%Geometry%GetType()).EQ.BIANGLESPHERICAL) THEN
!!NEC$ outerloop_unroll(8)
!      DO k=Mesh%KMIN,Mesh%KMAX
!        DO j=Mesh%JMIN,Mesh%JMAX
!!NEC$ ivdep
!          DO i=Mesh%IMIN,Mesh%IMAX
!            this%accel%data4d(i,j,k,1) = Mesh%OMEGA*(Mesh%OMEGA*&
!                   this%centproj%data4d(i,j,k,1) + 2.0*this%cos1%data4d(i,j,k,1)*&
!                   pvar%data4d(i,j,k,Physics%YVELOCITY))
!            this%accel%data4d(i,j,k,2) = -Mesh%OMEGA*2.0*this%cos1%data4d(i,j,k,1)*&
!                   pvar%data4d(i,j,k,Physics%XVELOCITY)
!          END DO
!        END DO
!      END DO
!    ! 2. OMEGA is always perpendicular to other curvilinear coordinates
!    ELSE
!NEC$ outerloop_unroll(8)
      DO k=mesh%KMIN,mesh%KMAX
        DO j=mesh%JMIN,mesh%JMAX
  !NEC$ ivdep
          DO i=mesh%IMIN,mesh%IMAX
            ! components of centrifugal and coriolis acceleration
            this%accel%data4d(i,j,k,1) = mesh%OMEGA*(mesh%OMEGA*this%cent%data4d(i,j,k,1) &
                 + 2.0*pvar%data4d(i,j,k,physics%YVELOCITY))
            this%accel%data4d(i,j,k,2) = mesh%OMEGA*(mesh%OMEGA*this%cent%data4d(i,j,k,2) &
                 - 2.0*pvar%data4d(i,j,k,physics%XVELOCITY))
          END DO
        END DO
      END DO
!    END IF

    ! inertial forces source terms
    CALL physics%ExternalSources(this%accel,pvar,cvar,sterm)
  END SUBROUTINE externalsources_single

  SUBROUTINE convert2rotatingframe(this,Mesh,Physics,pvar)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(sources_rotframe), INTENT(IN)    :: this
    CLASS(mesh_base),        INTENT(IN)    :: Mesh
    CLASS(physics_base),     INTENT(IN)    :: Physics
    REAL, DIMENSION(Mesh%IGMIN:Mesh%IGMAX,Mesh%JGMIN:Mesh%JGMAX,Mesh%KGMIN:Mesh%KGMAX,Physics%VNUM), &
                             INTENT(INOUT) :: pvar
    !------------------------------------------------------------------------!
    ! Convert velocities to the rotating frame
!    IF ((Mesh%Geometry%GetType()).EQ.BIANGLESPHERICAL) THEN
!      pvar(:,:,:,Physics%XVELOCITY) = pvar(:,:,:,Physics%XVELOCITY)
!      pvar(:,:,:,Physics%YVELOCITY) = pvar(:,:,:,Physics%YVELOCITY) - &
!                            Mesh%OMEGA*SIN(Mesh%bcenter(:,:,:,1))*this%gparam
!    ELSE
      pvar(:,:,:,physics%XVELOCITY) = pvar(:,:,:,physics%XVELOCITY) + &
                                      mesh%OMEGA * this%cent%data4d(:,:,:,2)
      pvar(:,:,:,physics%YVELOCITY) = pvar(:,:,:,physics%YVELOCITY) - &
                                      mesh%OMEGA * this%cent%data4d(:,:,:,1)
!    END IF
  END SUBROUTINE convert2rotatingframe

  SUBROUTINE finalize(this)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(sources_rotframe), INTENT(INOUT) :: this
    !------------------------------------------------------------------------!
    CALL this%accel%Destroy()
    CALL this%cent%Destroy()
    CALL this%centproj%Destroy()
    CALL this%cos1%Destroy()

    CALL this%sources_c_accel%Finalize()
  END SUBROUTINE finalize

END MODULE sources_rotframe_mod