Streaming filter

[c2xu]

In the MOM_streaming_filter module, a system of two coupled ODEs is configured as a streaming band-pass filter that takes the broadband model output as the input and returns the filtered signal at a user specified tidal frequency as the output. It is capable of detecting the instantaneous tidal signals while the model is running, and can be used to impose frequency-dependent wave drag or to de-tide model output. An example of filtering the M2 and K1 barotropic velocities is provided in the MOM_barotropic module.

[Hallberg-NOAA]

This commit is a step toward an interesting new capability for on-line detection of periodic signals. I definitely appreciate the direction that this is going, but I think that there are number of issues that should be addressed before this can be added to the MOM6 code base.

1. Whenever we add new diagnostic capabilities, we always make sure that they do not have an adverse impact on simulations where they are not used. Usually this is done by adding a logical flag to avoid the extra calculations related to the diagnostic when they are not being used. Please add a run-time logical flag (perhaps using a diagnostic id) to determine whether Filt_accum() is called from btstep().

2. The 4 pointers um2, uk1, vm2, and vk1 are temporary local variables in btstep(), but I do not see where they are actually being used, either in subsequent model calculations or in any diagnostics. What is their purpose? Should this PR be revised so that these fields are optionally output somewhere?

3. Inside of Filt_accum(), the code is working on the entire array size, including the halo regions. However, we do not generally guarantee that there are valid values throughout the halo regions, and there could be NaNs, which could be problematic. Elsewhere in the MOM6 code, we only work on the computational domain plus a specified (via an argument) halo size, and the same should perhaps be done here. This would require more extensive integration of this code with the MOM6 infrastructure (e.g., in the use of the contents of a hor_index_type argument) and it might require different versions of Filt_accum() for various horizontal locations (e.g., tracer points vs. u-points vs. v-points), but I think that it would be safer.

4. Rather than having a streaming_filter_CS that consists of a linked list of the Filt_node types for each of the filtered variables, have you considered having the calling routine set up (via a call to a revised filt_register) and store a separate Filt_node type for each of the fields that it is working on? I don't see how that approach would lose any generality, and it would simplify a lot of the complexity of this code that tries to figure out which node to work on with each call to Filt_accum(). There may be something that I am missing with this suggestion, but if so that might be useful to layout in the conversation surrounding this PR.

I may raise other issues after these first 4 have been addressed, but for now these seem to me like the most prominent issues raised by this very promising pull request.

[c2xu]

Thanks for the comments. Please see below for my point-by-point reply.

This commit is a step toward an interesting new capability for on-line detection of periodic signals. I definitely appreciate the direction that this is going, but I think that there are number of issues that should be addressed before this can be added to the MOM6 code base.

1. Whenever we add new diagnostic capabilities, we always make sure that they do not have an adverse impact on simulations where they are not used. Usually this is done by adding a logical flag to avoid the extra calculations related to the diagnostic when they are not being used. Please add a run-time logical flag (perhaps using a diagnostic id) to determine whether Filt_accum() is called from btstep().

Run-time flags have been added for the M2 and K1 filters separately.

2. The 4 pointers um2, uk1, vm2, and vk1 are temporary local variables in btstep(), but I do not see where they are actually being used, either in subsequent model calculations or in any diagnostics. What is their purpose? Should this PR be revised so that these fields are optionally output somewhere?

My intention is to use them for implementing the frequency-dependent internal wave drag, but this involves some major modification of the way that the linear wave drag is implemented in the MOM_barotropic module, and I think it would be better to do so in a separate pull request (which will follow shortly). For now, they are not being used anywhere in the code.

3. Inside of Filt_accum(), the code is working on the entire array size, including the halo regions. However, we do not generally guarantee that there are valid values throughout the halo regions, and there could be NaNs, which could be problematic.

I thought it would be okay for Filt_accum to work with NaNs for the halo regions. It should not affect the calculation in non-halo regions.

Elsewhere in the MOM6 code, we only work on the computational domain plus a specified (via an argument) halo size, and the same should perhaps be done here. This would require more extensive integration of this code with the MOM6 infrastructure (e.g., in the use of the contents of a hor_index_type argument)

I'm not very familiar of how this works. Is there an example that I can follow?

and it might require different versions of Filt_accum() for various horizontal locations (e.g., tracer points vs. u-points vs. v-points), but I think that it would be safer.

Actually, one reason that I use UBOUND and LBOUND to determine the array size is to avoid the implementation of different versions of Filt_accum for different horizontal grids.

4. Rather than having a streaming_filter_CS that consists of a linked list of the Filt_node types for each of the filtered variables, have you considered having the calling routine set up (via a call to a revised filt_register) and store a separate Filt_node type for each of the fields that it is working on? I don't see how that approach would lose any generality, and it would simplify a lot of the complexity of this code that tries to figure out which node to work on with each call to Filt_accum(). There may be something that I am missing with this suggestion, but if so that might be useful to layout in the conversation surrounding this PR.

Yes, you're right. I simplified the code by removing the master control structure, and now each filter has its own control structure.

I may raise other issues after these first 4 have been addressed, but for now these seem to me like the most prominent issues raised by this very promising pull request.

[c2xu]

This commit is a step toward an interesting new capability for on-line detection of periodic signals. I definitely appreciate the direction that this is going, but I think that there are number of issues that should be addressed before this can be added to the MOM6 code base.

1. Whenever we add new diagnostic capabilities, we always make sure that they do not have an adverse impact on simulations where they are not used. Usually this is done by adding a logical flag to avoid the extra calculations related to the diagnostic when they are not being used. Please add a run-time logical flag (perhaps using a diagnostic id) to determine whether Filt_accum() is called from btstep().
2. The 4 pointers um2, uk1, vm2, and vk1 are temporary local variables in btstep(), but I do not see where they are actually being used, either in subsequent model calculations or in any diagnostics. What is their purpose? Should this PR be revised so that these fields are optionally output somewhere?
3. Inside of Filt_accum(), the code is working on the entire array size, including the halo regions. However, we do not generally guarantee that there are valid values throughout the halo regions, and there could be NaNs, which could be problematic. Elsewhere in the MOM6 code, we only work on the computational domain plus a specified (via an argument) halo size, and the same should perhaps be done here. This would require more extensive integration of this code with the MOM6 infrastructure (e.g., in the use of the contents of a hor_index_type argument) and it might require different versions of Filt_accum() for various horizontal locations (e.g., tracer points vs. u-points vs. v-points), but I think that it would be safer.
4. Rather than having a streaming_filter_CS that consists of a linked list of the Filt_node types for each of the filtered variables, have you considered having the calling routine set up (via a call to a revised filt_register) and store a separate Filt_node type for each of the fields that it is working on? I don't see how that approach would lose any generality, and it would simplify a lot of the complexity of this code that tries to figure out which node to work on with each call to Filt_accum(). There may be something that I am missing with this suggestion, but if so that might be useful to layout in the conversation surrounding this PR.

I may raise other issues after these first 4 have been addressed, but for now these seem to me like the most prominent issues raised by this very promising pull request.

The filtering is a point-wise operation. It does not involve any spatial differentiation or integration. Results at a grid point will not be affected by the results of its neighbouring points in any way, if this was your concern when you raised the issue of halo regions, @Hallberg-NOAA.

[codecov]

Codecov Report

All modified and coverable lines are covered by tests ✅

Project coverage is 37.01%. Comparing base (1b39d07) to head (7106b09).
Report is 19 commits behind head on dev/gfdl.

Additional details and impacted files

@@             Coverage Diff              @@
##           dev/gfdl     #675      +/-   ##
============================================
+ Coverage     36.99%   37.01%   +0.02%     
============================================
  Files           272      274       +2     
  Lines         82403    82591     +188     
  Branches      15413    15450      +37     
============================================
+ Hits          30487    30575      +88     
- Misses        46226    46295      +69     
- Partials       5690     5721      +31     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

[Hallberg-NOAA]

Thank you for these revisions, which mostly address the concerns I raised with the review of the previous version.

I do see from the code that this is a point-wise filter, and so what happens at one point will not impact others. The one concern is that if there are NaNs (from uninitialized values) in the halo regions, doing any calculations with them can lead to an exception that will cause the model to abort (depending on compiler settings), even if they don't actually impact any subsequent values. This is most common when the model is compiled for debugging, but having the model come down due to a fake problem can be very annoying when trying to debug a real problem. If we knew that the fields that are being filtered are always being initialized (e.g, to 0), this would not be an issue, but there is no way to guarantee that this will always be true.

Examples where we use a hor_index_type argument to select only the computational domain can be found in the various routines in MOM6/src/framework/MOM_checksums.F90.

[c2xu]

Thank you for these revisions, which mostly address the concerns I raised with the review of the previous version.

I do see from the code that this is a point-wise filter, and so what happens at one point will not impact others. The one concern is that if there are NaNs (from uninitialized values) in the halo regions, doing any calculations with them can lead to an exception that will cause the model to abort (depending on compiler settings), even if they don't actually impact any subsequent values. This is most common when the model is compiled for debugging, but having the model come down due to a fake problem can be very annoying when trying to debug a real problem. If we knew that the fields that are being filtered are always being initialized (e.g, to 0), this would not be an issue, but there is no way to guarantee that this will always be true.

Examples where we use a hor_index_type argument to select only the computational domain can be found in the various routines in MOM6/src/framework/MOM_checksums.F90.

Hi @Hallberg-NOAA, I followed what's been done in register_barotropic_restarts and tried to make the dimensions of s1 and u1 consistent with variables defined in barotropic_CS, but I still don't fully understand how hor_index_type works and am not confident if my code is correct (though it compiles and runs). Please let me know if further changes are needed. Thanks!

[Hallberg-NOAA]

With these changes, this is all looking good to me. Thank you for this contribution, @c2xu, and your patience as we sorted out the details of this PR.

[marshallward]

Gaea regression: https://gitlab.gfdl.noaa.gov/ogrp/mom6ci/MOM6/-/jobs/137854 ✔️ 🟡

Squash-merging these 9 commits.