Merging devel into master.

Several bugs fixes included:
* Duplicate quantity codes corrected and output tables have been
  modified accordingly.   Affected code ranges include:
   1900--2000 (kinetic energy equation codes)
   1400--1600 (thermal energy equation codes)
   1600--1800 (magnetic energy equation codes)

* Several logical issues that produced seg faults at output time
  for certain quantity-code combinations have been corrected.

* Configure script now accepts the following options:
  --with-custom
  -devel

  (run ./configure --help for details)

commit 351a04f
Author: Nick Featherstone <[email protected]>
Date:   Sat Apr 28 13:58:31 2018 -0600

    Squashed commit of the following:

    commit 0307b3f
    Author: Nick Featherstone <[email protected]>
    Date:   Sat Apr 28 13:57:38 2018 -0600

        Updated copy command for custom directory in Makefile

    commit 983ea62
    Author: Nick Featherstone <[email protected]>
    Date:   Sat Apr 28 09:21:30 2018 -0600

        removed production_shear_pke from DiagnosticsTurbKE_Budget.F90.
        It is already calculated in diagnostics_KE_Flux.F90

    commit e67ef25
    Author: Nick Featherstone <[email protected]>
    Date:   Sat Apr 28 09:06:33 2018 -0600

        Fixed logical issue affecting terms involving v' dot grad v'

    commit 17c5e3b
    Author: Nick Featherstone <[email protected]>
    Date:   Tue Apr 24 15:34:30 2018 -0600

        Using indices.F in the main loop now

    commit 3f75296
    Author: Nick Featherstone <[email protected]>
    Date:   Tue Apr 24 15:30:24 2018 -0600

        All Diagnostics now include indices.F for their preprocessor-macro definitions.
        Indices.F is stored in new directory src/Includes
        Makefile has been updated accordingly.
        Create dev_build script

    commit 8b3d556
    Author: Nick Featherstone <[email protected]>
    Date:   Tue Apr 24 08:29:19 2018 -0600

        created indices.f (not used yet)

commit c5e1212
Author: Nick Featherstone <[email protected]>
Date:   Mon Apr 23 15:26:39 2018 -0600

    Small  modification to Makefile so that build does not fail in the event of specified,
    but empty custom directory.

commit dfc5729
Author: Ryan Orvedahl <[email protected]>
Date:   Mon Apr 23 15:14:38 2018 -0600

    add Custom directory option to configure with the --with-custom=* flag (#9)

    Incorporating --with-custom flag to configure script (work by Ryan Orvedahl).

commit 5bb736b
Merge: 7541ab1 2a80aa7
Author: Nick Featherstone <[email protected]>
Date:   Mon Apr 9 13:52:47 2018 -0600

    Merge branch 'devel' of github.com:geodynamics/Rayleigh into devel

commit 7541ab1
Author: Nick Featherstone <[email protected]>
Date:   Mon Apr 9 13:52:28 2018 -0600

    Configure now accepts the -devel option.

commit 2a80aa7
Author: Nick Featherstone <[email protected]>
Date:   Thu Apr 5 20:14:27 2018 +0000

    harmless edit (two carriage returns)

commit 1e237e0
Merge: fa0e45d c141780
Author: Nick Featherstone <[email protected]>
Date:   Wed Apr 4 15:41:15 2018 -0600

    Merge pull request #8 from geodynamics/unstable

    Bug Fixes

commit c141780
Author: Nick Featherstone <[email protected]>
Date:   Wed Apr 4 14:56:03 2018 -0600

    make distclean now removes Makefile and object_list

commit 2c74bf1
Author: Nick Featherstone <[email protected]>
Date:   Wed Apr 4 14:47:29 2018 -0600

    viscous fluxes have been corrected.
    Rayleigh output quantities updated

commit 15cb7e1
Author: Nick Featherstone <[email protected]>
Date:   Wed Apr 4 13:33:05 2018 -0600

    removed nr_old write

commit 6329a71
Author: Nick Featherstone <[email protected]>
Date:   Wed Apr 4 13:26:13 2018 -0600

    --Fixed mira instructions in INSTALL
    --viscous flux code is updated, but quite working

commit e5c45d8
Author: Nick Featherstone <[email protected]>
Date:   Tue Mar 27 11:03:14 2018 -0600

    -Cleared several double-booked quantity codes.
    -Added consistency check to build_menu_pdf.py to help catch inconsistencies in the future.
    -"make distclean" now clears .o and .mod files as well

commit 8f5ab10
Author: Nick Featherstone <[email protected]>
Date:   Tue Mar 27 07:57:44 2018 -0600

    The maxq parameter is:
    (1) now module-wide in rayleigh_diagnostics.py
    (2) set to 4000 to maintain consistency with Spherical_IO.F90

commit 5570b39
Author: Nick Featherstone <[email protected]>
Date:   Tue Mar 27 07:17:23 2018 -0600

    testing branch synchronization

commit 3a35bdf
Author: Nick Featherstone <[email protected]>
Date:   Tue Mar 27 07:12:46 2018 -0600

    removed .unstable file

commit 1e4e80c
Merge: 961a16e 8449e98
Author: Nick Featherstone <[email protected]>
Date:   Tue Mar 27 07:12:10 2018 -0600

    Merge branch 'unstable' of github.com:geodynamics/Rayleigh into unstable

commit 961a16e
Author: Nick Featherstone <[email protected]>
Date:   Tue Mar 27 07:00:02 2018 -0600

    added .unstable file

commit fa0e45d
Merge: 84d5e37 8449e98
Author: Nick Featherstone <[email protected]>
Date:   Mon Feb 5 14:35:23 2018 -0700

    Merging unstable changes into devel:
    -FFT Renormalization (optimization)
    -configure script now supports AVX 512

commit 8449e98
Author: Nick Featherstone <[email protected]>
Date:   Mon Feb 5 14:33:35 2018 -0700

    Configure script now supports AVX-512 vectorization for Intel processors

commit 72f33b7
Author: Nick Featherstone <[email protected]>
Date:   Tue Jan 30 12:07:21 2018 -0700

    FFT renormalization weights are wrapped into Legendre integration weights

commit 84d5e37
Author: Nick Featherstone <[email protected]>
Date:   Sat Jan 27 14:49:35 2018 -0700

    Devel branch -- initial commit

INSTALL

@@ -187,13 +187,13 @@ The example below shows how to configure the code for Mira, an IBM Blue Gene/Q.
 export BLAS_LIB=/soft/libraries/alcf/current/xl/BLAS/lib
 export LAPACK_LIB=/soft/libraries/alcf/current/xl/LAPACK/lib
 export FFTW_LIB=/soft/libraries/alcf/current/xl/FFTW3/lib
-export FFTW_INC=/soft/libraries/alcf/current/xl/FFTW3/inc
+export FFTW_INC=/soft/libraries/alcf/current/xl/FFTW3/include
 
 ./configure \
 --FC=mpixlf2003_r \
 --CC=gcc \
 --LIBFLAGS="-L$FFTW_LIB -lfftw3 -L$LAPACK_LIB -llapack -L$BLAS_LIB -lblas" \
---INCLUDE="-I$FFTW_INC" \
+--INCLUDE="-I . -I$FFTW_INC" \
 --FFLAGS_OPT="-qfree=f90 -O2" \
 --FFLAGS_DBG="-qfree=f90 -O0 -g -traceback  -qcheck  -C" \
 --nodirs

Makefile

@@ -1,6 +1,10 @@
 include make.inc
 SRC=src
 BUILD=$(SRC)/build
+
+# make the CUSTOMROOT variable available to sub-make processes
+export CUSTOMROOT
+
 rayleigh:
 	@mkdir -p $(BUILD)/compiled
 	@cp $(SRC)/parallel_framework/*.F90 $(BUILD)/.
@@ -14,10 +18,16 @@ rayleigh:
 	@cp $(SRC)/Diagnostics/*.F $(BUILD)/.
 	@cp $(SRC)/Utility/*.F90 $(BUILD)/.
 	@cp $(SRC)/Utility/*.c $(BUILD)/.
+	@cp $(SRC)/Include/*.F $(BUILD)/.
 	@cp $(SRC)/Makefile $(BUILD)/.
 	@cp $(SRC)/object_list $(BUILD)/.
 ifeq ($(NODIRS),1)
 	cp $(SRC)/Utility/MakeDir.F90_IBM $(BUILD)/MakeDir.F90
+endif
+ifdef CUSTOMROOT
+	@echo Custom directory specified.
+	@echo Any files from $(CUSTOMROOT) will overwrite standard Rayleigh source files.
+	@cp $(CUSTOMROOT)/* $(BUILD)/. 2>/dev/null || :
 endif
 	@$(MAKE) --no-print-directory --directory=$(BUILD) clean_exec
 	@$(MAKE) --no-print-directory --directory=$(BUILD) all
@@ -42,6 +52,10 @@ distclean:
 	rm -f $(BUILD)/*.F
 	rm -f $(BUILD)/*.F90
 	rm -f $(BUILD)/*.c
+	rm -f $(BUILD)/*.o
+	rm -f $(BUILD)/*.mod
+	rm -f $(BUILD)/object_list
+	rm -f $(BUILD)/Makefile
 	rm -f $(BUILD)/Machine_Definitions
 	rm -f $(BUILD)/machine.no_comments
 	@rm -f $(PREFIX)/rayleigh.*

README.md

@@ -5,6 +5,8 @@
 
 Rayleigh is a 3-D convection code designed for the study of dynamo behavior in spherical geometry.  It evolves the incompressible and anelastic MHD equations in spherical geometry using a pseudo-spectral approach.  Rayleigh employs spherical harmonics in the horizontal direction and Chebyshev polynomials in the radial direction.  The code has undergone extensive accuracy testing using the [Christensen et al. (2001)](http://adsabs.harvard.edu/abs/2001PEPI..128...25C) Boussinesq benchmarks and the [Jones et al. (2011)](http://adsabs.harvard.edu/abs/2011Icar..216..120J) anelastic benchmarks.   Rayleigh has been developed with NSF support through the Computational Infrastructure for Geodynamics ([CIG](https://geodynamics.org/cig/news/newsletters/may-2016/)).
 
+
+
 Parallelization
 -----------------------------
 The pseudo-spectral nature of Rayleigh means that its parallelization necessarily relies heavily on *global communication* patterns.  That said, Rayleigh's parallelization is based around a 2-D domain decomposition and large-message-size all-to-alls.  These features allow the code to overcome many of the obstacles that traditionally limit the scalability of spectral methods.   The end result is a pseudo-spectral code optimized for petascale machines.  Rayleigh's pure-MPI mode has demonstrated highly efficient strong scaling on  131,000 cores of the Mira Blue Gene/Q supercomputer for problems with approximately 2048^3 grid points (2048 spherical harmonics).  Performance numbers from Mira are shown below.  A summary of Rayleigh's performance and how it compares against other popular dynamo codes (albeit at at smaller process counts) may be found in the recent performance benchmark results of [Matsui et al. (2016)](http://onlinelibrary.wiley.com/doi/10.1002/2015GC006159/full).

configure

@@ -96,6 +96,11 @@ function displayhelp {
     echo "      The C-compiler command DOES NOT need to be MPI-capable."
     echo ""
 
+    echo "  --devel / -devel "
+    echo "      Forces Rayleigh to build with debugging flags and no."
+    echo "      vectorization (useful for quick development builds)."
+    echo ""
+
     echo "  --FC=<fortran compiler name> "
     echo "      Specifies the Fortran compiler command to be invoked by Make."
     echo "      FC must be in your PATH, or else configure will fail.  FC must"
@@ -233,6 +238,11 @@ function displayhelp {
     echo "          <LAPACK ROOT>/lib/liblapack.a"
     echo ""
 
+    echo "  --with-custom=<CUSTOM ROOT>"
+    echo "      Set this flag to specify the location of custom Rayleigh source files."
+    echo "      These files will be copied to the build directory and override any of"
+    echo "      the standard Rayleigh source files."
+    echo ""
 
 }
 
@@ -245,6 +255,7 @@ RAYLEIGHROOT=$PWD
 PREFIX=$RAYLEIGHROOT
 NEED_HELP=no
 USE_MKL="FALSE"
+CUSTOMROOT="none"
 
 
 ###############################################################################################
@@ -260,6 +271,10 @@ do
       exit 0;
       ;;
 
+    --with-custom=*)
+      CUSTOMROOT=`expr "x$option" : "x-*with-custom=\(.*\)"`;
+      ;;
+
     -prefix=* | --prefix=*)
       PREFIX=`expr "x$option" : "x-*prefix=\(.*\)"`
       ;;
@@ -298,6 +313,10 @@ do
       NODIRS="TRUE";
       ;;
 
+    --devel | -devel)
+      DEVELMODE="TRUE";
+      ;;
+
     -mkl)
       if [ ! -z "$MKLROOT" ]
       then
@@ -372,6 +391,12 @@ then
       then
         echo "NODIRS = 1" >> make.inc
       fi
+      if [[ $CUSTOMROOT == "none" ]]
+      then
+        echo "CUSTOMROOT = " >> make.inc
+      else
+        echo "CUSTOMROOT = "$CUSTOMROOT >> make.inc
+      fi
       exit 0
 
 fi
@@ -508,7 +533,18 @@ FFLAGS[1]=$DBG_FLAGS
 VECLEVEL[1]=""
 nvers=2
 MULTIVEC=FALSE
-if [[ $FVERSION == "INTEL" ]]
+
+if [[ $DEVELMODE == "TRUE" ]]
+then
+  echo "************************************************************************"
+  echo " Development mode has been activated.  Both rayleigh.opt"
+  echo " and rayleigh.dbg will be compiled with debugging flags "
+  echo " and no vectorization."
+  FFLAGS[0]=$DBG_FLAGS
+fi
+
+
+if [[ $FVERSION == "INTEL" ]] &&  [[ $DEVELMODE != "TRUE" ]]
 then
 
   echo "***********************************************************************"
@@ -518,8 +554,9 @@ then
   echo " -- Westmere cores                  :  SSE4.2"
   echo " -- Sandy Bridge / Ivy Bridge cores :  AVX"
   echo " -- Haswell /Broadwell              :  AVX2 "
+  echo " -- SkyLake / Knights Landing       :  AVX512 "
   echo ""
-    OPTIONS="SSE4.2 AVX AVX2 ALL NONE"
+    OPTIONS="SSE4.2 AVX AVX2 AVX512 ALL NONE"
     select opt in $OPTIONS; do
        if [ "$opt" = "SSE4.2" ]; then
         FFLAGS[0]=$OPT_FLAGS" -xSSE4.2"
@@ -530,18 +567,24 @@ then
        elif [ "$opt" = "AVX2" ]; then
         FFLAGS[0]=$OPT_FLAGS" -xCORE-AVX2"
         break            
+       elif [ "$opt" = "AVX512" ]; then
+        FFLAGS[0]=$OPT_FLAGS" -xCORE-AVX512"
+        break      
        elif [ "$opt" = "ALL" ]; then
          DEFAULT='avx2'
-         nvers=4
+         nvers=5
          MULTIVEC=TRUE
          FFLAGS[0]=$OPT_FLAGS" -xCORE-AVX2"
          FFLAGS[2]=$OPT_FLAGS" -xAVX"
          FFLAGS[3]=$OPT_FLAGS" -xSSE4.2"
+         FFLAGS[4]=$OPT_FLAGS" -xCORE-AVX512"
          RVALUES[0]=$DEFAULT 
          RVALUES[2]=avx
          RVALUES[3]=sse
+         RVAKYES[4]=avx512
          RVTAGS[2]="AVX"
          RVTAGS[3]="SSE"
+         RVTAGS[4]="AVX512"
          break
        elif [ "$opt" == "NONE" ]; then
          FFLAGS[0]=$OPT_FLAGS
@@ -552,6 +595,7 @@ then
        fi
     done
   echo "***********************************************************************"
+
 fi
 
 
@@ -736,6 +780,13 @@ if [[ $NODIRS == "TRUE" ]]
 then
   echo "NODIRS = 1" >> make.inc
 fi
+if [[ $CUSTOMROOT == "none" ]]
+then
+  echo "CUSTOMROOT = " >> make.inc
+else
+  echo "CUSTOMROOT = "$CUSTOMROOT >> make.inc
+fi
+
 echo " ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~"
 echo " Configuration is complete."
 echo " A Machine_Definitions files has been placed in $RAYLEIGHROOT/src/build."

dev_build

@@ -0,0 +1,3 @@
+make distclean
+./configure -devel
+make

etc/analysis/rayleigh_diagnostics.py

@@ -1,7 +1,7 @@
 from __future__ import print_function
 import numpy as np
 import os
-
+maxq = 4000
 class RayleighTiming:
 
     def __init__(self,filename,byteswap=True):
@@ -203,7 +203,7 @@ def __init__(self,filename='none',path='G_Avgs/'):
             self.vals[i,:] = tmp
             self.time[i] = swapread(fd,dtype='float64',count=1,swap=bs)
             self.iters[i] = swapread(fd,dtype='int32',count=1,swap=bs)
-        maxq = 801
+
         lut = np.zeros(maxq)+int(1000)
         self.lut = lut.astype('int32')
         for i,q in enumerate(self.qv):
@@ -285,7 +285,7 @@ def __init__(self,filename='none',path='Shell_Avgs/',ntheta=0):
                 print('This ShellAverage file is version 2, and ntheta was not provided.')
                 print('The 2nd, 3rd and 4th moments are set to zero')   
                 self.vals[:,1,:,:] = 0.0            
-        maxq = 3000
+
         lut = np.zeros(maxq)+int(1000)
         self.lut = lut.astype('int32')
         for i,q in enumerate(self.qv):
@@ -348,7 +348,7 @@ def __init__(self,filename='none',path='AZ_Avgs/'):
             self.vals[:,:,:,i] = tmp
             self.time[i] = swapread(fd,dtype='float64',count=1,swap=bs)
             self.iters[i] = swapread(fd,dtype='int32',count=1,swap=bs)
-        maxq = 801
+
         lut = np.zeros(maxq)+int(1000)
         self.lut = lut.astype('int32')
         for i,q in enumerate(self.qv):
@@ -445,7 +445,7 @@ def __init__(self,filename='none',path='Point_Probes/'):
 
         #print 'iters: ', self.iters
         #print 'times: ', self.time
-        maxq = 4000
+
         lut = np.zeros(maxq)+int(1000)
         self.lut = lut.astype('int32')
         for i,q in enumerate(self.qv):
@@ -521,7 +521,7 @@ def __init__(self,filename='none',path='Meridional_Slices/'):
             self.vals[:,:,:,:,i] = tmp
             self.time[i] = swapread(fd,dtype='float64',count=1,swap=bs)
             self.iters[i] = swapread(fd,dtype='int32',count=1,swap=bs)
-        maxq = 801
+
         lut = np.zeros(maxq)+int(1000)
         self.lut = lut.astype('int32')
         for i,q in enumerate(self.qv):
@@ -588,7 +588,7 @@ def __init__(self,filename='none',path='Equatorial_Slices/'):
             self.vals[:,:,:,i] = tmp
             self.time[i] = swapread(fd,dtype='float64',count=1,swap=bs)
             self.iters[i] = swapread(fd,dtype='int32',count=1,swap=bs)
-        maxq = 801
+
         lut = np.zeros(maxq)+int(1000)
         self.lut = lut.astype('int32')
         for i,q in enumerate(self.qv):
@@ -673,7 +673,7 @@ def __init__(self,filename='none',path='Shell_Slices/',slice_spec = [], rec0 = F
 
 
         qv = np.reshape(swapread(fd,dtype='int32',count=nq,swap=bs),(nq), order = 'F')
-        maxq = 801
+
         lut_max = 1000
         lut = np.zeros(maxq)+int(lut_max)
         self.lut = lut.astype('int32')
@@ -879,7 +879,7 @@ def __init__(self,filename='none',path='SPH_Modes/'):
                         for m in range(1,nm):
                             self.vals[m,:,j,q,k] = self.vals[m,:,j,q,k]/sqrttwo
 
-        maxq = 801
+
         lut = np.zeros(maxq)+int(1000)
         self.lut = lut.astype('int32')
         for i,q in enumerate(self.qv):
@@ -993,7 +993,7 @@ def __init__(self,filename='none',path='Shell_Spectra/'):
                         for m in range(1,nm):
                             self.vals[:,m,j,q,k] = self.vals[:,m,j,q,k]/sqrttwo
 
-        maxq = 801
+
         lut = np.zeros(maxq)+int(1000)
         self.lut = lut.astype('int32')
         for i,q in enumerate(self.qv):

etc/diagnostic_codes/Rayleigh_Output_Variables.log

@@ -1,4 +1,4 @@
-This is pdfTeX, Version 3.14159265-2.6-1.40.16 (TeX Live 2015/Debian) (preloaded format=pdflatex 2017.8.29)  23 JAN 2018 15:17
+This is pdfTeX, Version 3.14159265-2.6-1.40.16 (TeX Live 2015/Debian) (preloaded format=pdflatex 2017.8.29)  4 APR 2018 14:21
 entering extended mode
  restricted \write18 enabled.
  %&-line parsing enabled.
@@ -377,7 +377,7 @@ sfonts/cm/cmsy10.pfb></usr/share/texlive/texmf-dist/fonts/type1/public/amsfonts
 /cm/cmsy5.pfb></usr/share/texlive/texmf-dist/fonts/type1/public/amsfonts/cm/cms
 y7.pfb></usr/share/texlive/texmf-dist/fonts/type1/public/amsfonts/cm/cmti10.pfb
 >
-Output written on Rayleigh_Output_Variables.pdf (33 pages, 323670 bytes).
+Output written on Rayleigh_Output_Variables.pdf (33 pages, 323636 bytes).
 PDF statistics:
  187 PDF objects out of 1000 (max. 8388607)
  131 compressed objects within 2 object streams

etc/diagnostic_codes/build_menu_pdf.py

@@ -1,3 +1,5 @@
+nqmax = 4000  # The maximum possible quantity code.  Should be consistent with the number in SphericalIO.F90
+glob_echecks = (nqmax+1)*[0]  # used to check for duplicates across all quantity codes
 def classify_line(codeline):
     """Description:  Classifies a line of code as offset-code, quantity-code, or other
        Inputs:  
@@ -108,6 +110,7 @@ def parse_offset_line(codeline):
     qcodes = []
     texcodes=[]
 
+    loc_echecks = (nqmax+1)*[0]  # Used to check for duplicated indices within each .F file 
     for aline in lines:
         ltype=classify_line(aline)
         if (ltype == 1):
@@ -119,14 +122,20 @@ def parse_offset_line(codeline):
             varnames.append(vname)
             qcodes.append(qcode+offset)
             texcodes.append(tcode)
-
+            loc_echecks[int(qcode+offset)]+=1
+            glob_echecks[int(qcode+offset)]+=1
         if (ltype == 2):
             pieces=parse_offset_line(aline)
             vname=pieces[0]
             roffset=pieces[1]
             relvname=pieces[2]
             mydict[vname]=mydict[relvname]+roffset
             offset=mydict[vname]
+
+    for ecvi, ecv in enumerate(loc_echecks):
+        if (ecv > 1):
+            print("Error in ", infile, ' qcode duplicated: ', ecvi) 
+
     tfcnt=0
     ofile=inpref[0]+'_table'+str(tfcnt)+'.tex'
     myfile= open(ofile,"w")
@@ -172,3 +181,8 @@ def parse_offset_line(codeline):
                 myfile= open(ofile,"w")
     #print(len(varnames))
     myfile.close()
+
+#Finally, perform a consistency check across all quantity codes parsed
+for ecvi, ecv in enumerate(glob_echecks):
+    if (ecv > 1):
+        print('Global Error -- qcode duplicated: ', ecvi)