diff --git a/Apps/Regrid_Util.F90 b/Apps/Regrid_Util.F90 index 8b4810c4d073..3e472d62d345 100644 --- a/Apps/Regrid_Util.F90 +++ b/Apps/Regrid_Util.F90 @@ -96,7 +96,7 @@ subroutine process_command_line(this,rc) this%lat_range=uninit this%shave=64 this%deflate=0 - this%quantize_algorithm=1 + this%quantize_algorithm=0 this%quantize_level=0 this%use_weights = .false. nargs = command_argument_count() diff --git a/CHANGELOG.md b/CHANGELOG.md index c7737c04e4fb..2b007a0a7842 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -37,6 +37,9 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 - Update `esma_add_fortran_submodules` function - Move MPI detection out of FindBaselibs - Add SMOD to submodule generator +- Add support for preliminary CF Conventions quantization properties + - Add new quantization keyword `granular_bitround` to History. This will be the preferred keyword for quantization in the future + replacing `GranularBR` ### Fixed @@ -47,6 +50,9 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ### Deprecated +- Deprecate `GranularBR` as a quantization method keyword in History. We will prefer `granular_bitround` in the future to match + draft CF conventions. This will be removed in MAPL 3. + ## [2.47.1] - 2024-07-17 ### Fixed diff --git a/gridcomps/History/MAPL_HistoryGridComp.F90 b/gridcomps/History/MAPL_HistoryGridComp.F90 index 4b048912bd70..7aa0afff2a1c 100644 --- a/gridcomps/History/MAPL_HistoryGridComp.F90 +++ b/gridcomps/History/MAPL_HistoryGridComp.F90 @@ -425,6 +425,7 @@ subroutine Initialize ( gc, import, dumexport, clock, rc ) logical :: has_conservative_keyword, has_regrid_keyword integer :: create_mode character(len=:), allocatable :: uppercase_algorithm + character(len=2) :: tmpchar ! Begin !------ @@ -838,40 +839,51 @@ subroutine Initialize ( gc, import, dumexport, clock, rc ) call ESMF_ConfigGetAttribute ( cfg, list(n)%quantize_algorithm_string, default='NONE', & label=trim(string) // 'quantize_algorithm:' ,_RC ) + call ESMF_ConfigGetAttribute ( cfg, list(n)%quantize_level, default=0, & + label=trim(string) // 'quantize_level:' ,_RC ) + ! Uppercase the algorithm string just to allow for any case + ! CF Conventions will prefer 'bitgroom', 'bitround', and 'granular_bitround' + ! but we will allow 'GranularBR' in MAPL2, deprecate it, and remove it in MAPL3 uppercase_algorithm = ESMF_UtilStringUpperCase(list(n)%quantize_algorithm_string,_RC) select case (trim(uppercase_algorithm)) case ('NONE') - list(n)%quantize_algorithm = MAPL_Quantize_Disabled + list(n)%quantize_algorithm = MAPL_NOQUANTIZE + ! If quantize_algorithm is 0, then quantize_level must be 0 + _ASSERT( list(n)%quantize_level == 0, 'quantize_algorithm is none, so quantize_level must be "none"') case ('BITGROOM') - list(n)%quantize_algorithm = MAPL_Quantize_BitGroom - case ('GRANULARBR') - list(n)%quantize_algorithm = MAPL_Quantize_GranularBR + list(n)%quantize_algorithm = MAPL_QUANTIZE_BITGROOM + case ('GRANULARBR', 'GRANULAR_BITROUND') + list(n)%quantize_algorithm = MAPL_QUANTIZE_GRANULAR_BITROUND case ('BITROUND') - list(n)%quantize_algorithm = MAPL_Quantize_BitRound + list(n)%quantize_algorithm = MAPL_QUANTIZE_BITROUND case default - _FAIL('Invalid quantize_algorithm. Allowed values are NONE, BitGroom, GranularBR, BitRound') + _FAIL('Invalid quantize_algorithm. Allowed values are none, bitgroom, granular_bitround, granularbr (deprecated), and bitround') end select - call ESMF_ConfigGetAttribute ( cfg, list(n)%quantize_level, default=0, & - label=trim(string) // 'quantize_level:' ,_RC ) - ! If nbits_to_keep < MAPL_NBITS_UPPER_LIMIT (24) and quantize_algorithm greater than 0, then a user might be doing different ! shaving algorithms. We do not allow this - _ASSERT( .not. ( (list(n)%nbits_to_keep < MAPL_NBITS_UPPER_LIMIT) .and. (list(n)%quantize_algorithm > 0) ), 'nbits < 24 and quantize_algorithm > 0 is not allowed. Choose one bit grooming method.') - - ! quantize_algorithm must be between 0 and 3 where 0 means not enabled - _ASSERT( (list(n)%quantize_algorithm >= 0) .and. (list(n)%quantize_algorithm <= 3), 'quantize_algorithm must be between 0 and 3, where 0 means not enabled') + _ASSERT( .not. ( (list(n)%nbits_to_keep < MAPL_NBITS_UPPER_LIMIT) .and. (list(n)%quantize_algorithm > MAPL_NOQUANTIZE) ), 'nbits < 24 and quantize_algorithm not "none" is not allowed. Choose a supported quantization method.') ! Now we test in the case that a valid quantize algorithm is chosen - if (list(n)%quantize_algorithm == 0) then - ! If quantize_algorithm is 0, then quantize_level must be 0 - _ASSERT( list(n)%quantize_level == 0, 'quantize_algorithm is 0, so quantize_level must be 0') - else + if (list(n)%quantize_algorithm /= MAPL_NOQUANTIZE) then ! If quantize_algorithm is greater than 0, then quantize_level must be greater than or equal to 0 _ASSERT( list(n)%quantize_level >= 0, 'netCDF quantize has been enabled, so quantize_level must be greater than or equal to 0') end if + ! If a user has chosen MAPL_QUANTIZE_BITROUND, then we allow a maximum of 23 bits to be kept + if (list(n)%quantize_algorithm == MAPL_QUANTIZE_BITROUND) then + write(tmpchar, '(I2)') MAPL_QUANTIZE_MAX_NSB + _ASSERT( list(n)%quantize_level <= MAPL_QUANTIZE_MAX_NSB, 'netCDF bitround has been enabled, so number of significant bits (quantize_level) must be less than or equal to ' // trim(tmpchar)) + end if + + ! For MAPL_QUANTIZE_GRANULAR_BITROUND and MAPL_QUANTIZE_BITGROOM, these use number of + ! significant digits, so for single precision, we allow a maximum of 7 digits to be kept + if (list(n)%quantize_algorithm == MAPL_QUANTIZE_GRANULAR_BITROUND .or. list(n)%quantize_algorithm == MAPL_QUANTIZE_BITGROOM) then + write(tmpchar, '(I2)') MAPL_QUANTIZE_MAX_NSD + _ASSERT( list(n)%quantize_level <= MAPL_QUANTIZE_MAX_NSD, 'netCDF granular bitround or bitgroom has been enabled, so number of significant digits (quantize_level) must be less than or equal to ' // trim(tmpchar)) + end if + tm_default = -1 call ESMF_ConfigGetAttribute ( cfg, list(n)%tm, default=tm_default, & label=trim(string) // 'tm:', _RC ) diff --git a/griddedio/GriddedIO.F90 b/griddedio/GriddedIO.F90 index fe2e6fb45be5..7eb3cf53c7f7 100644 --- a/griddedio/GriddedIO.F90 +++ b/griddedio/GriddedIO.F90 @@ -26,6 +26,7 @@ module MAPL_GriddedIOMod use, intrinsic :: ISO_C_BINDING use, intrinsic :: iso_fortran_env, only: REAL64 use ieee_arithmetic, only: isnan => ieee_is_nan + use netcdf, only: nf90_inq_libvers implicit none private @@ -51,7 +52,7 @@ module MAPL_GriddedIOMod type(VerticalData) :: vdata type(GriddedIOitemVector) :: items integer :: deflateLevel = 0 - integer :: quantizeAlgorithm = 1 + integer :: quantizeAlgorithm = MAPL_NOQUANTIZE integer :: quantizeLevel = 0 integer, allocatable :: chunking(:) logical :: itemOrderAlphabetical = .true. @@ -60,6 +61,7 @@ module MAPL_GriddedIOMod contains procedure :: CreateFileMetaData procedure :: CreateVariable + procedure :: CreateQuantizationInfo procedure :: modifyTime procedure :: modifyTimeIncrement procedure :: bundlePost @@ -190,7 +192,6 @@ subroutine CreateFileMetaData(this,items,bundle,timeInfo,vdata,ogrid,global_attr call this%timeInfo%add_time_to_metadata(this%metadata,rc=status) _VERIFY(status) - iter = this%items%begin() if (.not.allocated(this%chunking)) then call this%set_default_chunking(rc=status) _VERIFY(status) @@ -198,11 +199,16 @@ subroutine CreateFileMetaData(this,items,bundle,timeInfo,vdata,ogrid,global_attr call this%check_chunking(this%vdata%lm,_RC) end if - order = this%metadata%get_order(rc=status) _VERIFY(status) metadataVarsSize = order%size() + ! If quantize algorithm is set, create a quantization_info variable + if (this%quantizeAlgorithm /= MAPL_NOQUANTIZE) then + call this%CreateQuantizationInfo(_RC) + end if + + iter = this%items%begin() do while (iter /= this%items%end()) item => iter%get() if (item%itemType == ItemTypeScalar) then @@ -423,6 +429,31 @@ subroutine CreateVariable(this,itemName,rc) #else call v%add_attribute('regrid_method', regrid_method_int_to_string(this%regrid_method)) #endif + ! The CF Convention will soon support quantization. This requires three new attributes + ! if enabled: + ! 1. quantization --> Will point to a quantization_info container with the quantization algorithm + ! (NOTE: this will need to be programmatic when per-variable quantization is enabled) + ! 2a. quantization_nsb --> Number of significant bits (only for bitround) + ! 2b. quantization_nsd --> Number of significant digits (only for bitgroom and granular_bitround) + ! 3. quantization_maximum_relative_error --> Maximum relative error (defined as 2^(-nsb) for bitround, and UNDEFINED? for bitgroom and granular_bitround) + + ! Bitround + if (this%quantizeAlgorithm == MAPL_QUANTIZE_BITROUND) then + call v%add_attribute('quantization', 'quantization_info') + call v%add_attribute('quantization_nsb', this%quantizeLevel) + call v%add_attribute('quantization_maximum_relative_error', 0.5 * 2.0**(-this%quantizeLevel)) + end if + ! granular_bitround and bitgroom + if (this%quantizeAlgorithm == MAPL_QUANTIZE_BITGROOM .or. this%quantizeAlgorithm == MAPL_QUANTIZE_GRANULAR_BITROUND) then + call v%add_attribute('quantization', 'quantization_info') + call v%add_attribute('quantization_nsd', this%quantizeLevel) + ! Per czender, these have maximum_absolute_error. We use the calculate_mae function below + ! which replicates a table in doi 10.5194/gmd-12-4099-2019 + ! NOTE: This might not be the right formula. As the CF Convention draft is updated, + ! we will update this code. + call v%add_attribute('quantization_maximum_absolute_error', calculate_mae(this%quantizeLevel)) + end if + call factory%append_variable_metadata(v) call this%metadata%add_variable(trim(varName),v,rc=status) _VERIFY(status) @@ -442,6 +473,81 @@ subroutine CreateVariable(this,itemName,rc) end subroutine CreateVariable + function calculate_mae(nsd) result(mae) + + ! This function is based on Table 3 of doi 10.5194/gmd-12-4099-2019 + ! The algorithm is weird, but it does duplicate the table + + integer, intent(in) :: nsd + real(kind=REAL32) :: mae + real(kind=REAL32) :: mae_base + integer :: correction + + mae_base = 4.0 * (1.0/16.0)**floor(real(nsd)/2.0) * (1.0/8.0)**ceiling(real(nsd)/2.0) + + correction = 1 + if ( (nsd > 2 .and. mod(nsd, 2) == 0) .or. nsd == 7 ) then + correction = 2 + end if + + mae = mae_base * correction + end function calculate_mae + + subroutine CreateQuantizationInfo(this,rc) + class (MAPL_GriddedIO), intent(inout) :: this + integer, optional, intent(out) :: rc + + integer :: status + + class (AbstractGridFactory), pointer :: factory + character(len=:), allocatable :: varName, netcdf_version + type(Variable) :: v + + factory => get_factory(this%output_grid,_RC) + + v = Variable(type=PFIO_CHAR) + + ! In the future when we can do per variable quantization, we will need + ! to do things like quantization_info1, quantization_info2, etc. + ! For now, we will just use quantization_info as it is per collection + varName = "quantization_info" + + ! We need to convert the quantization algorithm to a string + select case (this%quantizeAlgorithm) + case (MAPL_QUANTIZE_BITGROOM) + call v%add_attribute('algorithm', 'bitgroom') + case (MAPL_QUANTIZE_BITROUND) + call v%add_attribute('algorithm', 'bitround') + case (MAPL_QUANTIZE_GRANULAR_BITROUND) + call v%add_attribute('algorithm', 'granular_bitround') + case default + _FAIL('Unknown quantization algorithm') + end select + + ! Next add the implementation details + ! 3. implementation: This property contains free-form text + ! that concisely conveys the algorithm provenance, including the + ! name of the library or client that performed the quantization, + ! the software version, and the name of the author(s) if deemed + ! relevant. + ! + ! In the current case, all algorithms are from libnetcdf + ! we make a string using nf90_inq_libvers() + + netcdf_version = 'libnetcdf ' // nf90_inq_libvers() + call v%add_attribute('implementation', netcdf_version) + + ! NOTE: In the future if we add the MAPL bit-shaving + ! to use the quantization parts of the code, it will + ! need a different implementation string + + call factory%append_variable_metadata(v) + call this%metadata%add_variable(trim(varName),v,_RC) + + _RETURN(_SUCCESS) + + end subroutine CreateQuantizationInfo + subroutine modifyTime(this, oClients, rc) class(MAPL_GriddedIO), intent(inout) :: this type (ClientManager), optional, intent(inout) :: oClients diff --git a/pfio/Variable.F90 b/pfio/Variable.F90 index 9d42cf97f7f2..624954bfd162 100644 --- a/pfio/Variable.F90 +++ b/pfio/Variable.F90 @@ -28,7 +28,7 @@ module pFIO_VariableMod type (StringAttributeMap) :: attributes type (UnlimitedEntity) :: const_value integer :: deflation = 0 ! default no compression - integer :: quantize_algorithm = 1 ! default bitgroom + integer :: quantize_algorithm = 0 ! default no quantization integer :: quantize_level = 0 ! default no quantize_level integer, allocatable :: chunksizes(:) contains @@ -85,7 +85,7 @@ function new_Variable(unusable, type, dimensions, chunksizes,const_value, deflat var%type = -1 var%deflation = 0 - var%quantize_algorithm = 1 + var%quantize_algorithm = 0 var%quantize_level = 0 var%chunksizes = empty var%dimensions = StringVector() diff --git a/shared/Constants/InternalConstants.F90 b/shared/Constants/InternalConstants.F90 index ac2935ea9911..3cad2914e86f 100644 --- a/shared/Constants/InternalConstants.F90 +++ b/shared/Constants/InternalConstants.F90 @@ -166,13 +166,17 @@ module MAPL_InternalConstantsMod integer, parameter :: MAPL_NBITS_NOT_SET = 1000 integer, parameter :: MAPL_NBITS_UPPER_LIMIT = 24 - ! Constants for netCDF quantize + ! Constants for netCDF quantize (these echo the values in the netcdf-fortran library) enum, bind(c) - enumerator MAPL_Quantize_Disabled - enumerator MAPL_Quantize_BitGroom - enumerator MAPL_Quantize_GranularBR - enumerator MAPL_Quantize_BitRound + enumerator MAPL_NOQUANTIZE + enumerator MAPL_QUANTIZE_BITGROOM + enumerator MAPL_QUANTIZE_GRANULAR_BITROUND + enumerator MAPL_QUANTIZE_BITROUND endenum + ! Maximum number of significant digits for quantization (bitgroom, granular_bitround) + integer, parameter :: MAPL_QUANTIZE_MAX_NSD = 7 + ! Maximum number of significant bits for quantization (bitround) + integer, parameter :: MAPL_QUANTIZE_MAX_NSB = 23 ! Constant masking enum, bind(c) enumerator MAPL_MASK_OUT