diff --git a/src/atlas/field/FieldSet.cc b/src/atlas/field/FieldSet.cc index 0627ed7bb..01a0db916 100644 --- a/src/atlas/field/FieldSet.cc +++ b/src/atlas/field/FieldSet.cc @@ -11,6 +11,7 @@ #include #include "atlas/field/Field.h" +#include "atlas/field/detail/FieldInterface.h" #include "atlas/field/FieldSet.h" #include "atlas/grid/Grid.h" #include "atlas/runtime/Exception.h" @@ -86,6 +87,39 @@ std::vector FieldSetImpl::field_names() const { // C wrapper interfaces to C++ routines extern "C" { +void atlas__FieldSet__data_int_specf(FieldSetImpl* This, char* name, int*& data, int& rank, int*& shapef, int*& stridesf) { + atlas__Field__data_int_specf(This->field(name).get(), data, rank, shapef, stridesf); +} + +void atlas__FieldSet__data_long_specf(FieldSetImpl* This, char* name, long*& data, int& rank, int*& shapef, int*& stridesf) { + atlas__Field__data_long_specf(This->field(name).get(), data, rank, shapef, stridesf); +} + +void atlas__FieldSet__data_float_specf(FieldSetImpl* This, char* name, float*& data, int& rank, int*& shapef, int*& stridesf) { + atlas__Field__data_float_specf(This->field(name).get(), data, rank, shapef, stridesf); +} + +void atlas__FieldSet__data_double_specf(FieldSetImpl* This, char* name, double*& data, int& rank, int*& shapef, int*& stridesf) { + atlas__Field__data_double_specf(This->field(name).get(), data, rank, shapef, stridesf); +} + +void atlas__FieldSet__data_int_specf_by_idx(FieldSetImpl* This, int& idx, int*& data, int& rank, int*& shapef, int*& stridesf) { + atlas__Field__data_int_specf(This->operator[](idx).get(), data, rank, shapef, stridesf); +} + +void atlas__FieldSet__data_long_specf_by_idx(FieldSetImpl* This, int& idx, long*& data, int& rank, int*& shapef, int*& stridesf) { + atlas__Field__data_long_specf(This->operator[](idx).get(), data, rank, shapef, stridesf); +} + +void atlas__FieldSet__data_float_specf_by_idx(FieldSetImpl* This, int& idx, float*& data, int& rank, int*& shapef, int*& stridesf) { + atlas__Field__data_float_specf(This->operator[](idx).get(), data, rank, shapef, stridesf); +} + +void atlas__FieldSet__data_double_specf_by_idx(FieldSetImpl* This, int& idx, double*& data, int& rank, int*& shapef, int*& stridesf) { + atlas__Field__data_double_specf(This->operator[](idx).get(), data, rank, shapef, stridesf); +} + + FieldSetImpl* atlas__FieldSet__new(char* name) { FieldSetImpl* fset = new FieldSetImpl(std::string(name)); fset->name() = name; diff --git a/src/atlas/field/FieldSet.h b/src/atlas/field/FieldSet.h index 4bf5e1ff3..83bfc1f15 100644 --- a/src/atlas/field/FieldSet.h +++ b/src/atlas/field/FieldSet.h @@ -24,6 +24,7 @@ #include "eckit/deprecated.h" +#include "atlas/array_fwd.h" #include "atlas/field/Field.h" #include "atlas/library/config.h" #include "atlas/runtime/Exception.h" @@ -139,6 +140,22 @@ const char* atlas__FieldSet__name(FieldSetImpl* This); idx_t atlas__FieldSet__size(const FieldSetImpl* This); FieldImpl* atlas__FieldSet__field_by_name(FieldSetImpl* This, char* name); FieldImpl* atlas__FieldSet__field_by_idx(FieldSetImpl* This, idx_t idx); +void atlas__FieldSet__data_int_specf(FieldSetImpl* This, char* name, int*& field_data, int& rank, int*& field_shapef, + int*& field_stridesf); +void atlas__FieldSet__data_long_specf(FieldSetImpl* This, char* name, long*& field_data, int& rank, int*& field_shapef, + int*& field_stridesf); +void atlas__FieldSet__data_float_specf(FieldSetImpl* This, char* name, float*& field_data, int& rank, int*& field_shapef, + int*& field_stridesf); +void atlas__FieldSet__data_double_specf(FieldSetImpl* This, char* name, double*& field_data, int& rank, int*& field_shapef, + int*& field_stridesf); +void atlas__FieldSet__data_int_specf_by_idx(FieldSetImpl* This, int& idx, int*& field_data, int& rank, int*& field_shapef, + int*& field_stridesf); +void atlas__FieldSet__data_long_specf_by_idx(FieldSetImpl* This, int& idx, long*& field_data, int& rank, int*& field_shapef, + int*& field_stridesf); +void atlas__FieldSet__data_float_specf_by_idx(FieldSetImpl* This, int& idx, float*& field_data, int& rank, int*& field_shapef, + int*& field_stridesf); +void atlas__FieldSet__data_double_specf_by_idx(FieldSetImpl* This, int& idx, double*& field_data, int& rank, int*& field_shapef, + int*& field_stridesf); void atlas__FieldSet__set_dirty(FieldSetImpl* This, int value); void atlas__FieldSet__halo_exchange(FieldSetImpl* This, int on_device); } diff --git a/src/atlas/functionspace/BlockStructuredColumns.h b/src/atlas/functionspace/BlockStructuredColumns.h index f4a2f3c93..6f30ae314 100644 --- a/src/atlas/functionspace/BlockStructuredColumns.h +++ b/src/atlas/functionspace/BlockStructuredColumns.h @@ -42,8 +42,6 @@ class BlockStructuredColumns : public FunctionSpace { bool valid() const { return functionspace_; } idx_t size() const { return functionspace_->size(); } -// idx_t sizeOwned() const { return functionspace_->sizeOwned(); } -// idx_t sizeHalo() const { return functionspace_->sizeHalo(); } idx_t levels() const { return functionspace_->levels(); } const Vertical& vertical() const { return functionspace_->vertical(); } @@ -54,10 +52,6 @@ class BlockStructuredColumns : public FunctionSpace { std::string checksum(const Field&) const; idx_t index(idx_t blk, idx_t rof) const { return functionspace_->index(blk, rof); } -// idx_t i_begin(idx_t j) const { return functionspace_->i_begin(j); } -// idx_t i_end(idx_t j) const { return functionspace_->i_end(j); } -// idx_t j_begin() const { return functionspace_->j_begin(); } -// idx_t j_end() const { return functionspace_->j_end(); } idx_t k_begin() const { return functionspace_->k_begin(); } idx_t k_end() const { return functionspace_->k_end(); } idx_t nproma() const { return functionspace_->nproma(); } diff --git a/src/atlas/functionspace/detail/BlockStructuredColumns.cc b/src/atlas/functionspace/detail/BlockStructuredColumns.cc index 6041ca1cc..c18b1756d 100644 --- a/src/atlas/functionspace/detail/BlockStructuredColumns.cc +++ b/src/atlas/functionspace/detail/BlockStructuredColumns.cc @@ -156,7 +156,6 @@ void rev_block_copy(const Field loc, Field sloc, const functionspace::detail::Bl } } - void transpose_nonblocked_to_blocked(const Field& nonblocked, Field& blocked, const functionspace::detail::BlockStructuredColumns& fs) { auto kind = nonblocked.datatype().kind(); if (kind == array::DataType::kind()) { diff --git a/src/atlas/functionspace/detail/StructuredColumns.h b/src/atlas/functionspace/detail/StructuredColumns.h index 392ed577c..8385a5777 100644 --- a/src/atlas/functionspace/detail/StructuredColumns.h +++ b/src/atlas/functionspace/detail/StructuredColumns.h @@ -317,6 +317,7 @@ class StructuredColumns : public FunctionSpaceImpl { friend struct BlockStructuredColumnsFortranAccess; Map2to1 ij2gp_; + friend class BlockStructuredColumns; void setup(const grid::Distribution& distribution, const eckit::Configuration& config); friend class BlockStructuredColumns; diff --git a/src/atlas_f/CMakeLists.txt b/src/atlas_f/CMakeLists.txt index f015676a5..29d55ab99 100644 --- a/src/atlas_f/CMakeLists.txt +++ b/src/atlas_f/CMakeLists.txt @@ -175,7 +175,7 @@ ecbuild_add_library( TARGET atlas_f functionspace/atlas_functionspace_BlockStructuredColumns_module.F90 functionspace/atlas_functionspace_Spectral_module.F90 functionspace/atlas_functionspace_PointCloud_module.F90 - field/atlas_FieldSet_module.F90 + field/atlas_FieldSet_module.fypp field/atlas_State_module.F90 field/atlas_Field_module.fypp grid/atlas_Grid_module.F90 diff --git a/src/atlas_f/field/atlas_FieldSet_module.F90 b/src/atlas_f/field/atlas_FieldSet_module.fypp similarity index 65% rename from src/atlas_f/field/atlas_FieldSet_module.F90 rename to src/atlas_f/field/atlas_FieldSet_module.fypp index a8a09457c..e337c958c 100644 --- a/src/atlas_f/field/atlas_FieldSet_module.F90 +++ b/src/atlas_f/field/atlas_FieldSet_module.fypp @@ -8,9 +8,13 @@ #include "atlas/atlas_f.h" +#:include "atlas/atlas_f.fypp" +#:include "internals/atlas_generics.fypp" + module atlas_FieldSet_module use fckit_owned_object_module, only: fckit_owned_object +use atlas_field_module, only: atlas_field, array_c_to_f use atlas_kinds_module, only : ATLAS_KIND_IDX implicit none @@ -53,10 +57,33 @@ module atlas_FieldSet_module procedure, public :: set_dirty procedure, public :: halo_exchange +#:for rank in ranks +#:for dtype in dtypes + procedure, private :: access_data_${dtype}$_r${rank}$_by_name + procedure, private :: access_data_${dtype}$_r${rank}$_by_idx + procedure, private :: access_data_${dtype}$_r${rank}$_slice_by_name + procedure, private :: access_data_${dtype}$_r${rank}$_slice_by_idx +#:endfor +#:endfor + + generic, public :: data => & +#:for rank in ranks +#:for dtype in dtypes + & access_data_${dtype}$_r${rank}$_by_name, & + & access_data_${dtype}$_r${rank}$_by_idx, & + & access_data_${dtype}$_r${rank}$_slice_by_name, & + & access_data_${dtype}$_r${rank}$_slice_by_idx, & +#:endfor +#:endfor + & dummy + + procedure, private :: dummy + #if FCKIT_FINAL_NOT_INHERITING final :: atlas_FieldSet__final_auto #endif procedure, public :: has_field => has ! deprecated ! + END TYPE atlas_FieldSet !------------------------------------------------------------------------------ @@ -74,6 +101,86 @@ module atlas_FieldSet_module ! ----------------------------------------------------------------------------- ! FieldSet routines +#:for rank in ranks +#:for dtype,ftype,ctype in types +subroutine access_data_${dtype}$_r${rank}$_by_name(this, name, field) + use fckit_c_interop_module, only: c_str + use atlas_fieldset_c_binding + use, intrinsic :: iso_c_binding, only : c_ptr, c_int, c_long, c_float, c_double + class(atlas_FieldSet), intent(in) :: this + character(len=*), intent(in) :: name + ${ftype}$, pointer, intent(inout) :: field(${dim[rank]}$) + type(c_ptr) :: field_cptr + type(c_ptr) :: shape_cptr + type(c_ptr) :: strides_cptr + integer(c_int) :: rank + call atlas__FieldSet__data_${ctype}$_specf(this%CPTR_PGIBUG_A, c_str(name), field_cptr, rank, shape_cptr, strides_cptr) + call array_c_to_f(field_cptr, rank, shape_cptr, strides_cptr, field) +end subroutine +subroutine access_data_${dtype}$_r${rank}$_by_idx(this, idx, field) + use fckit_c_interop_module, only: c_str + use atlas_fieldset_c_binding + use, intrinsic :: iso_c_binding, only : c_ptr, c_int, c_long, c_float, c_double + class(atlas_FieldSet), intent(in) :: this + integer, intent(in) :: idx + ${ftype}$, pointer, intent(inout) :: field(${dim[rank]}$) + type(c_ptr) :: field_cptr + type(c_ptr) :: shape_cptr + type(c_ptr) :: strides_cptr + integer(c_int) :: rank + call atlas__FieldSet__data_${ctype}$_specf_by_idx(this%CPTR_PGIBUG_A, idx-1, field_cptr, rank, shape_cptr, strides_cptr) + call array_c_to_f(field_cptr, rank, shape_cptr, strides_cptr, field) +end subroutine +subroutine access_data_${dtype}$_r${rank}$_slice_by_name(this, name, slice, iblk) + use fckit_c_interop_module, only: c_str + use atlas_fieldset_c_binding + use, intrinsic :: iso_c_binding, only : c_ptr, c_int, c_long, c_float, c_double + class(atlas_FieldSet), intent(in) :: this + character(len=*), intent(in) :: name +#:if rank > 1 + ${ftype}$, pointer, intent(inout) :: slice(${dimr[rank]}$) +#:else + ${ftype}$, pointer, intent(inout) :: slice +#:endif + integer, intent(in) :: iblk + ${ftype}$, pointer :: field(${dim[rank]}$) + call access_data_${dtype}$_r${rank}$_by_name(this, c_str(name), field) +#:if rank > 1 + slice => field(${dimr[rank]}$, iblk) +#:else + slice => field(iblk) +#:endif +end subroutine +subroutine access_data_${dtype}$_r${rank}$_slice_by_idx(this, idx, slice, iblk) + use fckit_c_interop_module, only: c_str + use atlas_fieldset_c_binding + use, intrinsic :: iso_c_binding, only : c_ptr, c_int, c_long, c_float, c_double + class(atlas_FieldSet), intent(in) :: this + integer, intent(in) :: idx +#:if rank > 1 + ${ftype}$, pointer, intent(inout) :: slice(${dimr[rank]}$) +#:else + ${ftype}$, pointer, intent(inout) :: slice +#:endif + integer, intent(in) :: iblk + ${ftype}$, pointer :: field(${dim[rank]}$) + call access_data_${dtype}$_r${rank}$_by_idx(this, idx, field) +#:if rank > 1 + slice => field(${dimr[rank]}$, iblk) +#:else + slice => field(iblk) +#:endif +end subroutine +!------------------------------------------------------------------------------- + +#:endfor +#:endfor +subroutine dummy(this) + use atlas_fieldset_c_binding + class(atlas_FieldSet), intent(in) :: this + FCKIT_SUPPRESS_UNUSED(this) +end subroutine + function atlas_FieldSet__cptr(cptr) result(fieldset) use, intrinsic :: iso_c_binding, only: c_ptr type(atlas_FieldSet) :: fieldset diff --git a/src/atlas_f/field/atlas_Field_module.fypp b/src/atlas_f/field/atlas_Field_module.fypp index 6f859a5b2..68cc5f04a 100644 --- a/src/atlas_f/field/atlas_Field_module.fypp +++ b/src/atlas_f/field/atlas_Field_module.fypp @@ -23,6 +23,7 @@ public :: atlas_real public :: atlas_integer public :: atlas_logical public :: atlas_data_type +public :: array_c_to_f private @@ -75,6 +76,7 @@ contains #:for dtype in dtypes procedure, private :: access_data_${dtype}$_r${rank}$ procedure, private :: access_data_${dtype}$_r${rank}$_shape + procedure, private :: access_data_${dtype}$_r${rank}$_slice #:endfor #:endfor @@ -83,6 +85,7 @@ contains #:for dtype in dtypes & access_data_${dtype}$_r${rank}$, & & access_data_${dtype}$_r${rank}$_shape, & + & access_data_${dtype}$_r${rank}$_slice, & #:endfor #:endfor & dummy @@ -135,7 +138,6 @@ interface array_c_to_f end interface !------------------------------------------------------------------------------- - private :: fckit_owned_object private :: atlas_Config @@ -161,6 +163,7 @@ subroutine array_c_to_f_${dtype}$_r${rank}$(array_cptr,rank,shape_cptr,strides_c integer :: j if( rank /= ${rank}$ ) then + write(0,*) rank, "!=", ${rank}$ @:ATLAS_ABORT("Rank mismatch") endif @@ -199,6 +202,29 @@ subroutine access_data_${dtype}$_r${rank}$(this, field) call atlas__Field__data_${ctype}$_specf(this%CPTR_PGIBUG_A, field_cptr, rank, shape_cptr, strides_cptr) call array_c_to_f(field_cptr,rank,shape_cptr,strides_cptr, field) end subroutine +subroutine access_data_${dtype}$_r${rank}$_slice(this, slice, iblk) + use atlas_field_c_binding + use, intrinsic :: iso_c_binding, only : c_ptr, c_int, c_long, c_float, c_double + class(atlas_Field), intent(in) :: this + ${ftype}$, pointer :: field(${dim[rank]}$) +#:if rank > 1 + ${ftype}$, pointer, intent(inout) :: slice(${dimr[rank]}$) +#:else + ${ftype}$, pointer, intent(inout) :: slice +#:endif + integer, intent(in) :: iblk + type(c_ptr) :: field_cptr + type(c_ptr) :: shape_cptr + type(c_ptr) :: strides_cptr + integer(c_int) :: rank + call atlas__Field__data_${ctype}$_specf(this%CPTR_PGIBUG_A, field_cptr, rank, shape_cptr, strides_cptr) + call array_c_to_f(field_cptr, rank, shape_cptr, strides_cptr, field) +#:if rank > 1 + slice => field(${dimr[rank]}$, iblk) +#:else + slice => field(iblk) +#:endif +end subroutine !------------------------------------------------------------------------------- diff --git a/src/atlas_f/internals/atlas_generics.fypp b/src/atlas_f/internals/atlas_generics.fypp index 78d78c353..d99c7f565 100644 --- a/src/atlas_f/internals/atlas_generics.fypp +++ b/src/atlas_f/internals/atlas_generics.fypp @@ -10,6 +10,7 @@ #:set ranks = [1,2,3,4] #:set dim = ['',':',':,:',':,:,:',':,:,:,:',':,:,:,:,:'] +#:set dimr = ['','',':',':,:',':,:,:',':,:,:,:'] #:set ftypes = ['integer(c_int)','integer(c_long)','real(c_float)','real(c_double)', 'logical'] #:set ctypes = ['int','long','float','double', 'int'] diff --git a/src/tests/field/fctest_field.F90 b/src/tests/field/fctest_field.F90 index a0264f7e4..8ea98e09d 100644 --- a/src/tests/field/fctest_field.F90 +++ b/src/tests/field/fctest_field.F90 @@ -172,7 +172,7 @@ module fcta_Field_fixture implicit none type(atlas_Field) :: field real(c_double), pointer :: view(:,:,:) - field = atlas_Field("field_2",atlas_real(c_double),(/3,5,10/),alignment=4) + field = atlas_Field("field_3",atlas_real(c_double),(/3,5,10/),alignment=4) call field%data(view) FCTEST_CHECK_EQUAL( size(view,1) , 3 ) FCTEST_CHECK_EQUAL( size(view,2) , 5 ) @@ -183,6 +183,72 @@ module fcta_Field_fixture FCTEST_CHECK_EQUAL( field%stride(1), 1 ) FCTEST_CHECK_EQUAL( field%stride(2), 4 ) FCTEST_CHECK_EQUAL( field%stride(3), 4*5 ) + call field%final() +END_TEST + +TEST( test_fieldset_slice ) + implicit none + type(atlas_FieldSet) :: fieldset + type(atlas_Field) :: field + integer, pointer :: view1d(:), view3d(:,:,:) + integer, pointer :: slice0d, slice2d(:,:) + + ! slicing of a three-dimensional field + field = atlas_Field("field_4",atlas_integer(),(/3,5,10/)) + call field%data(view3d) + view3d(1,2,3) = 123 + call field%data(slice2d,3) + FCTEST_CHECK_EQUAL( size(slice2d,1) , 3 ) + FCTEST_CHECK_EQUAL( size(slice2d,2) , 5 ) + FCTEST_CHECK_EQUAL( slice2d(1,2), 123 ) + slice2d(1,2) = slice2d(1,2) * 2 + call field%data(view3d) + FCTEST_CHECK_EQUAL( view3d(1,2,3), 246 ) + + ! slicing of a one-dimensional field + field = atlas_Field("field_5",atlas_integer(),(/3/)) + call field%data(view1d) + view1d(2) = 123 + call field%data(slice0d,2) + FCTEST_CHECK_EQUAL( slice0d, 123 ) + slice0d = slice0d * 2 + call field%data(view1d) + FCTEST_CHECK_EQUAL( view1d(2), 246 ) + call field%final() + + ! slicing of a three-dimensional field through a fieldset by name and by idx + fieldset = atlas_FieldSet() + field = atlas_Field("field_6",atlas_integer(),(/3,5,10/)) + call fieldset%add(field) + field = atlas_Field("field_7",atlas_integer(),(/3/)) + call fieldset%add(field) + call fieldset%data("field_6", view3d) + view3d(1,2,3) = 122 + call fieldset%data(1, view3d) + view3d(1,2,3) = 123 + call fieldset%data("field_6", slice2d, 3) + slice2d(1,2) = slice2d(1,2) + 1 + call fieldset%data(1, slice2d, 3) + slice2d(1,2) = slice2d(1,2) - 1 + FCTEST_CHECK_EQUAL( size(slice2d,1) , 3 ) + FCTEST_CHECK_EQUAL( size(slice2d,2) , 5 ) + FCTEST_CHECK_EQUAL( slice2d(1,2), 123 ) + slice2d(1,2) = slice2d(1,2) * 2 + call fieldset%data('field_6', view3d) + FCTEST_CHECK_EQUAL( view3d(1,2,3), 246 ) + + ! slicing of a one-dimensional field through a fieldset by name and by idx + call fieldset%data('field_7', view1d) + view1d(2) = 122 + call fieldset%data(2, view1d) + view1d(2) = 123 + call field%data(slice0d, 2) + FCTEST_CHECK_EQUAL( slice0d, 123 ) + slice0d = slice0d * 2 + call field%data(view1d) + FCTEST_CHECK_EQUAL( view1d(2), 246 ) + call field%final() + call fieldset%final() END_TEST ! ----------------------------------------------------------------------------- diff --git a/src/tests/functionspace/fctest_blockstructuredcolumns.F90 b/src/tests/functionspace/fctest_blockstructuredcolumns.F90 index 708f6e60e..6d69bef50 100644 --- a/src/tests/functionspace/fctest_blockstructuredcolumns.F90 +++ b/src/tests/functionspace/fctest_blockstructuredcolumns.F90 @@ -7,6 +7,7 @@ ! This File contains Unit Tests for testing the ! C++ / Fortran Interfaces to the State Datastructure +! ! @author Willem Deconinck ! @author Slavko Brdar