src/examples/itm_grid_example10_solpsdata.f90

Go to the documentation of this file.

program itm_grid_example10_solpsdata
  use itm_types         ! ITM standard Fortran data types, e.g. real(R8)
  use itm_constants
  use euITM_schemas     ! ITM data structure/CPO definitions
  use euITM_routines    ! ITM UAL routines
  use itm_grid_common   ! ITM grid service library: constant definitions like COORDTYPE_R, COORDTYPE_Z

  use itm_grid_2dpolygon
  use itm_grid_subgrid
  use itm_grid_structured

  implicit none

  ! parameters
  integer, parameter :: NPOINTR = 6 
  integer, parameter :: NPOINTZ = 5   
  integer, parameter :: SHOTNUM_IN = 63254
  integer, parameter :: RUNNUM_IN = 1

  integer, parameter :: SHOTNUM_OUT = SHOTNUM_IN
  integer, parameter :: RUNNUM_OUT = RUNNUM_IN + 1

  ! variables
  type (type_edge) :: edgecpo_in
  type (type_edge), pointer :: edgecpo_out(:)
  integer :: idx

  integer :: nodeSg, cellSg

  type(GridObject) ::  cell
  type(GridObject),dimension(:),allocatable :: nodes
  real(R8), dimension(:,:), allocatable :: coords

  real(R8) :: r(NPOINTR), z(NPOINTZ), nodePolygonIndex(NPOINTR, NPOINTZ), nodeData(NPOINTR, NPOINTZ)
  real(R8) :: rmin, rmax, zmin, zmax, p(2)
  integer :: rFrom, rTo, zFrom, zTo, iNode, ir, iz, i, iCell, iCellData, iNodeData

  ! === 1. Read CPO ===

  write (*,*) "Example 10: reading from shot ", SHOTNUM_IN, ", run ", RUNNUM_IN
  !call euitm_open_env( 'euitm', SHOTNUM_IN, RUNNUM_IN, idx, 'coster', 'aug', '4.09a')
  call euitm_open_env( 'euitm', SHOTNUM_IN, RUNNUM_IN, idx, 'klingshi', 'jet', '4.09a')
  call euitm_get_slice( idx, "edge", edgecpo_in, 2.5d0, UAL_PREVIOUS_SAMPLE ) 
  call euitm_close(idx)

  ! figure out subgrid ids for nodes and cells
  nodeSg = gridFindSubgridByName(edgecpo_in%grid, "Nodes")
  cellSg = gridFindSubgridByName(edgecpo_in%grid, "Cells")

  ! Find ne entry for cell subgrid
  do i = 1, size(edgecpo_in%fluid%ne%value)
      if (edgecpo_in%fluid%ne%value(i)%subgrid == cellSg) then
          iCellData = i
      end if
      if (edgecpo_in%fluid%ne%value(i)%subgrid == nodeSg) then
          iNodeData = i
      end if
  end do

  ! === 2. Set up an output CPO with a structured grid ===
  allocate(edgecpo_out(1))

  rmin = 2.7_R8
  rmax = 2.85_R8
  zmin = -1.6_R8
  zmax = -1.5_R8

      r = (/ ( rmin + ((rmax - rmin)/(NPOINTR-1)) * i, i=0,NPOINTR-1) /)
  z = (/ ( zmin + ((zmax - zmin)/(NPOINTZ-1)) * i, i=0,NPOINTZ-1) /)

  call gridSetupStructuredSep( edgecpo_out(1)%grid, &
      & ndim = 2, &
      & c1 = COORDTYPE_R, &
      & x1 = r, &
      & c2 = COORDTYPE_Z, &
      & x2 = z,&
      & id = '2d_structured' )

  ! === 3. Interpolate data to it ===

  if (.true.) then

  nodeData = 0.0_R8
  nodePolygonIndex = 0.0_R8


  ! Consider every polygon in the grid...
  do iCell = 1, gridSubGridSize( edgecpo_in%grid%subgrids(cellSg) )

      ! Get nodes bounding the polygon
      cell = subGridGetObject(edgecpo_in%grid%subgrids(cellSg), iCell)
      call get2dPolygonNodes( edgecpo_in%grid, cell, nodes )
      
      !call writeObjectList(nodes)

      ! get the coordinates of the nodes
      if (allocated(coords)) deallocate(coords)
      allocate( coords(size(nodes), gridNDim(edgecpo_in%grid)) )
      do iNode = 1, size(nodes)
          !write (*,*) gridNodeCoord( edgecpo_in%grid, nodes(iNode)%ind )
          coords(iNode,:) = gridNodeCoord( edgecpo_in%grid, nodes(iNode)%ind )

      end do

      ! find bounding rectangle for the polygon
      rmin = minval(coords(:,1))
      rmax = maxval(coords(:,1))
      zmin = minval(coords(:,2))
      zmax = maxval(coords(:,2))

      ! if no overlap with target grid, cycle
      if (rmin > r(NPOINTR) .or. rmax < r(1) .or. zmin > z(NPOINTZ) .or. zmax < z(1)) then
          cycle
      end if

      !write (*,*) "Cell ", icell, ": box ", rmin, rmax, zmin, zmax


      ! cut away parts outside of the target grid
      !if (rmin < r(1)) rmin = r(1)
      rmin = max(rmin, r(1))
      rmax = min(rmax, r(NPOINTR))
      zmin = max(zmin, z(1))
      zmax = min(zmax, z(NPOINTZ))
      !write (*,*) "Cell ", icell, ": overlap box ", rmin, rmax, zmin, zmax

      ! Find index ranges for candidate points. Extend it a bit to make sure we really cover the polygon
      rFrom = max( cellIndex( r, rmin ) - 1, 1 )
      rTo = min( cellIndex( r, rmax ) + 1, NPOINTR )
 zFrom = max( cellIndex( z, zmin ) - 1, 1 )
      zTo = min( cellIndex( z, zmax ) + 1, NPOINTZ )

      !write (*,*) "Cell ", icell, ": grid box ", r(rFrom), r(rTo), z(zFrom), z(zTo), rFrom, rTo, zFrom, zTo

      if (rFrom == 0 .or. rTo == 0 .or. zFrom == 0 .or. zTo == 0) then          
          write (*,*) "WARNING: did not find cell index. This should not happen."
          write (*,*) rmin, rmax, zmin, zmax
          write (*,*) rFrom, rTo, zFrom, zTo
          cycle
      end if

      ! check for every candidate point whether it's inside the polygon
      do ir = rFrom, rTo
          do iz = zFrom, zTo
              p(1) = r(ir)
              p(2) = z(iz)
              if ( pointInPolygon( coords, p ) ) then
                  if (nodePolygonIndex(ir, iz) > 0.0_R8) then
                      write (*,*) "WARNING: found multiple cells matching the same node"
                      write (*,*) "Node ", ir, iz, p(1), p(2)
                      write (*,*) "Cell: ", iCell
                      write (*,*) "Cell x: ", coords(:,1)
                      write (*,*) "Cell y: ", coords(:,2)
                  else

                      nodePolygonIndex(ir, iz) = iCell

                      ! just take cell value...
                      !nodeData(ir, iz) = edgecpo_in%fluid%ne%value(iCellData)%scalar(iCell)

                      ! ...or do arithmetic average of four corner nodes
                      nodeData(ir, iz) = 0.0_R8
                      do iNode = 1, 4
                          nodeData(ir, iz) = nodeData(ir, iz) + &
                              & edgecpo_in%fluid%ne%value(iNodeData)%scalar(&
                              & subGridGetIndexForObject(edgecpo_in%grid%subgrids(nodeSg), nodes(iNode)) )
                      end do
                      nodeData(ir, iz) = nodeData(ir, iz) / 4.0_R8



                  end if
              end if
          end do
      end do
      
      deallocate(nodes)
  end do

  write (*,*) "Of ", NPOINTZ * NPOINTR, " points, found for ", count(nodeData > 0.0_R8), " a polygon"

  endif

  ! === 4. Write output CPO ===
  allocate( edgecpo_out(1)%datainfo%dataprovider(1) )
  edgecpo_out(1)%datainfo%dataprovider="IMP3"
  allocate( edgecpo_out(1)%codeparam%codename(1) )
  edgecpo_out(1)%codeparam%codename(1)="itm_grid_example10"
  edgecpo_out(1)%time = edgecpo_in%time

  allocate(edgecpo_out(1)%fluid%ne%value(1))
  call gridStructWriteData2d( edgecpo_out(1)%grid, edgecpo_out(1)%fluid%ne%value(1), GRID_STRUCT_NODES, nodeData )
  allocate(edgecpo_out(1)%fluid%te%value(1))
  call gridStructWriteData2d( edgecpo_out(1)%grid, edgecpo_out(1)%fluid%te%value(1), GRID_STRUCT_NODES, nodePolygonIndex )


  write (*,*) "Example 10: writing to shot ", SHOTNUM_OUT, ", run ", RUNNUM_OUT
  call euitm_create( 'euitm', SHOTNUM_OUT, RUNNUM_OUT, 0, 0, idx )
  call euitm_put(idx, "edge", edgecpo_out)
  call euitm_close(idx)


contains

  integer function cellIndex( x, val )
    real(R8), intent(in) :: x(:), val

    ! internal
    integer :: i

    do i = 1, size(x) - 2
        if (x(i) <= val .and. x(i+1) > val) then
            cellIndex = i
            return
        end if
    end do

    ! special treatment for rightmost cell
    if (x(size(x)-1) <= val .and. x(size(x)) >= val) then        
        cellIndex = size(x)-1
    else
        ! no cell found
        cellIndex = 0
    end if

  end function cellIndex




end program itm_grid_example10_solpsdata