boundary_shearing.f90

Go to the documentation of this file.

  !#############################################################################
  !#                                                                           #
  !# fosite - 3D hydrodynamical simulation program                             #
  !# module: boundary_shearing.f90                                             #
  !#                                                                           #
  !# Copyright (C) 2006-2018                                                   #
  !# 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 boundary_shearing_mod
  USE boundary_base_mod
  USE boundary_periodic_mod
  USE marray_compound_mod
  USE mesh_base_mod
  USE physics_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
  CHARACTER(LEN=32), PARAMETER  :: boundcond_name = "shearing"

  TYPE, EXTENDS(boundary_periodic) :: boundary_shearing
    REAL                        :: velocity_offset
    REAL, DIMENSION(:), POINTER :: velocity_shift
  CONTAINS
    PROCEDURE :: initboundary_shearing
    PROCEDURE :: setboundarydata
    PROCEDURE :: finalize
  END TYPE
  !--------------------------------------------------------------------------!
  PUBLIC :: boundary_shearing
  !--------------------------------------------------------------------------!

CONTAINS

  SUBROUTINE initboundary_shearing(this,Mesh,Physics,dir,config)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(boundary_shearing), INTENT(INOUT) :: this
    CLASS(mesh_base),         INTENT(IN)    :: Mesh
    CLASS(physics_base),      INTENT(IN)    :: Physics
    TYPE(Dict_TYP), POINTER                 :: config
    INTEGER,                  INTENT(IN)    :: dir
    !------------------------------------------------------------------------!
    INTEGER            :: err, l
    !------------------------------------------------------------------------!
    CALL this%InitBoundary(mesh,physics,shearing,boundcond_name,dir,config)

    ALLOCATE( &
             this%velocity_shift(physics%VNUM+physics%PNUM), &
             stat = err)
    IF (err.NE.0) THEN
       CALL this%Error("InitBoundary_shearing", "Unable to allocate memory.")
    END IF

    DO l=1,physics%VNUM + physics%PNUM
      ! this part for WEST-EAST shear
      IF (mesh%shear_dir.EQ.2) THEN
        IF (l.EQ.physics%YVELOCITY) THEN
          IF (mesh%FARGO.EQ.0) THEN
            this%velocity_shift(l) = mesh%Q*mesh%OMEGA*(mesh%xmax-mesh%xmin)
          ELSE IF (mesh%FARGO.EQ.3) THEN
            this%velocity_shift(l) = 0.0
          ELSE
            CALL this%Error("InitTimedisc", &
            "Shearing boundaries are only compatible without Fargo or Fargo type 3 (shearing box).")
          END IF
        ELSE
            this%velocity_shift(l) = 0.0
        END IF
      ! this part for SOUTH-NORTH shear
      ELSE IF (mesh%shear_dir.EQ.1) THEN
        IF (l.EQ.physics%XVELOCITY) THEN
          IF (mesh%FARGO.EQ.0) THEN
            this%velocity_shift(l) = mesh%Q*mesh%OMEGA*(mesh%ymax-mesh%ymin)
          ELSE IF (mesh%FARGO.EQ.3) THEN
            this%velocity_shift(l) = 0.0
          ELSE
            CALL this%Error("InitTimedisc", &
            "Shearing boundaries are only compatible without Fargo or Fargo type 3 (shearing box).")
          END IF
        ELSE
            this%velocity_shift(l) = 0.0
        END IF
      ELSE
        CALL this%Error("InitBoundary", &
        "Shearing boundaries in top/south direction not allowed, yet.")
      END IF
    END DO
    this%velocity_offset = mesh%Q*mesh%OMEGA*(mesh%xmax-mesh%xmin)/mesh%dy

  END SUBROUTINE initboundary_shearing

  PURE SUBROUTINE setboundarydata(this,Mesh,Physics,time,pvar)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(boundary_shearing), INTENT(INOUT) :: this
    CLASS(mesh_base),         INTENT(IN)    :: mesh
    CLASS(physics_base),      INTENT(IN)    :: physics
    REAL,                     INTENT(IN)    :: time
    CLASS(marray_compound),   INTENT(INOUT) :: pvar
    !------------------------------------------------------------------------!
    INTEGER            :: i,j,k,l,intshift
    REAL               :: offremain,offset,offset_tmp
#ifdef PARALLEL
    INTEGER            :: status(mpi_status_size)
    INTEGER            :: ierror,req(4)
    CHARACTER(LEN=80)  :: str
    REAL               :: mpi_buf(2*mesh%gnum)
#endif
    !------------------------------------------------------------------------!
    ! the routine below expect that periodic boundaries are already applied
#ifndef PARALLEL
    CALL this%boundary_periodic%SetBoundaryData(mesh,physics,time,pvar)
#endif

    ! make sure all MPI processes use the same step if domain is decomposed
    ! along the y-direction (can be different due to round-off errors)
    offset = -this%velocity_offset*time
    DO WHILE(offset<0)
      offset = offset + mesh%JNUM
    END DO

    DO l=1,physics%VNUM+physics%PNUM
      ! chose velocity offset for fargo and considered pvar
      SELECT CASE(this%GetDirection())
      CASE(west)
        offset_tmp = offset
        intshift = floor(offset_tmp)
        offremain = offset_tmp - intshift
        DO i=1,mesh%GNUM
          !------- integral shift ---------------------------------------------!
          pvar%data4d(mesh%IMIN-i,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,l)  = &
            cshift(pvar%data4d(mesh%IMIN-i,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,l),intshift)
          !------- residual shift ---------------------------------------------!
          pvar%data4d(mesh%IMIN-i,mesh%JMAX+1,:,l) = pvar%data4d(mesh%IMIN-i,mesh%JMIN,:,l)
          DO k=mesh%KMIN,mesh%KMAX
            DO j=mesh%JMIN,mesh%JMAX
              pvar%data4d(mesh%IMIN-i,j,k,l) = (1.0 - offremain)*pvar%data4d(mesh%IMIN-i,j,k,l) + &
                offremain*pvar%data4d(mesh%IMIN-i,j+1,k,l) + this%velocity_shift(l)
            END DO
          END DO
        END DO
      CASE(east)
        offset_tmp = -offset
        intshift = floor(offset_tmp)
        offremain = offset_tmp - intshift
        DO i=1,mesh%GNUM
          !------- integral shift ---------------------------------------------!
          pvar%data4d(mesh%IMAX+i,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,l)  = &
            cshift(pvar%data4d(mesh%IMAX+i,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,l),intshift)
          !------- residual shift ---------------------------------------------!
          pvar%data4d(mesh%IMAX+i,mesh%JMAX+1,:,l) = pvar%data4d(mesh%IMAX+i,mesh%JMIN,:,l)
          DO k=mesh%KMIN,mesh%KMAX
            DO j=mesh%JMIN,mesh%JMAX
              pvar%data4d(mesh%IMAX+i,j,k,l) = (1.0 - offremain)*pvar%data4d(mesh%IMAX+i,j,k,l) + &
                offremain*pvar%data4d(mesh%IMAX+i,j+1,k,l) - this%velocity_shift(l)
            END DO
          END DO
        END DO
      CASE(south)
        offset_tmp = -offset
        intshift = floor(offset_tmp)
        offremain = offset_tmp - intshift
        DO j=1,mesh%GNUM
          !------- integral shift ---------------------------------------------!
          pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN-j,mesh%KMIN:mesh%KMAX,l)  = &
            cshift(pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN-j,mesh%KMIN:mesh%KMAX,l),intshift)
          !------- residual shift ---------------------------------------------!
          pvar%data4d(mesh%IMAX+1,mesh%JMIN-j,:,l) = pvar%data4d(mesh%IMIN,mesh%JMIN-j,:,l)
          DO k=mesh%KMIN,mesh%KMAX
            DO i=mesh%IMIN,mesh%IMAX
              pvar%data4d(i,mesh%JMIN-j,k,l) = (1.0 - offremain)*pvar%data4d(i,mesh%JMIN-j,k,l) + &
                offremain*pvar%data4d(i+1,mesh%JMIN-j,k,l) - this%velocity_shift(l)
            END DO
          END DO
        END DO
      CASE(north)
        offset_tmp = offset
        intshift = floor(offset_tmp)
        offremain = offset_tmp - intshift
        DO j=1,mesh%GNUM
          !------- integral shift ---------------------------------------------!
          pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMAX+j,mesh%KMIN:mesh%KMAX,l)  = &
            cshift(pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMAX+j,mesh%KMIN:mesh%KMAX,l),intshift)
          !------- residual shift ---------------------------------------------!
          pvar%data4d(mesh%IMAX+1,mesh%JMAX+j,:,l) = pvar%data4d(mesh%IMIN,mesh%JMAX+j,:,l)
          DO k=mesh%KMIN,mesh%KMAX
            DO i=mesh%IMIN,mesh%IMAX
              pvar%data4d(i,mesh%JMAX+j,k,l) = (1.0 - offremain)*pvar%data4d(i,mesh%JMAX+j,k,l) + &
                offremain*pvar%data4d(i+1,mesh%JMAX+j,k,l) + this%velocity_shift(l)
            END DO
          END DO
        END DO
      CASE(bottom)
       !CALL this%Error("SetBoundary", "Shearing not supported in BOTTOM direction.")
      CASE(top)
       !CALL this%Error("SetBoundary", "Shearing not supported in TOP direction.")
      END SELECT
    END DO
  END SUBROUTINE setboundarydata


  SUBROUTINE finalize(this)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(boundary_shearing), INTENT(INOUT) :: this
    !------------------------------------------------------------------------!
    DEALLOCATE(this%velocity_shift)
    CALL this%Finalize_base()
  END SUBROUTINE finalize



END MODULE boundary_shearing_mod