From bc5fb81b73be9908ea7274352fada9443f7335ef Mon Sep 17 00:00:00 2001 From: Jeff Curtis Date: Tue, 12 Jul 2022 10:34:42 -0500 Subject: [PATCH 01/15] first pass at deaveraging --- CMakeLists.txt | 7 +++ src/aero_state.F90 | 100 +++++++++++++++++++++++++++++++++++++ src/bin_deaverage_comp.F90 | 90 +++++++++++++++++++++++++++++++++ 3 files changed, 197 insertions(+) create mode 100644 src/bin_deaverage_comp.F90 diff --git a/CMakeLists.txt b/CMakeLists.txt index f623759f6..74065413d 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -265,6 +265,13 @@ add_executable(bin_average_comp src/bin_average_comp.F90) target_link_libraries(bin_average_comp partmclib) +###################################################################### +# bin_deaverage_comp + +add_executable(bin_deaverage_comp src/bin_deaverage_comp.F90) + +target_link_libraries(bin_deaverage_comp partmclib) + ###################################################################### # bin_average_size diff --git a/src/aero_state.F90 b/src/aero_state.F90 index 7c45482e6..588165082 100644 --- a/src/aero_state.F90 +++ b/src/aero_state.F90 @@ -1939,6 +1939,106 @@ subroutine aero_state_bin_average_comp(aero_state, bin_grid, aero_data) end subroutine aero_state_bin_average_comp + subroutine aero_state_bin_deaverage_comp(aero_state, bin_grid, aero_data) + + !> Aerosol state to average. + type(aero_state_t), intent(inout) :: aero_state + !> Bin grid to average within. + type(bin_grid_t), intent(in) :: bin_grid + !> Aerosol data. + type(aero_data_t), intent(in) :: aero_data + + real(kind=dp) :: species_volume_conc(aero_data_n_spec(aero_data)) + real(kind=dp) :: total_volume_conc, particle_volume, num_conc + integer :: i_bin, i_class, i_entry, i_part, i_spec + real(kind=dp) :: particle_fractions(aero_data_n_spec(aero_data)) + integer :: particle_index(aero_data_n_spec(aero_data)) + integer :: n_part_spec, start_val, end_val, i, n_parts + logical :: edge_case + integer, allocatable :: shuffle_particles(:) + + call aero_state_sort(aero_state, aero_data, bin_grid) + + do i_bin = 1,bin_grid_size(bin_grid) + species_volume_conc = 0d0 + total_volume_conc = 0d0 + do i_class = 1,size(aero_state%awa%weight, 2) + do i_entry = 1,integer_varray_n_entry( & + aero_state%aero_sorted%size_class%inverse(i_bin, i_class)) + i_part = aero_state%aero_sorted%size_class%inverse(i_bin, & + i_class)%entry(i_entry) + num_conc = aero_weight_array_num_conc(aero_state%awa, & + aero_state%apa%particle(i_part), aero_data) + particle_volume = aero_particle_volume( & + aero_state%apa%particle(i_part)) + species_volume_conc = species_volume_conc & + + num_conc * aero_state%apa%particle(i_part)%vol + total_volume_conc = total_volume_conc + num_conc * particle_volume + end do + end do + do i_class = 1,size(aero_state%awa%weight, 2) + + ! Get the total particles + n_parts = integer_varray_n_entry( & + aero_state%aero_sorted%size_class%inverse(i_bin, i_class)) + + particle_fractions = n_parts * species_volume_conc / total_volume_conc + print*, particle_fractions + print*, 'total particles', n_parts + + shuffle_particles=[(i,i=1,n_parts)] + call shuffle_array(shuffle_particles, n_parts) + start_val = 1 + edge_case = .false. + do i_spec = 1,aero_data_n_spec(aero_data) + if (species_volume_conc(i_spec) > 0.0d0) then + if (edge_case) then + n_part_spec = prob_round(particle_fractions(i_spec) + 1) + edge_case = .false. + else + n_part_spec = prob_round(particle_fractions(i_spec)) + if (n_part_spec < particle_fractions(i_spec)) then + edge_case = .true. + end if + end if + do i_entry = start_val,min(start_val + n_part_spec - 1, n_parts) + i_part = aero_state%aero_sorted%size_class%inverse(i_bin, & + i_class)%entry(shuffle_particles(i_entry)) + particle_volume = aero_particle_volume( & + aero_state%apa%particle(i_part)) + aero_state%apa%particle(i_part)%vol = 0.0d0 + aero_state%apa%particle(i_part)%vol(i_spec) & + = particle_volume !* species_volume_conc(i_spec) / total_volume_conc + end do + start_val = start_val + n_part_spec + end if + end do + end do + end do + + end subroutine aero_state_bin_deaverage_comp + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + subroutine shuffle_array(array, n_values) + integer :: array(n_values) + integer :: n_values + integer :: temp + integer :: m, i, j + real(kind=dp) :: u + + do i=1,n_values-1 + u = pmc_random() + j = i + FLOOR((n_values-i+1)*u) + temp=array(j) + array(j)=array(i) + array(i)=temp + end do + + print*, array + + end subroutine + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !> Set each aerosol particle to have its original species ratios, diff --git a/src/bin_deaverage_comp.F90 b/src/bin_deaverage_comp.F90 new file mode 100644 index 000000000..418fcdf9e --- /dev/null +++ b/src/bin_deaverage_comp.F90 @@ -0,0 +1,90 @@ +! Copyright (C) 2009-2013 Matthew West +! Licensed under the GNU General Public License version 2 or (at your +! option) any later version. See the file COPYING for details. + +!> \file +!> The bin_average_comp program. + +!> Read a NetCDF file, average the composition of all particles within +!> each bin, and write the data out as another NetCDF file. +program bin_deaverage_comp + + use pmc_aero_state + use pmc_gas_data + use pmc_gas_state + use pmc_env_state + use pmc_aero_data + use pmc_bin_grid + use pmc_output + use pmc_rand + use netcdf + + character(len=1000) :: in_filename, out_prefix + type(bin_grid_t) :: bin_grid + type(aero_data_t) :: aero_data + type(aero_state_t) :: aero_state + type(gas_data_t) :: gas_data + type(gas_state_t) :: gas_state + type(env_state_t) :: env_state + integer :: n_bin, index, i_repeat, output_type + real(kind=dp) :: d_min, d_max, time, del_t + character(len=1000) :: tmp_str + logical :: record_removals, dry_volume, record_optical + character(len=PMC_UUID_LEN) :: uuid + + ! process commandline arguments + if (command_argument_count() .ne. 6) then + write(6,*) 'Usage: bin_average_comp ' & + // '<"wet" or "dry"> ' + write(6,*) '' + write(6,*) ' d_min: minimum bin diameter (m)' + write(6,*) ' d_max: maximum bin diameter (m)' + write(6,*) ' n_bin: number of bins' + write(6,*) ' wet/dry: average wet or dry sizes' + write(6,*) ' input_filename: like scenario_0001_00000001.nc' + write(6,*) ' output_prefix: like scenario_comp_average' + stop 2 + endif + call get_command_argument(1, tmp_str) + d_min = string_to_real(tmp_str) + call get_command_argument(2, tmp_str) + d_max = string_to_real(tmp_str) + call get_command_argument(3, tmp_str) + n_bin = string_to_integer(tmp_str) + call get_command_argument(4, tmp_str) + if (trim(tmp_str) == "wet") then + dry_volume = .false. + elseif (trim(tmp_str) == "dry") then + dry_volume = .true. + else + write(6,*) 'Argument 4 must be "wet" or "dry", not ' & + // trim(tmp_str) + stop 1 + end if + call get_command_argument(5, in_filename) + call get_command_argument(6, out_prefix) + + call pmc_mpi_init() + + call bin_grid_make(bin_grid, BIN_GRID_TYPE_LOG, n_bin, diam2rad(d_min), & + diam2rad(d_max)) + + call input_state(in_filename, index, time, del_t, i_repeat, uuid, & + aero_data, aero_state, gas_data, gas_state, env_state) + + if (dry_volume) then + call aero_state_make_dry(aero_state, aero_data) + end if + + call aero_state_bin_deaverage_comp(aero_state, bin_grid, aero_data) + + output_type = OUTPUT_TYPE_SINGLE + record_removals = .false. + record_optical = .true. + call output_state(out_prefix, output_type, aero_data, aero_state, & + gas_data, gas_state, env_state, index, time, del_t, i_repeat, & + record_removals, record_optical, uuid) + + call pmc_mpi_finalize() + +end program bin_deaverage_comp From 461cccf752262a1db79c3cf7cbe984509e811114 Mon Sep 17 00:00:00 2001 From: Jeff Curtis Date: Tue, 12 Jul 2022 10:55:23 -0500 Subject: [PATCH 02/15] add passing in of groups for deaveraging --- src/aero_state.F90 | 102 ++++++++++++++++++++++++++++++--------------- 1 file changed, 68 insertions(+), 34 deletions(-) diff --git a/src/aero_state.F90 b/src/aero_state.F90 index 588165082..7343c7529 100644 --- a/src/aero_state.F90 +++ b/src/aero_state.F90 @@ -1939,14 +1939,19 @@ subroutine aero_state_bin_average_comp(aero_state, bin_grid, aero_data) end subroutine aero_state_bin_average_comp - subroutine aero_state_bin_deaverage_comp(aero_state, bin_grid, aero_data) +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !> Aerosol state to average. + subroutine aero_state_bin_deaverage_comp(aero_state, bin_grid, aero_data, & + groups) + + !> Aerosol state to de-average. type(aero_state_t), intent(inout) :: aero_state - !> Bin grid to average within. + !> Bin grid to de-average within. type(bin_grid_t), intent(in) :: bin_grid !> Aerosol data. type(aero_data_t), intent(in) :: aero_data + !> Sets of species names to group together. + character(len=*), optional :: groups(:,:) real(kind=dp) :: species_volume_conc(aero_data_n_spec(aero_data)) real(kind=dp) :: total_volume_conc, particle_volume, num_conc @@ -1956,9 +1961,37 @@ subroutine aero_state_bin_deaverage_comp(aero_state, bin_grid, aero_data) integer :: n_part_spec, start_val, end_val, i, n_parts logical :: edge_case integer, allocatable :: shuffle_particles(:) + integer :: species_group_numbers(aero_data_n_spec(aero_data)) + integer :: n_group, i_group, i_name, next_group_number call aero_state_sort(aero_state, aero_data, bin_grid) + if (present(groups)) then + n_group = size(groups, 1) + ! species_group_numbers(i_spec) will give the group number for + ! each species + species_group_numbers = 0 + do i_group = 1, n_group + do i_name = 1, size(groups, 2) + if (len_trim(groups(i_group, i_name)) > 0) then + i_spec = aero_data_spec_by_name(aero_data, & + groups(i_group, i_name)) + call assert_msg(926066862, i_spec > 0, & + "unknown species: " // trim(groups(i_group, i_name))) + species_group_numbers(i_spec) = i_group + end if + end do + end do + ! Assign left overs to their own groups + next_group_number = n_group + 1 + do i_spec = 1,aero_data_n_spec(aero_data) + if (species_group_numbers(i_spec) == 0) then + species_group_numbers(i_spec) = next_group_number + next_group_number = next_group_number + 1 + end if + end do + end if + do i_bin = 1,bin_grid_size(bin_grid) species_volume_conc = 0d0 total_volume_conc = 0d0 @@ -1976,42 +2009,43 @@ subroutine aero_state_bin_deaverage_comp(aero_state, bin_grid, aero_data) total_volume_conc = total_volume_conc + num_conc * particle_volume end do end do + do i_class = 1,size(aero_state%awa%weight, 2) - ! Get the total particles - n_parts = integer_varray_n_entry( & + ! Get the total particles + n_parts = integer_varray_n_entry( & aero_state%aero_sorted%size_class%inverse(i_bin, i_class)) - particle_fractions = n_parts * species_volume_conc / total_volume_conc - print*, particle_fractions - print*, 'total particles', n_parts - - shuffle_particles=[(i,i=1,n_parts)] - call shuffle_array(shuffle_particles, n_parts) - start_val = 1 - edge_case = .false. - do i_spec = 1,aero_data_n_spec(aero_data) - if (species_volume_conc(i_spec) > 0.0d0) then - if (edge_case) then - n_part_spec = prob_round(particle_fractions(i_spec) + 1) - edge_case = .false. - else - n_part_spec = prob_round(particle_fractions(i_spec)) - if (n_part_spec < particle_fractions(i_spec)) then - edge_case = .true. + particle_fractions = n_parts * species_volume_conc / total_volume_conc + print*, particle_fractions + print*, 'total particles', n_parts + + shuffle_particles=[(i,i=1,n_parts)] + call shuffle_array(shuffle_particles, n_parts) + start_val = 1 + edge_case = .false. + do i_spec = 1,aero_data_n_spec(aero_data) + if (species_volume_conc(i_spec) > 0.0d0) then + if (edge_case) then + n_part_spec = prob_round(particle_fractions(i_spec) + 1) + edge_case = .false. + else + n_part_spec = prob_round(particle_fractions(i_spec)) + if (n_part_spec < particle_fractions(i_spec)) then + edge_case = .true. + end if + end if + do i_entry = start_val,min(start_val + n_part_spec - 1, n_parts) + i_part = aero_state%aero_sorted%size_class%inverse(i_bin, & + i_class)%entry(shuffle_particles(i_entry)) + particle_volume = aero_particle_volume( & + aero_state%apa%particle(i_part)) + aero_state%apa%particle(i_part)%vol = 0.0d0 + aero_state%apa%particle(i_part)%vol(i_spec) & + = particle_volume + end do + start_val = start_val + n_part_spec end if - end if - do i_entry = start_val,min(start_val + n_part_spec - 1, n_parts) - i_part = aero_state%aero_sorted%size_class%inverse(i_bin, & - i_class)%entry(shuffle_particles(i_entry)) - particle_volume = aero_particle_volume( & - aero_state%apa%particle(i_part)) - aero_state%apa%particle(i_part)%vol = 0.0d0 - aero_state%apa%particle(i_part)%vol(i_spec) & - = particle_volume !* species_volume_conc(i_spec) / total_volume_conc - end do - start_val = start_val + n_part_spec - end if end do end do end do From e3f3346ee4b9d1485b0598da1f7ff47a2092d75e Mon Sep 17 00:00:00 2001 From: Jeff Curtis Date: Tue, 12 Jul 2022 11:07:37 -0500 Subject: [PATCH 03/15] minor fix --- src/aero_state.F90 | 15 +++++++++++---- 1 file changed, 11 insertions(+), 4 deletions(-) diff --git a/src/aero_state.F90 b/src/aero_state.F90 index 7343c7529..5f40f940f 100644 --- a/src/aero_state.F90 +++ b/src/aero_state.F90 @@ -1963,6 +1963,7 @@ subroutine aero_state_bin_deaverage_comp(aero_state, bin_grid, aero_data, & integer, allocatable :: shuffle_particles(:) integer :: species_group_numbers(aero_data_n_spec(aero_data)) integer :: n_group, i_group, i_name, next_group_number + real(kind=dp), allocatable :: group_masses(:) call aero_state_sort(aero_state, aero_data, bin_grid) @@ -1971,8 +1972,8 @@ subroutine aero_state_bin_deaverage_comp(aero_state, bin_grid, aero_data, & ! species_group_numbers(i_spec) will give the group number for ! each species species_group_numbers = 0 - do i_group = 1, n_group - do i_name = 1, size(groups, 2) + do i_group = 1,n_group + do i_name = 1,size(groups, 2) if (len_trim(groups(i_group, i_name)) > 0) then i_spec = aero_data_spec_by_name(aero_data, & groups(i_group, i_name)) @@ -1993,6 +1994,9 @@ subroutine aero_state_bin_deaverage_comp(aero_state, bin_grid, aero_data, & end if do i_bin = 1,bin_grid_size(bin_grid) +! if (present(groups)) then +! +! else species_volume_conc = 0d0 total_volume_conc = 0d0 do i_class = 1,size(aero_state%awa%weight, 2) @@ -2019,8 +2023,10 @@ subroutine aero_state_bin_deaverage_comp(aero_state, bin_grid, aero_data, & particle_fractions = n_parts * species_volume_conc / total_volume_conc print*, particle_fractions print*, 'total particles', n_parts - - shuffle_particles=[(i,i=1,n_parts)] + allocate(shuffle_particles(n_parts)) + do i_part = 1,n_parts + shuffle_particles(i_part) = i_part + end do call shuffle_array(shuffle_particles, n_parts) start_val = 1 edge_case = .false. @@ -2047,6 +2053,7 @@ subroutine aero_state_bin_deaverage_comp(aero_state, bin_grid, aero_data, & start_val = start_val + n_part_spec end if end do + deallocate(shuffle_particles) end do end do From afd32ebfd20202d03697939c116d78c93b18a26c Mon Sep 17 00:00:00 2001 From: Jeff Curtis Date: Tue, 12 Jul 2022 15:46:28 -0500 Subject: [PATCH 04/15] added groups to deaveraging --- src/aero_state.F90 | 74 +++++++++++++++++++++++++++++++++----- src/bin_deaverage_comp.F90 | 9 ++++- 2 files changed, 73 insertions(+), 10 deletions(-) diff --git a/src/aero_state.F90 b/src/aero_state.F90 index 5f40f940f..2c4688f3e 100644 --- a/src/aero_state.F90 +++ b/src/aero_state.F90 @@ -1956,14 +1956,16 @@ subroutine aero_state_bin_deaverage_comp(aero_state, bin_grid, aero_data, & real(kind=dp) :: species_volume_conc(aero_data_n_spec(aero_data)) real(kind=dp) :: total_volume_conc, particle_volume, num_conc integer :: i_bin, i_class, i_entry, i_part, i_spec - real(kind=dp) :: particle_fractions(aero_data_n_spec(aero_data)) + real(kind=dp) :: particle_fractions(aero_data_n_spec(aero_data)), & + factors(aero_data_n_spec(aero_data)) integer :: particle_index(aero_data_n_spec(aero_data)) integer :: n_part_spec, start_val, end_val, i, n_parts logical :: edge_case integer, allocatable :: shuffle_particles(:) integer :: species_group_numbers(aero_data_n_spec(aero_data)) integer :: n_group, i_group, i_name, next_group_number - real(kind=dp), allocatable :: group_masses(:) + real(kind=dp), allocatable :: group_fractions(:), & + group_volume_conc(:) call aero_state_sort(aero_state, aero_data, bin_grid) @@ -1984,19 +1986,18 @@ subroutine aero_state_bin_deaverage_comp(aero_state, bin_grid, aero_data, & end do end do ! Assign left overs to their own groups - next_group_number = n_group + 1 do i_spec = 1,aero_data_n_spec(aero_data) if (species_group_numbers(i_spec) == 0) then - species_group_numbers(i_spec) = next_group_number - next_group_number = next_group_number + 1 + n_group = n_group + 1 + species_group_numbers(i_spec) = n_group end if end do + print*, species_group_numbers + allocate(group_fractions(n_group)) + allocate(group_volume_conc(n_group)) end if do i_bin = 1,bin_grid_size(bin_grid) -! if (present(groups)) then -! -! else species_volume_conc = 0d0 total_volume_conc = 0d0 do i_class = 1,size(aero_state%awa%weight, 2) @@ -2019,7 +2020,61 @@ subroutine aero_state_bin_deaverage_comp(aero_state, bin_grid, aero_data, & ! Get the total particles n_parts = integer_varray_n_entry( & aero_state%aero_sorted%size_class%inverse(i_bin, i_class)) - + group_volume_conc = 0.0d0 + if (present(groups)) then + do i_spec = 1,aero_data_n_spec(aero_data) + group_volume_conc(species_group_numbers(i_spec)) & + = group_volume_conc(species_group_numbers(i_spec)) & + + species_volume_conc(i_spec) + end do + group_fractions = 0.0d0 + do i_group = 1,n_group + group_fractions(i_group) = n_parts * group_volume_conc(i_group) / total_volume_conc + end do + print*, group_fractions + allocate(shuffle_particles(n_parts)) + do i_part = 1,n_parts + shuffle_particles(i_part) = i_part + end do + call shuffle_array(shuffle_particles, n_parts) + start_val = 1 + edge_case = .false. + do i_group = 1,n_group + if (group_volume_conc(i_group) > 0.0d0) then + if (edge_case) then + n_part_spec = prob_round(group_fractions(i_group) + 1) + edge_case = .false. + else + n_part_spec = prob_round(group_fractions(i_group)) + if (n_part_spec < group_fractions(i_group)) then + edge_case = .true. + end if + end if + + factors = 0.0d0 + do i_spec = 1,aero_data_n_spec(aero_data) + if (species_group_numbers(i_spec) == i_group) then + factors(i_spec) = species_volume_conc(i_spec) / group_volume_conc(i_group) + end if + end do + print*, i_group,'factors', factors + do i_entry = start_val,min(start_val + n_part_spec - 1, & + n_parts) + i_part = aero_state%aero_sorted%size_class%inverse(i_bin, & + i_class)%entry(shuffle_particles(i_entry)) + particle_volume = aero_particle_volume( & + aero_state%apa%particle(i_part)) + aero_state%apa%particle(i_part)%vol = 0.0d0 + do i_spec = 1,aero_data_n_spec(aero_data) + aero_state%apa%particle(i_part)%vol(i_spec) & + = particle_volume * factors(i_spec) + end do + end do + start_val = start_val + n_part_spec + end if + end do + deallocate(shuffle_particles) + else particle_fractions = n_parts * species_volume_conc / total_volume_conc print*, particle_fractions print*, 'total particles', n_parts @@ -2054,6 +2109,7 @@ subroutine aero_state_bin_deaverage_comp(aero_state, bin_grid, aero_data, & end if end do deallocate(shuffle_particles) + end if end do end do diff --git a/src/bin_deaverage_comp.F90 b/src/bin_deaverage_comp.F90 index 418fcdf9e..5dbc3f029 100644 --- a/src/bin_deaverage_comp.F90 +++ b/src/bin_deaverage_comp.F90 @@ -31,6 +31,12 @@ program bin_deaverage_comp character(len=1000) :: tmp_str logical :: record_removals, dry_volume, record_optical character(len=PMC_UUID_LEN) :: uuid + character(len=AERO_NAME_LEN), allocatable :: external_groups(:,:) + + allocate(external_groups(3, 4)) ! 3 groups, max 4 species per group + external_groups(1,:) = ["OC ", "BC ", " ", " "] + external_groups(2,:) = ["API1 ", "API2 ", "LIM1 ", "LIM2 "] + external_groups(3,:) = ["SO4 ", "NO3 ", "NH4 ", " "] ! process commandline arguments if (command_argument_count() .ne. 6) then @@ -76,7 +82,8 @@ program bin_deaverage_comp call aero_state_make_dry(aero_state, aero_data) end if - call aero_state_bin_deaverage_comp(aero_state, bin_grid, aero_data) + call aero_state_bin_deaverage_comp(aero_state, bin_grid, aero_data, & + external_groups) output_type = OUTPUT_TYPE_SINGLE record_removals = .false. From fd97d68e2a06a83c778ae06d3ab958ff56cb6e8d Mon Sep 17 00:00:00 2001 From: Jeff Curtis Date: Tue, 12 Jul 2022 15:58:39 -0500 Subject: [PATCH 05/15] move allocates to avoid warnings --- src/aero_state.F90 | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/src/aero_state.F90 b/src/aero_state.F90 index 2c4688f3e..fd8e39187 100644 --- a/src/aero_state.F90 +++ b/src/aero_state.F90 @@ -1993,8 +1993,6 @@ subroutine aero_state_bin_deaverage_comp(aero_state, bin_grid, aero_data, & end if end do print*, species_group_numbers - allocate(group_fractions(n_group)) - allocate(group_volume_conc(n_group)) end if do i_bin = 1,bin_grid_size(bin_grid) @@ -2020,13 +2018,15 @@ subroutine aero_state_bin_deaverage_comp(aero_state, bin_grid, aero_data, & ! Get the total particles n_parts = integer_varray_n_entry( & aero_state%aero_sorted%size_class%inverse(i_bin, i_class)) - group_volume_conc = 0.0d0 if (present(groups)) then + allocate(group_volume_conc(n_group)) + group_volume_conc = 0.0d0 do i_spec = 1,aero_data_n_spec(aero_data) group_volume_conc(species_group_numbers(i_spec)) & = group_volume_conc(species_group_numbers(i_spec)) & + species_volume_conc(i_spec) end do + allocate(group_fractions(n_group)) group_fractions = 0.0d0 do i_group = 1,n_group group_fractions(i_group) = n_parts * group_volume_conc(i_group) / total_volume_conc @@ -2074,6 +2074,8 @@ subroutine aero_state_bin_deaverage_comp(aero_state, bin_grid, aero_data, & end if end do deallocate(shuffle_particles) + deallocate(group_fractions) + deallocate(group_volume_conc) else particle_fractions = n_parts * species_volume_conc / total_volume_conc print*, particle_fractions From 8470a00652de52c1f252717b8f2b71c5d7d3a47e Mon Sep 17 00:00:00 2001 From: Jeff Curtis Date: Fri, 15 Jul 2022 15:50:40 -0500 Subject: [PATCH 06/15] move shuffle subroutine and clean up code --- src/aero_state.F90 | 97 ++++++++++++++++++---------------------------- src/rand.F90 | 23 +++++++++++ 2 files changed, 60 insertions(+), 60 deletions(-) diff --git a/src/aero_state.F90 b/src/aero_state.F90 index fd8e39187..2e4458b96 100644 --- a/src/aero_state.F90 +++ b/src/aero_state.F90 @@ -1992,7 +1992,6 @@ subroutine aero_state_bin_deaverage_comp(aero_state, bin_grid, aero_data, & species_group_numbers(i_spec) = n_group end if end do - print*, species_group_numbers end if do i_bin = 1,bin_grid_size(bin_grid) @@ -2029,14 +2028,14 @@ subroutine aero_state_bin_deaverage_comp(aero_state, bin_grid, aero_data, & allocate(group_fractions(n_group)) group_fractions = 0.0d0 do i_group = 1,n_group - group_fractions(i_group) = n_parts * group_volume_conc(i_group) / total_volume_conc + group_fractions(i_group) = n_parts & + * group_volume_conc(i_group) / total_volume_conc end do - print*, group_fractions allocate(shuffle_particles(n_parts)) do i_part = 1,n_parts shuffle_particles(i_part) = i_part end do - call shuffle_array(shuffle_particles, n_parts) + call pmc_rand_shuffle_array(shuffle_particles, n_parts) start_val = 1 edge_case = .false. do i_group = 1,n_group @@ -2054,10 +2053,10 @@ subroutine aero_state_bin_deaverage_comp(aero_state, bin_grid, aero_data, & factors = 0.0d0 do i_spec = 1,aero_data_n_spec(aero_data) if (species_group_numbers(i_spec) == i_group) then - factors(i_spec) = species_volume_conc(i_spec) / group_volume_conc(i_group) + factors(i_spec) = species_volume_conc(i_spec) & + / group_volume_conc(i_group) end if end do - print*, i_group,'factors', factors do i_entry = start_val,min(start_val + n_part_spec - 1, & n_parts) i_part = aero_state%aero_sorted%size_class%inverse(i_bin, & @@ -2077,67 +2076,45 @@ subroutine aero_state_bin_deaverage_comp(aero_state, bin_grid, aero_data, & deallocate(group_fractions) deallocate(group_volume_conc) else - particle_fractions = n_parts * species_volume_conc / total_volume_conc - print*, particle_fractions - print*, 'total particles', n_parts - allocate(shuffle_particles(n_parts)) - do i_part = 1,n_parts - shuffle_particles(i_part) = i_part - end do - call shuffle_array(shuffle_particles, n_parts) - start_val = 1 - edge_case = .false. - do i_spec = 1,aero_data_n_spec(aero_data) - if (species_volume_conc(i_spec) > 0.0d0) then - if (edge_case) then - n_part_spec = prob_round(particle_fractions(i_spec) + 1) - edge_case = .false. - else - n_part_spec = prob_round(particle_fractions(i_spec)) - if (n_part_spec < particle_fractions(i_spec)) then - edge_case = .true. + particle_fractions = n_parts * species_volume_conc / total_volume_conc + allocate(shuffle_particles(n_parts)) + do i_part = 1,n_parts + shuffle_particles(i_part) = i_part + end do + call pmc_rand_shuffle_array(shuffle_particles, n_parts) + start_val = 1 + edge_case = .false. + do i_spec = 1,aero_data_n_spec(aero_data) + if (species_volume_conc(i_spec) > 0.0d0) then + if (edge_case) then + n_part_spec = prob_round(particle_fractions(i_spec) + 1) + edge_case = .false. + else + n_part_spec = prob_round(particle_fractions(i_spec)) + if (n_part_spec < particle_fractions(i_spec)) then + edge_case = .true. + end if end if - end if - do i_entry = start_val,min(start_val + n_part_spec - 1, n_parts) - i_part = aero_state%aero_sorted%size_class%inverse(i_bin, & - i_class)%entry(shuffle_particles(i_entry)) - particle_volume = aero_particle_volume( & - aero_state%apa%particle(i_part)) - aero_state%apa%particle(i_part)%vol = 0.0d0 - aero_state%apa%particle(i_part)%vol(i_spec) & - = particle_volume - end do - start_val = start_val + n_part_spec - end if - end do - deallocate(shuffle_particles) + do i_entry = start_val,min(start_val + n_part_spec - 1, & + n_parts) + i_part = aero_state%aero_sorted%size_class%inverse(i_bin, & + i_class)%entry(shuffle_particles(i_entry)) + particle_volume = aero_particle_volume( & + aero_state%apa%particle(i_part)) + aero_state%apa%particle(i_part)%vol = 0.0d0 + aero_state%apa%particle(i_part)%vol(i_spec) & + = particle_volume + end do + start_val = start_val + n_part_spec + end if + end do + deallocate(shuffle_particles) end if end do end do end subroutine aero_state_bin_deaverage_comp -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - subroutine shuffle_array(array, n_values) - integer :: array(n_values) - integer :: n_values - integer :: temp - integer :: m, i, j - real(kind=dp) :: u - - do i=1,n_values-1 - u = pmc_random() - j = i + FLOOR((n_values-i+1)*u) - temp=array(j) - array(j)=array(i) - array(i)=temp - end do - - print*, array - - end subroutine - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !> Set each aerosol particle to have its original species ratios, diff --git a/src/rand.F90 b/src/rand.F90 index e7640b75c..fb1d68b09 100644 --- a/src/rand.F90 +++ b/src/rand.F90 @@ -206,6 +206,29 @@ end function pmc_rand_int_gsl end function pmc_rand_int +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + !> Shuffles randomly an array of integer values. + subroutine pmc_rand_shuffle_array(array, n_values) + !> Array of values. + integer, intent(inout) :: array(n_values) + !> Number of values. + integer, intent(in) :: n_values + + integer :: temp + integer :: i, j + real(kind=dp) :: u + + do i=1,n_values-1 + u = pmc_random() + j = i + floor((n_values-i+1)*u) + temp=array(j) + array(j)=array(i) + array(i)=temp + end do + + end subroutine pmc_rand_shuffle_array + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !> Round val to \c floor(val) or \c ceiling(val) with probability From c78c549107e42ee0ce2b45d0b41647237d4c2acf Mon Sep 17 00:00:00 2001 From: Jeff Curtis Date: Mon, 18 Jul 2022 16:10:34 -0500 Subject: [PATCH 07/15] clean up code. add gsl shuffle --- src/aero_state.F90 | 153 ++++++++++++++++++------------------- src/bin_deaverage_comp.F90 | 2 + src/rand.F90 | 26 +++++++ src/rand_gsl.c | 19 +++++ 4 files changed, 121 insertions(+), 79 deletions(-) diff --git a/src/aero_state.F90 b/src/aero_state.F90 index 2e4458b96..8da0083e6 100644 --- a/src/aero_state.F90 +++ b/src/aero_state.F90 @@ -2017,20 +2017,7 @@ subroutine aero_state_bin_deaverage_comp(aero_state, bin_grid, aero_data, & ! Get the total particles n_parts = integer_varray_n_entry( & aero_state%aero_sorted%size_class%inverse(i_bin, i_class)) - if (present(groups)) then - allocate(group_volume_conc(n_group)) - group_volume_conc = 0.0d0 - do i_spec = 1,aero_data_n_spec(aero_data) - group_volume_conc(species_group_numbers(i_spec)) & - = group_volume_conc(species_group_numbers(i_spec)) & - + species_volume_conc(i_spec) - end do - allocate(group_fractions(n_group)) - group_fractions = 0.0d0 - do i_group = 1,n_group - group_fractions(i_group) = n_parts & - * group_volume_conc(i_group) / total_volume_conc - end do + if (n_parts > 0) then allocate(shuffle_particles(n_parts)) do i_part = 1,n_parts shuffle_particles(i_part) = i_part @@ -2038,76 +2025,84 @@ subroutine aero_state_bin_deaverage_comp(aero_state, bin_grid, aero_data, & call pmc_rand_shuffle_array(shuffle_particles, n_parts) start_val = 1 edge_case = .false. - do i_group = 1,n_group - if (group_volume_conc(i_group) > 0.0d0) then - if (edge_case) then - n_part_spec = prob_round(group_fractions(i_group) + 1) - edge_case = .false. - else - n_part_spec = prob_round(group_fractions(i_group)) - if (n_part_spec < group_fractions(i_group)) then - edge_case = .true. + if (present(groups)) then + allocate(group_volume_conc(n_group)) + group_volume_conc = 0.0d0 + do i_spec = 1,aero_data_n_spec(aero_data) + group_volume_conc(species_group_numbers(i_spec)) & + = group_volume_conc(species_group_numbers(i_spec)) & + + species_volume_conc(i_spec) + end do + allocate(group_fractions(n_group)) + group_fractions = 0.0d0 + do i_group = 1,n_group + group_fractions(i_group) = n_parts & + * group_volume_conc(i_group) / total_volume_conc + end do + do i_group = 1,n_group + if (group_volume_conc(i_group) > 0.0d0) then + if (edge_case) then + n_part_spec = prob_round(group_fractions(i_group) + 1) + edge_case = .false. + else + n_part_spec = prob_round(group_fractions(i_group)) + if (n_part_spec < group_fractions(i_group)) then + edge_case = .true. + end if end if - end if - factors = 0.0d0 - do i_spec = 1,aero_data_n_spec(aero_data) - if (species_group_numbers(i_spec) == i_group) then - factors(i_spec) = species_volume_conc(i_spec) & - / group_volume_conc(i_group) - end if - end do - do i_entry = start_val,min(start_val + n_part_spec - 1, & - n_parts) - i_part = aero_state%aero_sorted%size_class%inverse(i_bin, & - i_class)%entry(shuffle_particles(i_entry)) - particle_volume = aero_particle_volume( & - aero_state%apa%particle(i_part)) - aero_state%apa%particle(i_part)%vol = 0.0d0 + factors = 0.0d0 do i_spec = 1,aero_data_n_spec(aero_data) + if (species_group_numbers(i_spec) == i_group) then + factors(i_spec) = species_volume_conc(i_spec) & + / group_volume_conc(i_group) + end if + end do + do i_entry = start_val,min(start_val + n_part_spec - 1, & + n_parts) + i_part = aero_state%aero_sorted%size_class%inverse( & + i_bin, i_class)%entry(shuffle_particles(i_entry)) + particle_volume = aero_particle_volume( & + aero_state%apa%particle(i_part)) + aero_state%apa%particle(i_part)%vol = 0.0d0 + do i_spec = 1,aero_data_n_spec(aero_data) + aero_state%apa%particle(i_part)%vol(i_spec) & + = particle_volume * factors(i_spec) + end do + end do + start_val = start_val + n_part_spec + end if + end do + deallocate(group_fractions) + deallocate(group_volume_conc) + else + particle_fractions = n_parts * species_volume_conc & + / total_volume_conc + do i_spec = 1,aero_data_n_spec(aero_data) + if (species_volume_conc(i_spec) > 0.0d0) then + if (edge_case) then + n_part_spec = prob_round(particle_fractions(i_spec)+1) + edge_case = .false. + else + n_part_spec = prob_round(particle_fractions(i_spec)) + if (n_part_spec < particle_fractions(i_spec)) then + edge_case = .true. + end if + end if + do i_entry = start_val,min(start_val + n_part_spec - 1, & + n_parts) + i_part = aero_state%aero_sorted%size_class%inverse( & + i_bin, i_class)%entry(shuffle_particles(i_entry)) + particle_volume = aero_particle_volume( & + aero_state%apa%particle(i_part)) + aero_state%apa%particle(i_part)%vol = 0.0d0 aero_state%apa%particle(i_part)%vol(i_spec) & - = particle_volume * factors(i_spec) + = particle_volume end do - end do - start_val = start_val + n_part_spec - end if - end do - deallocate(shuffle_particles) - deallocate(group_fractions) - deallocate(group_volume_conc) - else - particle_fractions = n_parts * species_volume_conc / total_volume_conc - allocate(shuffle_particles(n_parts)) - do i_part = 1,n_parts - shuffle_particles(i_part) = i_part - end do - call pmc_rand_shuffle_array(shuffle_particles, n_parts) - start_val = 1 - edge_case = .false. - do i_spec = 1,aero_data_n_spec(aero_data) - if (species_volume_conc(i_spec) > 0.0d0) then - if (edge_case) then - n_part_spec = prob_round(particle_fractions(i_spec) + 1) - edge_case = .false. - else - n_part_spec = prob_round(particle_fractions(i_spec)) - if (n_part_spec < particle_fractions(i_spec)) then - edge_case = .true. - end if - end if - do i_entry = start_val,min(start_val + n_part_spec - 1, & - n_parts) - i_part = aero_state%aero_sorted%size_class%inverse(i_bin, & - i_class)%entry(shuffle_particles(i_entry)) - particle_volume = aero_particle_volume( & - aero_state%apa%particle(i_part)) - aero_state%apa%particle(i_part)%vol = 0.0d0 - aero_state%apa%particle(i_part)%vol(i_spec) & - = particle_volume - end do - start_val = start_val + n_part_spec - end if - end do + start_val = start_val + n_part_spec + end if + end do + end if deallocate(shuffle_particles) end if end do diff --git a/src/bin_deaverage_comp.F90 b/src/bin_deaverage_comp.F90 index 5dbc3f029..b65612d38 100644 --- a/src/bin_deaverage_comp.F90 +++ b/src/bin_deaverage_comp.F90 @@ -72,6 +72,8 @@ program bin_deaverage_comp call pmc_mpi_init() + call pmc_srand(0, pmc_mpi_rank()) + call bin_grid_make(bin_grid, BIN_GRID_TYPE_LOG, n_bin, diam2rad(d_min), & diam2rad(d_max)) diff --git a/src/rand.F90 b/src/rand.F90 index fb1d68b09..0b504e1ff 100644 --- a/src/rand.F90 +++ b/src/rand.F90 @@ -219,6 +219,31 @@ subroutine pmc_rand_shuffle_array(array, n_values) integer :: i, j real(kind=dp) :: u +#ifdef PMC_USE_GSL + integer(kind=c_int) :: n_c + integer(kind=c_int) :: array_c(n_values) +#endif + +#ifdef PMC_USE_GSL +#ifndef DOXYGEN_SKIP_DOC + interface + integer(kind=c_int) function pmc_rand_shuffle_gsl(array, n) bind(c) + use iso_c_binding + integer(kind=c_int), value :: n + integer(kind=c_int) :: array(n) + end function pmc_rand_shuffle_gsl + end interface +#endif +#endif + +#ifdef PMC_USE_GSL + n_c = int(n_values, kind=c_int) + array_c = int(array, kind=c_int) + call rand_check_gsl(388234845, pmc_rand_shuffle_gsl(array_c, n_c)) + array = int(array_c) + print*, 'in fortran after call' + print*, array +#else do i=1,n_values-1 u = pmc_random() j = i + floor((n_values-i+1)*u) @@ -226,6 +251,7 @@ subroutine pmc_rand_shuffle_array(array, n_values) array(j)=array(i) array(i)=temp end do +#endif end subroutine pmc_rand_shuffle_array diff --git a/src/rand_gsl.c b/src/rand_gsl.c index db200abbf..be866f373 100644 --- a/src/rand_gsl.c +++ b/src/rand_gsl.c @@ -150,3 +150,22 @@ int pmc_rand_binomial_gsl(int n, double p, int *harvest) *harvest = gsl_ran_binomial(pmc_rand_gsl_rng, p, u); return PMC_RAND_GSL_SUCCESS; } + +int pmc_rand_shuffle_gsl(int *array, int n) +{ + printf("array size %u \n", n); + for (int i = 0 ; i < n ; i++) + { + printf("%u\n", array[i]); + } + if (!pmc_rand_gsl_rng) { + return PMC_RAND_GSL_NOT_INIT; + } + gsl_ran_shuffle(pmc_rand_gsl_rng, array, n, sizeof (int)); + for (int i = 0 ; i < n ; i++) + { + printf("%u\n", array[i]); + } + + return PMC_RAND_GSL_SUCCESS; +} From a541d2ca0cbc9f96e13eefcaa80f3d2c4c85007d Mon Sep 17 00:00:00 2001 From: Jeff Curtis Date: Mon, 18 Jul 2022 21:44:45 -0500 Subject: [PATCH 08/15] fix bug for gsl compiled comp-averaging --- src/bin_average_comp.F90 | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/bin_average_comp.F90 b/src/bin_average_comp.F90 index 9d3d32902..0183a3fe7 100644 --- a/src/bin_average_comp.F90 +++ b/src/bin_average_comp.F90 @@ -66,6 +66,8 @@ program bin_average_comp call pmc_mpi_init() + call pmc_srand(0, pmc_mpi_rank()) + call bin_grid_make(bin_grid, BIN_GRID_TYPE_LOG, n_bin, diam2rad(d_min), & diam2rad(d_max)) From 1e3ac5cf3e03b829594d346bc778a01b59311a29 Mon Sep 17 00:00:00 2001 From: Jeff Curtis Date: Mon, 18 Jul 2022 22:17:59 -0500 Subject: [PATCH 09/15] refactor code by creating groups each with a single species if no groups specified --- src/aero_state.F90 | 136 +++++++++++++++++++-------------------------- 1 file changed, 57 insertions(+), 79 deletions(-) diff --git a/src/aero_state.F90 b/src/aero_state.F90 index 8da0083e6..a5a011131 100644 --- a/src/aero_state.F90 +++ b/src/aero_state.F90 @@ -1941,6 +1941,12 @@ end subroutine aero_state_bin_average_comp !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + !> Set each aerosol particle to have its original total volume, but + !> species volume ratios given by the group species volume ratio + !> within each bin. If no groups are specified, each particle contain a + !> single species (fully externally mixed). If groups are specified, + !> particles contain the volume fractions of the speices within the group. + !> This preserves per-particle total volumes. subroutine aero_state_bin_deaverage_comp(aero_state, bin_grid, aero_data, & groups) @@ -1956,14 +1962,12 @@ subroutine aero_state_bin_deaverage_comp(aero_state, bin_grid, aero_data, & real(kind=dp) :: species_volume_conc(aero_data_n_spec(aero_data)) real(kind=dp) :: total_volume_conc, particle_volume, num_conc integer :: i_bin, i_class, i_entry, i_part, i_spec - real(kind=dp) :: particle_fractions(aero_data_n_spec(aero_data)), & - factors(aero_data_n_spec(aero_data)) - integer :: particle_index(aero_data_n_spec(aero_data)) + real(kind=dp) :: factors(aero_data_n_spec(aero_data)) integer :: n_part_spec, start_val, end_val, i, n_parts logical :: edge_case integer, allocatable :: shuffle_particles(:) integer :: species_group_numbers(aero_data_n_spec(aero_data)) - integer :: n_group, i_group, i_name, next_group_number + integer :: n_group, i_group, i_name real(kind=dp), allocatable :: group_fractions(:), & group_volume_conc(:) @@ -1992,6 +1996,10 @@ subroutine aero_state_bin_deaverage_comp(aero_state, bin_grid, aero_data, & species_group_numbers(i_spec) = n_group end if end do + else + do i_spec = 1,aero_data_n_spec(aero_data) + species_group_numbers(i_spec) = i_spec + end do end if do i_bin = 1,bin_grid_size(bin_grid) @@ -2013,8 +2021,6 @@ subroutine aero_state_bin_deaverage_comp(aero_state, bin_grid, aero_data, & end do do i_class = 1,size(aero_state%awa%weight, 2) - - ! Get the total particles n_parts = integer_varray_n_entry( & aero_state%aero_sorted%size_class%inverse(i_bin, i_class)) if (n_parts > 0) then @@ -2025,84 +2031,56 @@ subroutine aero_state_bin_deaverage_comp(aero_state, bin_grid, aero_data, & call pmc_rand_shuffle_array(shuffle_particles, n_parts) start_val = 1 edge_case = .false. - if (present(groups)) then - allocate(group_volume_conc(n_group)) - group_volume_conc = 0.0d0 - do i_spec = 1,aero_data_n_spec(aero_data) - group_volume_conc(species_group_numbers(i_spec)) & - = group_volume_conc(species_group_numbers(i_spec)) & - + species_volume_conc(i_spec) - end do - allocate(group_fractions(n_group)) - group_fractions = 0.0d0 - do i_group = 1,n_group - group_fractions(i_group) = n_parts & - * group_volume_conc(i_group) / total_volume_conc - end do - do i_group = 1,n_group - if (group_volume_conc(i_group) > 0.0d0) then - if (edge_case) then - n_part_spec = prob_round(group_fractions(i_group) + 1) - edge_case = .false. - else - n_part_spec = prob_round(group_fractions(i_group)) - if (n_part_spec < group_fractions(i_group)) then - edge_case = .true. - end if + allocate(group_volume_conc(n_group)) + group_volume_conc = 0.0d0 + do i_spec = 1,aero_data_n_spec(aero_data) + group_volume_conc(species_group_numbers(i_spec)) & + = group_volume_conc(species_group_numbers(i_spec)) & + + species_volume_conc(i_spec) + end do + allocate(group_fractions(n_group)) + group_fractions = 0.0d0 + do i_group = 1,n_group + group_fractions(i_group) = n_parts & + * group_volume_conc(i_group) / total_volume_conc + end do + do i_group = 1,n_group + if (group_volume_conc(i_group) > 0.0d0) then + if (edge_case) then + n_part_spec = prob_round(group_fractions(i_group) + 1) + edge_case = .false. + else + n_part_spec = prob_round(group_fractions(i_group)) + if (n_part_spec < group_fractions(i_group)) then + edge_case = .true. end if - - factors = 0.0d0 - do i_spec = 1,aero_data_n_spec(aero_data) - if (species_group_numbers(i_spec) == i_group) then - factors(i_spec) = species_volume_conc(i_spec) & - / group_volume_conc(i_group) - end if - end do - do i_entry = start_val,min(start_val + n_part_spec - 1, & - n_parts) - i_part = aero_state%aero_sorted%size_class%inverse( & - i_bin, i_class)%entry(shuffle_particles(i_entry)) - particle_volume = aero_particle_volume( & - aero_state%apa%particle(i_part)) - aero_state%apa%particle(i_part)%vol = 0.0d0 - do i_spec = 1,aero_data_n_spec(aero_data) - aero_state%apa%particle(i_part)%vol(i_spec) & - = particle_volume * factors(i_spec) - end do - end do - start_val = start_val + n_part_spec end if - end do - deallocate(group_fractions) - deallocate(group_volume_conc) - else - particle_fractions = n_parts * species_volume_conc & - / total_volume_conc - do i_spec = 1,aero_data_n_spec(aero_data) - if (species_volume_conc(i_spec) > 0.0d0) then - if (edge_case) then - n_part_spec = prob_round(particle_fractions(i_spec)+1) - edge_case = .false. - else - n_part_spec = prob_round(particle_fractions(i_spec)) - if (n_part_spec < particle_fractions(i_spec)) then - edge_case = .true. - end if + + factors = 0.0d0 + do i_spec = 1,aero_data_n_spec(aero_data) + if (species_group_numbers(i_spec) == i_group) then + factors(i_spec) = species_volume_conc(i_spec) & + / group_volume_conc(i_group) end if - do i_entry = start_val,min(start_val + n_part_spec - 1, & - n_parts) - i_part = aero_state%aero_sorted%size_class%inverse( & - i_bin, i_class)%entry(shuffle_particles(i_entry)) - particle_volume = aero_particle_volume( & - aero_state%apa%particle(i_part)) - aero_state%apa%particle(i_part)%vol = 0.0d0 + end do + + do i_entry = start_val,min(start_val + n_part_spec - 1, & + n_parts) + i_part = aero_state%aero_sorted%size_class%inverse( & + i_bin, i_class)%entry(shuffle_particles(i_entry)) + particle_volume = aero_particle_volume( & + aero_state%apa%particle(i_part)) + aero_state%apa%particle(i_part)%vol = 0.0d0 + do i_spec = 1,aero_data_n_spec(aero_data) aero_state%apa%particle(i_part)%vol(i_spec) & - = particle_volume + = particle_volume * factors(i_spec) end do - start_val = start_val + n_part_spec - end if - end do - end if + end do + start_val = start_val + n_part_spec + end if + end do + deallocate(group_fractions) + deallocate(group_volume_conc) deallocate(shuffle_particles) end if end do From 1ea3df59d7f068e97f61d695b737d3ad012b9617 Mon Sep 17 00:00:00 2001 From: Jeff Curtis Date: Tue, 19 Jul 2022 16:33:40 -0500 Subject: [PATCH 10/15] added option for naive algorithm --- src/aero_state.F90 | 116 +++++++++++++++++++++++++++++++-------------- 1 file changed, 81 insertions(+), 35 deletions(-) diff --git a/src/aero_state.F90 b/src/aero_state.F90 index a5a011131..d9e4522c3 100644 --- a/src/aero_state.F90 +++ b/src/aero_state.F90 @@ -1970,6 +1970,12 @@ subroutine aero_state_bin_deaverage_comp(aero_state, bin_grid, aero_data, & integer :: n_group, i_group, i_name real(kind=dp), allocatable :: group_fractions(:), & group_volume_conc(:) + ! Naive variables + logical, parameter :: naive_algorithm = .false. + real(kind=dp), allocatable :: cumulative_values(:) + real(kind=dp), allocatable :: factors_2d(:,:) + integer, allocatable :: particle_group(:) + real(kind=dp) :: x call aero_state_sort(aero_state, aero_data, bin_grid) @@ -2024,13 +2030,6 @@ subroutine aero_state_bin_deaverage_comp(aero_state, bin_grid, aero_data, & n_parts = integer_varray_n_entry( & aero_state%aero_sorted%size_class%inverse(i_bin, i_class)) if (n_parts > 0) then - allocate(shuffle_particles(n_parts)) - do i_part = 1,n_parts - shuffle_particles(i_part) = i_part - end do - call pmc_rand_shuffle_array(shuffle_particles, n_parts) - start_val = 1 - edge_case = .false. allocate(group_volume_conc(n_group)) group_volume_conc = 0.0d0 do i_spec = 1,aero_data_n_spec(aero_data) @@ -2044,44 +2043,91 @@ subroutine aero_state_bin_deaverage_comp(aero_state, bin_grid, aero_data, & group_fractions(i_group) = n_parts & * group_volume_conc(i_group) / total_volume_conc end do - do i_group = 1,n_group - if (group_volume_conc(i_group) > 0.0d0) then - if (edge_case) then - n_part_spec = prob_round(group_fractions(i_group) + 1) - edge_case = .false. - else - n_part_spec = prob_round(group_fractions(i_group)) - if (n_part_spec < group_fractions(i_group)) then - edge_case = .true. - end if - end if - - factors = 0.0d0 + if (naive_algorithm) then + allocate(cumulative_values(n_group+1)) + cumulative_values = 0.0d0 + do i_group = 2,n_group + 1 + cumulative_values(i_group) = cumulative_values(i_group-1) & + + group_fractions(i_group-1) + end do + allocate(particle_group(n_parts)) + do i_part = 1,n_parts + x = pmc_random() * n_parts + particle_group(i_part) = find_1d(n_group+1, cumulative_values, x) + end do + allocate(factors_2d(n_group,aero_data_n_spec(aero_data))) + factors_2d = 0.0d0 + do i_group = 1,n_group do i_spec = 1,aero_data_n_spec(aero_data) if (species_group_numbers(i_spec) == i_group) then - factors(i_spec) = species_volume_conc(i_spec) & + factors_2d(i_group,i_spec) = species_volume_conc(i_spec) & / group_volume_conc(i_group) end if end do + end do + + do i_entry = 1,n_parts + i_part = aero_state%aero_sorted%size_class%inverse( & + i_bin, i_class)%entry(i_entry) + particle_volume = aero_particle_volume( & + aero_state%apa%particle(i_part)) + aero_state%apa%particle(i_part)%vol = 0.0d0 + i_group = particle_group(i_entry) + do i_spec = 1,aero_data_n_spec(aero_data) + aero_state%apa%particle(i_part)%vol(i_spec) & + = particle_volume * factors_2d(i_group,i_spec) + end do + end do + deallocate(particle_group) + deallocate(cumulative_values) + deallocate(factors_2d) + else + allocate(shuffle_particles(n_parts)) + do i_part = 1,n_parts + shuffle_particles(i_part) = i_part + end do + call pmc_rand_shuffle_array(shuffle_particles, n_parts) + start_val = 1 + edge_case = .false. + do i_group = 1,n_group + if (group_volume_conc(i_group) > 0.0d0) then + if (edge_case) then + n_part_spec = prob_round(group_fractions(i_group) + 1) + edge_case = .false. + else + n_part_spec = prob_round(group_fractions(i_group)) + if (n_part_spec < group_fractions(i_group)) then + edge_case = .true. + end if + end if - do i_entry = start_val,min(start_val + n_part_spec - 1, & - n_parts) - i_part = aero_state%aero_sorted%size_class%inverse( & - i_bin, i_class)%entry(shuffle_particles(i_entry)) - particle_volume = aero_particle_volume( & - aero_state%apa%particle(i_part)) - aero_state%apa%particle(i_part)%vol = 0.0d0 + factors = 0.0d0 do i_spec = 1,aero_data_n_spec(aero_data) - aero_state%apa%particle(i_part)%vol(i_spec) & - = particle_volume * factors(i_spec) + if (species_group_numbers(i_spec) == i_group) then + factors(i_spec) = species_volume_conc(i_spec) & + / group_volume_conc(i_group) + end if end do - end do - start_val = start_val + n_part_spec - end if - end do + + do i_entry = start_val,min(start_val + n_part_spec - 1, & + n_parts) + i_part = aero_state%aero_sorted%size_class%inverse( & + i_bin, i_class)%entry(shuffle_particles(i_entry)) + particle_volume = aero_particle_volume( & + aero_state%apa%particle(i_part)) + aero_state%apa%particle(i_part)%vol = 0.0d0 + do i_spec = 1,aero_data_n_spec(aero_data) + aero_state%apa%particle(i_part)%vol(i_spec) & + = particle_volume * factors(i_spec) + end do + end do + start_val = start_val + n_part_spec + end if + end do + deallocate(shuffle_particles) + end if deallocate(group_fractions) deallocate(group_volume_conc) - deallocate(shuffle_particles) end if end do end do From 95a010dd40a1a7eda6b73a6071a12711d3f965e2 Mon Sep 17 00:00:00 2001 From: Jeff Curtis Date: Tue, 19 Jul 2022 16:34:14 -0500 Subject: [PATCH 11/15] remove printing in random shuffle code --- src/rand.F90 | 2 -- src/rand_gsl.c | 10 ---------- 2 files changed, 12 deletions(-) diff --git a/src/rand.F90 b/src/rand.F90 index 0b504e1ff..1e4472ef9 100644 --- a/src/rand.F90 +++ b/src/rand.F90 @@ -241,8 +241,6 @@ end function pmc_rand_shuffle_gsl array_c = int(array, kind=c_int) call rand_check_gsl(388234845, pmc_rand_shuffle_gsl(array_c, n_c)) array = int(array_c) - print*, 'in fortran after call' - print*, array #else do i=1,n_values-1 u = pmc_random() diff --git a/src/rand_gsl.c b/src/rand_gsl.c index be866f373..c14bcd62a 100644 --- a/src/rand_gsl.c +++ b/src/rand_gsl.c @@ -153,19 +153,9 @@ int pmc_rand_binomial_gsl(int n, double p, int *harvest) int pmc_rand_shuffle_gsl(int *array, int n) { - printf("array size %u \n", n); - for (int i = 0 ; i < n ; i++) - { - printf("%u\n", array[i]); - } if (!pmc_rand_gsl_rng) { return PMC_RAND_GSL_NOT_INIT; } gsl_ran_shuffle(pmc_rand_gsl_rng, array, n, sizeof (int)); - for (int i = 0 ; i < n ; i++) - { - printf("%u\n", array[i]); - } - return PMC_RAND_GSL_SUCCESS; } From c7aec39a3aa8f09b9071cbd5e2617d164868645f Mon Sep 17 00:00:00 2001 From: Jeff Curtis Date: Thu, 28 Jul 2022 17:59:17 -0500 Subject: [PATCH 12/15] fix minor error in low variance algorithm. clean up code --- src/aero_state.F90 | 86 ++++++++++++++++++++++++++++------------------ 1 file changed, 52 insertions(+), 34 deletions(-) diff --git a/src/aero_state.F90 b/src/aero_state.F90 index d9e4522c3..16f354cb9 100644 --- a/src/aero_state.F90 +++ b/src/aero_state.F90 @@ -1963,19 +1963,19 @@ subroutine aero_state_bin_deaverage_comp(aero_state, bin_grid, aero_data, & real(kind=dp) :: total_volume_conc, particle_volume, num_conc integer :: i_bin, i_class, i_entry, i_part, i_spec real(kind=dp) :: factors(aero_data_n_spec(aero_data)) - integer :: n_part_spec, start_val, end_val, i, n_parts - logical :: edge_case + integer :: n_part_spec, start_val, i, n_part integer, allocatable :: shuffle_particles(:) integer :: species_group_numbers(aero_data_n_spec(aero_data)) integer :: n_group, i_group, i_name real(kind=dp), allocatable :: group_fractions(:), & group_volume_conc(:) - ! Naive variables - logical, parameter :: naive_algorithm = .false. - real(kind=dp), allocatable :: cumulative_values(:) + logical, parameter :: do_naive_algorithm = .true. + real(kind=dp), allocatable :: cumulative_vals(:) real(kind=dp), allocatable :: factors_2d(:,:) integer, allocatable :: particle_group(:) real(kind=dp) :: x + integer, allocatable :: leftovers(:) + real(kind=dp) :: tot_leftovers call aero_state_sort(aero_state, aero_data, bin_grid) @@ -2003,9 +2003,11 @@ subroutine aero_state_bin_deaverage_comp(aero_state, bin_grid, aero_data, & end if end do else + ! Assign each species to a group do i_spec = 1,aero_data_n_spec(aero_data) species_group_numbers(i_spec) = i_spec end do + n_group = aero_data_n_spec(aero_data) end if do i_bin = 1,bin_grid_size(bin_grid) @@ -2027,9 +2029,9 @@ subroutine aero_state_bin_deaverage_comp(aero_state, bin_grid, aero_data, & end do do i_class = 1,size(aero_state%awa%weight, 2) - n_parts = integer_varray_n_entry( & + n_part = integer_varray_n_entry( & aero_state%aero_sorted%size_class%inverse(i_bin, i_class)) - if (n_parts > 0) then + if (n_part > 0) then allocate(group_volume_conc(n_group)) group_volume_conc = 0.0d0 do i_spec = 1,aero_data_n_spec(aero_data) @@ -2040,33 +2042,35 @@ subroutine aero_state_bin_deaverage_comp(aero_state, bin_grid, aero_data, & allocate(group_fractions(n_group)) group_fractions = 0.0d0 do i_group = 1,n_group - group_fractions(i_group) = n_parts & + group_fractions(i_group) = n_part & * group_volume_conc(i_group) / total_volume_conc end do - if (naive_algorithm) then - allocate(cumulative_values(n_group+1)) - cumulative_values = 0.0d0 + if (do_naive_algorithm) then + allocate(cumulative_vals(n_group+1)) + cumulative_vals = 0.0d0 do i_group = 2,n_group + 1 - cumulative_values(i_group) = cumulative_values(i_group-1) & + cumulative_vals(i_group) = cumulative_vals(i_group-1) & + group_fractions(i_group-1) end do - allocate(particle_group(n_parts)) - do i_part = 1,n_parts - x = pmc_random() * n_parts - particle_group(i_part) = find_1d(n_group+1, cumulative_values, x) + allocate(particle_group(n_part)) + do i_part = 1,n_part + x = pmc_random() * n_part + particle_group(i_part) = find_1d(n_group+1, & + cumulative_vals, x) end do allocate(factors_2d(n_group,aero_data_n_spec(aero_data))) factors_2d = 0.0d0 do i_group = 1,n_group do i_spec = 1,aero_data_n_spec(aero_data) if (species_group_numbers(i_spec) == i_group) then - factors_2d(i_group,i_spec) = species_volume_conc(i_spec) & + factors_2d(i_group,i_spec) = & + species_volume_conc(i_spec) & / group_volume_conc(i_group) end if end do end do - do i_entry = 1,n_parts + do i_entry = 1,n_part i_part = aero_state%aero_sorted%size_class%inverse( & i_bin, i_class)%entry(i_entry) particle_volume = aero_particle_volume( & @@ -2079,27 +2083,40 @@ subroutine aero_state_bin_deaverage_comp(aero_state, bin_grid, aero_data, & end do end do deallocate(particle_group) - deallocate(cumulative_values) + deallocate(cumulative_vals) deallocate(factors_2d) else - allocate(shuffle_particles(n_parts)) - do i_part = 1,n_parts + allocate(shuffle_particles(n_part)) + do i_part = 1,n_part shuffle_particles(i_part) = i_part end do - call pmc_rand_shuffle_array(shuffle_particles, n_parts) + call pmc_rand_shuffle_array(shuffle_particles, n_part) + + allocate(leftovers(n_group)) + allocate(cumulative_vals(n_group+1)) + cumulative_vals = 0.0d0 + leftovers = 0 + tot_leftovers = 0 + do i_group = 1,n_group + tot_leftovers = tot_leftovers + group_fractions(i_group) & + - floor(group_fractions(i_group)) + end do + do i_group = 1,n_group + cumulative_vals(i_group+1) = cumulative_vals(i_group) + & + (group_fractions(i_group) - floor( & + group_fractions(i_group))) / tot_leftovers + end do + do i_part = 1,nint(tot_leftovers) + x = pmc_random() + i_group = find_1d(n_group+1, cumulative_vals, x) + leftovers(i_group) = leftovers(i_group) + 1 + end do + start_val = 1 - edge_case = .false. do i_group = 1,n_group if (group_volume_conc(i_group) > 0.0d0) then - if (edge_case) then - n_part_spec = prob_round(group_fractions(i_group) + 1) - edge_case = .false. - else - n_part_spec = prob_round(group_fractions(i_group)) - if (n_part_spec < group_fractions(i_group)) then - edge_case = .true. - end if - end if + n_part_spec = floor(group_fractions(i_group)) & + + leftovers(i_group) factors = 0.0d0 do i_spec = 1,aero_data_n_spec(aero_data) @@ -2109,8 +2126,7 @@ subroutine aero_state_bin_deaverage_comp(aero_state, bin_grid, aero_data, & end if end do - do i_entry = start_val,min(start_val + n_part_spec - 1, & - n_parts) + do i_entry = start_val,(start_val + n_part_spec - 1) i_part = aero_state%aero_sorted%size_class%inverse( & i_bin, i_class)%entry(shuffle_particles(i_entry)) particle_volume = aero_particle_volume( & @@ -2125,6 +2141,8 @@ subroutine aero_state_bin_deaverage_comp(aero_state, bin_grid, aero_data, & end if end do deallocate(shuffle_particles) + deallocate(leftovers) + deallocate(cumulative_vals) end if deallocate(group_fractions) deallocate(group_volume_conc) From 529a580085dcb4b8d2dc2838e003709d8d48c043 Mon Sep 17 00:00:00 2001 From: Jeff Curtis Date: Sat, 30 Jul 2022 22:12:53 -0500 Subject: [PATCH 13/15] changed low variance algorithm --- src/aero_state.F90 | 178 +++++++++++++++++++++------------------------ 1 file changed, 84 insertions(+), 94 deletions(-) diff --git a/src/aero_state.F90 b/src/aero_state.F90 index 16f354cb9..7a597af28 100644 --- a/src/aero_state.F90 +++ b/src/aero_state.F90 @@ -1969,13 +1969,13 @@ subroutine aero_state_bin_deaverage_comp(aero_state, bin_grid, aero_data, & integer :: n_group, i_group, i_name real(kind=dp), allocatable :: group_fractions(:), & group_volume_conc(:) - logical, parameter :: do_naive_algorithm = .true. + logical, parameter :: do_naive_algorithm = .false. real(kind=dp), allocatable :: cumulative_vals(:) real(kind=dp), allocatable :: factors_2d(:,:) integer, allocatable :: particle_group(:) real(kind=dp) :: x - integer, allocatable :: leftovers(:) - real(kind=dp) :: tot_leftovers + + real(kind=dp) :: v_aerosol, v_bulk, p, v_group call aero_state_sort(aero_state, aero_data, bin_grid) @@ -2028,24 +2028,37 @@ subroutine aero_state_bin_deaverage_comp(aero_state, bin_grid, aero_data, & end do end do - do i_class = 1,size(aero_state%awa%weight, 2) - n_part = integer_varray_n_entry( & - aero_state%aero_sorted%size_class%inverse(i_bin, i_class)) - if (n_part > 0) then - allocate(group_volume_conc(n_group)) - group_volume_conc = 0.0d0 - do i_spec = 1,aero_data_n_spec(aero_data) - group_volume_conc(species_group_numbers(i_spec)) & - = group_volume_conc(species_group_numbers(i_spec)) & - + species_volume_conc(i_spec) - end do - allocate(group_fractions(n_group)) - group_fractions = 0.0d0 - do i_group = 1,n_group - group_fractions(i_group) = n_part & - * group_volume_conc(i_group) / total_volume_conc - end do - if (do_naive_algorithm) then + allocate(group_volume_conc(n_group)) + group_volume_conc = 0.0d0 + do i_spec = 1,aero_data_n_spec(aero_data) + group_volume_conc(species_group_numbers(i_spec)) & + = group_volume_conc(species_group_numbers(i_spec)) & + + species_volume_conc(i_spec) + end do + + allocate(factors_2d(n_group,aero_data_n_spec(aero_data))) + factors_2d = 0.0d0 + do i_group = 1,n_group + do i_spec = 1,aero_data_n_spec(aero_data) + if (species_group_numbers(i_spec) == i_group) then + factors_2d(i_group,i_spec) = & + species_volume_conc(i_spec) & + / group_volume_conc(i_group) + end if + end do + end do + + if (do_naive_algorithm) then + do i_class = 1,size(aero_state%awa%weight, 2) + n_part = integer_varray_n_entry( & + aero_state%aero_sorted%size_class%inverse(i_bin, i_class)) + if (n_part > 0) then + allocate(group_fractions(n_group)) + group_fractions = 0.0d0 + do i_group = 1,n_group + group_fractions(i_group) = n_part & + * group_volume_conc(i_group) / total_volume_conc + end do allocate(cumulative_vals(n_group+1)) cumulative_vals = 0.0d0 do i_group = 2,n_group + 1 @@ -2058,17 +2071,6 @@ subroutine aero_state_bin_deaverage_comp(aero_state, bin_grid, aero_data, & particle_group(i_part) = find_1d(n_group+1, & cumulative_vals, x) end do - allocate(factors_2d(n_group,aero_data_n_spec(aero_data))) - factors_2d = 0.0d0 - do i_group = 1,n_group - do i_spec = 1,aero_data_n_spec(aero_data) - if (species_group_numbers(i_spec) == i_group) then - factors_2d(i_group,i_spec) = & - species_volume_conc(i_spec) & - / group_volume_conc(i_group) - end if - end do - end do do i_entry = 1,n_part i_part = aero_state%aero_sorted%size_class%inverse( & @@ -2084,70 +2086,58 @@ subroutine aero_state_bin_deaverage_comp(aero_state, bin_grid, aero_data, & end do deallocate(particle_group) deallocate(cumulative_vals) - deallocate(factors_2d) - else - allocate(shuffle_particles(n_part)) - do i_part = 1,n_part - shuffle_particles(i_part) = i_part - end do - call pmc_rand_shuffle_array(shuffle_particles, n_part) - - allocate(leftovers(n_group)) - allocate(cumulative_vals(n_group+1)) - cumulative_vals = 0.0d0 - leftovers = 0 - tot_leftovers = 0 - do i_group = 1,n_group - tot_leftovers = tot_leftovers + group_fractions(i_group) & - - floor(group_fractions(i_group)) - end do - do i_group = 1,n_group - cumulative_vals(i_group+1) = cumulative_vals(i_group) + & - (group_fractions(i_group) - floor( & - group_fractions(i_group))) / tot_leftovers - end do - do i_part = 1,nint(tot_leftovers) - x = pmc_random() - i_group = find_1d(n_group+1, cumulative_vals, x) - leftovers(i_group) = leftovers(i_group) + 1 - end do - - start_val = 1 - do i_group = 1,n_group - if (group_volume_conc(i_group) > 0.0d0) then - n_part_spec = floor(group_fractions(i_group)) & - + leftovers(i_group) - - factors = 0.0d0 - do i_spec = 1,aero_data_n_spec(aero_data) - if (species_group_numbers(i_spec) == i_group) then - factors(i_spec) = species_volume_conc(i_spec) & - / group_volume_conc(i_group) - end if - end do - - do i_entry = start_val,(start_val + n_part_spec - 1) - i_part = aero_state%aero_sorted%size_class%inverse( & - i_bin, i_class)%entry(shuffle_particles(i_entry)) - particle_volume = aero_particle_volume( & - aero_state%apa%particle(i_part)) - aero_state%apa%particle(i_part)%vol = 0.0d0 - do i_spec = 1,aero_data_n_spec(aero_data) - aero_state%apa%particle(i_part)%vol(i_spec) & - = particle_volume * factors(i_spec) - end do - end do - start_val = start_val + n_part_spec - end if + deallocate(group_fractions) + end if + end do + else + n_part = aero_state_n_part(aero_state) + allocate(shuffle_particles(n_part)) + do i_part = 1,n_part + shuffle_particles(i_part) = i_part + end do + call pmc_rand_shuffle_array(shuffle_particles, n_part) + i_entry = 1 + v_aerosol = 0.0d0 + v_bulk = 0.0d0 + v_group = 0.0d0 + i_group = 0 + do while (i_entry <= n_part) + if (v_aerosol >= v_bulk) then + i_group = i_group + 1 + v_group = group_volume_conc(i_group) + v_bulk = v_bulk + v_group + end if + i_part = shuffle_particles(i_entry) + particle_volume = aero_particle_volume( & + aero_state%apa%particle(i_part)) + num_conc = aero_weight_array_num_conc(aero_state%awa, & + aero_state%apa%particle(i_part), aero_data) + if (particle_volume * num_conc + v_aerosol < v_bulk) then + aero_state%apa%particle(i_part)%vol = 0.0d0 + do i_spec = 1,aero_data_n_spec(aero_data) + aero_state%apa%particle(i_part)%vol(i_spec) = & + particle_volume * factors_2d(i_group,i_spec) end do - deallocate(shuffle_particles) - deallocate(leftovers) - deallocate(cumulative_vals) + v_aerosol = v_aerosol + particle_volume * num_conc + i_entry = i_entry + 1 + else + p = (v_bulk - max(v_aerosol, v_bulk - v_group)) & + / (particle_volume * num_conc) + if (pmc_random() < p) then + aero_state%apa%particle(i_part)%vol = 0.0d0 + do i_spec = 1,aero_data_n_spec(aero_data) + aero_state%apa%particle(i_part)%vol(i_spec) = & + particle_volume * factors_2d(i_group,i_spec) + end do + v_aerosol = v_aerosol + particle_volume * num_conc + i_entry = i_entry + 1 + end if end if - deallocate(group_fractions) - deallocate(group_volume_conc) - end if - end do + end do + deallocate(shuffle_particles) + end if + deallocate(group_volume_conc) + deallocate(factors_2d) end do end subroutine aero_state_bin_deaverage_comp From 594261febc0ad8b571c0ace4373988035b58bbc0 Mon Sep 17 00:00:00 2001 From: Jeff Curtis Date: Wed, 3 Aug 2022 11:28:16 -0500 Subject: [PATCH 14/15] fixed low variance method to work with bins --- src/aero_state.F90 | 102 +++++++++++++++++++++++++++------------------ 1 file changed, 62 insertions(+), 40 deletions(-) diff --git a/src/aero_state.F90 b/src/aero_state.F90 index 7a597af28..a07a1e146 100644 --- a/src/aero_state.F90 +++ b/src/aero_state.F90 @@ -1964,7 +1964,6 @@ subroutine aero_state_bin_deaverage_comp(aero_state, bin_grid, aero_data, & integer :: i_bin, i_class, i_entry, i_part, i_spec real(kind=dp) :: factors(aero_data_n_spec(aero_data)) integer :: n_part_spec, start_val, i, n_part - integer, allocatable :: shuffle_particles(:) integer :: species_group_numbers(aero_data_n_spec(aero_data)) integer :: n_group, i_group, i_name real(kind=dp), allocatable :: group_fractions(:), & @@ -1974,8 +1973,8 @@ subroutine aero_state_bin_deaverage_comp(aero_state, bin_grid, aero_data, & real(kind=dp), allocatable :: factors_2d(:,:) integer, allocatable :: particle_group(:) real(kind=dp) :: x - - real(kind=dp) :: v_aerosol, v_bulk, p, v_group + integer, allocatable :: particle_index(:) + real(kind=dp) :: v_aerosol, v_bulk, p, v_group, v_particle, v_group_prev call aero_state_sort(aero_state, aero_data, bin_grid) @@ -2090,51 +2089,74 @@ subroutine aero_state_bin_deaverage_comp(aero_state, bin_grid, aero_data, & end if end do else - n_part = aero_state_n_part(aero_state) - allocate(shuffle_particles(n_part)) - do i_part = 1,n_part - shuffle_particles(i_part) = i_part + n_part = 0 + do i_class = 1,size(aero_state%awa%weight, 2) + n_part = n_part + integer_varray_n_entry( & + aero_state%aero_sorted%size_class%inverse(i_bin, i_class)) end do - call pmc_rand_shuffle_array(shuffle_particles, n_part) - i_entry = 1 - v_aerosol = 0.0d0 - v_bulk = 0.0d0 - v_group = 0.0d0 - i_group = 0 - do while (i_entry <= n_part) - if (v_aerosol >= v_bulk) then - i_group = i_group + 1 - v_group = group_volume_conc(i_group) - v_bulk = v_bulk + v_group - end if - i_part = shuffle_particles(i_entry) - particle_volume = aero_particle_volume( & - aero_state%apa%particle(i_part)) - num_conc = aero_weight_array_num_conc(aero_state%awa, & - aero_state%apa%particle(i_part), aero_data) - if (particle_volume * num_conc + v_aerosol < v_bulk) then - aero_state%apa%particle(i_part)%vol = 0.0d0 - do i_spec = 1,aero_data_n_spec(aero_data) - aero_state%apa%particle(i_part)%vol(i_spec) = & - particle_volume * factors_2d(i_group,i_spec) + + if (n_part > 0) then + allocate(particle_index(n_part)) + i_part = 1 + do i_class = 1,size(aero_state%awa%weight, 2) + do i_entry = 1, integer_varray_n_entry( & + aero_state%aero_sorted%size_class%inverse(i_bin, i_class)) + particle_index(i_part) = & + aero_state%aero_sorted%size_class%inverse( & + i_bin, i_class)%entry(i_entry) + i_part = i_part + 1 end do - v_aerosol = v_aerosol + particle_volume * num_conc - i_entry = i_entry + 1 - else - p = (v_bulk - max(v_aerosol, v_bulk - v_group)) & - / (particle_volume * num_conc) - if (pmc_random() < p) then + end do + + call pmc_rand_shuffle_array(particle_index, n_part) + + i_entry = 1 + v_aerosol = 0.0d0 + v_bulk = 0.0d0 + v_particle = 0.0d0 + v_group = 0.0d0 + v_group_prev = 0.0d0 + i_group = 0 + do while (i_entry <= n_part) + if (v_aerosol + v_particle >= v_bulk) then + i_group = i_group + 1 + v_group = group_volume_conc(i_group) + v_group_prev = v_bulk + v_bulk = v_bulk + v_group + end if + i_part = particle_index(i_entry) + particle_volume = aero_particle_volume( & + aero_state%apa%particle(i_part)) + num_conc = aero_weight_array_num_conc(aero_state%awa, & + aero_state%apa%particle(i_part), aero_data) + v_particle = particle_volume * num_conc + if (min(v_particle + v_aerosol, sum(group_volume_conc)) & + <= v_bulk) then aero_state%apa%particle(i_part)%vol = 0.0d0 do i_spec = 1,aero_data_n_spec(aero_data) aero_state%apa%particle(i_part)%vol(i_spec) = & - particle_volume * factors_2d(i_group,i_spec) + particle_volume * factors_2d(i_group,i_spec) end do - v_aerosol = v_aerosol + particle_volume * num_conc + v_aerosol = v_aerosol + v_particle i_entry = i_entry + 1 + v_particle = 0.0d0 + else + p = (v_bulk - max(v_aerosol, v_group_prev)) & + / (v_particle) + if (pmc_random() < p) then + aero_state%apa%particle(i_part)%vol = 0.0d0 + do i_spec = 1,aero_data_n_spec(aero_data) + aero_state%apa%particle(i_part)%vol(i_spec) = & + particle_volume * factors_2d(i_group,i_spec) + end do + v_aerosol = v_aerosol + v_particle + i_entry = i_entry + 1 + v_particle = 0.0d0 + end if end if - end if - end do - deallocate(shuffle_particles) + end do + deallocate(particle_index) + end if end if deallocate(group_volume_conc) deallocate(factors_2d) From e463d99d0b305d9a30460f7d25d6574b7417b077 Mon Sep 17 00:00:00 2001 From: Jeffrey Curtis Date: Mon, 29 Apr 2024 10:33:09 -0500 Subject: [PATCH 15/15] add comment about algorithm selection --- src/aero_state.F90 | 1 + 1 file changed, 1 insertion(+) diff --git a/src/aero_state.F90 b/src/aero_state.F90 index a07a1e146..165343369 100644 --- a/src/aero_state.F90 +++ b/src/aero_state.F90 @@ -1968,6 +1968,7 @@ subroutine aero_state_bin_deaverage_comp(aero_state, bin_grid, aero_data, & integer :: n_group, i_group, i_name real(kind=dp), allocatable :: group_fractions(:), & group_volume_conc(:) + ! Flag to control algorithm selection - false uses low-variance method. logical, parameter :: do_naive_algorithm = .false. real(kind=dp), allocatable :: cumulative_vals(:) real(kind=dp), allocatable :: factors_2d(:,:)