geometry_logspherical.f90

Go to the documentation of this file.

!#############################################################################
!#                                                                           #
!# fosite - 3D hydrodynamical simulation program                             #
!# module: geometry_logspherical.f90                                         #
!#                                                                           #
!# Copyright (C) 2007-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 geometry_logspherical_mod
  USE geometry_base_mod
  USE geometry_spherical_mod, ONLY: geometry_spherical
  USE common_dict
  IMPLICIT NONE
  !--------------------------------------------------------------------------!
  TYPE, EXTENDS(geometry_spherical) :: geometry_logspherical
  CONTAINS
    PROCEDURE :: initgeometry_logspherical
    PROCEDURE :: scalefactors_1
    PROCEDURE :: scalefactors_2
    PROCEDURE :: scalefactors_3
    PROCEDURE :: scalefactors_4
    PROCEDURE :: radius_1
    PROCEDURE :: radius_2
    PROCEDURE :: radius_3
    PROCEDURE :: radius_4
    PROCEDURE :: positionvector_1
    PROCEDURE :: positionvector_2
    PROCEDURE :: positionvector_3
    PROCEDURE :: positionvector_4
    PROCEDURE :: convert2cartesian_coords_1
    PROCEDURE :: convert2cartesian_coords_2
    PROCEDURE :: convert2cartesian_coords_3
    PROCEDURE :: convert2cartesian_coords_4
    PROCEDURE :: convert2curvilinear_coords_1
    PROCEDURE :: convert2curvilinear_coords_2
    PROCEDURE :: convert2curvilinear_coords_3
    PROCEDURE :: convert2curvilinear_coords_4
    PROCEDURE :: finalize
  END TYPE
  PRIVATE
  CHARACTER(LEN=32), PARAMETER :: geometry_name = "logspherical"
  !--------------------------------------------------------------------------!
  PUBLIC :: geometry_logspherical
  !--------------------------------------------------------------------------!

CONTAINS

  SUBROUTINE initgeometry_logspherical(this,config)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(geometry_logspherical), INTENT(INOUT) :: this
    TYPE(DICT_TYP),POINTER                   :: config
    !------------------------------------------------------------------------!
    REAL                                     :: gparam
    !------------------------------------------------------------------------!
    CALL getattr(config, "gparam", gparam, 1.0)

    CALL this%SetScale(gparam)
    CALL this%InitGeometry(logspherical,geometry_name,config)
    CALL this%SetAzimuthIndex(3)
  END SUBROUTINE initgeometry_logspherical

  PURE SUBROUTINE scalefactors_1(this,coords,hx,hy,hz)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(geometry_logspherical), INTENT(IN) :: this
    REAL, INTENT(IN),  DIMENSION(:,:) :: coords
    REAL, INTENT(OUT), DIMENSION(:)   :: hx,hy,hz
    !------------------------------------------------------------------------!
    CALL scalefactors(this%geoparam(1),coords(:,1),coords(:,2),coords(:,3), &
                      hx(:),hy(:),hz(:))
  END SUBROUTINE scalefactors_1

  PURE SUBROUTINE scalefactors_2(this,coords,hx,hy,hz)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(geometry_logspherical), INTENT(IN) :: this
    REAL, INTENT(IN),  DIMENSION(:,:,:) :: coords
    REAL, INTENT(OUT), DIMENSION(:,:)   :: hx,hy,hz
    !------------------------------------------------------------------------!
    CALL scalefactors(this%geoparam(1),coords(:,:,1),coords(:,:,2), &
                      coords(:,:,3),hx(:,:),hy(:,:),hz(:,:))
  END SUBROUTINE scalefactors_2

  PURE SUBROUTINE scalefactors_3(this,coords,hx,hy,hz)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(geometry_logspherical), INTENT(IN) :: this
    REAL, INTENT(IN),  DIMENSION(:,:,:,:) :: coords
    REAL, INTENT(OUT), DIMENSION(:,:,:)   :: hx,hy,hz
    !------------------------------------------------------------------------!
    CALL scalefactors(this%geoparam(1),coords(:,:,:,1),coords(:,:,:,2), &
                      coords(:,:,:,3),hx(:,:,:),hy(:,:,:),hz(:,:,:))
  END SUBROUTINE scalefactors_3

  PURE SUBROUTINE scalefactors_4(this,coords,hx,hy,hz)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(geometry_logspherical), INTENT(IN) :: this
    REAL, INTENT(IN),  DIMENSION(:,:,:,:,:) :: coords
    REAL, INTENT(OUT), DIMENSION(:,:,:,:)   :: hx,hy,hz
    !------------------------------------------------------------------------!
    CALL scalefactors(this%geoparam(1),coords(:,:,:,:,1),coords(:,:,:,:,2), &
                      coords(:,:,:,:,3),hx(:,:,:,:),hy(:,:,:,:),hz(:,:,:,:))
  END SUBROUTINE scalefactors_4

  PURE SUBROUTINE radius_1(this,coords,r)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(geometry_logspherical), INTENT(IN) :: this
    REAL, DIMENSION(:,:), INTENT(IN) :: coords
    REAL, DIMENSION(:),  INTENT(OUT) :: r
    !------------------------------------------------------------------------!
    r = radius(this%geoparam(1),coords(:,1))
  END SUBROUTINE radius_1

  PURE SUBROUTINE radius_2(this,coords,r)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(geometry_logspherical), INTENT(IN) :: this
    REAL, DIMENSION(:,:,:), INTENT(IN) :: coords
    REAL, DIMENSION(:,:),  INTENT(OUT) :: r
    !------------------------------------------------------------------------!
    r = radius(this%geoparam(1),coords(:,:,1))
  END SUBROUTINE radius_2

  PURE SUBROUTINE radius_3(this,coords,r)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(geometry_logspherical), INTENT(IN) :: this
    REAL, DIMENSION(:,:,:,:), INTENT(IN) :: coords
    REAL, DIMENSION(:,:,:),  INTENT(OUT) :: r
    !------------------------------------------------------------------------!
    r = radius(this%geoparam(1),coords(:,:,:,1))
  END SUBROUTINE radius_3

  PURE SUBROUTINE radius_4(this,coords,r)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(geometry_logspherical), INTENT(IN) :: this
    REAL, DIMENSION(:,:,:,:,:), INTENT(IN) :: coords
    REAL, DIMENSION(:,:,:,:),  INTENT(OUT) :: r
    !------------------------------------------------------------------------!
    r = radius(this%geoparam(1),coords(:,:,:,:,1))
  END SUBROUTINE radius_4

  PURE SUBROUTINE positionvector_1(this,coords,posvec)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(geometry_logspherical), INTENT(IN)  :: this
    REAL, DIMENSION(:,:), INTENT(IN)  :: coords
    REAL, DIMENSION(:,:), INTENT(OUT) :: posvec
    !------------------------------------------------------------------------!
    CALL positionvector(this%geoparam(1),coords(:,1),posvec(:,1), &
                        posvec(:,2),posvec(:,3))
  END SUBROUTINE positionvector_1

  PURE SUBROUTINE positionvector_2(this,coords,posvec)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(geometry_logspherical), INTENT(IN)  :: this
    REAL, DIMENSION(:,:,:), INTENT(IN)  :: coords
    REAL, DIMENSION(:,:,:), INTENT(OUT) :: posvec
    !------------------------------------------------------------------------!
    CALL positionvector(this%geoparam(1),coords(:,:,1),posvec(:,:,1), &
                        posvec(:,:,2),posvec(:,:,3))
  END SUBROUTINE positionvector_2

  PURE SUBROUTINE positionvector_3(this,coords,posvec)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(geometry_logspherical), INTENT(IN)  :: this
    REAL, DIMENSION(:,:,:,:), INTENT(IN)  :: coords
    REAL, DIMENSION(:,:,:,:), INTENT(OUT) :: posvec
    !------------------------------------------------------------------------!
    CALL positionvector(this%geoparam(1),coords(:,:,:,1),posvec(:,:,:,1), &
                        posvec(:,:,:,2),posvec(:,:,:,3))
  END SUBROUTINE positionvector_3

  PURE SUBROUTINE positionvector_4(this,coords,posvec)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(geometry_logspherical), INTENT(IN)  :: this
    REAL, DIMENSION(:,:,:,:,:), INTENT(IN)  :: coords
    REAL, DIMENSION(:,:,:,:,:), INTENT(OUT) :: posvec
    !------------------------------------------------------------------------!
    CALL positionvector(this%geoparam(1),coords(:,:,:,:,1),posvec(:,:,:,:,1), &
                        posvec(:,:,:,:,2),posvec(:,:,:,:,3))
  END SUBROUTINE positionvector_4

  PURE SUBROUTINE convert2cartesian_coords_1(this,curv,cart)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(geometry_logspherical), INTENT(IN) :: this
    REAL, DIMENSION(:,:), INTENT(IN)  :: curv
    REAL, DIMENSION(:,:), INTENT(OUT) :: cart
    !------------------------------------------------------------------------!
    CALL convert2cartesian_coords(this%geoparam(1),curv(:,1),curv(:,2),curv(:,3), &
                                  cart(:,1),cart(:,2),cart(:,3))
  END SUBROUTINE convert2cartesian_coords_1

  PURE SUBROUTINE convert2cartesian_coords_2(this,curv,cart)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(geometry_logspherical), INTENT(IN) :: this
    REAL, DIMENSION(:,:,:), INTENT(IN)  :: curv
    REAL, DIMENSION(:,:,:), INTENT(OUT) :: cart
    !------------------------------------------------------------------------!
    CALL convert2cartesian_coords(this%geoparam(1),curv(:,:,1),curv(:,:,2),curv(:,:,3), &
                                  cart(:,:,1),cart(:,:,2),cart(:,:,3))
  END SUBROUTINE convert2cartesian_coords_2

  PURE SUBROUTINE convert2cartesian_coords_3(this,curv,cart)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(geometry_logspherical), INTENT(IN) :: this
    REAL, DIMENSION(:,:,:,:), INTENT(IN)  :: curv
    REAL, DIMENSION(:,:,:,:), INTENT(OUT) :: cart
    !------------------------------------------------------------------------!
    CALL convert2cartesian_coords(this%geoparam(1),curv(:,:,:,1),curv(:,:,:,2),curv(:,:,:,3), &
                                  cart(:,:,:,1),cart(:,:,:,2),cart(:,:,:,3))
  END SUBROUTINE convert2cartesian_coords_3

  PURE SUBROUTINE convert2cartesian_coords_4(this,curv,cart)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(geometry_logspherical), INTENT(IN) :: this
    REAL, DIMENSION(:,:,:,:,:), INTENT(IN)  :: curv
    REAL, DIMENSION(:,:,:,:,:), INTENT(OUT) :: cart
    !------------------------------------------------------------------------!
    CALL convert2cartesian_coords(this%geoparam(1),curv(:,:,:,:,1),curv(:,:,:,:,2),curv(:,:,:,:,3), &
                                  cart(:,:,:,:,1),cart(:,:,:,:,2),cart(:,:,:,:,3))
  END SUBROUTINE convert2cartesian_coords_4

  PURE SUBROUTINE convert2curvilinear_coords_1(this,cart,curv)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(geometry_logspherical), INTENT(IN) :: this
    REAL, DIMENSION(:,:), INTENT(IN)  :: cart
    REAL, DIMENSION(:,:), INTENT(OUT) :: curv
    !------------------------------------------------------------------------!
    CALL convert2curvilinear_coords(this%geoparam(1),cart(:,1),cart(:,2),cart(:,3), &
                                    curv(:,1),curv(:,2),curv(:,3))
  END SUBROUTINE convert2curvilinear_coords_1

  PURE SUBROUTINE convert2curvilinear_coords_2(this,cart,curv)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(geometry_logspherical), INTENT(IN) :: this
    REAL, DIMENSION(:,:,:), INTENT(IN)  :: cart
    REAL, DIMENSION(:,:,:), INTENT(OUT) :: curv
    !------------------------------------------------------------------------!
    CALL convert2curvilinear_coords(this%geoparam(1),cart(:,:,1),cart(:,:,2),cart(:,:,3), &
                                    curv(:,:,1),curv(:,:,2),curv(:,:,3))
  END SUBROUTINE convert2curvilinear_coords_2

  PURE SUBROUTINE convert2curvilinear_coords_3(this,cart,curv)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(geometry_logspherical), INTENT(IN) :: this
    REAL, DIMENSION(:,:,:,:), INTENT(IN)  :: cart
    REAL, DIMENSION(:,:,:,:), INTENT(OUT) :: curv
    !------------------------------------------------------------------------!
    CALL convert2curvilinear_coords(this%geoparam(1),cart(:,:,:,1),cart(:,:,:,2),cart(:,:,:,3), &
                                    curv(:,:,:,1),curv(:,:,:,2),curv(:,:,:,3))
  END SUBROUTINE convert2curvilinear_coords_3

  PURE SUBROUTINE convert2curvilinear_coords_4(this,cart,curv)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(geometry_logspherical), INTENT(IN) :: this
    REAL, DIMENSION(:,:,:,:,:), INTENT(IN)  :: cart
    REAL, DIMENSION(:,:,:,:,:), INTENT(OUT) :: curv
    !------------------------------------------------------------------------!
    CALL convert2curvilinear_coords(this%geoparam(1),cart(:,:,:,:,1),cart(:,:,:,:,2),cart(:,:,:,:,3), &
                                    curv(:,:,:,:,1),curv(:,:,:,:,2),curv(:,:,:,:,3))
  END SUBROUTINE convert2curvilinear_coords_4

  ELEMENTAL SUBROUTINE convert2cartesian_vectors_0(this,xi,eta,phi,vxi,veta,vphi,vx,vy,vz)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(geometry_logspherical), INTENT(IN) :: this
    REAL, INTENT(IN)                      :: xi,eta,phi,vxi,veta,vphi
    REAL, INTENT(OUT)                     :: vx,vy,vz
    !------------------------------------------------------------------------!
    vx = vxi*sin(eta)*cos(phi) + veta*cos(eta)*cos(phi) - vphi*sin(phi)
    vy = vxi*sin(eta)*sin(phi) + veta*cos(eta)*sin(phi) + vphi*cos(phi)
    vz = vxi*cos(eta) - veta*sin(eta)
  END SUBROUTINE convert2cartesian_vectors_0

  ELEMENTAL SUBROUTINE convert2curvilinear_vectors_0(this,xi,eta,phi,vx,vy,vz,vxi,veta,vphi)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(geometry_logspherical), INTENT(IN) :: this
    REAL, INTENT(IN)                      :: xi,eta,phi,vx,vy,vz
    REAL, INTENT(OUT)                     :: vxi,veta,vphi
    !------------------------------------------------------------------------!
    vxi  = vx*sin(eta)*cos(phi) + vy*sin(eta)*sin(phi) + vz*cos(eta)
    veta = vx*cos(eta)*cos(phi) + vy*cos(eta)*sin(phi) - vz*sin(eta)
    vphi = -vx*sin(phi) + vy*cos(phi)
  END SUBROUTINE convert2curvilinear_vectors_0

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

  ELEMENTAL SUBROUTINE scalefactors(gp,logr,theta,phi,hlogr,htheta,hphi)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    REAL, INTENT(IN)  :: gp,logr,theta,phi
    REAL, INTENT(OUT) :: hlogr,htheta,hphi
    !------------------------------------------------------------------------!
    hlogr  = gp*exp(logr)
    htheta = hlogr
    hphi   = hlogr * sin(theta)
  END SUBROUTINE scalefactors

  ELEMENTAL FUNCTION radius(gp,logr)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    REAL, INTENT(IN)  :: gp,logr
    REAL :: radius
    !------------------------------------------------------------------------!
    radius = gp*exp(logr)
  END FUNCTION radius

  ELEMENTAL SUBROUTINE positionvector(gp,logr,rx,ry,rz)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    REAL, INTENT(IN)  :: gp,logr
    REAL, INTENT(OUT) :: rx,ry,rz
    !------------------------------------------------------------------------!
    rx = radius(gp,logr)
    ry = 0.0
    rz = 0.0
  END SUBROUTINE positionvector

  ELEMENTAL SUBROUTINE convert2cartesian_coords(gp,logr,theta,phi,x,y,z)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    REAL, INTENT(IN)  :: gp,logr,theta,phi
    REAL, INTENT(OUT) :: x,y,z
    !------------------------------------------------------------------------!
    REAL :: r
    !------------------------------------------------------------------------!
    r = gp*exp(logr)
    x = r*sin(theta)*cos(phi)
    y = r*sin(theta)*sin(phi)
    z = r*cos(theta)
  END SUBROUTINE convert2cartesian_coords

  ELEMENTAL SUBROUTINE convert2curvilinear_coords(gp,x,y,z,logr,theta,phi)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    REAL, INTENT(IN)  :: gp,x,y,z
    REAL, INTENT(OUT) :: logr,theta,phi
    !------------------------------------------------------------------------!
    REAL :: r
    !------------------------------------------------------------------------!
    r = sqrt(x*x+y*y+z*z)
    logr = log(r/gp)
    theta = acos(z/r)
    phi = atan2(y,x)
    IF(phi.LT.0.0) THEN
        phi = phi + 2.0*pi
    END IF
  END SUBROUTINE convert2curvilinear_coords

END MODULE geometry_logspherical_mod