test_carre.f

Go to the documentation of this file.

C> @file tests/test_carre.f
C> Unit tests of subroutines from formats/carre.f and convert/carre.f
C GRIDMAN, grid managment library. Author: Vladislav Kotov, v.kotov@fz-juelich.de

!    Copyright (c) 2017  Forschungszentrum Juelich GmbH
!                        Vladislav Kotov 
!
!    This file is part of GRIDMAN.
!
!    GRIDMAN 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 3 of the License, or
!    (at your option) any later version.
!
!    GRIDMAN 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.  See the
!    GNU General Public License for more details.
!
!    You should have received a copy of the GNU General Public License
!    along with GRIDMAN.  If not, see <http://www.gnu.org/licenses/>.

      PROGRAM test_carre
      USE gridman
      USE gridman_lib
      IMPLICIT NONE

      INTERFACE
       SUBROUTINE check_mappings(LX,LY,NX,NY,BR,BZ,LXB,LYB,NXB,NYB,GRID)
        USE gridman
        IMPLICIT NONE
        INTEGER(GRIDMAN_SP),INTENT(IN) :: lx,ly,nx,ny
        INTEGER(GRIDMAN_SP),INTENT(IN) :: lxb,lyb,nxb,nyb
        REAL(GRIDMAN_DP) :: br(lxb:nxb,lyb:nyb,4),bz(lxb:nxb,lyb:nyb,4)
        TYPE(gridman_grid) :: grid
       END SUBROUTINE
      END INTERFACE

      TYPE(gridman_grid) :: fort30_grid,carre_grid,grid
      CHARACTER(256) :: stitle
      INTEGER(GRIDMAN_SP) :: nx,ny,nncut,nniso 
      INTEGER :: ierr,res
      INTEGER(GRIDMAN_SP),ALLOCATABLE :: 
     i                    nxcut30(:,:),nycut30(:,:),
     i                    nxiso30(:,:),nyiso30(:,:),
     i                    nxcut(:),nycut(:),nxiso(:),nyiso(:)
      REAL(GRIDMAN_DP),ALLOCATABLE :: br(:,:,:),bz(:,:,:),pit(:,:),
     r                                bc(:,:,:)
      REAL(GRIDMAN_DP) :: rbt

      gridman_dbg=.false. 
      gridman_check=.true.
C
C     GRID FROM FORT30 FILE
C

      CALL gridman_carre_read30('./input/fort.30',stitle,nx,ny,
     c                          nncut,nxcut30,nycut30,
     c                          nniso,nxiso30,nyiso30,br,bz,ierr)
      IF(ierr.NE.0) stop "TEST_CARRE TERMINATED"

      CALL gridman_carre2grid_fort30(fort30_grid,
     c                               1_gridman_sp,1_gridman_sp,
     c                               nx,ny,nncut,nxcut30,nycut30,
     c                               nniso,nxiso30,nyiso30,br,bz,
     c                               1_gridman_sp,1_gridman_sp,
     c                               nx,ny,ierr)
      IF(ierr.NE.0) stop "TEST_CARRE TERMINATED"

      CALL gridman_carre_read30(fort30_grid,'./input/fort.30',ierr)
      IF(ierr.NE.0) stop "TEST_CARRE TERMINATED"

      CALL gridman_grid_write(fort30_grid,'carre.grd',ierr)
      IF(ierr.NE.0) stop "TEST_CARRE TERMINATED"

      CALL gridman_grid_read(grid,'carre.grd',ierr)
      IF(ierr.NE.0) stop "TEST_CARRE TERMINATED"

      CALL gridman_grid_compare(fort30_grid,grid,res,ierr)
      IF(res.NE.0.OR.ierr.NE.0) stop "TEST_CARRE TERMINATED"

      CALL check_mappings(1_gridman_sp,1_gridman_sp,
     c                    nx,ny,br,bz,1_gridman_sp,1_gridman_sp,
     c                    nx,ny,fort30_grid)

C
C     GRID FROM SONNET FILE    
C
      CALL gridman_carre_readsonnet('./input/iterm.carre.105',nx,ny,
     s                         nncut,nxcut,nycut,nniso,nxiso,nyiso,
     s                         rbt,br,bz,pit,bc,ierr)
      IF(ierr.NE.0) stop "TEST_CARRE TERMINATED"

      CALL gridman_carre2grid(carre_grid,
     s                        0_gridman_sp,0_gridman_sp,
     s                        nx+1,ny+1,nncut,nxcut,nycut,
     s                        nniso,nxiso,nyiso,br,bz,
     s                        1_gridman_sp,1_gridman_sp,
     s                        nx,ny,ierr)
      IF(ierr.NE.0) stop "TEST_CARRE TERMINATED"

      CALL gridman_carre_readsonnet(carre_grid,
     s                             './input/iterm.carre.105',ierr)
      IF(ierr.NE.0) stop "TEST_CARRE TERMINATED"

      gridman_tol=1e-5
      carre_grid%DESCRIPTION=stitle
      CALL gridman_grid_compare(fort30_grid,carre_grid,res,ierr)
      IF(res.NE.0.OR.ierr.NE.0) stop "TEST_CARRE TERMINATED"
C
C    CHECK DIFFERENT TYPES OF INDEXING
C
      WRITE(*,*) "CHECK 1"
      CALL gridman_carre2grid(carre_grid,
     s                        0_gridman_sp,0_gridman_sp,
     s                        nx+1,ny+1,nncut,nxcut,nycut,
     s                        nniso,nxiso,nyiso,br,bz,
     s                        0_gridman_sp,0_gridman_sp,
     s                        nx+1,ny+1,ierr)
      IF(ierr.NE.0) stop "TEST_CARRE TERMINATED"
      CALL check_mappings(0_gridman_sp,0_gridman_sp,nx+1,ny+1,br,bz,
     c                  0_gridman_sp,0_gridman_sp, nx+1,ny+1,carre_grid)

      WRITE(*,*) "CHECK 2"
      CALL gridman_carre2grid(carre_grid,
     s                        0_gridman_sp,0_gridman_sp,
     s                        nx+1,ny+1,nncut,nxcut,nycut,
     s                        nniso,nxiso,nyiso,br,bz,
     s                        1_gridman_sp,0_gridman_sp,
     s                        nx,ny+1,ierr)
      IF(ierr.NE.0) stop "TEST_CARRE TERMINATED"
      CALL check_mappings(1_gridman_sp,0_gridman_sp,nx,ny+1,br,bz,
     c                   0_gridman_sp,0_gridman_sp,nx+1,ny+1,carre_grid)

      WRITE(*,*) "CHECK 3"
      CALL gridman_carre2grid(carre_grid,
     s                        0_gridman_sp,0_gridman_sp,
     s                        nx+1,ny+1,nncut,nxcut,nycut,
     s                        nniso,nxiso,nyiso,br,bz,
     s                        0_gridman_sp,1_gridman_sp,
     s                        nx+1,ny,ierr)
      IF(ierr.NE.0) stop "TEST_CARRE TERMINATED"
      CALL check_mappings(0_gridman_sp,1_gridman_sp,nx+1,ny,br,bz,
     c                   0_gridman_sp,0_gridman_sp,nx+1,ny+1,carre_grid)

      WRITE(*,*) "CHECK 4"
      CALL gridman_carre2grid(carre_grid,
     s                        0_gridman_sp,0_gridman_sp,
     s                        nx+1,ny+1,nncut,nxcut,nycut,
     s                        nniso,nxiso,nyiso,br,bz,
     s                        0_gridman_sp,2_gridman_sp,
     s                        nx+1,ny,ierr)
      IF(ierr.NE.0) stop "TEST_CARRE TERMINATED"
      CALL check_mappings(0_gridman_sp,2_gridman_sp,nx+1,ny,br,bz,
     c                   0_gridman_sp,0_gridman_sp,nx+1,ny+1,carre_grid)

      WRITE(*,*) "CHECK 5"
      CALL gridman_carre2grid(carre_grid,
     s                        0_gridman_sp,0_gridman_sp,
     s                        nx+1,ny+1,nncut,nxcut,nycut,
     s                        nniso,nxiso,nyiso,br,bz,
     s                        1_gridman_sp,2_gridman_sp,
     s                        nx-1,ny,ierr)
      IF(ierr.NE.0) stop "TEST_CARRE TERMINATED"
      CALL check_mappings(1_gridman_sp,2_gridman_sp,nx-1,ny,br,bz,
     c                   0_gridman_sp,0_gridman_sp,nx+1,ny+1,carre_grid)
    
      WRITE(*,*) "CHECK 6"
C
C     GRID FROM FORT30 FILE PRODUCED BY SOLPS5
C     
      CALL gridman_carre_read30(fort30_grid,'./input/fort.30.b25',ierr)
      IF(ierr.NE.0) stop "TEST_CARRE TERMINATED"

C     DEALLOCATION
      CALL gridman_grid_deallocate(carre_grid,ierr)
      IF(ierr.NE.0) stop "TEST_CARRE TERMINATED"
      CALL gridman_grid_deallocate(grid,ierr)
      IF(ierr.NE.0) stop "TEST_CARRE TERMINATED"

      WRITE(*,*) "TEST_CARRE COMPLETED"

      END PROGRAM test_carre

C
C
C
      SUBROUTINE check_mappings(LX,LY,NX,NY,BR,BZ,LXB,LYB,NXB,NYB,GRID)
      USE gridman
      USE gridman_lib
      IMPLICIT NONE
      INTEGER(GRIDMAN_SP),INTENT(IN) :: lx,ly,nx,ny
      INTEGER(GRIDMAN_SP),INTENT(IN) :: lxb,lyb,nxb,nyb
      REAL(GRIDMAN_DP) :: br(lxb:nxb,lyb:nyb,4),bz(lxb:nxb,lyb:nyb,4)
      TYPE(gridman_grid) :: grid
      REAL(GRIDMAN_DP) :: volb2(lx:nx,ly:ny),vol(grid%ncells),
     r                    sx(lx-1:nx,ly:ny),sy(lx:nx,ly-1:ny),
     r                    sarea(grid%NEDGES)

      INTEGER(GRIDMAN_SP) :: ix,iy,ie,ix1,ix2,iy1,iy2,icell,iedge,m
      INTEGER :: ierr
      REAL(GRIDMAN_DP) :: s213,s143,r213,r143,x1,x2,y1,y2
            
C     REGION WHERE VOLUMES AND SURFACES WILL BE COMPARED
      ix1=lx+3
      ix2=nx-3
      iy1=ly+3
      iy2=ny-3

      ix1=lx
      ix2=nx
      iy1=ly
      iy2=ny

C     VOLUME OF B2 GRID
      DO iy=ly,ny
       DO ix=lx,nx
        s213=br(ix,iy,2)*(bz(ix,iy,1)-bz(ix,iy,3))+
     +       br(ix,iy,1)*(bz(ix,iy,3)-bz(ix,iy,2))+
     +       br(ix,iy,3)*(bz(ix,iy,2)-bz(ix,iy,1))
        s143=br(ix,iy,1)*(bz(ix,iy,4)-bz(ix,iy,3))+
     +       br(ix,iy,4)*(bz(ix,iy,3)-bz(ix,iy,1))+
     +       br(ix,iy,3)*(bz(ix,iy,1)-bz(ix,iy,4))
        r213=br(ix,iy,2)+br(ix,iy,1)+br(ix,iy,3)
        r143=br(ix,iy,1)+br(ix,iy,4)+br(ix,iy,3)
       volb2(ix,iy)=gridman_pi*abs((r213*s213+r143*s143))/3.  
       END DO
      END DO
     
C     VOLUME OF THE GRID CELLS
      CALL gridman_grid2d_cylvolumes(grid,vol,ierr)
      IF(ierr.NE.0) stop "TEST_CARRE TERMINATED"
     
      IF(grid%NCELLINDEX.NE.1) THEN
       WRITE(*,*) "ERORR: wrong number of indices"
       WRITE(*,*) " NCELLINDEX ",grid%NCELLINDEX
      END IF

      IF(grid%CELLINDEX(1)%NELEMENTS.NE.grid%NCELLS) THEN
       WRITE(*,*) "ERORR: wrong cell index"
       WRITE(*,*) " NELEMENTS, NCELLS ",
     w             grid%CELLINDEX(1)%NELEMENTS,grid%NCELLS
       stop "TEST_CARRE TERMINATED"
      END IF

      DO ie=1,grid%CELLINDEX(1)%NELEMENTS
       icell=grid%CELLINDEX(1)%INDEXES(0,ie)
       IF(ie.NE.icell) THEN 
        WRITE(*,*) "ERORR: wrong cell index"
        WRITE(*,*) " index mismatch"
        WRITE(*,*) " IE, ICELL ",ie, icell
        stop "TEST_CARRE TERMINATED"
       END IF
       ix=grid%CELLINDEX(1)%INDEXES(1,ie)
       iy=grid%CELLINDEX(1)%INDEXES(2,ie)
       IF(ix.LT.lx.OR.ix.GT.nx.OR.
     c    iy.LT.ly.OR.iy.GT.ny) THEN
        WRITE(*,*) "ERORR: cell index out of range"
        WRITE(*,*) " ICELL ",icell
        WRITE(*,*) " IX, LX, NX ",ix,lx,nx
        WRITE(*,*) " IY, LY, NY ",iy,ly,ny
        stop "TEST_CARRE TERMINATED"
       END IF 
       IF(ix.GE.ix1.AND.ix.LE.ix2.AND.
     .    iy.GE.iy1.AND.iy.LE.iy2) THEN
        IF(abs(vol(icell)-volb2(ix,iy)).GT.
     .     gridman_tol*(abs(vol(icell))+abs(volb2(ix,iy)))) THEN 
         WRITE(*,*) "ERORR: wrong cell index"
         WRITE(*,*) " volume mismatch between ",
     w              "original and mapped cells"
         WRITE(*,*) "ICELL, IX, IY, VOL, VOLB2"
         WRITE(*,*)  icell,ix,iy,vol(icell),volb2(ix,iy)
         stop "TEST_CARRE TERMINATED"
        END IF
       END IF 
      END DO
         
C     SURFACE AREAS OF B2 GRID
      DO iy=ly,ny
       DO ix=lx,nx
        x1=br(ix,iy,1)
        y1=bz(ix,iy,1)
        x2=br(ix,iy,2)
        y2=bz(ix,iy,2)
        sx(ix,iy)=gridman_pi*(x1+x2)*sqrt((x1-x2)**2+(y1-y2)**2)
        x1=br(ix,iy,3)
        y1=bz(ix,iy,3)
        x2=br(ix,iy,2)
        y2=bz(ix,iy,2)
        sy(ix,iy)=gridman_pi*(x1+x2)*sqrt((x1-x2)**2+(y1-y2)**2)
       END DO
      END DO
      ix=lx
      DO iy=ly,ny
       x1=br(ix,iy,3)
       y1=bz(ix,iy,3)
       x2=br(ix,iy,4)
       y2=bz(ix,iy,4)
       sx(ix-1,iy)=gridman_pi*(x1+x2)*sqrt((x1-x2)**2+(y1-y2)**2)
      END DO
      iy=ly
      DO ix=lx,nx
       x1=br(ix,iy,4)
       y1=bz(ix,iy,4)
       x2=br(ix,iy,1)
       y2=bz(ix,iy,1)
       sy(ix,iy-1)=gridman_pi*(x1+x2)*sqrt((x1-x2)**2+(y1-y2)**2)
      END DO

C     TOTAL SURFACE OF GRID
      CALL gridman_grid2d_cylareas(grid,sarea,ierr)
      IF(ierr.NE.0) stop "TEST_CARRE TERMINATED"

C     COMPARING ORIGINAL AND MAPPED SURFACE AREAS
      IF(grid%NEDGEINDEX.NE.3) THEN
       WRITE(*,*) "ERORR: wrong number of indices"
       WRITE(*,*) " NEDGEINDEX ",grid%NEDGEINDEX
       stop "TEST_CARRE TERMINATED"
      END IF

      m=2*(nx-lx+1)+2*(ny-ly+1)
      IF(grid%EDGEINDEX(1)%NINDEX.NE.1.OR.
     f   grid%EDGEINDEX(1)%NELEMENTS.NE.m) THEN 
       WRITE(*,*) "ERORR: wrong 1st edge index"
       WRITE(*,*) " NINDEX, NELEMENTS, M ",
     w              grid%EDGEINDEX(1)%NINDEX,
     w              grid%EDGEINDEX(1)%NELEMENTS,m
       stop "TEST_CARRE TERMINATED"
      END IF

      m=(nx-lx+1)*(ny-ly+2)
C     Poloidal edges 
      IF(grid%EDGEINDEX(2)%NELEMENTS.NE.m) THEN
       WRITE(*,*) "ERORR: wrong 1st edge index"
       WRITE(*,*) " NELEMENTS, NX*(NY+1) ", 
     w              grid%EDGEINDEX(1)%NELEMENTS,m
       stop "TEST_CARRE TERMINATED"
      END IF
      DO ie=1,grid%EDGEINDEX(2)%NELEMENTS
       iedge=grid%EDGEINDEX(2)%INDEXES(0,ie)
       ix=grid%EDGEINDEX(2)%INDEXES(1,ie)
       iy=grid%EDGEINDEX(2)%INDEXES(2,ie)
       IF(iedge.LT.1.OR.iedge.GT.grid%NEDGES.OR.
     f    ix.LT.lx.OR.ix.GT.nx.OR.iy.LT.ly-1.OR.iy.GT.ny) THEN
        WRITE(*,*) "ERORR: poloidal edge index is out of range"
        WRITE(*,*) " IEDGE, NEDGES ",iedge,grid%NEDGES
        WRITE(*,*) " IX, LX, NX ",ix,lx,nx
        WRITE(*,*) " IY, LY-1, NY ",iy,ly-1,ny
        stop "TEST_CARRE TERMINATED"
       END IF        
       IF(ix.GE.ix1.AND.ix.LE.ix2.AND.
     .    iy.GE.ly-1.AND.iy.LE.ny) THEN
        IF(abs(sarea(iedge)-sy(ix,iy)).GT.
     .     gridman_tol*(abs(sarea(iedge))+abs(sy(ix,iy)))) THEN 
         WRITE(*,*) "ERORR: surface area mismatch between ",
     w              "original and mapped Y-surfaces"
         WRITE(*,*) "IEDGE, IX, IY, SAREA, SY"
         WRITE(*,*)  iedge,ix,iy,sarea(iedge),sy(ix,iy)
         stop "TEST_CARRE TERMINATED"
        END IF
       END IF 
      END DO
C     Radial edges
      m=(nx-lx+2)*(ny-ly+1)
      IF(grid%EDGEINDEX(3)%NELEMENTS.NE.m) THEN
       WRITE(*,*) "ERORR: wrong 2nd edge index"
       WRITE(*,*) " NELEMENTS, (NX-1)*NY ", 
     w              grid%EDGEINDEX(2)%NELEMENTS,m
       stop "TEST_CARRE TERMINATED"
      END IF
      DO ie=1,grid%EDGEINDEX(3)%NELEMENTS
       iedge=grid%EDGEINDEX(3)%INDEXES(0,ie)
       ix=grid%EDGEINDEX(3)%INDEXES(1,ie)
       iy=grid%EDGEINDEX(3)%INDEXES(2,ie)
       IF(iedge.LT.1.OR.iedge.GT.grid%NEDGES.OR.
     f    ix.LT.lx-1.OR.ix.GT.nx.OR.iy.LT.ly.OR.iy.GT.ny) THEN
        WRITE(*,*) "ERORR: poloidal edge index is out of range"
        WRITE(*,*) " IEDGE, NEDGES ",iedge,grid%NEDGES
        WRITE(*,*) " IX, LX-1, NX ",ix,lx-1,nx
        WRITE(*,*) " IY, LY, NY ",iy,ly,ny
        stop "TEST_CARRE TERMINATED"
       END IF      
       IF(ix.GE.lx-1.AND.ix.LE.nx.AND.
     .    iy.GE.iy1.AND.iy.LE.iy2) THEN
        IF(abs(sarea(iedge)-sx(ix,iy)).GT.
     .     gridman_tol*(abs(sarea(iedge))+abs(sx(ix,iy)))) THEN 
         WRITE(*,*) "ERORR: surface area mismatch between ",
     w              "original and mapped X-surfaces"
         WRITE(*,*) "IEDGE, IX, IY, SAREA, SX"
         WRITE(*,*)  iedge,ix,iy,sarea(iedge),sx(ix,iy)
         stop "TEST_CARRE TERMINATED"
        END IF
       END IF 
      END DO

      END SUBROUTINE check_mappings