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

Update MAPL3 time averaging to use a pointwise counter #3189

Merged
Merged
Show file tree
Hide file tree
Changes from 13 commits
Commits
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
75 changes: 71 additions & 4 deletions field/FieldPointerUtilities.F90
Original file line number Diff line number Diff line change
Expand Up @@ -33,6 +33,8 @@ module MAPL_FieldPointerUtilities
module procedure assign_fptr_r8_rank2
module procedure assign_fptr_r4_rank3
module procedure assign_fptr_r8_rank3
module procedure assign_fptr_i4_rank1
module procedure assign_fptr_i8_rank1
end interface assign_fptr

interface FieldGetCptr
Expand Down Expand Up @@ -846,8 +848,11 @@ subroutine MAPL_FieldGetLocalElementCount(field,local_count,rc)

real(kind=ESMF_KIND_R4), pointer :: r4_1d(:),r4_2d(:,:),r4_3d(:,:,:),r4_4d(:,:,:,:)
real(kind=ESMF_KIND_R8), pointer :: r8_1d(:),r8_2d(:,:),r8_3d(:,:,:),r8_4d(:,:,:,:)
integer(kind=ESMF_KIND_I4), pointer :: i4_1d(:),i4_2d(:,:),i4_3d(:,:,:),i4_4d(:,:,:,:)
integer(kind=ESMF_KIND_I8), pointer :: i8_1d(:),i8_2d(:,:),i8_3d(:,:,:),i8_4d(:,:,:,:)

call ESMF_FieldGet(field,rank=rank,typekind=tk,_RC)
_ASSERT(rank > 0 .and. rank < 5, "Unsupported rank")
if (tk == ESMF_TypeKind_R4) then
if (rank==1) then
call ESMF_FieldGet(field,0,farrayptr=r4_1d,_RC)
Expand All @@ -861,8 +866,6 @@ subroutine MAPL_FieldGetLocalElementCount(field,local_count,rc)
else if (rank ==4) then
call ESMF_FieldGet(field,0,farrayptr=r4_4d,_RC)
local_count = shape(r4_4d)
else
_FAIL("Unsupported rank")
end if
else if (tk == ESMF_TypeKind_R8) then
if (rank==1) then
Expand All @@ -877,8 +880,34 @@ subroutine MAPL_FieldGetLocalElementCount(field,local_count,rc)
else if (rank ==4) then
call ESMF_FieldGet(field,0,farrayptr=r8_4d,_RC)
local_count = shape(r8_4d)
else
_FAIL("Unsupported rank")
end if
else if (tk == ESMF_TypeKind_I4) then
if (rank==1) then
call ESMF_FieldGet(field,0,farrayptr=i4_1d,_RC)
local_count = shape(i4_1d)
else if (rank ==2) then
call ESMF_FieldGet(field,0,farrayptr=i4_2d,_RC)
local_count = shape(i4_2d)
else if (rank ==3) then
call ESMF_FieldGet(field,0,farrayptr=i4_3d,_RC)
local_count = shape(i4_3d)
else if (rank ==4) then
call ESMF_FieldGet(field,0,farrayptr=i4_4d,_RC)
local_count = shape(i4_4d)
end if
darianboggs marked this conversation as resolved.
Show resolved Hide resolved
else if (tk == ESMF_TypeKind_I8) then
if (rank==1) then
call ESMF_FieldGet(field,0,farrayptr=i8_1d,_RC)
local_count = shape(i8_1d)
else if (rank ==2) then
call ESMF_FieldGet(field,0,farrayptr=i8_2d,_RC)
local_count = shape(i8_2d)
else if (rank ==3) then
call ESMF_FieldGet(field,0,farrayptr=i8_3d,_RC)
local_count = shape(i8_3d)
else if (rank ==4) then
call ESMF_FieldGet(field,0,farrayptr=i8_4d,_RC)
local_count = shape(i8_4d)
end if
else
_FAIL("Unsupported type")
Expand Down Expand Up @@ -990,4 +1019,42 @@ subroutine Destroy(Field,RC)

end subroutine Destroy

subroutine assign_fptr_i4_rank1(x, fptr, rc)
type(ESMF_Field), intent(inout) :: x
integer(kind=ESMF_KIND_I4), pointer, intent(out) :: fptr(:)
integer, optional, intent(out) :: rc

! local declarations
type(c_ptr) :: cptr
integer(ESMF_KIND_I8), allocatable :: fp_shape(:)
integer(ESMF_KIND_I8) :: local_size
integer :: status

local_size = FieldGetLocalSize(x, _RC)
fp_shape = [ local_size ]
call FieldGetCptr(x, cptr, _RC)
call c_f_pointer(cptr, fptr, fp_shape)

_RETURN(_SUCCESS)
end subroutine assign_fptr_i4_rank1

subroutine assign_fptr_i8_rank1(x, fptr, rc)
type(ESMF_Field), intent(inout) :: x
integer(kind=ESMF_KIND_I8), pointer, intent(out) :: fptr(:)
integer, optional, intent(out) :: rc

! local declarations
type(c_ptr) :: cptr
integer(ESMF_KIND_I8), allocatable :: fp_shape(:)
integer(ESMF_KIND_I8) :: local_size
integer :: status

local_size = FieldGetLocalSize(x, _RC)
fp_shape = [ local_size ]
call FieldGetCptr(x, cptr, _RC)
call c_f_pointer(cptr, fptr, fp_shape)

_RETURN(_SUCCESS)
end subroutine assign_fptr_i8_rank1

end module MAPL_FieldPointerUtilities
86 changes: 60 additions & 26 deletions generic3g/actions/AccumulatorAction.F90
Original file line number Diff line number Diff line change
Expand Up @@ -15,16 +15,19 @@ module mapl3g_AccumulatorAction
type(ESMF_Field) :: result_field
real(kind=ESMF_KIND_R4) :: CLEAR_VALUE_R4 = 0.0_ESMF_KIND_R4
logical :: update_calculated = .FALSE.
type(ESMF_TypeKind_Flag) :: typekind = ESMF_TYPEKIND_R4
contains
! Implementations of deferred procedures
procedure :: invalidate
procedure :: initialize
procedure :: update
! Helpers
procedure :: accumulate
procedure :: initialized
procedure :: clear_accumulator
procedure :: accumulate
procedure :: accumulate_R4
procedure :: clear
procedure :: create_fields
procedure :: update_result
end type AccumulatorAction

contains
Expand All @@ -36,22 +39,20 @@ logical function initialized(this) result(lval)

end function initialized

subroutine clear_accumulator(this, rc)
subroutine clear(this, rc)
class(AccumulatorAction), intent(inout) :: this
integer, optional, intent(out) :: rc

integer :: status
type(ESMF_TypeKind_Flag) :: tk

call ESMF_FieldGet(this%accumulation_field, typekind=tk, _RC)
if(tk == ESMF_TYPEKIND_R4) then
if(this%typekind == ESMF_TYPEKIND_R4) then
call FieldSet(this%accumulation_field, this%CLEAR_VALUE_R4, _RC)
else
_FAIL('Unsupported typekind')
end if
_RETURN(_SUCCESS)

end subroutine clear_accumulator
end subroutine clear

subroutine initialize(this, importState, exportState, clock, rc)
class(AccumulatorAction), intent(inout) :: this
Expand All @@ -62,23 +63,45 @@ subroutine initialize(this, importState, exportState, clock, rc)

integer :: status
type(ESMF_Field) :: import_field, export_field
logical :: fields_are_conformable
type(ESMF_TypeKind_Flag) :: typekind
logical :: conformable
logical :: same_typekind

conformable = .FALSE.
same_typekind = .FALSE.
call get_field(importState, import_field, _RC)
call ESMF_FieldGet(import_field, typekind=typekind, _RC)
_ASSERT(typekind==ESMF_TYPEKIND_R4, 'Only ESMF_TYPEKIND_R4 is supported.')
call get_field(exportState, export_field, _RC)
conformable = FieldsAreConformable(import_field, export_field, _RC)
_ASSERT(conformable, 'Import and export fields are not conformable.')
same_typekind = FieldsAreSameTypeKind(import_field, export_field, _RC)
_ASSERT(same_typekind, 'Import and export fields are different typekinds.')
this%typekind = typekind
darianboggs marked this conversation as resolved.
Show resolved Hide resolved
call this%create_fields(import_field, export_field, _RC)
call this%clear(_RC)
_RETURN(_SUCCESS)
_UNUSED_DUMMY(clock)

end subroutine initialize

subroutine create_fields(this, import_field, export_field, rc)
class(AccumulatorAction), intent(inout) :: this
type(ESMF_Field), intent(inout) :: import_field
type(ESMF_Field), intent(inout) :: export_field
integer, optional, intent(out) :: rc

integer :: status

if(this%initialized()) then
call ESMF_FieldDestroy(this%accumulation_field, _RC)
call ESMF_FieldDestroy(this%result_field, _RC)
end if
this%accumulation_field = ESMF_FieldCreate(import_field, _RC)
this%result_field = ESMF_FieldCreate(export_field, _RC)

call this%clear_accumulator(_RC)
_RETURN(_SUCCESS)
_UNUSED_DUMMY(clock)

end subroutine initialize
end subroutine create_fields

subroutine update(this, importState, exportState, clock, rc)
class(AccumulatorAction), intent(inout) :: this
Expand All @@ -92,19 +115,30 @@ subroutine update(this, importState, exportState, clock, rc)

_ASSERT(this%initialized(), 'Accumulator has not been initialized.')
if(.not. this%update_calculated) then
call FieldCopy(this%accumulation_field, this%result_field, _RC)
this%update_calculated = .TRUE.
call this%update_result(_RC)
end if
call get_field(exportState, export_field, _RC)
call FieldCopy(this%result_field, export_field, _RC)

call this%clear_accumulator(_RC)
call this%clear(_RC)
_RETURN(_SUCCESS)
_UNUSED_DUMMY(clock)
_UNUSED_DUMMY(importState)
_RETURN(_SUCCESS)

end subroutine update

subroutine update_result(this, rc)
class(AccumulatorAction), intent(inout) :: this
integer, optional, intent(out) :: rc

integer :: status

call FieldCopy(this%accumulation_field, this%result_field, _RC)
this%update_calculated = .TRUE.
tclune marked this conversation as resolved.
Show resolved Hide resolved
_RETURN(_SUCCESS)

end subroutine update_result

subroutine invalidate(this, importState, exportState, clock, rc)
class(AccumulatorAction), intent(inout) :: this
type(ESMF_State) :: importState
Expand All @@ -119,9 +153,9 @@ subroutine invalidate(this, importState, exportState, clock, rc)
this%update_calculated = .FALSE.
call get_field(importState, import_field, _RC)
call this%accumulate(import_field, _RC)
_RETURN(_SUCCESS)
_UNUSED_DUMMY(clock)
_UNUSED_DUMMY(exportState)
_RETURN(_SUCCESS)

end subroutine invalidate

Expand Down Expand Up @@ -151,12 +185,11 @@ subroutine accumulate(this, update_field, rc)
integer, optional, intent(out) :: rc

integer :: status
type(ESMF_TypeKind_Flag) :: tk, tk_field
type(ESMF_TypeKind_Flag) :: tk_field

call ESMF_FieldGet(this%accumulation_field, typekind=tk, _RC)
call ESMF_FieldGet(update_field, typekind=tk_field, _RC)
_ASSERT(tk == tk_field, 'Update field must be the same typekind as the accumulation field.')
if(tk == ESMF_TYPEKIND_R4) then
_ASSERT(this%typekind == tk_field, 'Update field must be the same typekind as the accumulation field.')
if(this%typekind == ESMF_TYPEKIND_R4) then
call this%accumulate_R4(update_field, _RC)
else
_FAIL('Unsupported typekind value')
Expand All @@ -174,15 +207,16 @@ subroutine accumulate_R4(this, update_field, rc)
integer :: status
real(kind=ESMF_KIND_R4), pointer :: current(:)
real(kind=ESMF_KIND_R4), pointer :: latest(:)
real(kind=ESMF_KIND_R4) :: undef
real(kind=ESMF_KIND_R4), parameter :: UNDEF = MAPL_UNDEFINED_REAL

undef = MAPL_UNDEFINED_REAL
current => null()
latest => null()
call assign_fptr(this%accumulation_field, current, _RC)
call assign_fptr(update_field, latest, _RC)
where(current /= undef .and. latest /= undef)
where(current /= UNDEF .and. latest /= UNDEF)
current = current + latest
elsewhere(latest == undef)
current = undef
elsewhere(latest == UNDEF)
current = UNDEF
end where
_RETURN(_SUCCESS)

Expand Down
Loading