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.

Sign up for GitHub

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

Qceff inflation fix; Explicitly handling BNRHF transform failures. #794

Open
wants to merge 5 commits into
base: main
Choose a base branch
from
Open

Qceff inflation fix; Explicitly handling BNRHF transform failures. #794

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
6731e28
Existing code did not take appropriate action when the bnrh transform…
jlaucar Dec 9, 2024
844163f
Corrected use of 1 in a loop over i that was checking to see if each
jlaucar Dec 9, 2024
7bc2724
Bringing the qceff_inflation_fix branch up-to-date with main for review.
jlaucar Dec 9, 2024
1f2b7e1
Explicitly copy the input file into the output if the transform fails
jlaucar Dec 23, 2024
f0e7f99
Corrected syntax error in probit_tranform_mod.f90. Reran tests with
jlaucar Dec 23, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
49 changes: 34 additions & 15 deletions assimilation_code/modules/assimilation/assim_tools_mod.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -346,7 +346,7 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, &
integer(i8), allocatable :: my_state_indx(:)
integer(i8), allocatable :: my_obs_indx(:)

integer :: my_num_obs, i, j, owner, owners_index, my_num_state
integer :: my_num_obs, i, j, owner, owners_index, my_num_state, ierr
integer :: obs_mean_index, obs_var_index
integer :: grp_beg(num_groups), grp_end(num_groups), grp_size, grp_bot, grp_top, group
integer :: num_close_obs, obs_index, num_close_states
Expand All @@ -373,6 +373,7 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, &
type(obs_def_type) :: obs_def
type(time_type) :: obs_time

logical, allocatable :: obs_probit_trans_ok(:), state_probit_trans_ok(:)
logical :: allow_missing_in_state
logical :: local_single_ss_inflate
logical :: local_varying_ss_inflate
Expand All @@ -395,13 +396,15 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, &
my_obs_indx( obs_ens_handle%my_num_vars), &
my_obs_kind( obs_ens_handle%my_num_vars), &
my_obs_type( obs_ens_handle%my_num_vars), &
my_obs_loc( obs_ens_handle%my_num_vars))
my_obs_loc( obs_ens_handle%my_num_vars), &
obs_probit_trans_ok(obs_ens_handle%my_num_vars))

allocate(close_state_dist( ens_handle%my_num_vars), &
close_state_ind( ens_handle%my_num_vars), &
my_state_indx( ens_handle%my_num_vars), &
my_state_kind( ens_handle%my_num_vars), &
my_state_loc( ens_handle%my_num_vars))
my_state_loc( ens_handle%my_num_vars), &
state_probit_trans_ok(ens_handle%my_num_vars))
! end alloc

! Initialize assim_tools_module if needed
Expand Down Expand Up @@ -503,8 +506,9 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, &
! Convert all my state variables to appropriate probit space
call transform_to_probit(ens_size, ens_handle%copies(1:ens_size, i), dist_for_state, &
state_dist_params(i), probit_ens, .false., &
bounded_below, bounded_above, lower_bound, upper_bound)
ens_handle%copies(1:ens_size, i) = probit_ens
bounded_below, bounded_above, lower_bound, upper_bound, ierr)
state_probit_trans_ok(i) = (ierr == 0)
Comment on lines +509 to +510
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is related to the question about the error code being an integer, if you had the error be a logical you could do:
call transform_to_probit(........., state_probit_trans_ok(i))

that way you there is less chance of a future person, putting something in the middle of these two statements.

call transform_to_probit(....., ierr)
call something_else_mathsie(.....,ierr)
state_probit_trans_ok(i) = (ierr == 0)

if(state_probit_trans_ok(i)) ens_handle%copies(1:ens_size, i) = probit_ens
end do

!> optionally convert all state location verticals
Expand Down Expand Up @@ -541,8 +545,9 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, &
! Convert all my obs (extended state) variables to appropriate probit space
call transform_to_probit(ens_size, obs_ens_handle%copies(1:ens_size, i), dist_for_obs, &
obs_dist_params(i), probit_ens, .false., &
bounded_below, bounded_above, lower_bound, upper_bound)
obs_ens_handle%copies(1:ens_size, i) = probit_ens
bounded_below, bounded_above, lower_bound, upper_bound, ierr)
obs_probit_trans_ok(i) = (ierr == 0)
if(obs_probit_trans_ok(i)) obs_ens_handle%copies(1:ens_size, i) = probit_ens
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The added if(obs_probit_trans_ok(i)) & at line 511 gives you the same answer as before, because of the bnrh

fail -> ! Just return the original ensemble
I guess this is for algorithm readability?

endif
end do

Expand Down Expand Up @@ -635,9 +640,11 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, &
orig_obs_prior_var = obs_ens_handle%copies(OBS_PRIOR_VAR_START: &
OBS_PRIOR_VAR_END, owners_index)

! If QC is okay, convert this observation ensemble from probit to regular space
call transform_from_probit(ens_size, obs_ens_handle%copies(1:ens_size, owners_index) , &
obs_dist_params(owners_index), obs_ens_handle%copies(1:ens_size, owners_index))
! Convert this observation ensemble from probit back to regular space if to_probit was successful
if(obs_probit_trans_ok(owners_index)) then
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If the to_probit was not ok, do you want to set the qc to failed fwd op?

The other processes don't have access to obs_probit_trans_ok(owners_index) and so they will go ahead and calculate obs_increments because the qc is 0, then transform_to_probit which may or may not work?

call transform_from_probit(ens_size, obs_ens_handle%copies(1:ens_size, owners_index) , &
obs_dist_params(owners_index), obs_ens_handle%copies(1:ens_size, owners_index))
Comment on lines +645 to +646
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

this is still dangerous since you have the same input for intent(in) and intent(out), but the original code does that so yippee ki yay I guess.

DART/assimilation_code/modules/assimilation/probit_transform_mod.f90

Lines 596 to 601 in f0e7f99

subroutine transform_from_probit(ens_size, probit_ens, p, state_ens)
integer, intent(in) :: ens_size
real(r8), intent(in) :: probit_ens(ens_size)
type(distribution_params_type), intent(inout) :: p
real(r8), intent(out) :: state_ens(ens_size)

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

since you return the original if the transform fails, might be a good time to change this to inout

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

i see that calling a routine with the same array with both an intent(in) and intent(out) goes several levels deep in this code so it would be pretty disruptive to change it. but in the future i'd avoid passing in the same array multiple times with different intents. it's leaving a trap for future code maintainers because someone looking only at the called subroutines could decide to set the outgoing state_ens(:) to 0 or MISSING_R8 as an initial value for debugging, which would result in the probit_ens(:) mysteriously getting overwritten without any protection from the intent(in). you're requiring some behavior in the called subroutine because of how the calling code is written.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I don't think it is standard conforming to alias in this way so not just people changing the code, its undefined what the compiler will do.

endif

obs_prior = obs_ens_handle%copies(1:ens_size, owners_index)
endif IF_QC_IS_OKAY
Expand Down Expand Up @@ -698,10 +705,14 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, &
! Convert the prior and posterior for this observation to probit space
call transform_to_probit(grp_size, obs_prior(grp_bot:grp_top), dist_for_obs, &
temp_dist_params, probit_obs_prior(grp_bot:grp_top), .false., &
bounded_below, bounded_above, lower_bound, upper_bound)
bounded_below, bounded_above, lower_bound, upper_bound, ierr)
! If probit transform fails for any group, this observation cannot be assimimlated
if(ierr /= 0) cycle SEQUENTIAL_OBS
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think it could be helpful to add a CAUTION comment here.
Because of the new error check, the interaction between the transformations and the groups can become messy. If it happens and the probit transform fails, then only part of the ensemble may be updated for a certain obs. Because of the cycle SEQUENTIAL_OBS the remaining groups, which could be transformed without issues, are disregarded.

call transform_to_probit(grp_size, obs_post(grp_bot:grp_top), dist_for_obs, &
temp_dist_params, probit_obs_post(grp_bot:grp_top), .true., &
bounded_below, bounded_above, lower_bound, upper_bound)
bounded_below, bounded_above, lower_bound, upper_bound, ierr)
! If probit transform fails for any group, this observation cannot be assimimlated
if(ierr /= 0) cycle SEQUENTIAL_OBS
! Free up the storage used for this transform
call deallocate_distribution_params(temp_dist_params)

Expand Down Expand Up @@ -761,6 +772,9 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, &
STATE_UPDATE: do j = 1, num_close_states
state_index = close_state_ind(j)

! If transform to probit failed for this state variable, don't update it
if(.not. state_probit_trans_ok(state_index)) cycle STATE_UPDATE

if ( allow_missing_in_state ) then
! Don't allow update of state ensemble with any missing values
if (any(ens_handle%copies(1:ens_size, state_index) == MISSING_R8)) cycle STATE_UPDATE
Expand Down Expand Up @@ -796,6 +810,9 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, &
OBS_UPDATE: do j = 1, num_close_obs
obs_index = close_obs_ind(j)

! If transform to probit failed for this observed variable, don't update it
if(.not. obs_probit_trans_ok(obs_index)) cycle OBS_UPDATE

! Only have to update obs that have not yet been used
if(my_obs_indx(obs_index) > i) then

Expand All @@ -820,7 +837,7 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, &

! Do the inverse probit transform for state variables
call transform_all_from_probit(ens_size, ens_handle%my_num_vars, ens_handle%copies, &
state_dist_params, ens_handle%copies)
state_dist_params, ens_handle%copies, state_probit_trans_ok)

! Every pe needs to get the current my_inflate and my_inflate_sd back
if(local_single_ss_inflate) then
Expand Down Expand Up @@ -876,13 +893,15 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, &
my_obs_type, &
close_obs_ind, &
vstatus, &
my_obs_loc)
my_obs_loc, &
state_probit_trans_ok)

deallocate(close_state_dist, &
my_state_indx, &
close_state_ind, &
my_state_kind, &
my_state_loc)
my_state_loc, &
obs_probit_trans_ok)

! end dealloc

Expand Down
26 changes: 15 additions & 11 deletions assimilation_code/modules/assimilation/filter_mod.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1561,7 +1561,7 @@ subroutine filter_ensemble_inflate(ens_handle, inflate_copy, inflate, ENS_MEAN_C
real(r8) :: probit_ens(ens_size), probit_ens_mean
logical :: bounded_below, bounded_above
real(r8) :: lower_bound, upper_bound
integer :: dist_type
integer :: dist_type, ierr

! Inflate each group separately; Divide ensemble into num_groups groups
grp_size = ens_size / num_groups
Expand Down Expand Up @@ -1603,16 +1603,20 @@ subroutine filter_ensemble_inflate(ens_handle, inflate_copy, inflate, ENS_MEAN_C

call transform_to_probit(grp_size, ens_handle%copies(grp_bot:grp_top, j), &
dist_type, dist_params, probit_ens(1:grp_size), .false., &
bounded_below, bounded_above, lower_bound, upper_bound)

! Compute the ensemble mean in transformed space
probit_ens_mean = sum(probit_ens(1:grp_size)) / grp_size
! Inflate in probit space
call inflate_ens(inflate, probit_ens(1:grp_size), probit_ens_mean, &
ens_handle%copies(inflate_copy, j))
! Transform back from probit space
call transform_from_probit(grp_size, probit_ens(1:grp_size), &
dist_params, ens_handle%copies(grp_bot:grp_top, j))
bounded_below, bounded_above, lower_bound, upper_bound, ierr)

! Only inflate if successful return from the probit transform
if(ierr == 0) then
! Compute the ensemble mean in transformed space
probit_ens_mean = sum(probit_ens(1:grp_size)) / grp_size
! Inflate in probit space
call inflate_ens(inflate, probit_ens(1:grp_size), probit_ens_mean, &
ens_handle%copies(inflate_copy, j))
! Transform back from probit space
call transform_from_probit(grp_size, probit_ens(1:grp_size), &
dist_params, ens_handle%copies(grp_bot:grp_top, j))
endif

end do
endif
end do
Expand Down
Loading
Loading