Fix ITP & improve performance

[iangrooms]

Description:

This PR does two main things. First, the ITP algorithm in the rootfinding mod is fixed. The original version was based on Eq (15) from their paper, which was a typo. The buggy version still converged but was slow.

Second, several changes have been made to improve performance.

• In the rootfinding mod, for bounded cases we don't immediately go to ITP, which was inefficient.
• In the kde distribution mod we now have the option to use fifth, seventh, or (the original) ninth-order Gaussian quadrature. Lower orders could degrade accuracy for small ensembles, but significantly improve performance (nearly 40% speedup for fifth order compared to ninth, since most of the time is spent in this function). This is controlled with an optional namelise kde_nml which has a single entry quadrature_order equal to 5, 7, or 9. Documentation has been updated to reflect this option.
• In the kde distribution mod the logic associated with the boundary correction has been updated to reduce unnecessary computation. This doesn't change answers but significantly improves performance.
• I removed an unused function integrate_pdf. Since this function was part of the test suite I also changed the test suite to test the quadrature approximation to the cdf. The expected number of 'PASS' results is the same.

Types of changes

• Bug fix (non-breaking change which fixes an issue)
• New feature (non-breaking change which adds functionality)
• Breaking change (fix or feature that would cause existing functionality to not work as expected)
• Documentation update
  Performance improvement isn't on this list, but that's a big part of this PR.

Documentation changes needed?

• My change requires a change to the documentation.
  • I have updated the documentation accordingly.

Tests

I ran the kde distribution test, which passes with both the rootfinding mod version of inv_cdf and the normal distribution mod version of inv_cdf. I tested performance by just looping that test 4000 times.

Checklist for merging

• Updated changelog entry
• Documentation updated
• Update conf.py

Checklist for release

• Merge into main
• Create release from the main branch with appropriate tag
• Delete feature-branch

Testing Datasets

• Dataset needed for testing available upon request
• Dataset download instructions included
• No dataset needed

[iangrooms]

Right now I've switched the quadrature from fifth to third order, but left the fifth order code in there, inactive. It would be nice to control this using a namelist parameter, but I'm not sure how to do that for a module. Is there a way to read a namelist variable directly into this module? If so, then I can set up the logic to use third or fifth order based on that variable.

[hkershaw-brown]

For a namelist in a module, you'll need to add a check for the namelist being read for any routine that needs the namlist option. DART follows this pattern, where there is an initializing routine (e.g. "rootfinding_init") called to read the namelist on the first call. Pseudo code:

! module storage to indicate module has been initialized or not
logical :: module_initialized = .false.

namelist / rootfinding_nml / your_option_name

---
subroutine rootfinding_init

module_initialized = .true.

! Read the namelist entry
call find_namelist_in_file("input.nml", "rootfinding_nml", iunit) ! <- this will kill the program if rootfinding_nml is not in the input.nml
read(iunit, nml = probit_transform_nml, iostat = io)
call check_namelist_read(iunit, io, "rootfinding_nml")

! These log the namelist values
if (do_nml_file()) write(nmlfileunit,nml=rootfiniding_nml)
if (do_nml_term()) write(     *     ,nml=rootfinding_nml)


end subroutine  rootfinding_init
---

---
function integrate_pdf(x, p) result(q)

! If not initialized, read in the namelist - for routines that depend on the namelist
if(.not. module_initialized) call initialize_rootfiniding

! use "your_option_name" however you need in the routine
...
---

You can take a look at

DART/assimilation_code/modules/assimilation/probit_transform_mod.f90

Line 57 in d908222

namelist /probit_transform_nml/ fix_bound_violations, &

or your favorite module for an example.

[iangrooms]

Thanks! I will give this a shot.

[iangrooms]

I'm working on using a namelist to set the quadrature order, and I think it's working. Is there a way to set it up so that if rootfinding_nml is not present in input.nml it will just use default values, rather than killing the program? Or maybe there's a reason to prefer for it to kill the program?

[hkershaw-brown]

I'm working on using a namelist to set the quadrature order, and I think it's working. Is there a way to set it up so that if rootfinding_nml is not present in input.nml it will just use default values, rather than killing the program? Or maybe there's a reason to prefer for it to kill the program?

My opinion is that optional namelists are a good idea, particularly for "advanced" user features or even experimental features, because otherwise all input.nmls have to have every namelist.

The change would be here, e.g. having an optional argument
find_namelist_in_file(.... missing_fine)
if missing_is_fine skip this error

DART/assimilation_code/modules/utilities/utilities_mod.f90

Lines 425 to 430 in d908222

	if(io /= 0) then
	! Reached end of file and didn't find this namelist
	write(msgstring1, *) 'Namelist entry &', trim(nml_name), &
	' must exist in file ', trim(namelist_file_name)
	call error_handler(E_ERR, 'find_namelist_in_file', msgstring1, source)
	else

Do you want to be the first optional namelist in dart?

[iangrooms]

I tried skipping the error message, but then it just hangs. The subroutine returns iunit and if I just skip the error then I'm guessing it returns a value for iunit that is not useful to the calling code. Maybe I should just leave it as a required namelist. I guess one benefit of making it required is that users will know about these options

[hkershaw-brown]

I tried skipping the error message, but then it just hangs. The subroutine returns iunit and if I just skip the error then I'm guessing it returns a value for iunit that is not useful to the calling code. Maybe I should just leave it as a required namelist. I guess one benefit of making it required is that users will know about these options

ok thanks, I'll take a look.

[hkershaw-brown]

Hi Ian, I put an optional namelist in for kde in https://github.com/iangrooms/DART/pull/1/files

The namelist values still get logged with the default values if the namelist is not found in input.nml

Up to you whether you want to make this namelist required or optional.
kde_distribution_mod is only called if someone puts kde in their qceff.csv, so kde_nml being required will mean people running the kde will need to update their input.nmls with the change. So you're kind of forcing them to think about it if they are using kde.

[iangrooms]

Thanks! Rather than forcing users to add this namelist, I think I'd rather update the docs to let them know that there is an option to change the quadrature order via a namelist and thus reduce computational cost, potentially at the cost of reduced accuracy in situations where the observation error is much smaller than the distance between ensemble members in observation space.