[Dart-dev] [5633] DART/branches/development/adaptive_inflate/adaptive_inflate_mod.f90: old improvement, never committed.
nancy at ucar.edu
nancy at ucar.edu
Fri Mar 30 11:13:56 MDT 2012
Revision: 5633
Author: nancy
Date: 2012-03-30 11:13:55 -0600 (Fri, 30 Mar 2012)
Log Message:
-----------
old improvement, never committed. if initial inflation mean/sd values were set
by namelist (not read in from a restart file), the original code did 2 transposes
to set the values in the ensemble handle. with large processor counts and
large state vectors this was slow. the new code finds which tasks are the
owners of the inflation mean and sd copies, and reads in the restart values
only on those tasks directly into the ens handle without having to leave the
'all vars' configuration. has been tested with different processor counts
and gives results identical to original code. no change in results; not a
bug fix per-se but an optimization. also moved the log file messages to
the end of the init routine; if you read in the inf mean/sd from a file
it now prints out the actual min/max values for each field into the log.
this change should probably be added to the trunk and kodiak branches after
a code review by tim.
Modified Paths:
--------------
DART/branches/development/adaptive_inflate/adaptive_inflate_mod.f90
-------------- next part --------------
Modified: DART/branches/development/adaptive_inflate/adaptive_inflate_mod.f90
===================================================================
--- DART/branches/development/adaptive_inflate/adaptive_inflate_mod.f90 2012-03-30 15:00:05 UTC (rev 5632)
+++ DART/branches/development/adaptive_inflate/adaptive_inflate_mod.f90 2012-03-30 17:13:55 UTC (rev 5633)
@@ -18,8 +18,8 @@
error_handler, E_ERR, E_MSG
use random_seq_mod, only : random_seq_type, random_gaussian, init_random_seq
use ensemble_manager_mod, only : ensemble_type, all_vars_to_all_copies, all_copies_to_all_vars, &
- read_ensemble_restart, write_ensemble_restart
-use mpi_utilities_mod, only : my_task_id
+ read_ensemble_restart, write_ensemble_restart, get_copy_owner_index
+use mpi_utilities_mod, only : my_task_id, send_to, receive_from
implicit none
private
@@ -92,7 +92,8 @@
character(len = *), intent(in) :: label
character(len = 128) :: det, tadapt, sadapt, akind, rsread, nmread
-integer :: restart_unit, io
+integer :: restart_unit, io, owner, owners_index
+real(r8) :: minmax_mean(2), minmax_sd(2)
! Record the module version if this is first initialize call
if(.not. initialized) then
@@ -105,49 +106,6 @@
if(.not. deterministic) call init_random_seq(inflate_handle%ran_seq, my_task_id()+1)
endif
-! Write to log file what kind of inflation is being used.
-if(deterministic) then
- det = 'deterministic,'
-else
- det = 'random-noise,'
-endif
-if (sd_initial > inf_lower_bound) then
- det = trim(det) // ' covariance adaptive,'
-endif
-if (inf_lower_bound < 1.0_r8) then
- det = trim(det) // ' deflation permitted,'
-endif
-if (sd_initial > 0.0_r8) then
- tadapt = ' time-adaptive,'
-else
- tadapt = ' time-constant,'
-endif
-
-select case(inf_flavor)
- case (0)
- det = ''
- tadapt = ''
- sadapt = ''
- akind = 'None '
- case (1)
- sadapt = ''
- akind = ' observation-space'
- case (2)
- sadapt = ' spatially-varying,'
- akind = ' state-space '
- case (3)
- sadapt = ' spatially-constant,'
- akind = ' state-space'
- case default
- write(errstring, *) 'Illegal inflation value for ', label
- call error_handler(E_ERR, 'adaptive_inflate_init', errstring, source, revision, revdate)
-end select
-
-! say in plain english what kind of inflation was selected.
-write(errstring, '(4A)') &
- trim(det), trim(tadapt), trim(sadapt), trim(akind)
-call error_handler(E_MSG, trim(label) // ' inflation:', errstring, source, revision, revdate)
-
! more information for users to document what they selected in the nml:
! if flavor > 0, look at read_from_restart for both mean and sd.
! print the actual value used if from namelist, or say
@@ -201,6 +159,11 @@
call error_handler(E_ERR, 'adaptive_inflate_init', errstring, source, revision, revdate)
endif
+! give this an initial value which defaults to time-constant inflation
+! change it below if the value(s) are > 0
+minmax_mean = 0.0_r8
+minmax_sd = 0.0_r8
+
!------ Block for state space inflation initialization ------
! Types 2 and 3 are state space inflation types
@@ -216,6 +179,10 @@
errstring, source, revision, revdate)
endif
+ ! set this, and then below if sd values are from a file, compute them
+ ! and send the value to PE0 if not already there.
+ minmax_sd = sd_initial
+
! Read in initial values from file OR get from subroutine arguments
! If either mean, sd, or both are to be read from the restart file, read them in.
! Then test to see if either are to be overwritten, and do it.
@@ -226,21 +193,89 @@
call read_ensemble_restart(ens_handle, ss_inflate_index, ss_inflate_sd_index, &
.true., in_file_name, force_single_file = .true.)
endif
- ! If one or both are false, we need to transpose and then set the values
+ ! If one or both are false, we need to set the array values from the namelist item
if (.not. mean_from_restart .or. .not. sd_from_restart) then
- ! Get initial values from higher level; requires pe's to have all copies of some vars
- call all_vars_to_all_copies(ens_handle)
- if (.not. mean_from_restart) ens_handle%copies(ss_inflate_index, :) = inf_initial
- if (.not. sd_from_restart) ens_handle%copies(ss_inflate_sd_index, :) = sd_initial
- call all_copies_to_all_vars(ens_handle)
+! original code requires an expensive transpose which is not necessary.
+! ! Get initial values from higher level; requires pe's to have all copies of some vars
+! call all_vars_to_all_copies(ens_handle)
+! if (.not. mean_from_restart) ens_handle%copies(ss_inflate_index, :) = inf_initial
+! if (.not. sd_from_restart) ens_handle%copies(ss_inflate_sd_index, :) = sd_initial
+! call all_copies_to_all_vars(ens_handle)
+! proposed alternate:
+ ! if setting initial values from the namelist, find out which task has the
+ ! inflation and inf sd values and set them only on that task. this saves us
+ ! a transpose.
+ if (.not. mean_from_restart) then
+ call get_copy_owner_index(ss_inflate_index, owner, owners_index)
+ if (owner == my_task_id()) ens_handle%vars(:, owners_index) = inf_initial
+ endif
+ if (.not. sd_from_restart) then
+ call get_copy_owner_index(ss_inflate_sd_index, owner, owners_index)
+ if (owner == my_task_id()) ens_handle%vars(:, owners_index) = sd_initial
+ endif
endif
+
+ ! this block figures out what the min/max value of the mean/sd is
+ ! if we are reading in the values from a restart file. it is used
+ ! in diagnostic output so it needs to get to PE0. we also could check it
+ ! against the limits in the namelist to be sure the file values aren't
+ ! already outside the requested limits. (not sure that's necessary -
+ ! depends on when the code that changes the values imposes the limits.)
+ if (mean_from_restart) then
+ call get_copy_owner_index(ss_inflate_index, owner, owners_index)
+ ! if inflation array is already on PE0, just figure out the
+ ! largest value in the array and we're done.
+ if (owner == 0) then
+ minmax_mean(1) = minval(ens_handle%vars(:, owners_index))
+ minmax_mean(2) = maxval(ens_handle%vars(:, owners_index))
+ else
+ ! someone else has the inf array. have the owner send the min/max
+ ! values to PE0. after this point only PE0 has the right value
+ ! in minmax_mean, but it is the only one who is going to print below.
+ if (my_task_id() == 0) then
+ call receive_from(owner, minmax_mean)
+ else if (my_task_id() == owner) then
+ minmax_mean(1) = minval(ens_handle%vars(:, owners_index))
+ minmax_mean(2) = maxval(ens_handle%vars(:, owners_index))
+ call send_to(0, minmax_mean)
+ endif
+ endif
+ endif
+ if (sd_from_restart) then
+ call get_copy_owner_index(ss_inflate_sd_index, owner, owners_index)
+ ! if inflation sd array is already on PE0, just figure out the
+ ! largest value in the array and we're done.
+ if (owner == 0) then
+ minmax_sd(1) = minval(ens_handle%vars(:, owners_index))
+ minmax_sd(2) = maxval(ens_handle%vars(:, owners_index))
+ else
+ ! someone else has the sd array. have the owner send the min/max
+ ! values to PE0. after this point only PE0 has the right value
+ ! in minmax_sd, but it is the only one who is going to print below.
+ if (my_task_id() == 0) then
+ call receive_from(owner, minmax_sd)
+ else if (my_task_id() == owner) then
+ minmax_sd(1) = minval(ens_handle%vars(:, owners_index))
+ minmax_sd(2) = maxval(ens_handle%vars(:, owners_index))
+ call send_to(0, minmax_sd)
+ endif
+ endif
+ endif
! Inflation type 3 is spatially-constant. Make sure the entire array is set to that
! value; the computation only uses index 1, but the diagnostics write out the entire
- ! array and it will be misleading if not constant.
+ ! array and it will be misleading if not constant. the inf values were set above.
+ ! if they were set by namelist, this code changes nothing. but if they were read in
+ ! from a file, then it is possible the values vary across the array. these lines
+ ! ensure the entire array contains a single constant value to match what the code uses.
if(inf_flavor == 3) then
- ens_handle%copies(ss_inflate_index, :) = ens_handle%copies(ss_inflate_index, 1)
- ens_handle%copies(ss_inflate_sd_index, :) = ens_handle%copies(ss_inflate_sd_index, 1)
+! this was wrong - we are var complete at this point...
+! ens_handle%copies(ss_inflate_index, :) = ens_handle%copies(ss_inflate_index, 1)
+! ens_handle%copies(ss_inflate_sd_index, :) = ens_handle%copies(ss_inflate_sd_index, 1)
+ call get_copy_owner_index(ss_inflate_index, owner, owners_index)
+ if (owner == my_task_id()) ens_handle%vars(:, owners_index) = ens_handle%vars(1, owners_index)
+ call get_copy_owner_index(ss_inflate_sd_index, owner, owners_index)
+ if (owner == my_task_id()) ens_handle%vars(:, owners_index) = ens_handle%vars(1, owners_index)
endif
!------ Block for obs. space inflation initialization ------
@@ -249,7 +284,7 @@
else if(inf_flavor == 1) then
! Initialize observation space inflation values from restart files
- ! Only values are inflation, inflation_sd
+ ! Only single values for inflation, inflation_sd (not arrays)
if(mean_from_restart .or. sd_from_restart) then
! Open the file
restart_unit = get_unit()
@@ -276,6 +311,64 @@
endif
+! Write to log file what kind of inflation is being used.
+! This used to be at the start, but if you are starting from a restart
+! file then you can't tell what sd values are being read in until here.
+if(deterministic) then
+ det = 'deterministic,'
+else
+ det = 'random-noise,'
+endif
+if (inflate_handle%sd_lower_bound > inf_lower_bound) then
+ det = trim(det) // ' covariance adaptive,'
+endif
+if (inflate_handle%inf_lower_bound < 1.0_r8) then
+ det = trim(det) // ' deflation permitted,'
+endif
+if (minmax_sd(2) > 0.0_r8) then
+ tadapt = ' time-adaptive,'
+else
+ tadapt = ' time-constant,'
+endif
+
+select case(inf_flavor)
+ case (0)
+ det = ''
+ tadapt = ''
+ sadapt = ''
+ akind = 'None '
+ case (1)
+ sadapt = ''
+ akind = ' observation-space'
+ case (2)
+ sadapt = ' spatially-varying,'
+ akind = ' state-space '
+ case (3)
+ sadapt = ' spatially-constant,'
+ akind = ' state-space'
+ case default
+ write(errstring, *) 'Illegal inflation value for ', label
+ call error_handler(E_ERR, 'adaptive_inflate_init', errstring, source, revision, revdate)
+end select
+
+! say in plain english what kind of inflation was selected.
+write(errstring, '(4A)') &
+ trim(det), trim(tadapt), trim(sadapt), trim(akind)
+call error_handler(E_MSG, trim(label) // ' inflation:', errstring, source, revision, revdate)
+
+if (inf_flavor > 0) then
+ if (mean_from_restart) then
+ write(errstring, '(A, F8.3, A, F8.3)') &
+ 'inf mean from restart file: min value: ', minmax_mean(1), ' max value: ', minmax_mean(2)
+ call error_handler(E_MSG, trim(label) // ' inflation:', errstring, source, revision, revdate)
+ endif
+ if (sd_from_restart) then
+ write(errstring, '(A, F8.3, A, F8.3)') &
+ 'inf stddev from restart file: min value: ', minmax_sd(1), ' max value: ', minmax_sd(2)
+ call error_handler(E_MSG, trim(label) // ' inflation:', errstring, source, revision, revdate)
+ endif
+endif
+
end subroutine adaptive_inflate_init
!------------------------------------------------------------------
More information about the Dart-dev
mailing list