search

E3SM-Project / e3sm_diags

E3SM Diagnostics package
https://e3sm-project.github.io/e3sm_diags
BSD 3-Clause "New" or "Revised" License
42 stars 32 forks source link

Update QBO amplitude calculation #816

Closed chengzhuzhang closed 2 months ago

chengzhuzhang commented 6 months ago

Change the NFFT0, frequency, and interpolation of period calculations

Description

@jjbenedict and @whannah1 identified an odd behavior that the QBO amplitude peak (bottom figure) gets weaker with longer years of simulation being analysed.

@chengzhuzhang tested with multiple year lengths (1-10, 1-20, …1-100) and can reproduce this behavior. QBO_diags_tests.pptx

This PR is contributed by @justin-richling and Jack Chen (cchen@ucar.edu) for coming up with a fix.

Checklist

If applicable:

chengzhuzhang commented 6 months ago

@justin-richling Thanks for the PR! I still won't be able to figure out why it is needed to add the default=str argument when json.dump (json.dump(json_dict, outfile, default=str)). But it seems to be a fine workaround for saving the files in json format.

I ran a test with the new algorithm and compared to the original one, the power spectrum density changed a lot, and I'm not sure if it is the intended behavior. Could you provide an explanation of the code change here? Thank you!

Before:

image

After: image

Command line to reproduce on Perlmutter:e3sm_diags qbo --no_viewer --reference_data_path '/global/cfs/cdirs/e3sm/diagnostics/observations/Atm/time-series' --test_data_path '/global/cfs/cdirs/e3sm/chengzhu/tutorial2024/v3.LR.historical_0101/post/atm/180x360_aave/ts/monthly/15yr' --results_dir '/global/cfs/cdirs/e3sm/www/chengzhu/tests/qbofix' --case_id 'QBO-ERA-Interim' --ref_timeseries_input --test_timeseries_input --run_type 'model_vs_obs' --sets 'qbo' --variables 'U' --seasons 'ANN' 'DJF' 'MAM' 'JJA' 'SON' --regions 'global' --regrid_tool 'esmf' --regrid_method 'conservative' --multiprocessing --num_workers '24' --main_title 'QBO index, amplitude, and power spectral density for U' --backend 'mpl' --output_format 'png' --canvas_size_w '1212' --canvas_size_h '1628' --figsize '8.5' '11.0' --dpi '150' --arrows --test_name 'extendedOutput.v3.LR.historical_0101' --short_test_name 'v3.LR.historical_0101' --test_colormap 'cet_rainbow.rgb' --ref_name 'ERA-Interim' --reference_name 'ERA-Interim' --reference_colormap 'cet_rainbow.rgb' --diff_title 'Model - Observations' --diff_colormap 'diverging_bwr.rgb' --granulate 'variables' 'plevs' 'regions' --selectors 'sets' 'seasons' --test_start_yr '2000' --test_end_yr '2014' --ref_start_yr '2000' --ref_end_yr '2014' --start_yr '2000' --end_yr '2014'

justin-richling commented 5 months ago

Sure thing, @chengzhuzhang!

From Jack:

Fast Fourier transform (FFT) is used for the QBO spectra calculator in which power spectra in frequency space is carried out in a discrete form.

The original code used an interpolation approach to obtain power spectra in the period space (inverse of frequency).

However, this approach could not produce consistent results under simulations of different durations. For a longer simulation, more frequencies will be resolved through FFT which could lower the peak amplitude in the frequency space.

The revised code uses a bin approach which sums all the power within each period bin (2 months). Consequently, this approach produces robust power spectra.

jjbenedict commented 5 months ago

@justin-richling -- so it looks like in this 15-year demonstration case the higher frequency activity is better resolved and the lower frequency activity is poorly resolved. I assume amplitude line in the bottom panel would become smoother if a longer time window were used, correct?

justin-richling commented 5 months ago

@jjbenedict Correct, for something like a typical 30-year run (ERA Interim), the result will be smoother.

chengzhuzhang commented 4 months ago

Thank you for the explanation @justin-richling! My understanding is that the original one does provide smoother power spectrum, but the amplitude of the peak can get weaker when more frequency is resolved given longer time series. The new one could be difficult to use for simulations shorter than i.e. 30 years. @jjbenedict Could you provide some recommendation how we should move forward?

jjbenedict commented 4 months ago

@chengzhuzhang I talked with @justin-richling and Jack earlier this week. Justin is actively working with Walter Hannah to implement the wavelet analysis (Walter had earlier written a python script to compute the wavelets). We're thinking it might be good to include in E3SM Diags both the traditional Fourier method (as implemented now, and appropriate for time windows of 30 yr or longer) -and- the wavelet method (which may be better for shorter time windows). We'll have to do some tests with the wavelet method to see how effective it is at shorter time windows, but the expectation is that is will be smoother than Jack's existing Fourier method. Maintaining the Fourier method would also allow greater consistency with QBO plots from previous E3SM versions. We realize the urgency of implementing the wavelet analysis, given that v3 simulations are now being run... we'll do our best to get this set up and tested as quickly as possible.

FYI -- I will be on vacation July 6-22 (though interrupted by a meeting July 16-17, when I will be checking email).

chengzhuzhang commented 4 months ago

@jjbenedict thank you for the updates, and the heads-up! It seems a good idea to include both methods, and we can discuss on the implementation details after getting testing results from the wavelet analysis.

minxu74 commented 4 months ago

Besides the smoothness, I wonder why the phase were different among the old and revised ones. The old result showed the phase shifted simulated by E3SM, but the new one looked perfectly good.

chengzhuzhang commented 2 months ago

The wavelet analysis PR #836 replaces this one.