sources_shearbox.f90

Go to the documentation of this file.

!#############################################################################
!#                                                                           #
!# fosite - 3D hydrodynamical simulation program                             #
!# module: sources_shearbox.f90                                              #
!#                                                                           #
!# Copyright (C) 2010-2024                                                   #
!# Jannes Klee <jklee@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 sources_shearbox_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
#ifdef PARALLEL
#ifdef HAVE_MPI_MOD
  USE mpi
#endif
#endif
  IMPLICIT NONE
#ifdef PARALLEL
#ifdef HAVE_MPIF_H
  include 'mpif.h'
#endif
#endif
  !--------------------------------------------------------------------------!
  PRIVATE
  CHARACTER(LEN=32) :: source_name = "forces in shearing-box"
  !--------------------------------------------------------------------------!
  TYPE, EXTENDS(sources_c_accel) :: sources_shearbox
    REAL      :: sign1, sign2
    INTEGER   :: vel1, vel2
    INTEGER   :: i1, i2
    INTEGER   :: momentum1, momentum2
  CONTAINS
    PROCEDURE :: initsources
    PROCEDURE :: infosources
    PROCEDURE :: externalsources
  END TYPE
 
  PUBLIC :: &
       ! classes
       sources_shearbox
 
CONTAINS
 
  SUBROUTINE initsources(this,Mesh,Physics,Fluxes,config,IO)
    USE physics_euler_mod, ONLY: physics_euler
    USE physics_eulerisotherm_mod, ONLY: physics_eulerisotherm
    USE geometry_cartesian_mod, ONLY: geometry_cartesian
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(sources_shearbox),INTENT(INOUT) :: 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%InitSources_base(stype,source_name)
 
    SELECT TYPE(physics)
    TYPE IS(physics_euler)
      ! do nothing
    TYPE IS(physics_eulerisotherm)
      ! do nothing
    CLASS DEFAULT
      CALL this%Error("InitSources_shearbox","physics not supported")
    END SELECT
 
    SELECT TYPE(geo=>mesh%Geometry)
    TYPE IS(geometry_cartesian)
      ! do nothing
    CLASS DEFAULT
      CALL this%Error("InitSources_shearbox","mesh not supported")
    END SELECT
 
    ALLOCATE(this%accel)
    this%accel = marray_base(physics%VDIM)
    this%accel%data1d(:) = 0.
 
    IF(mesh%shear_dir.EQ.2) THEN
      this%Vel1=physics%YVELOCITY
      this%Vel2=physics%XVELOCITY
      this%SIGN1 = 1.0
      this%SIGN2 = -1.0
      this%I1 = 1
      this%I2 = 2
      this%MOMENTUM1 = physics%XMOMENTUM
      this%MOMENTUM2 = physics%YMOMENTUM
    ELSE IF(mesh%shear_dir.EQ.1) THEN
      this%Vel1=physics%XVELOCITY
      this%Vel2=physics%YVELOCITY
      this%SIGN1 = -1.0
      this%SIGN2 = 1.0
      this%I1 = 2
      this%I2 = 1
      this%MOMENTUM1 = physics%YMOMENTUM
      this%MOMENTUM2 = physics%XMOMENTUM
    END IF
 
    ! set vertical gravitational acceleration for 3D shearing
    ! box simulations using bary center cartesian z-coordinate
    IF (mesh%KNUM.GT.1) THEN
      this%accel%data2d(:,3) = -mesh%OMEGA**2 * mesh%cart%data3d(:,2,3)
    END IF
 
    ! print some information
    CALL this%InfoSources(mesh)
  END SUBROUTINE initsources
 
  SUBROUTINE infosources(this,Mesh)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(sources_shearbox), INTENT(IN) :: this
    CLASS(mesh_base),        INTENT(IN) :: Mesh
    !------------------------------------------------------------------------!
    CHARACTER(LEN=32)       :: omega_str,q_str
    !------------------------------------------------------------------------!
    WRITE (omega_str,'(ES10.2)') mesh%OMEGA
    WRITE (q_str,'(ES10.2)')     mesh%Q
    CALL this%Info("            angular velocity:  " // trim(omega_str))
    CALL this%Info("            shearing parameter:" // trim(q_str))
  END SUBROUTINE infosources
 
  SUBROUTINE externalsources(this,Mesh,Physics,Fluxes,Sources,time,dt,pvar,cvar,sterm)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(sources_shearbox), 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
    !------------------------------------------------------------------------!
    SELECT CASE(mesh%fargo%GetType())
    CASE(0)
      ! fargo transport disabled
!NEC$ IVDEP
      this%accel%data4d(:,:,:,this%I1) = 2*mesh%OMEGA &
              * (mesh%Q*mesh%OMEGA*mesh%bcenter(:,:,:,this%I1) &
               + this%SIGN1*pvar%data4d(:,:,:,this%VEL1))
      this%accel%data4d(:,:,:,this%I2) = 2*mesh%OMEGA*this%SIGN2 &
              * pvar%data4d(:,:,:,this%VEL2)
      ! shearingsheet inertial forces source terms
      CALL physics%ExternalSources(this%accel,pvar,cvar,sterm)
    CASE(3)
      ! fargo transport type 3 enabled
      sterm%data2d(:,physics%DENSITY) = 0.0
!NEC$ IVDEP
      sterm%data2d(:,this%MOMENTUM1) = pvar%data2d(:,physics%DENSITY) &
        *mesh%OMEGA*2.0*this%SIGN1*pvar%data2d(:,this%VEL1)
      sterm%data2d(:,this%MOMENTUM2) = pvar%data2d(:,physics%DENSITY) &
        *mesh%OMEGA*(2.0-mesh%Q)*this%SIGN2*pvar%data2d(:,this%VEL2)
      IF (mesh%KNUM.GT.1) &
        sterm%data2d(:,physics%ZMOMENTUM) = pvar%data2d(:,physics%DENSITY) &
          *this%accel%data2d(:,3)
      IF (physics%PRESSURE .GT. 0) THEN
        sterm%data2d(:,physics%ENERGY) = &
             this%SIGN1*pvar%data2d(:,physics%DENSITY)*mesh%Q*mesh%OMEGA* &
             pvar%data2d(:,this%VEL2)*pvar%data2d(:,this%VEL1)
        IF (mesh%KNUM.GT.1) &
          sterm%data2d(:,physics%ENERGY) = sterm%data2d(:,physics%ENERGY) &
            + cvar%data2d(:,physics%ZMOMENTUM)*this%accel%data2d(:,3)
      END IF
    CASE DEFAULT
      ! other fargo transport schemes are not supported
      CALL this%Error("sources_shearbox::ExternalSources_single", &
                      "currently only Fargo transport type 3 is supported")
    END SELECT
  END SUBROUTINE externalsources
 
END MODULE sources_shearbox_mod