Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add the Multigrid Beta Filter (MGBF) for ensemble localization in preparation for the 3DRTMAv1 implementation #699

Closed
ManuelPondeca-NOAA opened this issue Feb 15, 2024 · 0 comments · Fixed by #700

Comments

@ManuelPondeca-NOAA
Copy link
Contributor

We plan to use the Multigrid Beta Filter (MGBF) for static background error modeling and ensemble localization in 3DRTMAv1, which is scheduled for implementation in 2025. The work here is intended to add to the GSI an initial version of the MGBF, which has been successfully tested for ensemble localization.

ShunLiu-NOAA pushed a commit that referenced this issue Mar 26, 2024
…#700)

**DUE DATE for merger of this PR into `develop` is 3/29/2024 (six weeks
after PR creation).**

**Description**

Resolves #699

This PR is to add the option to apply Multigrid Beta Filter (MGBF;
[Purser et al. 2022](https://doi.org/10.1175/MWR-D-20-0405.1)) for
ensemble localization instead of Recursive Filter (RF). This work
includes to add an initial version of the MGBF as a subdirectory in GSI.

To apply the MGBF, set "l_mgbf_loc=true" in the namelist and
additionally input "mgbf_loc01.nml". (In Scale/Variable-Dependent
Localization, input also "mgbf_locXX.nml" (XX=02,03,...) with the same
number of grid points.)

<details>

**<summary>How to set MGBF parameters in mgbf_locXX.nml</summary>**

An example of mgbf_locXX.nml is as follows:

```
&PARAMETERS_MGBETA
  mg_ampl01=1.125,            ! length of vertical beta filter (standard deviation; filter grid unit)
  mg_ampl02=2.4,              ! length of horizontal beta filter (standard deviation; filter grid unit)
  mg_ampl03=0.85,             ! length of 3D beta filter (standard deviation; filter grid unit)
  mg_weig1=0.,                ! weight of generation 1
  mg_weig2=0.,                ! weight of generation 2
  mg_weig3=0.,                ! weight of generation 3
  mg_weig4=1.,                ! weight of generation 4
  hx=5,                       ! number of halo grid points in x-direction
  hy=5,                       ! number of halo grid points in y-direction
  hz=3,                       ! number of halo grid points in z-direction
  p=2,                        ! beta filter exponent
  mgbf_line=.false.,          ! set false except for mgbf_proc=2,4,7
  mgbf_proc=8,                ! 1-2: 3D filter; 3-5: 2D filter for static B; 6-8: 2D filter for localization (1,3,6: radial filter; 2,4,7: line filter; 5,8: isotropic line filter)
  lm_a=65,                    ! number of vertical layers in analysis grid
  lm=33,                      ! number of vertical layers in filter grid
  km2=0,                      ! number of 2D variables (set 0 for localization)
  km3=1,                      ! number of 3D variables (set 1 for localization)
  n_ens=30,                   ! ensemble size
  l_loc=.true.,               ! set true in localization
  l_filt_g1=.false.,          ! set false in skipping generation 1
  l_lin_vertical=.true.,      ! set true in applying linear vertical interpolation for analysis-filter mapping
  l_lin_horizontal=.true.,    ! set true in applying linear horizontal interpolation for analysis-filter mapping
  l_quad_horizontal=.false.,  ! set true in applying quadratic horizontal interpolation for analysis-filter mapping
  l_new_map=.true.,           ! set true in applying efficient vertical interpolation for analysis-filter mapping
  l_vertical_filter=.true.,   ! set true in applying vertical beta filter outside 2D filter
  ldelta=.false.,             ! (not used)
  lquart=.false.,             ! set true in applying quadratic horizontal interpolation for up/down-sending
  lhelm=.false.,              ! set true in applying Helmholtz differential operator for weighting
  nm0=1975,                   ! number of analysis grid points in x-direction
  mm0=1350,                   ! number of analysis grid points in y-direction
  gm_max=4,                   ! highest generation (max: 4)
  nxPE=79,                    ! number of MPI processors in x-direction
  nyPE=54,                    ! number of MPI processors in y-direction
  im_filt=8,                  ! number of filter grid points in each MPI processor in x-direction
  jm_filt=8,                  ! number of filter grid points in each MPI processor in y-direction
 /
```

Here, to make the result of MGBF-based localization similar to RF-based
one, we can set the beta filter length ( mg_ampl0[12] ) from the
recursive filter length in the GSI namelist ( s_ens_[vh] ) as:

- $\text{mg\\_ampl01} = \left[\text{s\\_ens\\_v (grid unit)} *
\frac{1}{\sqrt{2}} * \frac{\text{lm}-1}{\text{lm\\_a}-1} \right]^2$
- $\text{mg\\_ampl02} = \left[\frac{\text{s\\_ens\\_h
(km)}}{\text{analysis grid interval (km)}} * \frac{1}{\sqrt{2}} *
\frac{\text{im\\_filt} * \text{nxPE}}{\text{nm0}} * \frac{1}{2} *
\frac{1}{2} * \frac{1}{2} \right]^2$ (in case mg_weig[1-4]=[0,0,0,1])

Please note there are some limitations for the other MGBF parameters
such as:

- The number of MPI processors input in GSI should be nxPE x nyPE 
- (nm0, mm0, lm_a) should be the same as the GSI analysis grid
- nm0 should be divisible by nxPE
- mm0 should be divisible by nyPE
- nm0 / nxPE = mm0 / nyPE

</details>

<details>

**<summary>How to run RRFS regression tests with MGBF-based
localization</summary>**

Change settings in regression/ as follows, and run Test#3
(rrfs_3denvar_glbens)

```diff
diff --git a/regression/regression_namelists.sh b/regression/regression_namelists.sh
index 7ca183ef3..671d028ff 100755
--- a/regression/regression_namelists.sh
+++ b/regression/regression_namelists.sh
@@ -457,7 +457,7 @@ OBS_INPUT::
    beta_s0=0.15,s_ens_h=110,s_ens_v=3,
    regional_ensemble_option=1,
    pseudo_hybens = .false.,
-   grid_ratio_ens = 3,
+   grid_ratio_ens = 5.1,
    l_ens_in_diff_time=.true.,
    ensemble_path='',
    i_en_perts_io=1,
@@ -465,6 +465,7 @@ OBS_INPUT::
    fv3sar_bg_opt=0,
    readin_localization=.true.,
    ens_fast_read=.false.,
+   l_mgbf_loc=.true.,
  /
  &RAPIDREFRESH_CLDSURF
    dfi_radar_latent_heat_time_period=20.0,
```

```diff
diff --git a/regression/regression_namelists_db.sh b/regression/regression_namelists_db.sh
index e03917e88..36b8b6a22 100755
--- a/regression/regression_namelists_db.sh
+++ b/regression/regression_namelists_db.sh
@@ -438,7 +438,7 @@ OBS_INPUT::
    beta_s0=0.15,s_ens_h=110,s_ens_v=3,
    regional_ensemble_option=1,
    pseudo_hybens = .false.,
-   grid_ratio_ens = 3,
+   grid_ratio_ens = 5.1,
    l_ens_in_diff_time=.true.,
    ensemble_path='',
    i_en_perts_io=1,
@@ -446,6 +446,7 @@ OBS_INPUT::
    fv3sar_bg_opt=0,
    readin_localization=.true.,
    ens_fast_read=.false.,
+   l_mgbf_loc=.true.,
  /
  &RAPIDREFRESH_CLDSURF
    dfi_radar_latent_heat_time_period=20.0,
```

```diff
diff --git a/regression/regression_param.sh b/regression/regression_param.sh
index 2ac615fc4..6186acdbb 100755
--- a/regression/regression_param.sh
+++ b/regression/regression_param.sh
@@ -87,23 +87,23 @@ case $regtest in
     rrfs_3denvar_glbens)

         if [[ "$machine" = "Hera" ]]; then
-           topts[1]="0:15:00" ; popts[1]="20/1/"  ; ropts[1]="/1"
-           topts[2]="0:15:00" ; popts[2]="20/2/"  ; ropts[2]="/1"
+           topts[1]="0:15:00" ; popts[1]="30/2/"  ; ropts[1]="/1"
+           topts[2]="0:15:00" ; popts[2]="15/4/"  ; ropts[2]="/1"
         elif [[ "$machine" = "Orion" ]]; then
-           topts[1]="0:15:00" ; popts[1]="20/1/" ; ropts[1]="/1"
-           topts[2]="0:15:00" ; popts[2]="20/2/" ; ropts[2]="/2"
+           topts[1]="0:15:00" ; popts[1]="30/2/"  ; ropts[1]="/1"
+           topts[2]="0:15:00" ; popts[2]="15/4/"  ; ropts[2]="/2"
         elif [[ "$machine" = "Hercules" ]]; then
-           topts[1]="0:15:00" ; popts[1]="20/1/" ; ropts[1]="/1"
-           topts[2]="0:15:00" ; popts[2]="20/2/" ; ropts[2]="/2"
+           topts[1]="0:15:00" ; popts[1]="30/2/"  ; ropts[1]="/1"
+           topts[2]="0:15:00" ; popts[2]="15/4/"  ; ropts[2]="/2"
         elif [[ "$machine" = "Jet" ]]; then
-           topts[1]="0:15:00" ; popts[1]="20/1/"  ; ropts[1]="/1"
-           topts[2]="0:15:00" ; popts[2]="20/2/"  ; ropts[2]="/1"
+           topts[1]="0:15:00" ; popts[1]="15/4/"  ; ropts[1]="/1"
+           topts[2]="0:15:00" ; popts[2]="10/6/"  ; ropts[2]="/1"
         elif [[ "$machine" = "Gaea" ]]; then
-           topts[1]="0:15:00" ; popts[1]="18/1/"  ; ropts[1]="/1"
-           topts[2]="0:15:00" ; popts[2]="18/2/"  ; ropts[2]="/1"
+           topts[1]="0:15:00" ; popts[1]="15/4/"  ; ropts[1]="/1"
+           topts[2]="0:15:00" ; popts[2]="10/6/"  ; ropts[2]="/1"
         elif [[ "$machine" = "wcoss2" ]]; then
-           topts[1]="0:15:00" ; popts[1]="64/1/" ; ropts[1]="/1"
-           topts[2]="0:15:00" ; popts[2]="128/2/" ; ropts[2]="/1"
+           topts[1]="0:15:00" ; popts[1]="60/1/"  ; ropts[1]="/1"
+           topts[2]="0:15:00" ; popts[2]="30/2/"  ; ropts[2]="/1"
         fi

         if [ "$debug" = ".true." ] ; then
```

```diff
diff --git a/regression/rrfs_3denvar_glbens.sh b/regression/rrfs_3denvar_glbens.sh
index af5da5117..04fd73d57 100755
--- a/regression/rrfs_3denvar_glbens.sh
+++ b/regression/rrfs_3denvar_glbens.sh
@@ -272,6 +272,46 @@ $gsi_namelist

 EOF

+cat << EOF > mgbf_loc01.nml
+&PARAMETERS_MGBETA
+  mg_ampl01=1.125,            ! length of vertical beta filter (standard deviation; filter grid unit)
+  mg_ampl02=1.615,            ! length of horizontal beta filter (standard deviation; filter grid unit)
+  mg_ampl03=0.85,             ! length of 3D beta filter (standard deviation; filter grid unit)
+  mg_weig1=0.,                ! weight of generation 1
+  mg_weig2=1.,                ! weight of generation 2
+  mg_weig3=0.,                ! weight of generation 3
+  mg_weig4=0.,                ! weight of generation 4
+  hx=4,                       ! number of halo grid points in x-direction
+  hy=4,                       ! number of halo grid points in y-direction
+  hz=3,                       ! number of halo grid points in z-direction
+  p=2,                        ! beta filter exponent
+  mgbf_line=.false.,          ! set false except for mgbf_proc=2,4,7
+  mgbf_proc=8,                ! 1-2: 3D filter; 3-5: 2D filter for static B; 6-8: 2D filter for localization (1,3,6: radial filter; 2,4,7: line filter; 5,8: isotropic line filter)
+  lm_a=65,                    ! number of vertical layers in analysis grid
+  lm=33,                      ! number of vertical layers in filter grid
+  km2=0,                      ! number of 2D variables (set 0 for localization)
+  km3=1,                      ! number of 3D variables (set 1 for localization)
+  n_ens=10,                   ! ensemble size
+  l_loc=.true.,               ! set true in localization
+  l_filt_g1=.false.,          ! set false in skipping generation 1
+  l_lin_vertical=.true.,      ! set true in applying linear vertical interpolation for analysis-filter mapping
+  l_lin_horizontal=.true.,    ! set true in applying linear horizontal interpolation for analysis-filter mapping
+  l_quad_horizontal=.false.,  ! set true in applying quadratic horizontal interpolation for analysis-filter mapping
+  l_new_map=.true.,           ! set true in applying efficient vertical interpolation for analysis-filter mapping
+  l_vertical_filter=.true.,   ! set true in applying vertical beta filter outside 2D filter
+  ldelta=.false.,             ! (not used)
+  lquart=.false.,             ! set true in applying quadratic horizontal interpolation for up/down-sending
+  lhelm=.false.,              ! set true in applying Helmholtz differential operator for weighting
+  nm0=40,                     ! number of analysis grid points in x-direction
+  mm0=24,                     ! number of analysis grid points in y-direction
+  gm_max=2,                   ! highest generation (max: 4)
+  nxPE=10,                    ! number of MPI processors in x-direction
+  nyPE=6,                     ! number of MPI processors in y-direction
+  im_filt=4,                  ! number of filter grid points in each MPI processor in x-direction
+  jm_filt=4,                  ! number of filter grid points in each MPI processor in y-direction
+ /
+EOF
+
 # Copy executable and fixed files to $tmpdir
 if [[ $exp == *"updat"* ]]; then
    $ncp $gsiexec_updat  ./gsi.x
```

</details>

**Type of change**

Please delete options that are not relevant.

- [ ] Bug fix (non-breaking change which fixes an issue)
- [x] New feature (non-breaking change which adds functionality)
- [ ] Breaking change (fix or feature that would cause existing
functionality to not work as expected)
- [ ] This change requires a documentation update

**How Has This Been Tested?**

EnVar for NA-domain RRFS was tested with "mgbf_locXX.nml" (XX=01) shown
above on Orion. The resulting analysis increment was similar to the
original and the computation time for localization became short.
  
**Checklist**

- [x] My code follows the style guidelines of this project
- [x] I have performed a self-review of my own code
- [x] I have commented my code, particularly in hard-to-understand areas
- [x] New and existing tests pass with my changes
- [x] Any dependent changes have been merged and published

Co-authored-by: Sho Yokota <syokota@Orion-login-1.HPC.MsState.Edu>
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging a pull request may close this issue.

1 participant