Skip to content
New issue

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

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

Already on GitHub? Sign in to your account

Simpol multiple layers #28

Draft
wants to merge 27 commits into
base: main
Choose a base branch
from
Draft
Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
27 commits
Select commit Hold shift + click to select a range
c3f6ad9
addition of is_at_surface flag to num_phase_instances function
jtgasparik Aug 25, 2024
f04db9a
testing partmc coupling
jtgasparik Aug 25, 2024
85e781f
testing
jtgasparik Aug 25, 2024
122c1e3
fixed num_phase_instances function, now passes
jtgasparik Aug 26, 2024
11f8275
fixed function and added test
jtgasparik Aug 27, 2024
98ad87d
one more test
jtgasparik Aug 27, 2024
3e3c6d2
testing
jtgasparik Aug 27, 2024
a0a319b
reverting change
jtgasparik Aug 27, 2024
361f877
changes to phase_id function
jtgasparik Sep 3, 2024
992ea23
deleted duplicate num_instances variable
jtgasparik Sep 10, 2024
c7d5074
testing SIMPOL reaction with multiple layers added to the json - curr…
jtgasparik Sep 23, 2024
9a0c89b
update unique names function to have is_at_surface flag
mattldawson Oct 12, 2024
eb8ace3
remote diagnostic output
mattldawson Oct 12, 2024
b2e9b2a
Merge pull request #8 from open-atmos/develop-debug-SIMPOL
jtgasparik Oct 15, 2024
f91e6ac
added test for unique_name function
jtgasparik Oct 22, 2024
c2b608f
ive broken the SIMPOL rxn again by adding layers
jtgasparik Oct 23, 2024
178950d
adding mass
jtgasparik Oct 28, 2024
785e5d2
minor changes to SIMPOl test
jtgasparik Nov 10, 2024
26fbaea
fix to write statement - still not passing tests
jtgasparik Nov 10, 2024
0627631
added tests for phase_ids function to single particle
jtgasparik Dec 9, 2024
dee7c7f
test added for phase id for modal/binned rep
jtgasparik Dec 10, 2024
5c49679
added second particle to test suite for single particle c test - in p…
jtgasparik Dec 26, 2024
9be64bc
adding unit tests to single particle c functions with multiple partic…
jtgasparik Dec 27, 2024
e74abdb
adding second particle to single particle c tests -- in progress and…
jtgasparik Dec 27, 2024
fa093b1
adding another test particle and phase to single particle test aero r…
jtgasparik Dec 31, 2024
0b39752
testing the same phase in a different particle works. testing a diffe…
jtgasparik Jan 1, 2025
de1bf73
testing differnt particles and different layers
jtgasparik Jan 6, 2025
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
23 changes: 20 additions & 3 deletions src/aero_rep_data.F90
Original file line number Diff line number Diff line change
Expand Up @@ -276,7 +276,8 @@ end function get_size
!> Get a list of unique names for each element on the
!! \c camp_camp_state::camp_state_t::state_var array for this aerosol
!! representation.
function unique_names(this, phase_name, tracer_type, spec_name)
function unique_names(this, phase_name, tracer_type, spec_name, &
phase_is_at_surface)
use camp_util, only : string_t, i_kind
import :: aero_rep_data_t

Expand All @@ -290,6 +291,8 @@ function unique_names(this, phase_name, tracer_type, spec_name)
integer(kind=i_kind), optional, intent(in) :: tracer_type
!> Aerosol-phase species name
character(len=*), optional, intent(in) :: spec_name
!> Indicates if aerosol phase is at the surface of particle
logical, optional, intent(in) :: phase_is_at_surface

end function unique_names

Expand Down Expand Up @@ -342,7 +345,7 @@ end function spec_name
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!> Get the number of instances of a specified aerosol phase
function num_phase_instances(this, phase_name)
function num_phase_instances(this, phase_name, is_at_surface)
use camp_util, only : i_kind
import :: aero_rep_data_t

Expand All @@ -352,6 +355,8 @@ function num_phase_instances(this, phase_name)
class(aero_rep_data_t), intent(in) :: this
!> Aerosol phase name
character(len=*), intent(in) :: phase_name
!> Indicates if aerosol phase is at the surface of particle
logical, intent(in), optional :: is_at_surface

end function num_phase_instances

Expand Down Expand Up @@ -506,7 +511,18 @@ function phase_ids(this, phase_name, is_at_surface)

integer(kind=i_kind) :: num_instances, i_instance, i_phase

num_instances = this%num_phase_instances(phase_name)
if (present(is_at_surface)) then
if (is_at_surface) then
num_instances = this%num_phase_instances(phase_name, &
is_at_surface = .true.)
else
num_instances = this%num_phase_instances(phase_name, &
is_at_surface = .false.)
end if
else
num_instances = this%num_phase_instances(phase_name)
end if

allocate(phase_ids(num_instances))
if (present(is_at_surface)) then
if (is_at_surface) then
Expand Down Expand Up @@ -537,6 +553,7 @@ function phase_ids(this, phase_name, is_at_surface)
end if
end do
end if
call assert(642387392, num_instances == i_instance-1)

end function phase_ids

Expand Down
52 changes: 46 additions & 6 deletions src/aero_reps/aero_rep_modal_binned_mass.F90
Original file line number Diff line number Diff line change
Expand Up @@ -734,7 +734,8 @@ end function get_size
!!
!! ... and for modes are:
!! - "mode name.phase name.species name"
function unique_names(this, phase_name, tracer_type, spec_name)
function unique_names(this, phase_name, tracer_type, spec_name, &
phase_is_at_surface)

!> List of unique names
type(string_t), allocatable :: unique_names(:)
Expand All @@ -746,6 +747,8 @@ function unique_names(this, phase_name, tracer_type, spec_name)
integer(kind=i_kind), optional, intent(in) :: tracer_type
!> Aerosol-phase species name
character(len=*), optional, intent(in) :: spec_name
!> Aerosol-phase species is at surface
logical, optional, intent(in) :: phase_is_at_surface

integer(kind=i_kind) :: num_spec, i_spec, j_spec, i_phase, j_phase, &
i_section, i_bin
Expand All @@ -764,6 +767,12 @@ function unique_names(this, phase_name, tracer_type, spec_name)
if (phase_name.ne.curr_phase_name) cycle
end if

! Filter by phase is at surface
if (present(phase_is_at_surface)) then
if (phase_is_at_surface .neqv. &
this%aero_phase_is_at_surface(i_phase)) cycle
end if

! Filter by spec name and/or tracer type
if (present(spec_name).or.present(tracer_type)) then
spec_names = this%aero_phase(i_phase)%val%get_species_names()
Expand Down Expand Up @@ -810,6 +819,15 @@ function unique_names(this, phase_name, tracer_type, spec_name)
end if
end if

! Filter by phase is at surface
if (present(phase_is_at_surface)) then
if (phase_is_at_surface .neqv. &
this%aero_phase_is_at_surface(i_phase)) then
i_phase = i_phase + NUM_BINS_(i_section)
cycle
end if
end if

! Get the species names in this phase
spec_names = this%aero_phase(i_phase)%val%get_species_names()

Expand Down Expand Up @@ -941,23 +959,45 @@ end function spec_name
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!> Get the number of instances of a specified aerosol phase.
function num_phase_instances(this, phase_name)
function num_phase_instances(this, phase_name, is_at_surface)

!> Number of instances of the aerosol phase
integer(kind=i_kind) :: num_phase_instances
!> Aerosol representation data
class(aero_rep_modal_binned_mass_t), intent(in) :: this
!> Aerosol phase name
character(len=*), intent(in) :: phase_name
!> Indicates if aerosol phase is at the surface of particle
logical, intent(in), optional :: is_at_surface

integer(kind=i_kind) :: i_phase

num_phase_instances = 0
do i_phase = 1, size(this%aero_phase)
if (this%aero_phase(i_phase)%val%name().eq.phase_name) then
num_phase_instances = num_phase_instances + 1
if (present(is_at_surface)) then
if (is_at_surface) then
do i_phase = 1, size(this%aero_phase)
if (this%aero_phase(i_phase)%val%name().eq.phase_name .and. &
this%aero_phase_is_at_surface(i_phase)) then
num_phase_instances = num_phase_instances + 1
end if
end do
else
do i_phase = 1, size(this%aero_phase)
if (this%aero_phase(i_phase)%val%name().eq.phase_name .and. &
.not. this%aero_phase_is_at_surface(i_phase)) then
call die_msg(507144607, "Species must exist at surface "// &
"in modal/binned mass aerosol representation '"// &
this%rep_name//"'")
end if
end do
end if
end do
else
do i_phase = 1, size(this%aero_phase)
if (this%aero_phase(i_phase)%val%name().eq.phase_name) then
num_phase_instances = num_phase_instances + 1
end if
end do
end if

end function num_phase_instances

Expand Down
63 changes: 50 additions & 13 deletions src/aero_reps/aero_rep_single_particle.F90
Original file line number Diff line number Diff line change
Expand Up @@ -375,7 +375,7 @@ subroutine initialize(this, aero_phase_set, spec_state_id)
curr_phase = 1
do i_layer = 1, size(ordered_layer_id)
j_layer = ordered_layer_id(i_layer)

print *, "layer number ", i_layer
! Loop through the phases and make sure they exist
call phases(j_layer)%val_%iter_reset()
do i_phase = 1, phases(j_layer)%val_%size()
Expand All @@ -387,11 +387,11 @@ subroutine initialize(this, aero_phase_set, spec_state_id)
layer_names_unordered(j_layer)%string// &
"' in single-particle layer aerosol representation '"// &
this%rep_name//"'")

! find phase and set pointer and indices
found = .false.
do j_phase = 1, size(aero_phase_set)
if (aero_phase_set(j_phase)%val%name() .eq. phase_name) then
print *, "phase name : ", phase_name
found = .true.
do i_particle = 0, num_particles-1
this%aero_phase(i_particle*num_phases + curr_phase) = &
Expand All @@ -404,6 +404,7 @@ subroutine initialize(this, aero_phase_set, spec_state_id)
.false.
end if
end do
!print *, "aero_phase_is_at_surface", this%aero_phase_is_at_surface
PHASE_STATE_ID_(i_layer,i_phase) = curr_id
PHASE_MODEL_DATA_ID_(i_layer,i_phase) = j_phase
curr_id = curr_id + aero_phase_set(j_phase)%val%size()
Expand All @@ -415,6 +416,7 @@ subroutine initialize(this, aero_phase_set, spec_state_id)
call phases(j_layer)%val_%iter_next()
end do
end do
print *, "aero_phase_is_at_surface", this%aero_phase_is_at_surface
PARTICLE_STATE_SIZE_ = curr_id - spec_state_id

! Initialize the aerosol representation id
Expand Down Expand Up @@ -485,7 +487,8 @@ end function per_particle_size
!! For a single particle representation, the unique names will be a 'P'
!! followed by the computational particle number, a '.', the phase name,
!! another '.', and the species name.
function unique_names(this, phase_name, tracer_type, spec_name)
function unique_names(this, phase_name, tracer_type, spec_name, &
phase_is_at_surface)

use camp_util, only : integer_to_string
!> List of unique names
Expand All @@ -498,6 +501,8 @@ function unique_names(this, phase_name, tracer_type, spec_name)
integer(kind=i_kind), optional, intent(in) :: tracer_type
!> Aerosol-phase species name
character(len=*), optional, intent(in) :: spec_name
!> Indicates if aerosol phase is at the surface of particle
logical, optional, intent(in) :: phase_is_at_surface

integer :: i_particle, i_layer, i_phase, i_spec, j_spec
integer :: num_spec, curr_tracer_type
Expand All @@ -518,6 +523,10 @@ function unique_names(this, phase_name, tracer_type, spec_name)
if (present(phase_name)) then
if(phase_name .ne. this%aero_phase(i_phase)%val%name()) cycle
end if
if (present(phase_is_at_surface)) then
if (phase_is_at_surface .neqv. &
this%aero_phase_is_at_surface(i_phase)) cycle
end if
if (present(spec_name) .or. present(tracer_type)) then
spec_names = this%aero_phase(i_phase)%val%get_species_names()
do j_spec = 1, size(spec_names)
Expand Down Expand Up @@ -547,6 +556,10 @@ function unique_names(this, phase_name, tracer_type, spec_name)
if (present(phase_name)) then
if(phase_name .ne. this%aero_phase(i_phase)%val%name()) cycle
end if
if (present(phase_is_at_surface)) then
if (phase_is_at_surface .neqv. &
this%aero_phase_is_at_surface(i_phase)) cycle
end if
spec_names = this%aero_phase(i_phase)%val%get_species_names()
do j_spec = 1, this%aero_phase(i_phase)%val%size()
if (present(spec_name)) then
Expand Down Expand Up @@ -638,25 +651,49 @@ end function spec_name
!> Get the number of instances of a specified aerosol phase. In the single
!! particle representation with layers, a phase can exist in multiple layers
!! in one particle.
integer(kind=i_kind) function num_phase_instances(this, phase_name)
integer(kind=i_kind) function num_phase_instances(this, phase_name, &
is_at_surface)

!> Aerosol representation data
class(aero_rep_single_particle_t), intent(in) :: this
!> Aerosol phase name
character(len=*), intent(in) :: phase_name
!> Indicates if aerosol phase is at the surface of particle
logical, intent(in), optional :: is_at_surface

integer(kind=i_kind) :: i_phase, i_layer, phase_index

num_phase_instances = 0
phase_index = 0
do i_layer = 1, NUM_LAYERS_
do i_phase = LAYER_PHASE_START_(i_layer), LAYER_PHASE_END_(i_layer)
if (this%aero_phase(i_phase)%val%name() .eq. phase_name) then
phase_index = phase_index + 1
end if
if (present(is_at_surface)) then
if (is_at_surface) then
do i_layer = 1, NUM_LAYERS_
do i_phase = LAYER_PHASE_START_(i_layer), LAYER_PHASE_END_(i_layer)
if (this%aero_phase(i_phase)%val%name() .eq. phase_name .and. &
this%aero_phase_is_at_surface(i_phase)) then
num_phase_instances = num_phase_instances + 1
end if
end do
end do
else
do i_layer = 1, NUM_LAYERS_
do i_phase = LAYER_PHASE_START_(i_layer), LAYER_PHASE_END_(i_layer)
if (this%aero_phase(i_phase)%val%name() .eq. phase_name .and. &
.not. this%aero_phase_is_at_surface(i_phase)) then
num_phase_instances = num_phase_instances + 1
end if
end do
end do
end if
else
do i_layer = 1, NUM_LAYERS_
do i_phase = LAYER_PHASE_START_(i_layer), LAYER_PHASE_END_(i_layer)
if (this%aero_phase(i_phase)%val%name() .eq. phase_name) then
num_phase_instances = num_phase_instances + 1
end if
end do
end do
end do
num_phase_instances = phase_index * MAX_PARTICLES_
end if
num_phase_instances = num_phase_instances * MAX_PARTICLES_

end function num_phase_instances

Expand All @@ -675,7 +712,7 @@ function num_jac_elem(this, phase_id)

integer(kind=i_kind) :: i_phase

call assert_msg(927040495, phase_id .ge. 1 .and. &
call assert_msg(927040495, phase_id .ge. 0 .and. &
phase_id .le. size( this%aero_phase ), &
"Aerosol phase index out of range. Got "// &
trim( integer_to_string( phase_id ) )//", expected 1:"// &
Expand Down
19 changes: 15 additions & 4 deletions src/rxns/rxn_SIMPOL_phase_transfer.F90
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
! Copyright (C) 2021 Barcelona Supercomputing Center and University of
! Copyright (C) 2021 Barcelona :/repSupercomputing Center and University of
! Illinois at Urbana-Champaign
! SPDX-License-Identifier: MIT

Expand Down Expand Up @@ -221,6 +221,7 @@ subroutine initialize(this, chem_spec_data, aero_rep, n_cells)
this%property_set%get_string(key_name, aero_spec_name), &
"Missing aerosol-phase species in SIMPOL.1 phase-transfer "// &
"reaction")
print *, "species name ", aero_spec_name

! Set up a general error message
error_msg = " for SIMPOL.1 phase transfer of gas species '"// &
Expand All @@ -245,7 +246,8 @@ subroutine initialize(this, chem_spec_data, aero_rep, n_cells)
! Get the unique names in this aerosol representation for the
! partitioning species
unique_spec_names = aero_rep(i_aero_rep)%val%unique_names( &
phase_name = phase_name, spec_name = aero_spec_name)
phase_name = phase_name, spec_name = aero_spec_name, &
phase_is_at_surface = .true.)

! Skip aerosol representations that do not contain this phase
if (.not.allocated(unique_spec_names)) cycle
Expand All @@ -256,6 +258,7 @@ subroutine initialize(this, chem_spec_data, aero_rep, n_cells)
! Get the number of Jacobian elements for calculations of mass, volume,
! number, etc. for this partitioning into this phase
phase_ids = aero_rep(i_aero_rep)%val%phase_ids(phase_name, is_at_surface=.true.)
print *, "phase_ids", phase_ids
do i_phase = 1, size(phase_ids)
n_aero_jac_elem = n_aero_jac_elem + &
aero_rep(i_aero_rep)%val%num_jac_elem(phase_ids(i_phase))
Expand Down Expand Up @@ -306,12 +309,14 @@ subroutine initialize(this, chem_spec_data, aero_rep, n_cells)
! Get the unique names in this aerosol representation for the
! partitioning species
unique_spec_names = aero_rep(i_aero_rep)%val%unique_names( &
phase_name = phase_name, spec_name = aero_spec_name)
phase_name = phase_name, spec_name = aero_spec_name, &
phase_is_at_surface = .true.)

! Find the corresponding activity coefficients, if specified
if (has_act_coeff) then
unique_act_names = aero_rep(i_aero_rep)%val%unique_names( &
phase_name = phase_name, spec_name = act_name)
phase_name = phase_name, spec_name = act_name, &
phase_is_at_surface = .true.)
call assert_msg(236251734, size(unique_act_names).eq. &
size(unique_spec_names), &
"Mismatch of SIMPOL species and activity coeffs"// &
Expand Down Expand Up @@ -420,12 +425,18 @@ subroutine initialize(this, chem_spec_data, aero_rep, n_cells)
! Check the sizes of the data arrays
tmp_size = PHASE_INT_LOC_(i_aero_id - 1) + 1 + &
2*NUM_AERO_PHASE_JAC_ELEM_(i_aero_id - 1) - 1
print *, "PHASE_INT_LOC_", PHASE_INT_LOC_(i_aero_id - 1)
print *, "NUM_AERO_PHASE_JAC_ELEM_", NUM_AERO_PHASE_JAC_ELEM_(i_aero_id - 1)
print *, "tmp_size_1", tmp_size
print *, "size condensed data int",size(this%condensed_data_int)
call assert_msg(625802519, size(this%condensed_data_int) .eq. tmp_size, &
"int array size mismatch"//error_msg)
tmp_size = PHASE_REAL_LOC_(i_aero_id - 1) + &
4*NUM_AERO_PHASE_JAC_ELEM_(i_aero_id - 1) - 1
call assert_msg(391089510, size(this%condensed_data_real) .eq. tmp_size, &
"real array size mismatch"//error_msg)
print *, "tmp_size_2", tmp_size
print *, "size condensed data real",size(this%condensed_data_real)

end subroutine initialize

Expand Down
Loading
Loading