From 35cae670bcd8b312343cc6dd89ef0bba565c08ba Mon Sep 17 00:00:00 2001 From: Slavko Brdar Date: Tue, 14 Nov 2023 14:20:28 +0000 Subject: [PATCH] improve fortran sandbox/atlas-filter to actually do a filtering --- src/sandbox/interpolation/atlas-filter.F90 | 95 +++++++++++++++------- 1 file changed, 67 insertions(+), 28 deletions(-) diff --git a/src/sandbox/interpolation/atlas-filter.F90 b/src/sandbox/interpolation/atlas-filter.F90 index d749f56f7..852f37b08 100644 --- a/src/sandbox/interpolation/atlas-filter.F90 +++ b/src/sandbox/interpolation/atlas-filter.F90 @@ -27,38 +27,63 @@ module filter_module implicit none -public :: atlas_Filter +public :: atlas_Filter, JPRB private - INTEGER, PARAMETER :: JPRB = SELECTED_REAL_KIND(13,300) +INTEGER, PARAMETER:: JPRB = SELECTED_REAL_KIND(13,300) +! -------------------------------------------------------------------------------------------- type :: atlas_Filter - TYPE(atlas_Redistribution) :: src_redist, tgt_redist - TYPE(atlas_Interpolation) :: interpolation_st, interpolation_ts - TYPE(atlas_Field) :: tgt_field, src_field_tgtpart, tgt_field_srcpart +private + type(atlas_FunctionSpace) :: src_fs, tgt_fs + type(atlas_Redistribution):: src_redist, tgt_redist + type(atlas_Interpolation) :: interpolation_st, interpolation_ts + type(atlas_Field) :: tgt_field, src_field_tgtpart, tgt_field_srcpart contains - procedure, public :: setup => filter_setup procedure, public :: execute => filter_execute + procedure, public :: source => filter_src_fs + procedure, public :: target => filter_tgt_fs final :: filter_finalise end type atlas_Filter +! -------------------------------------------------------------------------------------------- + +interface atlas_Filter + module procedure atlas_Filter__create +end interface CONTAINS ! -------------------------------------------------------------------------------------------- -SUBROUTINE FILTER_SETUP(this) - class(atlas_Filter), intent(inout) :: this - type(atlas_Grid) :: src_grid, tgt_grid - type(atlas_MeshGenerator) :: meshgen +function filter_src_fs(this) result(fs) + class(atlas_Filter), intent(in) :: this + type(atlas_FunctionSpace) :: fs + fs = this%src_fs +end function filter_src_fs + +! -------------------------------------------------------------------------------------------- + +function filter_tgt_fs(this) result(fs) + class(atlas_Filter), intent(in) :: this + type(atlas_FunctionSpace) :: fs + fs = this%tgt_fs +end function filter_tgt_fs + +! -------------------------------------------------------------------------------------------- + +function atlas_Filter__create() result(this) + type(atlas_Filter) :: this + type(atlas_Grid) :: src_grid, tgt_grid + type(atlas_MeshGenerator) :: meshgen type(atlas_GridDistribution) :: griddist - type(atlas_Mesh) :: src_mesh, tgt_mesh - type(atlas_Mesh) :: src_mesh_tgtpart, tgt_mesh_srcpart - type(atlas_Config) :: interpolation_config - type(atlas_FunctionSpace) :: src_fs, tgt_fs - type(atlas_FunctionSpace) :: src_fs_tgtpart, tgt_fs_srcpart - type(atlas_Redistribution) :: src_redist, tgt_redist - real(kind=JPRB), pointer :: src_val(:) + type(atlas_Redistribution) :: src_redist, tgt_redist + type(atlas_Mesh) :: src_mesh, tgt_mesh + type(atlas_Mesh) :: src_mesh_tgtpart, tgt_mesh_srcpart + type(atlas_Config) :: interpolation_config + type(atlas_FunctionSpace) :: src_fs, tgt_fs + type(atlas_FunctionSpace) :: src_fs_tgtpart, tgt_fs_srcpart + real(kind=JPRB), pointer :: src_val(:) src_grid = atlas_StructuredGrid("O80") tgt_grid = atlas_StructuredGrid("O40") @@ -69,6 +94,8 @@ SUBROUTINE FILTER_SETUP(this) tgt_mesh = meshgen%generate(tgt_grid, griddist) src_fs = atlas_functionspace_NodeColumns(src_mesh, halo=4) tgt_fs = atlas_functionspace_NodeColumns(tgt_mesh, halo=2) + this%src_fs = src_fs + this%tgt_fs = tgt_fs ! // redistribution setup griddist = atlas_GridDistribution(src_grid, atlas_MatchingPartitioner(tgt_mesh)) @@ -77,8 +104,9 @@ SUBROUTINE FILTER_SETUP(this) tgt_mesh_srcpart = meshgen%generate(tgt_grid, griddist) src_fs_tgtpart = atlas_functionspace_NodeColumns(src_mesh_tgtpart) tgt_fs_srcpart = atlas_functionspace_NodeColumns(tgt_mesh_srcpart) - src_redist = atlas_Redistribution(src_fs, src_fs_tgtpart) - tgt_redist = atlas_Redistribution(tgt_fs, tgt_fs_srcpart) + + this%src_redist = atlas_Redistribution(src_fs, src_fs_tgtpart) + this%tgt_redist = atlas_Redistribution(tgt_fs, tgt_fs_srcpart) ! // interpolation setup interpolation_config = atlas_Config() @@ -98,9 +126,7 @@ SUBROUTINE FILTER_SETUP(this) call tgt_mesh_srcpart%final() call src_fs_tgtpart%final() call tgt_fs_srcpart%final() - call src_fs%final() - call tgt_fs%final() -END SUBROUTINE FILTER_SETUP +end function atlas_Filter__create ! -------------------------------------------------------------------------------------------- @@ -155,15 +181,28 @@ end module filter_module program filtering use atlas_module -use filter_module, only: atlas_Filter +use filter_module, only: JPRB, atlas_Filter implicit none -type(atlas_Grid) :: grid -type(atlas_Filter) :: filter + type(atlas_Grid) :: grid + type(atlas_Filter) :: filter + type(atlas_FunctionSpace) :: fs, ts + type(atlas_Field) :: sfield + real(kind=JPRB), pointer :: sfield_v(:) + + call atlas_library%initialise() + + filter = atlas_Filter() + fs = filter%source() + ts = filter%target() + + sfield = fs%create_field(name="sfield", kind=atlas_real(JPRB)) + call sfield%data(sfield_v) + sfield_v = 1._JPRB + + call filter.execute(sfield) -call filter%setup() -grid = atlas_Grid("O40") -print *, "size ", grid%size() + call atlas_library%finalise() end program filtering