XIOS reduction bug

check.py

#module load python-compute pc-numpy pc-netcdf4-python

import numpy as np
from netCDF4 import Dataset
import glob

class TwoFile(object):
    def __init__(self, basename):
        self.filenames = [basename+'_0.nc', basename+'_1.nc']
        self.datasets = [Dataset(fn) for fn in self.filenames]
        
    def variable(self, name):
        vars = [ds.variables[name] for ds in self.datasets]
        dims = vars[0].dimensions
        try:
            concat_dim = dims.index(u'lat')
        except ValueError:
            # not a fancy dimensions
            return vars[0]
        
        return np.concatenate(vars, axis=concat_dim)

data = TwoFile('cube')
# first time index only...
temp = data.variable('temp')[0]

# Recall axes are ordered time, z, lat lon
x = 3600*np.arange(60)[np.newaxis, np.newaxis, :]
y = 60*np.arange(60)[np.newaxis, :, np.newaxis]
z = np.arange(60)[:, np.newaxis, np.newaxis]
comp = x+y+z

print "Full output OK?", np.all(comp==temp)

# first time index only...
zmin_file = data.variable('zmin')[0]
zmean_file = data.variable('zmean')[0]
zmax_file = data.variable('zmax')[0]

zmin_computed = temp.min(axis=0)
zmean_computed = temp.mean(axis=0)
zmax_computed = temp.max(axis=0)

print "Min OK?", np.all(zmin_computed == zmin_file)
print "Mean OK?", np.all(zmean_computed == zmean_file)
print "Max OK?", np.all(zmax_computed == zmax_file)

cube.f90

PROGRAM cube
  USE xios

  IMPLICIT NONE


  INCLUDE "mpif.h"
  INTEGER :: comm, rank, size, ierr
  INTEGER :: cart_comm, ndims = 3, size3d(3) = 2, rank3d(3), periodic(3) = 0
  
  INTEGER :: ts
  TYPE(xios_duration) :: dtime

  ! global domain sizes
  INTEGER :: gloDomSize(3) = 60
  ! local domain size and offset
  INTEGER :: locDomSize(3), locDomStart(3), locDomEnd(3)
  
  ! local data
  DOUBLE PRECISION, ALLOCATABLE :: l_x(:), l_y(:), l_z(:), l_temp(:,:,:)
  ! Indexing
  INTEGER :: li, lj, lk, gi, gj, gk
  
  ! XIOS initialization
  CALL xios_initialize("cube", return_comm=comm)
  
  ! Figure out our decomposition etc
  CALL MPI_Comm_rank(comm, rank, ierr)
  CALL MPI_Comm_size(comm, size, ierr)
  ! Assume this will run on 8 ranks
  IF (size3d(1)*size3d(2)*size3d(3) .ne. size) THEN
     CALL MPI_Abort(comm,1)
  END IF

  CALL MPI_Cart_create(comm, ndims, size3d, periodic, 1, cart_comm, ierr)
  CALL MPI_Cart_get(cart_comm, ndims, size3d, periodic, rank3d, ierr)

  locDomSize = gloDomSize / size3d
  locDomStart = rank3d*locDomSize
  locDomEnd = locDomStart + locDomSize
  
  ALLOCATE(l_x(0:locDomSize(1)-1))
  ALLOCATE(l_y(0:locDomSize(2)-1))
  ALLOCATE(l_z(0:locDomSize(3)-1))
  ALLOCATE(l_temp(0:locDomSize(1)-1, 0:locDomSize(2)-1, 0:locDomSize(3)-1))

  DO li = 0, locDomSize(1)
     gi = locDomStart(1) + li
     l_x(li) = gi
  END DO
  DO lj = 0, locDomSize(2)
     gj = locDomStart(2) + lj
     l_y(lj) = gj
  END DO
  DO lk = 0, locDomSize(3)
     gk = locDomStart(3) + lk
     l_z(lk) = gk
  END DO

  DO li = 0, locDomSize(1)-1
     gi = locDomStart(1) + li
     DO lj = 0, locDomSize(2)-1
        gj = locDomStart(2) + lj
        DO lk = 0, locDomSize(3)-1
           gk = locDomStart(3) + lk
           l_temp(li, lj, lk) = (gi*gloDomSize(2) + gj) * gloDomSize(3) + gk
        END DO
     END DO
  END DO
  
  CALL xios_context_initialize("cube", comm)
  ! Set the sizes of the grid
  CALL xios_set_domain_attr("xy", ni_glo=gloDomSize(1), ibegin=locDomStart(1), ni=locDomSize(1))
  CALL xios_set_domain_attr("xy", nj_glo=gloDomSize(2), jbegin=locDomStart(2), nj=locDomSize(2))
  CALL xios_set_axis_attr("z", n_glo=gloDomSize(3), begin=locDomStart(3), n=locDomSize(3))

  ! Setting time step
  dtime%second=3600
  CALL xios_set_timestep(dtime)
  ! Closing definition
  CALL xios_close_context_definition()
  
  ! Entering time loop
  DO ts=1,96
     CALL xios_update_calendar(ts)
     CALL xios_send_field("temp", l_temp)
  ENDDO
  
  ! XIOS finalization
  CALL xios_context_finalize()
  CALL xios_finalize()
  
END PROGRAM cube

iodef.xml

<?xml version="1.0"?>
<simulation>
  <context id="cube">
    <calendar type="Gregorian"
	      start_date="2012-03-01 15:00:00"
	      time_origin="2012-02-28 15:00:00 + 1d" />
    <axis_definition>
      <axis id="z" />
    </axis_definition>

    <domain_definition>
      <domain id="xy" type="rectilinear"/>
    </domain_definition>

    <grid_definition>
      <grid id="threeD">
	<domain domain_ref="xy" />
	<axis axis_ref="z" />
      </grid>

      <grid id="zmin">
	<domain domain_ref="xy" />
	<scalar id="zmin">
	  <reduce_axis operation="min" />
	</scalar>
      </grid>
      <grid id="zmean">
	<domain domain_ref="xy" />
	<scalar id="zmean">
	  <reduce_axis operation="average" />
	</scalar>
      </grid>
      <grid id="zmax">
	<domain domain_ref="xy" />
	<scalar id="zmax">
	  <reduce_axis operation="max" />
	</scalar>
      </grid>
    </grid_definition>
    
    <field_definition>
      <field id="temp"
	     grid_ref="threeD"
	     operation="instant" />
    </field_definition>
    
    <file_definition type="multiple_file" output_freq="1d">
      <file id="output" name="cube" output_freq="1d">
	<field field_ref="temp" />
	<field name="zmin" field_ref="temp" grid_ref="zmin"
	       />
	<field name="zmean" field_ref="temp" grid_ref="zmean"
	       />
	<field name="zmax" field_ref="temp" grid_ref="zmax"
	       />
      </file>
    </file_definition>
    
  </context>
    
  <context id="xios">
    <variable_definition>
      <variable_group id="buffer">
        <variable id="optimal_buffer_size" type="string">performance</variable>
        <variable id="buffer_size_factor" type="double">1.0</variable>
      </variable_group>

      <variable_group id="parameters" >
        <variable id="using_server" type="bool">false</variable>
        <variable id="info_level" type="int">50</variable>
        <variable id="print_file" type="bool">true</variable>
      </variable_group>
    </variable_definition>
  </context>

</simulation>

Makefile

XIOS_BASE=$(HOME)/xios
XIOS_VERSION=trunk
XIOS_DIR = $(XIOS_BASE)/$(XIOS_VERSION)

NAME = cube

FC = ftn
FFLAGS = -I$(XIOS_DIR)/inc -D__NONE__ -O0 -g -traceback -I /opt/cray/netcdf-hdf5parallel/4.4.1.1/INTEL/15.0/include

$(NAME).exe : $(NAME).o
	$(FC) -nofor-main -o $@ $^ -L$(XIOS_DIR)/lib -lxios -Wl,"--allow-multiple-definition" -Wl,"-Bstatic" -L /opt/cray/netcdf-hdf5parallel/4.4.1.1/INTEL/15.0/lib -L /opt/cray/netcdf-hdf5parallel/4.4.1.1/INTEL/15.0/lib  -lnetcdf -lnetcdff -lhdf5_hl -lhdf5 -lz -lstdc++

$(NAME).o : $(NAME).f90
	$(FC) -o $@ $(FFLAGS) -c $<


.PHONY: clean
clean :
	-rm $(NAME).o $(NAME).exe

test.sh

#!/bin/bash --login
#
#PBS -l select=2
#PBS -q short
#PBS -l walltime=00:10:00
#PBS -A n02-cms
#PBS -N test_cube
#PBS -W umask=0022
#PBS -j oe

export PBS_O_WORKDIR=$(readlink -f $PBS_O_WORKDIR)
cd $PBS_O_WORKDIR

export OMP_NUM_THREADS=1

XIOS_BUILD_DIR=/home/n02/n02/rnashn02/xios/trunk
export PATH=$XIOS_BUILD_DIR/bin:$PATH

rm -f xios_server.exe
ln -s $XIOS_BUILD_DIR/bin/xios_server.exe

aprun  -n 8 -N 8 ./cube.exe  : -n 2 -N 2 ./xios_server.exe