diff --git a/Makefile b/Makefile index a5d1adb5f..101501fd0 100644 --- a/Makefile +++ b/Makefile @@ -22,6 +22,7 @@ OPERATE_ON ?= infrastructure \ mesh_tools \ applications/skeleton \ applications/simple_diffusion \ + applications/coupled \ applications/lbc_demo \ applications/io_demo diff --git a/applications/coupled/example/configuration_glo.nml b/applications/coupled/example/configuration_glo.nml index e212a5e2a..b2d2bb1f5 100644 --- a/applications/coupled/example/configuration_glo.nml +++ b/applications/coupled/example/configuration_glo.nml @@ -58,6 +58,9 @@ coord_system='native' &partitioning partitioner='cubedsphere', panel_decomposition = 'auto', + tile_size_x=1, + tile_size_y=1, + inner_halo_tiles=.false. / &planet diff --git a/applications/coupled/example/configuration_lam.nml b/applications/coupled/example/configuration_lam.nml index 5e92d9fe9..a584cfa73 100644 --- a/applications/coupled/example/configuration_lam.nml +++ b/applications/coupled/example/configuration_lam.nml @@ -57,6 +57,9 @@ coord_system='native' &partitioning partitioner='planar', panel_decomposition = 'auto', + tile_size_x=1, + tile_size_y=1, + inner_halo_tiles=.false. / &planet diff --git a/applications/coupled/source/coupled.f90 b/applications/coupled/source/coupled.f90 index a74455853..00abd995c 100644 --- a/applications/coupled/source/coupled.f90 +++ b/applications/coupled/source/coupled.f90 @@ -50,7 +50,8 @@ program coupled call init_config( filename, coupled_required_namelists, & config=modeldb%config ) - call init_logger( modeldb%mpi%get_comm(), & + call init_logger( modeldb%config, & + modeldb%mpi%get_comm(), & program_name//"_"//cpl_component_name ) write(log_scratch_space,'(A)') & diff --git a/applications/coupled/source/driver/coupled_driver_mod.f90 b/applications/coupled/source/driver/coupled_driver_mod.f90 index ee3fda2aa..205990e8b 100644 --- a/applications/coupled/source/driver/coupled_driver_mod.f90 +++ b/applications/coupled/source/driver/coupled_driver_mod.f90 @@ -75,11 +75,17 @@ subroutine initialise( program_name, modeldb, calendar ) integer(i_def) :: geometry integer(i_def) :: method integer(i_def) :: number_of_layers + integer(i_def) :: tile_size_x + integer(i_def) :: tile_size_y real(r_def) :: domain_bottom real(r_def) :: domain_height real(r_def) :: scaled_radius logical :: check_partitions + logical :: inner_halo_tiles + logical :: prepartitioned + + integer(i_def), allocatable :: tile_size(:,:) integer(i_def) :: i integer(i_def), parameter :: one_layer = 1_i_def @@ -91,6 +97,17 @@ subroutine initialise( program_name, modeldb, calendar ) domain_height = modeldb%config%extrusion%domain_height() number_of_layers = modeldb%config%extrusion%number_of_layers() scaled_radius = modeldb%config%planet%scaled_radius() + prepartitioned = modeldb%config%base_mesh%prepartitioned() + + if (prepartitioned) then + tile_size_x = 1 + tile_size_y = 1 + inner_halo_tiles = .false. + else + tile_size_x = maxval([1,modeldb%config%partitioning%tile_size_x()]) + tile_size_y = maxval([1,modeldb%config%partitioning%tile_size_y()]) + inner_halo_tiles = modeldb%config%partitioning%inner_halo_tiles() + end if ! Initialise mesh ! Determine the required meshes @@ -120,23 +137,29 @@ subroutine initialise( program_name, modeldb, calendar ) ! Create the required meshes stencil_depth = 1 check_partitions = .false. + allocate(tile_size(2,size(base_mesh_names))) + tile_size(1,:) = tile_size_x + tile_size(2,:) = tile_size_y call init_mesh( modeldb%config, & modeldb%mpi%get_comm_rank(), & modeldb%mpi%get_comm_size(), & base_mesh_names, extrusion, & + inner_halo_tiles, tile_size, & stencil_depth, check_partitions ) allocate( twod_names, source=base_mesh_names ) do i=1, size(twod_names) twod_names(i) = trim(twod_names(i))//'_2d' end do + call create_mesh( base_mesh_names, extrusion_2d, & + inner_halo_tiles, tile_size, & alt_name=twod_names ) call assign_mesh_maps( twod_names ) ! Build the FEM function spaces and coordinate fields - call init_fem( mesh_collection, chi_inventory, panel_id_inventory ) + call init_fem( modeldb%config, chi_inventory, panel_id_inventory ) ! Create and initialise prognostic fields mesh => mesh_collection%get_mesh(prime_mesh_name) diff --git a/applications/coupled/source/driver/init_coupled_mod.X90 b/applications/coupled/source/driver/init_coupled_mod.X90 index 3028d5a06..8715a45a2 100644 --- a/applications/coupled/source/driver/init_coupled_mod.X90 +++ b/applications/coupled/source/driver/init_coupled_mod.X90 @@ -76,6 +76,11 @@ module init_coupled_mod procedure(write_interface), pointer :: tmp_ptr integer(i_def) :: order_h, order_v + integer(i_def) :: coord_system + integer(i_def) :: geometry + integer(i_def) :: topology + real(r_def) :: scaled_radius + type(function_space_type), pointer :: fs call log_event( 'coupled: Initialising app ...', LOG_LEVEL_INFO ) @@ -83,8 +88,12 @@ module init_coupled_mod ! Get the name of the coupling component call modeldb%values%get_value("cpl_name", cpl_component_name) - order_h = modeldb%config%finite_element%element_order_h() - order_v = modeldb%config%finite_element%element_order_v() + order_h = modeldb%config%finite_element%element_order_h() + order_v = modeldb%config%finite_element%element_order_v() + coord_system = modeldb%config%finite_element%coord_system() + geometry = modeldb%config%base_mesh%geometry() + topology = modeldb%config%base_mesh%topology() + scaled_radius = modeldb%config%planet%scaled_radius() fs => function_space_collection%get_fs(mesh, order_h, order_v, W3) @@ -122,8 +131,10 @@ module init_coupled_mod call depository%get_field( trim(name), field_1_ptr) ! Initialise the values in the field that will be sent to the coupler. ! Set them to the longitude of the cell-centre (converted to degrees) - call invoke(compute_latlon_kernel_type(field_2, field_1_ptr, & - chi, panel_id), & + call invoke(compute_latlon_kernel_type(field_2, field_1_ptr, & + chi, panel_id, geometry, & + topology, coord_system, & + scaled_radius), & inc_a_times_X(radians_to_degrees, field_1_ptr)) ! Add that field to the coupling 2d "send" field collection abs_field_ptr => field_1_ptr diff --git a/applications/io_demo/example/configuration.nml b/applications/io_demo/example/configuration.nml index 03aee6514..4770bf3a0 100644 --- a/applications/io_demo/example/configuration.nml +++ b/applications/io_demo/example/configuration.nml @@ -82,6 +82,9 @@ coord_system='native' &partitioning partitioner='cubedsphere', panel_decomposition = 'auto', + tile_size_x=1, + tile_size_y=1, + inner_halo_tiles=.false. / &planet diff --git a/applications/io_demo/rose-meta/lfric-io_demo/vn3.1/rose-meta.conf b/applications/io_demo/rose-meta/lfric-io_demo/vn3.1/rose-meta.conf index bd3389f88..2af56e5b0 100644 --- a/applications/io_demo/rose-meta/lfric-io_demo/vn3.1/rose-meta.conf +++ b/applications/io_demo/rose-meta/lfric-io_demo/vn3.1/rose-meta.conf @@ -1,6 +1,28 @@ import=lfric-driver/vn3.1 -[namelist:io=multifile_io] +[namelist:io_demo] +compulsory=true +description=Provides options for configuring the runtime behaviour of the IO_Demo app +ns=namelist/io_demo +sort-key=Section-A02 +title=IO_Demo + +[namelist:io_demo=benchmark_sleep_time] +compulsory=true +description=Number of seconds to sleep for each timestep in I/O benchmark mode +!kind=default +type=integer + +[namelist:io_demo=io_benchmark] +compulsory=true +description=Configure application to run as an I/O benchmarking tool +help=Configure application to run as an I/O benchmarking tool +!kind=default +trigger=namelist:io_demo=benchmark_sleep_time: .true. ; + =namelist:io_demo=n_benchmark_fields: .true. ; +type=logical + +[namelist:io_demo=multifile_io] compulsory=true description=Use multifile_io functionality help=This is used to turn the multifile_io functionality in the io_demo app @@ -8,6 +30,12 @@ help=This is used to turn the multifile_io functionality in the io_demo app !kind=default type=logical +[namelist:io_demo=n_benchmark_fields] +compulsory=true +description=Number of fields created in I/O benchmark +!kind=default +type=integer + [namelist:multifile_io] compulsory=false duplicate=true diff --git a/applications/io_demo/source/driver/io_demo_driver_mod.f90 b/applications/io_demo/source/driver/io_demo_driver_mod.f90 index 24eb9d8d2..a56c77c52 100644 --- a/applications/io_demo/source/driver/io_demo_driver_mod.f90 +++ b/applications/io_demo/source/driver/io_demo_driver_mod.f90 @@ -86,15 +86,22 @@ subroutine initialise(program_name, modeldb) integer(i_def) :: stencil_depth(1) integer(i_def) :: geometry + integer(i_def) :: topology integer(i_def) :: method integer(i_def) :: number_of_layers + integer(i_def) :: tile_size_x + integer(i_def) :: tile_size_y real(r_def) :: domain_bottom real(r_def) :: domain_height real(r_def) :: scaled_radius - logical :: check_partitions - logical :: multifile_io - logical :: io_benchmark + logical :: check_partitions + logical :: inner_halo_tiles + logical :: prepartitioned + logical :: multifile_io + logical :: io_benchmark + + integer(i_def), allocatable :: tile_size(:,:) integer(i_def), parameter :: one_layer = 1_i_def integer(i_def) :: i @@ -107,13 +114,25 @@ subroutine initialise(program_name, modeldb) !======================================================================= prime_mesh_name = modeldb%config%base_mesh%prime_mesh_name() geometry = modeldb%config%base_mesh%geometry() + topology = modeldb%config%base_mesh%topology() method = modeldb%config%extrusion%method() domain_height = modeldb%config%extrusion%domain_height() number_of_layers = modeldb%config%extrusion%number_of_layers() scaled_radius = modeldb%config%planet%scaled_radius() + prepartitioned = modeldb%config%base_mesh%prepartitioned() multifile_io = modeldb%config%io_demo%multifile_io() io_benchmark = modeldb%config%io_demo%io_benchmark() + if (prepartitioned) then + tile_size_x = 1 + tile_size_y = 1 + inner_halo_tiles = .false. + else + tile_size_x = maxval([1,modeldb%config%partitioning%tile_size_x()]) + tile_size_y = maxval([1,modeldb%config%partitioning%tile_size_y()]) + inner_halo_tiles = modeldb%config%partitioning%inner_halo_tiles() + end if + !======================================================================= ! Mesh !======================================================================= @@ -161,10 +180,14 @@ subroutine initialise(program_name, modeldb) ! --------------------------------------------------------- stencil_depth = 1 check_partitions = .false. + allocate(tile_size(2,size(base_mesh_names))) + tile_size(1,:) = tile_size_x + tile_size(2,:) = tile_size_y call init_mesh( modeldb%config, & modeldb%mpi%get_comm_rank(), & modeldb%mpi%get_comm_size(), & base_mesh_names, extrusion, & + inner_halo_tiles, tile_size, & stencil_depth, check_partitions ) allocate( twod_names, source=base_mesh_names ) @@ -172,13 +195,14 @@ subroutine initialise(program_name, modeldb) twod_names(i) = trim(twod_names(i))//'_2d' end do call create_mesh( base_mesh_names, extrusion_2d, & + inner_halo_tiles, tile_size, & alt_name=twod_names ) call assign_mesh_maps(twod_names) !======================================================================= ! Build the FEM function spaces and coordinate fields !======================================================================= - call init_fem( mesh_collection, chi_inventory, panel_id_inventory ) + call init_fem( modeldb%config, chi_inventory, panel_id_inventory ) !======================================================================= ! Setup multifile reading @@ -201,10 +225,11 @@ subroutine initialise(program_name, modeldb) if (associated(files_init_ptr)) then call init_io( program_name, prime_mesh_name, modeldb, & chi_inventory, panel_id_inventory, & - populate_filelist=files_init_ptr ) + geometry, topology, populate_filelist=files_init_ptr ) else call init_io( program_name, prime_mesh_name, modeldb, & - chi_inventory, panel_id_inventory ) + chi_inventory, panel_id_inventory, & + geometry, topology ) end if diff --git a/applications/io_demo/source/driver/multifile_io/multifile_io_mod.F90 b/applications/io_demo/source/driver/multifile_io/multifile_io_mod.F90 index f15a5e764..5f2acba77 100644 --- a/applications/io_demo/source/driver/multifile_io/multifile_io_mod.F90 +++ b/applications/io_demo/source/driver/multifile_io/multifile_io_mod.F90 @@ -9,7 +9,7 @@ module multifile_io_mod use calendar_mod, only: calendar_type - use constants_mod, only: str_def, i_def + use constants_mod, only: str_def, i_def, r_def use driver_model_data_mod, only: model_data_type use driver_modeldb_mod, only: modeldb_type use empty_io_context_mod, only: empty_io_context_type @@ -123,11 +123,21 @@ subroutine step_multifile_io(modeldb, chi_inventory, panel_id_inventory) procedure(event_action), pointer :: context_advance procedure(callback_clock_arg), pointer :: before_close + integer(i_def) :: geometry + integer(i_def) :: topology + integer(i_def) :: coord_system + real(r_def) :: scaled_radius + nullify(mesh) nullify(chi) nullify(panel_id) nullify(before_close) + geometry = modeldb%config%base_mesh%geometry() + topology = modeldb%config%base_mesh%topology() + coord_system = modeldb%config%finite_element%coord_system() + scaled_radius = modeldb%config%planet%scaled_radius() + call iter%initialise(modeldb%config%multifile_io) do while (iter%has_next()) @@ -163,6 +173,8 @@ subroutine step_multifile_io(modeldb, chi_inventory, panel_id_inventory) chi, panel_id, & modeldb%clock, tmp_calendar, & before_close, & + geometry, topology, & + coord_system, scaled_radius, & start_at_zero=.true. ) ! Attach context advancement to the model's clock diff --git a/applications/io_demo/source/io_demo.f90 b/applications/io_demo/source/io_demo.f90 index 4a5c69e73..cbd310409 100644 --- a/applications/io_demo/source/io_demo.f90 +++ b/applications/io_demo/source/io_demo.f90 @@ -60,7 +60,9 @@ program io_demo deallocate( filename ) - call init_logger(modeldb%mpi%get_comm(), program_name) + call init_logger( modeldb%config, & + modeldb%mpi%get_comm(), & + program_name ) subroutine_timers = modeldb%config%io%subroutine_timers() timer_output_path = modeldb%config%io%timer_output_path() diff --git a/applications/lbc_demo/rose-meta/lfric-lbc_demo/version22_30.py b/applications/lbc_demo/rose-meta/lfric-lbc_demo/version22_30.py index c2098bef6..864b94e95 100644 --- a/applications/lbc_demo/rose-meta/lfric-lbc_demo/version22_30.py +++ b/applications/lbc_demo/rose-meta/lfric-lbc_demo/version22_30.py @@ -39,9 +39,8 @@ class vn22_t4231(MacroUpgrade): def upgrade(self, config, meta_config=None): # Add settings return config, self.reports -""" - +""" class vn22_t34(MacroUpgrade): # Upgrade macro for 34 by jennifer hickson @@ -50,3 +49,4 @@ class vn22_t34(MacroUpgrade): def upgrade(self, config, meta_config=None): return config, self.reports +""" diff --git a/applications/lbc_demo/source/driver/lbc_demo_driver_mod.f90 b/applications/lbc_demo/source/driver/lbc_demo_driver_mod.f90 index ce9a4a29d..9a1fc448c 100644 --- a/applications/lbc_demo/source/driver/lbc_demo_driver_mod.f90 +++ b/applications/lbc_demo/source/driver/lbc_demo_driver_mod.f90 @@ -86,10 +86,15 @@ subroutine initialise( program_name, modeldb) integer(i_def) :: geometry integer(i_def) :: method integer(i_def) :: number_of_layers + integer(i_def) :: tile_size_x + integer(i_def) :: tile_size_y real(r_def) :: domain_bottom real(r_def) :: domain_height real(r_def) :: scaled_radius - logical :: check_partitions + + logical :: check_partitions + logical :: inner_halo_tiles + logical :: prepartitioned logical :: write_diag @@ -100,22 +105,32 @@ subroutine initialise( program_name, modeldb) integer :: i + integer(i_def), allocatable :: tile_size(:,:) integer(i_def), parameter :: one_layer = 1_i_def !======================================================================= ! Extract configuration variables !======================================================================= - prime_mesh_name = modeldb%config%base_mesh%prime_mesh_name() - geometry = modeldb%config%base_mesh%geometry() - topology = modeldb%config%base_mesh%topology() - + prime_mesh_name = modeldb%config%base_mesh%prime_mesh_name() + geometry = modeldb%config%base_mesh%geometry() + topology = modeldb%config%base_mesh%topology() method = modeldb%config%extrusion%method() domain_height = modeldb%config%extrusion%domain_height() number_of_layers = modeldb%config%extrusion%number_of_layers() + scaled_radius = modeldb%config%planet%scaled_radius() + prepartitioned = modeldb%config%base_mesh%prepartitioned() - scaled_radius = modeldb%config%planet%scaled_radius() - write_diag = modeldb%config%io%write_diag() + if (prepartitioned) then + tile_size_x = 1 + tile_size_y = 1 + inner_halo_tiles = .false. + else + tile_size_x = maxval([1,modeldb%config%partitioning%tile_size_x()]) + tile_size_y = maxval([1,modeldb%config%partitioning%tile_size_y()]) + inner_halo_tiles = modeldb%config%partitioning%inner_halo_tiles() + end if + write_diag = modeldb%config%io%write_diag() enable_lbc = modeldb%config%lbc_demo%enable_lbc() apply_lbc = modeldb%config%lbc_demo%apply_lbc() write_lbc = modeldb%config%lbc_demo%write_lbc() @@ -185,11 +200,15 @@ subroutine initialise( program_name, modeldb) !------------------------------------------------------------------------- stencil_depth = 1 check_partitions = .false. + allocate(tile_size(2,size(base_mesh_names))) + tile_size(1,:) = tile_size_x + tile_size(2,:) = tile_size_y call init_mesh( modeldb%config, & modeldb%mpi%get_comm_rank(), & modeldb%mpi%get_comm_size(), & base_mesh_names, extrusion, & + inner_halo_tiles, tile_size, & stencil_depth, check_partitions ) allocate( twod_names, source=base_mesh_names ) @@ -198,14 +217,16 @@ subroutine initialise( program_name, modeldb) twod_names(i) = trim(twod_names(i))//'_2d' end do - call create_mesh( base_mesh_names, extrusion_2d, alt_name=twod_names ) + call create_mesh( base_mesh_names, extrusion_2d, & + inner_halo_tiles, tile_size, & + alt_name=twod_names ) call assign_mesh_maps( twod_names ) !======================================================================= ! Build the FEM function spaces and coordinate fields !======================================================================= - call init_fem( mesh_collection, chi_inventory, panel_id_inventory ) + call init_fem( modeldb%config, chi_inventory, panel_id_inventory ) !======================================================================= ! Setup general I/O system. @@ -229,7 +250,8 @@ subroutine initialise( program_name, modeldb) call log_event(log_scratch_space, log_level_info) call init_io( context_name, output_mesh_name, modeldb, & - chi_inventory, panel_id_inventory ) + chi_inventory, panel_id_inventory, & + geometry, topology ) end if !======================================================================= diff --git a/applications/lbc_demo/source/lbc_demo.f90 b/applications/lbc_demo/source/lbc_demo.f90 index 1500ada59..0f53931d6 100644 --- a/applications/lbc_demo/source/lbc_demo.f90 +++ b/applications/lbc_demo/source/lbc_demo.f90 @@ -45,7 +45,9 @@ program lbc_demo call init_config(filename, required_namelists, & config=modeldb%config) - call init_logger( modeldb%mpi%get_comm(), program_name ) + call init_logger( modeldb%config, & + modeldb%mpi%get_comm(), & + program_name ) ! Before anything else, test that the mesh provided was a regional domain. ! This application is not intended for cubed-sphere meshes. diff --git a/applications/simple_diffusion/example/configuration.nml b/applications/simple_diffusion/example/configuration.nml index 4694f1750..e32c05b0a 100644 --- a/applications/simple_diffusion/example/configuration.nml +++ b/applications/simple_diffusion/example/configuration.nml @@ -58,6 +58,9 @@ coord_system='native' &partitioning partitioner='cubedsphere', panel_decomposition = 'auto', + tile_size_x=1, + tile_size_y=1, + inner_halo_tiles=.false. / &planet diff --git a/applications/simple_diffusion/source/driver/simple_diffusion_driver_mod.f90 b/applications/simple_diffusion/source/driver/simple_diffusion_driver_mod.f90 index 8a5506f3b..b7d295a33 100644 --- a/applications/simple_diffusion/source/driver/simple_diffusion_driver_mod.f90 +++ b/applications/simple_diffusion/source/driver/simple_diffusion_driver_mod.f90 @@ -82,13 +82,20 @@ subroutine initialise( program_name, modeldb) integer(i_def) :: stencil_depth(1) integer(i_def) :: geometry + integer(i_def) :: topology integer(i_def) :: method integer(i_def) :: number_of_layers + integer(i_def) :: tile_size_x + integer(i_def) :: tile_size_y real(r_def) :: domain_bottom real(r_def) :: domain_height real(r_def) :: scaled_radius - logical :: check_partitions + logical :: check_partitions + logical :: inner_halo_tiles + logical :: prepartitioned + + integer(i_def), allocatable :: tile_size(:,:) integer(i_def), parameter :: one_layer = 1_i_def integer(i_def) :: i @@ -97,10 +104,22 @@ subroutine initialise( program_name, modeldb) !======================================================================= prime_mesh_name = modeldb%config%base_mesh%prime_mesh_name() geometry = modeldb%config%base_mesh%geometry() + topology = modeldb%config%base_mesh%topology() method = modeldb%config%extrusion%method() domain_height = modeldb%config%extrusion%domain_height() number_of_layers = modeldb%config%extrusion%number_of_layers() scaled_radius = modeldb%config%planet%scaled_radius() + prepartitioned = modeldb%config%base_mesh%prepartitioned() + + if (prepartitioned) then + tile_size_x = 1 + tile_size_y = 1 + inner_halo_tiles = .false. + else + tile_size_x = maxval([1,modeldb%config%partitioning%tile_size_x()]) + tile_size_y = maxval([1,modeldb%config%partitioning%tile_size_y()]) + inner_halo_tiles = modeldb%config%partitioning%inner_halo_tiles() + end if !======================================================================= ! 1.0 Mesh @@ -149,10 +168,15 @@ subroutine initialise( program_name, modeldb) ! --------------------------------------------------------- stencil_depth = 1 check_partitions = .false. + allocate(tile_size(2,size(base_mesh_names))) + tile_size(1,:) = tile_size_x + tile_size(2,:) = tile_size_y + call init_mesh( modeldb%config, & modeldb%mpi%get_comm_rank(), & modeldb%mpi%get_comm_size(), & base_mesh_names, extrusion, & + inner_halo_tiles, tile_size, & stencil_depth, check_partitions ) allocate( twod_names, source=base_mesh_names ) @@ -160,21 +184,23 @@ subroutine initialise( program_name, modeldb) twod_names(i) = trim(twod_names(i))//'_2d' end do call create_mesh( base_mesh_names, extrusion_2d, & + inner_halo_tiles, tile_size, & alt_name=twod_names ) call assign_mesh_maps(twod_names) !======================================================================= ! 2.0 Build the FEM function spaces and coordinate fields !======================================================================= - call init_fem( mesh_collection, chi_inventory, panel_id_inventory ) + call init_fem( modeldb%config, chi_inventory, panel_id_inventory ) !======================================================================= ! 3.0 Setup I/O system. !======================================================================= ! Initialise I/O context - call init_io( program_name, prime_mesh_name, & - modeldb, chi_inventory, panel_id_inventory ) + call init_io( program_name, prime_mesh_name, modeldb, & + chi_inventory, panel_id_inventory, & + geometry, topology ) !======================================================================= diff --git a/applications/simple_diffusion/source/simple_diffusion.f90 b/applications/simple_diffusion/source/simple_diffusion.f90 index b990fe6cf..d6d01842f 100644 --- a/applications/simple_diffusion/source/simple_diffusion.f90 +++ b/applications/simple_diffusion/source/simple_diffusion.f90 @@ -45,7 +45,9 @@ program simple_diffusion simple_diffusion_required_namelists, & config=modeldb%config ) - call init_logger( modeldb%mpi%get_comm(), program_name ) + call init_logger( modeldb%config, & + modeldb%mpi%get_comm(), & + program_name ) write(log_scratch_space,& '("Application built with ", A, "-bit real numbers")') & diff --git a/applications/skeleton/example/configuration.nml b/applications/skeleton/example/configuration.nml index 1f1175e64..2041e3701 100644 --- a/applications/skeleton/example/configuration.nml +++ b/applications/skeleton/example/configuration.nml @@ -58,6 +58,9 @@ coord_system='native' &partitioning partitioner='cubedsphere', panel_decomposition = 'auto', + tile_size_x=1, + tile_size_y=1, + inner_halo_tiles=.false. / &planet diff --git a/applications/skeleton/source/driver/skeleton_driver_mod.f90 b/applications/skeleton/source/driver/skeleton_driver_mod.f90 index 7334ab132..f2fd62ad9 100644 --- a/applications/skeleton/source/driver/skeleton_driver_mod.f90 +++ b/applications/skeleton/source/driver/skeleton_driver_mod.f90 @@ -79,11 +79,17 @@ subroutine initialise(program_name, modeldb) integer(i_def) :: geometry integer(i_def) :: method integer(i_def) :: number_of_layers + integer(i_def) :: tile_size_x + integer(i_def) :: tile_size_y real(r_def) :: domain_bottom real(r_def) :: domain_height real(r_def) :: scaled_radius - logical :: check_partitions + logical :: check_partitions + logical :: inner_halo_tiles + logical :: prepartitioned + + integer(i_def), allocatable :: tile_size(:,:) integer(i_def) :: i integer(i_def), parameter :: one_layer = 1_i_def @@ -102,6 +108,17 @@ subroutine initialise(program_name, modeldb) domain_height = modeldb%config%extrusion%domain_height() number_of_layers = modeldb%config%extrusion%number_of_layers() scaled_radius = modeldb%config%planet%scaled_radius() + prepartitioned = modeldb%config%base_mesh%prepartitioned() + + if (prepartitioned) then + tile_size_x = 1 + tile_size_y = 1 + inner_halo_tiles = .false. + else + tile_size_x = maxval([1,modeldb%config%partitioning%tile_size_x()]) + tile_size_y = maxval([1,modeldb%config%partitioning%tile_size_y()]) + inner_halo_tiles = modeldb%config%partitioning%inner_halo_tiles() + end if !======================================================================= ! Mesh @@ -141,17 +158,24 @@ subroutine initialise(program_name, modeldb) !----------------------------------------------------------------------- stencil_depth = 1 check_partitions = .false. + allocate(tile_size(2,size(base_mesh_names))) + tile_size(1,:) = tile_size_x + tile_size(2,:) = tile_size_y + call init_mesh( modeldb%config, & modeldb%mpi%get_comm_rank(), & modeldb%mpi%get_comm_size(), & base_mesh_names, extrusion, & + inner_halo_tiles, tile_size, & stencil_depth, check_partitions ) allocate( twod_names, source=base_mesh_names ) do i=1, size(twod_names) twod_names(i) = trim(twod_names(i))//'_2d' end do + call create_mesh( base_mesh_names, extrusion_2d, & + inner_halo_tiles, tile_size, & alt_name=twod_names ) call assign_mesh_maps( twod_names ) @@ -160,7 +184,7 @@ subroutine initialise(program_name, modeldb) ! Build the FEM function spaces and coordinate fields !======================================================================= ! Create FEM specifics (function spaces and chi field) - call init_fem(mesh_collection, chi_inventory, panel_id_inventory) + call init_fem(modeldb%config, chi_inventory, panel_id_inventory) !======================================================================= ! Create and initialise prognostic fields diff --git a/applications/skeleton/source/skeleton.f90 b/applications/skeleton/source/skeleton.f90 index cd35a0ee8..03da201ab 100644 --- a/applications/skeleton/source/skeleton.f90 +++ b/applications/skeleton/source/skeleton.f90 @@ -45,7 +45,9 @@ program skeleton call init_config( filename, skeleton_required_namelists, & config=modeldb%config ) - call init_logger( modeldb%mpi%get_comm(), program_name ) + call init_logger( modeldb%config, & + modeldb%mpi%get_comm(), & + program_name ) write(log_scratch_space,'(A)') & 'Application built with '// trim(precision_real) // & diff --git a/components/driver/rose-meta/lfric-driver/HEAD/rose-meta.conf b/components/driver/rose-meta/lfric-driver/HEAD/rose-meta.conf index 980ecd18b..4e1fa8b1c 100644 --- a/components/driver/rose-meta/lfric-driver/HEAD/rose-meta.conf +++ b/components/driver/rose-meta/lfric-driver/HEAD/rose-meta.conf @@ -52,7 +52,6 @@ type=integer # PRIMARY GLOBAL MESH #============================================================================== [namelist:base_mesh] -compulsory=true description=Provides information to define the LFRic infrastructure principle mesh. help=Lfric must use at least one mesh (prime) to run. =This panel specifies details of the mesh and its @@ -186,7 +185,6 @@ values='fully_periodic', 'non_periodic' # 2D MESH EXTRUSION #============================================================================== [namelist:extrusion] -compulsory=true description=Settings for the selected vertical mesh extrusion method. help=Settings for the uniform, quadratic, geometric and DCMIP mesh extrusion =profiles to extrude 2D to 3D mesh using non-dimensional vertical coordinate. @@ -305,7 +303,6 @@ values='linear', 'smooth' # FINITE ELEMENT #============================================================================== [namelist:finite_element] -compulsory=true description=Settings to define the choice of finite elements used help=Settings to define which finite elements create the function spaces used =in the model @@ -416,7 +413,6 @@ type=logical # IO #============================================================================== [namelist:io] -compulsory=true description=Sets I/O options for diagnostic output, checkpointing and dumps help=?????? ns=namelist/Job/IO @@ -568,7 +564,6 @@ type=logical # SYSTEM LOGGING #============================================================================== [namelist:logging] -compulsory=true ns=namelist/Job/IO/System [namelist:logging=log_to_rank_zero_only] @@ -607,7 +602,6 @@ values='error','warning','info','debug','trace' # MULTIGRID #============================================================================== [namelist:multigrid] -compulsory=true description=?????? help=?????? =?????? @@ -629,14 +623,35 @@ sort-key=Panel-A04 !string_length=default type=character +[namelist:multigrid=coarsen_multigrid_tiles] +compulsory=true +description=Reduce x and y tile sizes by a factor of 2 in each multigrid level +help=Enables using larger tiles at higher resolution levels by automatically + =reducing tile sizes in coarser levels, which can improve performance. +!kind=default +sort-key=Panel-A10 +trigger=namelist:multigrid=max_tiled_multigrid_level: .true. ; +type=logical + +[namelist:multigrid=max_tiled_multigrid_level] +compulsory=true +description=Coarsest multigrid level to be tiled +help=Revert to 1x1 tiling (equivalent to colouring) for multigrid levels + =above this threshold (level 1 has highest resolution); tiling is + =typically more beneficial for higher resolutions. +!kind=default +range=1: +sort-key=Panel-A09 +type=integer + [namelist:multigrid=multigrid_chain_nitems] compulsory=true description=Number of items in multigrid function space chain -fail-if=this < 1 ; +fail-if=this < 0 ; help=?????? =?????? !kind=default -range=1: +range=0: sort-key=Panel-A02 type=integer @@ -660,7 +675,6 @@ type=real # GLOBAL MESH PARTITIONING #============================================================================== [namelist:partitioning] -compulsory=true description=Global mesh panel partitioning. help=For parallel computing, the 2D global mesh is divided up into partitions. =Each process rank runs an instance of the model on one partition. The @@ -669,14 +683,6 @@ help=For parallel computing, the 2D global mesh is divided up into partitions. ns=namelist/Model/Mesh/Partitioning sort-key=Section-A02 -[namelist:partitioning=coarsen_multigrid_tiles] -compulsory=false -description=Reduce x and y tile sizes by a factor of 2 in each multigrid level -help=Enables using larger tiles at higher resolution levels by automatically - =reducing tile sizes in coarser levels, which can improve performance. -sort-key=Panel-A10 -type=logical - [namelist:partitioning=generate_inner_halos] compulsory=true description=Generate inner halo regions @@ -687,24 +693,17 @@ sort-key=Panel-A05 type=logical [namelist:partitioning=inner_halo_tiles] -compulsory=false +compulsory=true description=Tile inner halos separately from partition interior. help=Tiling inner halos separately from the partition interior guarantees =that tiles never cross the boundary between interior and inner halo, =which can be useful when overlapping communication and computation. +!kind=default sort-key=Panel-A08 +trigger=namelist:partitioning=tile_size_x: .true. ; + =namelist:partitioning=tile_size_y: .true. ; type=logical -[namelist:partitioning=max_tiled_multigrid_level] -compulsory=false -description=Coarsest multigrid level to be tiled -help=Revert to 1x1 tiling (equivalent to colouring) for multigrid levels - =above this threshold (level 1 has highest resolution); tiling is - =typically more beneficial for higher resolutions. -range=1: -sort-key=Panel-A09 -type=integer - [namelist:partitioning=panel_decomposition] compulsory=true description=Panel partition decomposition @@ -769,25 +768,27 @@ value-titles=Planar, Cubedsphere values='planar', 'cubedsphere' [namelist:partitioning=tile_size_x] -compulsory=false +compulsory=true description=Tile size (number of cells) in x-direction help=Tiling reorders computation of cells in the horizontal mesh to maximise =processor cache reuse. It is currently only available for partitioned =meshes and where mesh partitions have a rectangular shape. Tiles sizes =along partition borders are automatically adjusted to fit, but sizes that =are larger than partition dimensions are not accepted. +!kind=default range=1: sort-key=Panel-A06 type=integer [namelist:partitioning=tile_size_y] -compulsory=false +compulsory=true description=Tile size (number of cells) in y-direction help=Tiling reorders computation of cells in the horizontal mesh to maximise =processor cache reuse. It is currently only available for partitioned =meshes and where mesh partitions have a rectangular shape. Tiles sizes =along partition borders are automatically adjusted to fit, but sizes that =are larger than partition dimensions are not accepted. +!kind=default range=1: sort-key=Panel-A07 type=integer @@ -796,7 +797,6 @@ type=integer # PLANET #============================================================================== [namelist:planet] -compulsory=true description=?????? help=?????? =?????? @@ -824,7 +824,6 @@ type=real # TIME CONTROL #============================================================================== [namelist:time] -compulsory=true description=Time options help=At the moment, this just sets the start and end timestep for the run ns=namelist/Job/Time @@ -878,7 +877,6 @@ type=character # TIMESTEPPING #============================================================================== [namelist:timestepping] -compulsory=true description=?????? help=?????? =?????? diff --git a/components/driver/rose-meta/lfric-driver/versions.py b/components/driver/rose-meta/lfric-driver/versions.py index 01798ad2b..62eead4b2 100644 --- a/components/driver/rose-meta/lfric-driver/versions.py +++ b/components/driver/rose-meta/lfric-driver/versions.py @@ -1,3 +1,4 @@ +import re import sys from metomi.rose.upgrade import MacroUpgrade # noqa: F401 @@ -20,14 +21,44 @@ def __repr__(self): """ Copy this template and complete to add your macro - class vnXX_txxx(MacroUpgrade): # Upgrade macro for by - BEFORE_TAG = "vnX.X" AFTER_TAG = "vnX.X_txxx" - def upgrade(self, config, meta_config=None): # Add settings return config, self.reports """ + + +class vn31_t324(MacroUpgrade): + """Upgrade macro for LFRic Core PR#324 by Ricky Wong.""" + + BEFORE_TAG = "vn3.1" + AFTER_TAG = "vn3.1_t324" + + def upgrade(self, config, meta_config=None): + # Commands From: rose-meta/lfric-driver + # Only add in new configuration settings if the namelists + # are already present + # + if config.get(["namelist:partitioning"]) is not None: + self.add_setting( + config, ["namelist:partitioning", "inner_halo_tiles"], ".false." + ) + self.add_setting( + config, ["namelist:partitioning", "tile_size_x"], "1" + ) + self.add_setting( + config, ["namelist:partitioning", "tile_size_y"], "1" + ) + if config.get(["namelist:multigrid"]) is not None: + self.add_setting( + config, + ["namelist:multigrid", "coarsen_multigrid_tiles"], ".false." + ) + self.add_setting( + config, ["namelist:multigrid", "max_tiled_multigrid_level"], "1" + ) + + return config, self.reports diff --git a/components/driver/rose-meta/lfric-driver/vn3.1/rose-meta.conf b/components/driver/rose-meta/lfric-driver/vn3.1/rose-meta.conf index 980ecd18b..754acdcb8 100644 --- a/components/driver/rose-meta/lfric-driver/vn3.1/rose-meta.conf +++ b/components/driver/rose-meta/lfric-driver/vn3.1/rose-meta.conf @@ -607,7 +607,6 @@ values='error','warning','info','debug','trace' # MULTIGRID #============================================================================== [namelist:multigrid] -compulsory=true description=?????? help=?????? =?????? diff --git a/components/driver/source/driver_coordinates_mod.F90 b/components/driver/source/driver_coordinates_mod.F90 index 686b196f2..04983b088 100644 --- a/components/driver/source/driver_coordinates_mod.F90 +++ b/components/driver/source/driver_coordinates_mod.F90 @@ -7,24 +7,22 @@ !> @brief Module to assign the values of the coordinates of the mesh to a field. module driver_coordinates_mod - use base_mesh_config_mod, only: geometry, & - geometry_planar, & - geometry_spherical, & - topology, & - topology_fully_periodic, & - topology_non_periodic - use constants_mod, only: r_def, i_def, l_def, & - radians_to_degrees, & - i_halo_index, eps, pi - use log_mod, only: log_event, log_scratch_space, & - log_level_error - use planet_config_mod, only: scaled_radius - use coord_transform_mod, only: xyz2llr, llr2xyz, identify_panel, & - xyz2alphabetar, alphabetar2xyz, & - schmidt_transform_xyz, & - inverse_schmidt_transform_xyz - use finite_element_config_mod, only: coord_system, & - coord_system_xyz + use constants_mod, only: r_def, i_def, l_def, & + radians_to_degrees, & + i_halo_index, eps, pi + use log_mod, only: log_event, log_scratch_space, & + log_level_error + use coord_transform_mod, only: xyz2llr, llr2xyz, identify_panel, & + xyz2alphabetar, alphabetar2xyz, & + schmidt_transform_xyz, & + inverse_schmidt_transform_xyz + + ! Configuration modules + use base_mesh_config_mod, only: geometry_planar, & + geometry_spherical, & + topology_fully_periodic, & + topology_non_periodic + use finite_element_config_mod, only: coord_system_xyz implicit none @@ -56,7 +54,13 @@ module driver_coordinates_mod !> @param[in,out] chi Model coordinate array of size 3 of fields !> @param[in] panel_id Field giving the ID of mesh panels !> @param[in] mesh Mesh on which this field is attached - subroutine assign_coordinate_field(chi, panel_id, mesh) + !> @param[in] geometry Mesh geometry enumeration value + !> @param[in] topology Mesh topology enumeration value + !> @param[in] coord_system Finite-element coordinate syatem enumeration value + !> @param[in] scaled_radius Scaled planet radius + subroutine assign_coordinate_field(chi, panel_id, mesh, & + geometry, topology, & + coord_system, scaled_radius ) use domain_mod, only: domain_type use field_mod, only: field_type, field_proxy_type @@ -73,10 +77,15 @@ subroutine assign_coordinate_field(chi, panel_id, mesh) type( field_type ), intent( inout ) :: panel_id type( mesh_type ), intent( in ), pointer :: mesh - integer(i_def), pointer :: map(:,:) => null() - integer(i_def), pointer :: map_pid(:,:) => null() - real(kind=r_def), pointer :: dof_coords(:,:) => null() - class(reference_element_type), pointer :: reference_element => null() + integer(i_def), intent(in) :: geometry + integer(i_def), intent(in) :: topology + integer(i_def), intent(in) :: coord_system + real(r_def), intent(in) :: scaled_radius + + integer(i_def), pointer :: map(:,:) + integer(i_def), pointer :: map_pid(:,:) + real(kind=r_def), pointer :: dof_coords(:,:) + class(reference_element_type), pointer :: reference_element type(field_proxy_type) :: chi_proxy(3) type(field_proxy_type) :: panel_id_proxy @@ -86,7 +95,7 @@ subroutine assign_coordinate_field(chi, panel_id, mesh) real(kind=r_def) :: domain_min_y real(kind=r_def), allocatable :: column_coords(:,:,:) - real(kind=r_def), allocatable :: dz(:) ! dz(nlayers) array + real(kind=r_def), allocatable :: dz(:) real(kind=r_def), allocatable :: vertex_coords(:,:) integer(i_def) :: cell @@ -106,6 +115,8 @@ subroutine assign_coordinate_field(chi, panel_id, mesh) real(kind=r_def) :: inverse_rot_matrix(3,3) real(kind=r_def) :: stretch_factor + nullify( map, map_pid, dof_coords, reference_element ) + ! Break encapsulation and get the proxy. chi_proxy(1) = chi(1)%get_proxy() chi_proxy(2) = chi(2)%get_proxy() @@ -176,12 +187,14 @@ subroutine assign_coordinate_field(chi, panel_id, mesh) do cell = 1,chi_proxy(1)%vspace%get_ncell() - call calc_panel_id( nlayers_pid, & - ndf_pid, undf_pid, & - map_pid(:,cell), & - panel_id_proxy%data, & - global_dof_id, & - panel_ncells ) + call calc_panel_id( nlayers_pid, & + ndf_pid, undf_pid, & + map_pid(:,cell), & + panel_id_proxy%data, & + geometry, & + topology, & + global_dof_id, & + panel_ncells ) call mesh%get_column_coords(cell,column_coords) @@ -200,6 +213,9 @@ subroutine assign_coordinate_field(chi, panel_id, mesh) domain_max_x, & domain_min_y, & panel_id_proxy%data, & + geometry, & + topology, & + scaled_radius, & ndf_pid, & undf_pid, & map_pid(:,cell) ) @@ -210,32 +226,34 @@ subroutine assign_coordinate_field(chi, panel_id, mesh) do cell = 1,chi_proxy(1)%vspace%get_ncell() - call calc_panel_id( nlayers_pid, & - ndf_pid, undf_pid, & - map_pid(:,cell), & - panel_id_proxy%data, & - global_dof_id, & - panel_ncells ) + call calc_panel_id( nlayers_pid, & + ndf_pid, undf_pid, & + map_pid(:,cell), & + panel_id_proxy%data, & + geometry, topology, & + global_dof_id, & + panel_ncells ) call mesh%get_column_coords(cell,column_coords) - call assign_coordinate_lonlatz( nlayers, & - ndf, & - nverts, & - undf, & - map(:,cell), & - chi_proxy(1)%data, & - chi_proxy(2)%data, & - chi_proxy(3)%data, & - column_coords, & - dof_coords, & - vertex_coords, & - to_rotate, & - inverse_rot_matrix, & - panel_id_proxy%data, & - ndf_pid, & - undf_pid, & - map_pid(:,cell) ) + call assign_coordinate_lonlatz( nlayers, & + ndf, & + nverts, & + undf, & + map(:,cell), & + chi_proxy(1)%data, & + chi_proxy(2)%data, & + chi_proxy(3)%data, & + column_coords, & + dof_coords, & + vertex_coords, & + to_rotate, & + inverse_rot_matrix, & + panel_id_proxy%data, & + scaled_radius, & + ndf_pid, & + undf_pid, & + map_pid(:,cell) ) end do else if ( geometry == geometry_spherical .and. & @@ -243,33 +261,35 @@ subroutine assign_coordinate_field(chi, panel_id, mesh) do cell = 1,chi_proxy(1)%vspace%get_ncell() - call calc_panel_id( nlayers_pid, & - ndf_pid, undf_pid, & - map_pid(:,cell), & - panel_id_proxy%data, & - global_dof_id, & - panel_ncells ) + call calc_panel_id( nlayers_pid, & + ndf_pid, undf_pid, & + map_pid(:,cell), & + panel_id_proxy%data, & + geometry, topology, & + global_dof_id, & + panel_ncells ) call mesh%get_column_coords(cell,column_coords) - call assign_coordinate_alphabetaz( nlayers, & - ndf, & - nverts, & - undf, & - map(:,cell), & - chi_proxy(1)%data, & - chi_proxy(2)%data, & - chi_proxy(3)%data, & - column_coords, & - dof_coords, & - vertex_coords, & - to_rotate, & - inverse_rot_matrix, & - stretch_factor, & - panel_id_proxy%data, & - ndf_pid, & - undf_pid, & - map_pid(:,cell) ) + call assign_coordinate_alphabetaz( nlayers, & + ndf, & + nverts, & + undf, & + map(:,cell), & + chi_proxy(1)%data, & + chi_proxy(2)%data, & + chi_proxy(3)%data, & + column_coords, & + dof_coords, & + vertex_coords, & + to_rotate, & + inverse_rot_matrix, & + stretch_factor, & + panel_id_proxy%data, & + scaled_radius, & + ndf_pid, & + undf_pid, & + map_pid(:,cell) ) end do else @@ -294,18 +314,21 @@ end subroutine assign_coordinate_field !! be the panel IDs which are calculated from the coordinates. !! For planar geometry the ID is just 1 everywhere. !> - !> @param[in] nlayers Number of layers for the panel_id field - !> @param[in] ndf_pid Number of DoFs per cell for the panel_id field - !> @param[in] undf_pid Universal number of DoFs for the panel_id field - !> @param[in] map_pid DoF map for the panel_id field - !> @param[out] panel_id Field (to be calculated) with the ID of cubed sphere panels - !> @param[in] global_dof_id Array of global id's - !> @param[in] panel_ncells Number of cells per cubed sphere panel + !> @param[in] nlayers Number of layers for the panel_id field + !> @param[in] ndf_pid Number of DoFs per cell for the panel_id field + !> @param[in] undf_pid Universal number of DoFs for the panel_id field + !> @param[in] map_pid DoF map for the panel_id field + !> @param[out] panel_id Field (to be calculated) with the ID of cubed sphere panels + !> @param[in] geometry Mesh geometry enumeration value + !> @param[in] topology Mesh topology enumeration value + !> @param[in] global_dof_id Array of global id's + !> @param[in] panel_ncells Number of cells per cubed sphere panel subroutine calc_panel_id( nlayers, & ndf_pid, & undf_pid, & map_pid, & panel_id, & + geometry, topology, & global_dof_id, & panel_ncells ) @@ -314,6 +337,8 @@ subroutine calc_panel_id( nlayers, & integer(kind=i_def), intent(in) :: nlayers, ndf_pid, undf_pid integer(kind=i_def), intent(in) :: map_pid(ndf_pid) real(kind=r_def), intent(out) :: panel_id(undf_pid) + integer(kind=i_def), intent(in) :: geometry + integer(kind=i_def), intent(in) :: topology integer(kind=i_def), intent(in) :: global_dof_id(undf_pid) integer(kind=i_def), intent(in) :: panel_ncells @@ -351,6 +376,9 @@ end subroutine calc_panel_id !> @param[in] domain_x Domain extent in x direction for planar mesh !> @param[in] domain_y Domain extent in y direction for planar mesh !> @param[in] panel_id Field giving IDs of mesh panels + !> @param[in] geometry Mesh geometry enumeration value + !> @param[in] topology Mesh topology enumeration value + !> @param[in] scaled_radius Scaled planet radius !> @param[in] ndf_pid Number of DoFs per cell for panel_id space !> @param[in] undf_pid Number of universal DoFs for panel_id space !> @param[in] map_pid DoF map for panel_id space @@ -369,6 +397,9 @@ subroutine assign_coordinate_xyz( nlayers, & domain_x, & domain_y, & panel_id, & + geometry, & + topology, & + scaled_radius, & ndf_pid, & undf_pid, & map_pid ) @@ -387,6 +418,8 @@ subroutine assign_coordinate_xyz( nlayers, & real(kind=r_def), intent(in) :: chi_hat_node(3,ndf), chi_hat_vert(nverts,3) real(kind=r_def), intent(in) :: domain_x, domain_y real(kind=r_def), intent(in) :: panel_id(undf_pid) + integer(i_def), intent(in) :: geometry, topology + real(r_def), intent(in) :: scaled_radius ! Internal variables integer(kind=i_def) :: k, df, dfk, vert @@ -471,6 +504,7 @@ end subroutine assign_coordinate_xyz !! Cartesian coordinates from physical ones !> @param[in] stretch_factor Stretch factor for Schmidt transform !> @param[in] panel_id Field giving IDs of mesh panels + !> @param[in] scaled_radius Scaled planet radius !> @param[in] ndf_pid Number of DoFs per cell for panel_id space !> @param[in] undf_pid Number of universal DoFs for panel_id space !> @param[in] map_pid DoF map for panel_id space @@ -489,6 +523,7 @@ subroutine assign_coordinate_alphabetaz( nlayers, & inverse_rot_matrix, & stretch_factor, & panel_id, & + scaled_radius, & ndf_pid, & undf_pid, & map_pid ) @@ -507,6 +542,7 @@ subroutine assign_coordinate_alphabetaz( nlayers, & integer(kind=i_def), intent(in) :: ndf_pid, undf_pid integer(kind=i_def), intent(in) :: map_pid(ndf_pid) real(kind=r_def), intent(in) :: panel_id(undf_pid) + real(kind=r_def), intent(in) :: scaled_radius ! Internal variables integer(kind=i_def) :: k, df, dfk, vert @@ -587,6 +623,7 @@ end subroutine assign_coordinate_alphabetaz !> @param[in] inverse_rot_matrix Rotation matrix to apply to obtain native !! Cartesian coordinates from physical ones !> @param[in] panel_id Field giving IDs of mesh panels + !> @param[in] scaled_radius Scaled planet radius !> @param[in] ndf_pid Number of DoFs per cell for panel_id space !> @param[in] undf_pid Number of universal DoFs for panel_id space !> @param[in] map_pid DoF map for panel_id space @@ -604,6 +641,7 @@ subroutine assign_coordinate_lonlatz( nlayers, & to_rotate, & inverse_rot_matrix, & panel_id, & + scaled_radius, & ndf_pid, & undf_pid, & map_pid ) @@ -621,6 +659,7 @@ subroutine assign_coordinate_lonlatz( nlayers, & integer(kind=i_def), intent(in) :: ndf_pid, undf_pid integer(kind=i_def), intent(in) :: map_pid(ndf_pid) real(kind=r_def), intent(in) :: panel_id(undf_pid) + real(kind=r_def), intent(in) :: scaled_radius ! Internal variables integer(kind=i_def) :: k, df, dfk, vert diff --git a/components/driver/source/driver_counter_mod.f90 b/components/driver/source/driver_counter_mod.f90 index c4cfb300a..21865e95f 100644 --- a/components/driver/source/driver_counter_mod.f90 +++ b/components/driver/source/driver_counter_mod.f90 @@ -7,10 +7,10 @@ !> module driver_counter_mod - use count_mod, only : count_type, halo_calls - use io_config_mod, only : subroutine_counters, & - counter_output_suffix - use timer_mod, only : timer, output_timer, init_timer + use config_mod, only: config_type + use constants_mod, only: str_max_filename, l_def + use count_mod, only: count_type, halo_calls + use timer_mod, only: timer, output_timer, init_timer implicit none @@ -24,15 +24,21 @@ module driver_counter_mod !> As well as initialising the system a "top level" counter is set up !? for tracking halo calls. !> + !> @param[in] config Configuration object !> @param[in] identifier Top level halo name. !> - subroutine init_counters( identifier ) + subroutine init_counters(config, identifier) implicit none - character(*), intent(in) :: identifier + type(config_type), intent(in) :: config + character(*), intent(in) :: identifier - if (subroutine_counters) then + logical(l_def) :: subroutine_counters + + subroutine_counters = config%io%subroutine_counters() + + if ( subroutine_counters ) then allocate( halo_calls, source=count_type('halo_calls') ) call halo_calls%counter( identifier ) end if @@ -48,15 +54,23 @@ end subroutine init_counters !> @todo Reconsider the existance of the simple counter system once the !> profiler is integrated. !> + !> @param[in] config Configuration object !> @param[in] identifier Top level counter name. !> - subroutine final_counters( identifier ) + subroutine final_counters(config, identifier) implicit none - character(*), intent(in) :: identifier + type(config_type), intent(in) :: config + character(*), intent(in) :: identifier - if ( subroutine_counters ) then + logical(l_def) :: subroutine_counters + character(str_max_filename) :: counter_output_suffix + + subroutine_counters = config%io%subroutine_counters() + counter_output_suffix = config%io%counter_output_suffix() + + if (subroutine_counters) then call halo_calls%counter( identifier ) call halo_calls%output_counters( counter_output_suffix ) end if diff --git a/components/driver/source/driver_fem_mod.f90 b/components/driver/source/driver_fem_mod.f90 index 7ba02690e..bce83645b 100644 --- a/components/driver/source/driver_fem_mod.f90 +++ b/components/driver/source/driver_fem_mod.f90 @@ -12,11 +12,9 @@ !> * Initialises function space chains for use by the model. module driver_fem_mod - use sci_chi_transform_mod, only: init_chi_transforms, & - final_chi_transforms - use constants_mod, only: i_def, l_def, str_def + use config_mod, only: config_type + use constants_mod, only: i_def, r_def, l_def, str_def, imdi use extrusion_mod, only: TWOD, PRIME_EXTRUSION - use finite_element_config_mod, only: coord_order use field_mod, only: field_type use fs_continuity_mod, only: W0, W2, W3, Wtheta, Wchi, W2v, W2h use function_space_mod, only: function_space_type @@ -35,9 +33,10 @@ module driver_fem_mod LOG_LEVEL_ERROR, & log_scratch_space use mesh_mod, only: mesh_type - use mesh_collection_mod, only: mesh_collection_type + use mesh_collection_mod, only: mesh_collection - use base_mesh_config_mod, only: geometry, topology + use sci_chi_transform_mod, only: init_chi_transforms, & + final_chi_transforms implicit none @@ -48,47 +47,64 @@ module driver_fem_mod !> @brief Initialises the coordinate fields (chi) and FEM components. !> - !> @param[in] mesh_collection Collection of all meshes to set up - !! coordinates for + !> @param[in] config Configuration object !> @param[in,out] chi_inventory Inventory object, containing all of !! the chi fields indexed by mesh !> @param[in,out] panel_id_inventory Inventory object, containing all of !! the fields with the ID of mesh panels - subroutine init_fem( mesh_collection, chi_inventory, panel_id_inventory ) + subroutine init_fem(config, chi_inventory, panel_id_inventory) implicit none ! Coordinate field - type(mesh_collection_type), intent(in) :: mesh_collection - type(inventory_by_mesh_type), intent(inout) :: chi_inventory - type(inventory_by_mesh_type), intent(inout) :: panel_id_inventory + type(config_type), intent(in) :: config + + type(inventory_by_mesh_type), intent(inout) :: chi_inventory + type(inventory_by_mesh_type), intent(inout) :: panel_id_inventory character(str_def), allocatable :: all_mesh_names(:) - type(mesh_type), pointer :: mesh => null() - type(mesh_type), pointer :: twod_mesh => null() + type(mesh_type), pointer :: mesh + type(mesh_type), pointer :: twod_mesh type(field_type) :: chi(3) type(field_type) :: panel_id - type(function_space_type), pointer :: fs => null() - integer(kind=i_def) :: chi_space, coord, i + type(function_space_type), pointer :: fs + integer(i_def) :: chi_space, coord, i character(str_def) :: mesh_name + integer(i_def) :: coord_order, geometry, topology, coord_system + real(r_def) :: scaled_radius + + call log_event( 'FEM specifics: creating function spaces...', & + log_level_info ) + + nullify(mesh, twod_mesh, fs) + + if (config%namelist_exists('base_mesh')) then + geometry = config%base_mesh%geometry() + topology = config%base_mesh%topology() + else + geometry = imdi + topology = imdi + end if - call log_event( 'FEM specifics: creating function spaces...', log_level_info ) + coord_system = config%finite_element%coord_system() + coord_order = config%finite_element%coord_order() + scaled_radius = config%planet%scaled_radius() ! ======================================================================== ! ! Initialise coordinates ! ======================================================================== ! ! Initialise coordinate transformations - call init_chi_transforms( geometry, topology, & - mesh_collection=mesh_collection ) + call init_chi_transforms( geometry, topology, mesh_collection=mesh_collection ) ! To loop through mesh collection, get all mesh names ! Then get mesh from collection using these names all_mesh_names = mesh_collection%get_mesh_names() - call chi_inventory%initialise(name="chi", table_len=SIZE(all_mesh_names)) - call panel_id_inventory%initialise(name="panel_id", table_len=SIZE(all_mesh_names)) + call chi_inventory%initialise(name="chi", table_len=size(all_mesh_names)) + call panel_id_inventory%initialise(name="panel_id", & + table_len=size(all_mesh_names)) ! ======================================================================== ! ! Loop through all 3D meshes @@ -104,28 +120,35 @@ subroutine init_fem( mesh_collection, chi_inventory, panel_id_inventory ) ! Initialise panel ID field object --------------------------------------- twod_mesh => mesh_collection%get_mesh(mesh, TWOD) fs => function_space_collection%get_fs(twod_mesh, 0, 0, W3) - call panel_id%initialise( vector_space = fs, halo_depth = twod_mesh%get_halo_depth() ) + call panel_id%initialise( vector_space=fs, & + halo_depth=twod_mesh%get_halo_depth() ) ! Initialise chi field object -------------------------------------------- if ( coord_order == 0 ) then chi_space = W0 write(log_scratch_space,'(A)') & - 'Computing W0 coordinate fields for ' // trim(mesh_name) // 'mesh' + 'Computing W0 coordinate fields for ' // & + trim(mesh_name) // 'mesh' call log_event( log_scratch_space, log_level_info ) else chi_space = Wchi write(log_scratch_space,'(A)') & - 'Computing Wchi coordinate fields for ' // trim(mesh_name) // 'mesh' + 'Computing Wchi coordinate fields for ' // & + trim(mesh_name) // 'mesh' call log_event( log_scratch_space, log_level_info ) end if - fs => function_space_collection%get_fs(mesh, coord_order, coord_order, chi_space) + fs => function_space_collection%get_fs( mesh, coord_order, & + coord_order, chi_space ) do coord = 1, size(chi) - call chi(coord)%initialise(vector_space = fs, halo_depth = twod_mesh%get_halo_depth() ) + call chi(coord)%initialise( vector_space=fs, & + halo_depth=twod_mesh%get_halo_depth() ) end do ! Set coordinate fields -------------------------------------------------- - call assign_coordinate_field(chi, panel_id, mesh) + call assign_coordinate_field( chi, panel_id, mesh, & + geometry, topology, & + coord_system, scaled_radius ) ! Add fields to inventory call chi_inventory%copy_field_array(chi, mesh) @@ -140,22 +163,24 @@ subroutine init_fem( mesh_collection, chi_inventory, panel_id_inventory ) end subroutine init_fem !> @brief Initialises the function space chains used in multigrid. - !> @param[in] mesh_collection Collection of all meshes to set up - !! coordinates for - !> @param[in] multigrid_mesh_names Names of the multigrid meshes - subroutine init_function_space_chains( mesh_collection, multigrid_mesh_names ) + !> @param[in] multigrid_mesh_names Names of the multigrid meshes + subroutine init_function_space_chains(multigrid_mesh_names) implicit none - type(mesh_collection_type), intent(in) :: mesh_collection - character(str_def), intent(in) :: multigrid_mesh_names(:) + character(str_def), intent(in) :: multigrid_mesh_names(:) - type(mesh_type), pointer :: mesh => null() - type(mesh_type), pointer :: twod_mesh => null() - type(function_space_type), pointer :: fs => null() - integer(kind=i_def) :: i + type(mesh_type), pointer :: mesh + type(mesh_type), pointer :: twod_mesh + + type(function_space_type), pointer :: fs + + integer(i_def) :: i + + nullify(mesh, twod_mesh, fs) - call log_event( 'FEM specifics: creating function space chains...', LOG_LEVEL_INFO ) + call log_event( 'FEM specifics: creating function space chains...', & + log_level_info ) ! ======================================================================== ! ! Create function space chains @@ -167,7 +192,7 @@ subroutine init_function_space_chains( mesh_collection, multigrid_mesh_names ) w2h_multigrid_function_space_chain = function_space_chain_type() wtheta_multigrid_function_space_chain = function_space_chain_type() - write(log_scratch_space,'(A,I1,A)') & + write(log_scratch_space,'(A,I1,A)') & 'Initialising MultiGrid ', size(multigrid_mesh_names), & '-level function space chain.' call log_event( log_scratch_space, LOG_LEVEL_INFO ) diff --git a/components/driver/source/driver_io_mod.F90 b/components/driver/source/driver_io_mod.F90 index 2beadda81..9fd4252aa 100644 --- a/components/driver/source/driver_io_mod.F90 +++ b/components/driver/source/driver_io_mod.F90 @@ -10,7 +10,7 @@ !> module driver_io_mod - use constants_mod, only: str_def, i_def, l_def + use constants_mod, only: str_def, i_def, l_def, r_def use driver_modeldb_mod, only: modeldb_type use driver_model_data_mod, only: model_data_type use empty_io_context_mod, only: empty_io_context_type @@ -56,6 +56,8 @@ end subroutine filelist_populator !> @param[inout] modeldb Model state !> @param[in] chi_inventory Contains the model's coordinate fields !> @param[in] panel_id_inventory Contains the model's panel ID fields + !> @param[in] geometry Mesh geometry enumeration value + !> @param[in] topology Mesh topology enumeration value !> @param[in] populate_filelist Optional procedure for creating a list of !! file descriptions used by the model I/O @@ -68,17 +70,22 @@ subroutine init_io( context_name, & modeldb, & chi_inventory, & panel_id_inventory, & + geometry, topology, & populate_filelist, & alt_mesh_names, & before_close ) + implicit none - character(*), intent(in) :: context_name - character(*), intent(in) :: mesh_name - class(modeldb_type), intent(inout) :: modeldb - type(inventory_by_mesh_type), intent(in) :: chi_inventory - type(inventory_by_mesh_type), intent(in) :: panel_id_inventory + character(*), intent(in) :: context_name + character(*), intent(in) :: mesh_name + class(modeldb_type), intent(inout) :: modeldb + type(inventory_by_mesh_type), intent(in) :: chi_inventory + type(inventory_by_mesh_type), intent(in) :: panel_id_inventory + integer(i_def), intent(in) :: geometry + integer(i_def), intent(in) :: topology + procedure(filelist_populator), & pointer, optional, intent(in) :: populate_filelist character(len=str_def), optional, intent(in) :: alt_mesh_names(:) @@ -105,6 +112,7 @@ subroutine init_io( context_name, & chi_inventory, & panel_id_inventory, & before_close_ptr, & + geometry, topology, & populate_filelist, & alt_mesh_names ) #else @@ -155,6 +163,8 @@ end subroutine init_empty_io_context !> @param[in] chi_inventory Contains the model's coordinate fields !> @param[in] panel_id_inventory Contains the model's panel ID fields !> @param[in] before_close Routine to be called before context closes + !> @param[in] geometry Mesh geometry enumeration value + !> @param[in] topology Mesh topology enumeration value !> @param[in] populate_filelist Optional procedure for creating a list of !! file descriptions used by the model I/O !> @param[in] alt_mesh_names Optional array of names for other meshes @@ -165,17 +175,21 @@ subroutine init_xios_io_context( context_name, & chi_inventory, & panel_id_inventory, & before_close, & + geometry, topology, & populate_filelist, & alt_mesh_names ) implicit none - character(*), intent(in) :: context_name - character(*), intent(in) :: mesh_name - class(modeldb_type), intent(inout) :: modeldb - type(inventory_by_mesh_type), intent(in) :: chi_inventory - type(inventory_by_mesh_type), intent(in) :: panel_id_inventory - procedure(callback_clock_arg), pointer, intent(in) :: before_close + character(*), intent(in) :: context_name + character(*), intent(in) :: mesh_name + class(modeldb_type), intent(inout) :: modeldb + type(inventory_by_mesh_type), intent(in) :: chi_inventory + type(inventory_by_mesh_type), intent(in) :: panel_id_inventory + procedure(callback_clock_arg), intent(in), pointer :: before_close + integer(i_def), intent(in) :: geometry + integer(i_def), intent(in) :: topology + procedure(filelist_populator), & pointer, optional, intent(in) :: populate_filelist character(len=str_def), optional, intent(in) :: alt_mesh_names(:) @@ -197,6 +211,9 @@ subroutine init_xios_io_context( context_name, & integer(i_def) :: num_meshes, i, j + integer(i_def) :: coord_system + real(r_def) :: scaled_radius + mesh => null() chi => null() panel_id => null() @@ -207,6 +224,9 @@ subroutine init_xios_io_context( context_name, & !============================================================== + coord_system = modeldb%config%finite_element%coord_system() + scaled_radius = modeldb%config%planet%scaled_radius() + call tmp_io_context%initialise(context_name) call modeldb%io_contexts%add_context(tmp_io_context) call modeldb%io_contexts%get_io_context(context_name, io_context) @@ -256,6 +276,9 @@ subroutine init_xios_io_context( context_name, & modeldb%clock, & modeldb%calendar, & before_close, & + geometry, topology, & + coord_system, & + scaled_radius, & alt_coords, & alt_panel_ids ) deallocate(alt_coords) @@ -265,7 +288,10 @@ subroutine init_xios_io_context( context_name, & chi, panel_id, & modeldb%clock, & modeldb%calendar, & - before_close ) + before_close, & + geometry, topology, & + coord_system, & + scaled_radius ) end if ! Attach context advancement to the model's clock diff --git a/components/driver/source/driver_log_mod.f90 b/components/driver/source/driver_log_mod.f90 index 28f416313..4c8187b1c 100644 --- a/components/driver/source/driver_log_mod.f90 +++ b/components/driver/source/driver_log_mod.f90 @@ -1,7 +1,8 @@ module driver_log_mod -use constants_mod, only: i_def +use constants_mod, only: i_def, l_def use convert_to_upper_mod, only: convert_to_upper +use config_mod, only: config_type use lfric_mpi_mod, only: lfric_comm_type use log_mod, only: log_event, & log_set_level, & @@ -14,14 +15,14 @@ module driver_log_mod LOG_LEVEL_INFO, & LOG_LEVEL_DEBUG, & LOG_LEVEL_TRACE -use logging_config_mod, only: run_log_level, & - key_from_run_log_level, & - RUN_LOG_LEVEL_ERROR, & - RUN_LOG_LEVEL_INFO, & - RUN_LOG_LEVEL_DEBUG, & - RUN_LOG_LEVEL_TRACE, & - RUN_LOG_LEVEL_WARNING, & - log_to_rank_zero_only + +use logging_config_mod, only: key_from_run_log_level, & + RUN_LOG_LEVEL_ERROR, & + RUN_LOG_LEVEL_INFO, & + RUN_LOG_LEVEL_DEBUG, & + RUN_LOG_LEVEL_TRACE, & + RUN_LOG_LEVEL_WARNING + implicit none @@ -32,17 +33,24 @@ module driver_log_mod !> @brief Initialises the logging system from a namelist. !> +!> @param[in] config Application configuration object. !> @param[in] communicator MPI communicator to use for logging. !> @param[in] program_name Identifies the running program. !> -subroutine init_logger(communicator, program_name) +subroutine init_logger(config, communicator, program_name) implicit none - character(len=*), intent(in) :: program_name - type(lfric_comm_type), intent(in) :: communicator + type(config_type), intent(in) :: config + type(lfric_comm_type), intent(in) :: communicator + character(len=*), intent(in) :: program_name integer(i_def) :: log_level + integer(i_def) :: run_log_level + logical(l_def) :: log_to_rank_zero_only + + run_log_level = config%logging%run_log_level() + log_to_rank_zero_only = config%logging%log_to_rank_zero_only() call initialise_logging( communicator%get_comm_mpi_val(), program_name, & log_to_rank_zero_only=log_to_rank_zero_only) diff --git a/components/driver/source/driver_mesh_mod.f90 b/components/driver/source/driver_mesh_mod.f90 index 79adb4407..40157459a 100644 --- a/components/driver/source/driver_mesh_mod.f90 +++ b/components/driver/source/driver_mesh_mod.f90 @@ -79,6 +79,8 @@ module driver_mesh_mod !> @param[in] total_ranks Total number of MPI ranks in this job. !> @param[in] mesh_names Mesh names to load from the mesh input file(s). !> @param[in] extrusion Extrusion object to be applied to meshes. +!> @param[in] inner_halo_tiles Apply tiling to inner halos. +!> @param[in] tile_size Tile sizes to apply to inner halos if applicable. !> @param[in] stencil_depths_in Required stencil depth for each mesh for !! the application. If this array is of size 1 then !! the single value is applied to all meshes. @@ -94,6 +96,8 @@ module driver_mesh_mod subroutine init_mesh( config, & local_rank, total_ranks, & mesh_names, extrusion, & + inner_halo_tiles, & + tile_size, & stencil_depths_in, & check_partitions, & alt_names ) @@ -107,8 +111,10 @@ subroutine init_mesh( config, & character(str_def), intent(in) :: mesh_names(:) class(extrusion_type), intent(in) :: extrusion - integer(i_def), intent(in) :: stencil_depths_in(:) - logical(l_def), intent(in) :: check_partitions + logical(l_def), intent(in) :: inner_halo_tiles + integer(i_def), intent(in) :: tile_size(:,:) + integer(i_def), intent(in) :: stencil_depths_in(:) + logical(l_def), intent(in) :: check_partitions character(str_def), optional, intent(in) :: alt_names(:) @@ -137,9 +143,10 @@ subroutine init_mesh( config, & class(panel_decomposition_type), allocatable :: decomposition - integer(i_def) :: i, n_digit character(str_def) :: fmt_str, number_str + integer(i_def) :: i, n_digit + !============================================================================ ! Extract configuration variables !============================================================================ @@ -183,7 +190,7 @@ subroutine init_mesh( config, & end do ! Currently only quad elements are fully functional - if (cellshape /= CELLSHAPE_QUADRILATERAL) then + if (cellshape /= cellshape_quadrilateral) then call log_event( "Reference_element must be QUAD for now...", & LOG_LEVEL_ERROR ) end if @@ -328,12 +335,13 @@ subroutine init_mesh( config, & end if ! prepartitioned - !============================================================================ ! 3.0 Extrude the specified meshes from local mesh objects into ! mesh objects on the given extrusion. !============================================================================ - call create_mesh( mesh_names, extrusion, alt_name=names ) + call create_mesh( mesh_names, extrusion, & + inner_halo_tiles, tile_size, & + alt_name=names ) !============================================================================ diff --git a/components/driver/source/driver_modeldb_mod.f90 b/components/driver/source/driver_modeldb_mod.f90 index d26db8254..3743bed6a 100644 --- a/components/driver/source/driver_modeldb_mod.f90 +++ b/components/driver/source/driver_modeldb_mod.f90 @@ -13,12 +13,12 @@ module driver_modeldb_mod use calendar_mod, only: calendar_type + use config_mod, only: config_type use driver_model_data_mod, only: model_data_type use key_value_collection_mod, only: key_value_collection_type use lfric_mpi_mod, only: lfric_mpi_type use model_clock_mod, only: model_clock_type use namelist_collection_mod, only: namelist_collection_type - use config_mod, only: config_type use io_context_collection_mod, only: io_context_collection_type implicit none diff --git a/components/driver/source/mesh/create_mesh_mod.f90 b/components/driver/source/mesh/create_mesh_mod.f90 index a91e6d5fb..5c4eec6ec 100644 --- a/components/driver/source/mesh/create_mesh_mod.f90 +++ b/components/driver/source/mesh/create_mesh_mod.f90 @@ -17,10 +17,7 @@ module create_mesh_mod use extrusion_mod, only: extrusion_type, & uniform_extrusion_type, & geometric_extrusion_type, & - quadratic_extrusion_type, & - PRIME_EXTRUSION, & - SHIFTED, & - DOUBLE_LEVEL + quadratic_extrusion_type use local_mesh_mod, only: local_mesh_type use mesh_mod, only: mesh_type use sci_query_mod, only: check_lbc @@ -34,14 +31,6 @@ module create_mesh_mod method_geometric, & method_quadratic - use multigrid_config_mod, only: chain_mesh_tags - - use partitioning_config_mod, only: tile_size_x, & - tile_size_y, & - inner_halo_tiles, & - max_tiled_multigrid_level, & - coarsen_multigrid_tiles - implicit none private @@ -77,15 +66,15 @@ function create_extrusion( extrusion_method, & select case (extrusion_method) case (method_uniform) - allocate( new, source=uniform_extrusion_type( & + allocate( new, source=uniform_extrusion_type( & domain_bottom, domain_height, & n_layers, extrusion_id ) ) case (method_quadratic) - allocate( new, source=quadratic_extrusion_type( & + allocate( new, source=quadratic_extrusion_type( & domain_bottom, domain_height, & n_layers, extrusion_id ) ) case (method_geometric) - allocate( new, source=geometric_extrusion_type( & + allocate( new, source=geometric_extrusion_type( & domain_bottom, domain_height, & n_layers, extrusion_id ) ) case default @@ -102,16 +91,22 @@ end function create_extrusion !! !> @param[in] local_mesh_names Names of the local_mesh_types to extrude. !> @param[in] extrusion Extrusion to employ. +!> @param[in] inner_halo_tiles +!> @param[in] tile_size !> @param[in] alt_name Optional, Alternative names for the !! extruded meshes, defaults to local_mesh_names !! if absent. -subroutine create_mesh_multiple( local_mesh_names, & - extrusion, & +subroutine create_mesh_multiple( local_mesh_names, extrusion, & + inner_halo_tiles, tile_size, & alt_name ) + implicit none character(str_def), intent(in) :: local_mesh_names(:) class(extrusion_type), intent(in) :: extrusion + logical(l_def), intent(in) :: inner_halo_tiles + integer(i_def), intent(in) :: tile_size(:,:) + character(str_def), intent(in), & optional :: alt_name(:) @@ -120,7 +115,6 @@ subroutine create_mesh_multiple( local_mesh_names, & character(str_def), allocatable :: names(:) if (present(alt_name)) then - if ( size(alt_name) /= size(local_mesh_names) ) then write(log_scratch_space, '(A)') & 'Number of alternative mesh names does not match '// & @@ -129,16 +123,13 @@ subroutine create_mesh_multiple( local_mesh_names, & end if allocate(names, source=alt_name) - else - allocate(names, source=local_mesh_names) - end if do i=1, size(local_mesh_names) - call create_mesh_single( local_mesh_names(i), & - extrusion, & + call create_mesh_single( local_mesh_names(i), extrusion, & + inner_halo_tiles, tile_size(:,i), & alt_name=names(i) ) end do @@ -154,30 +145,31 @@ end subroutine create_mesh_multiple !> @param[in] local_mesh_name Name of local_mesh_type object in !! application local_mesh_collection. !> @param[in] extrusion Extrusion to employ for this mesh_type object +!> @param[in] inner_halo_tiles +!> @param[in] tile_size !> @param[in] alt_name Optional, Alternative name for the !! extruded mesh, defaults to local_mesh_name !! if absent. -subroutine create_mesh_single( local_mesh_name, & - extrusion, & +subroutine create_mesh_single( local_mesh_name, extrusion, & + inner_halo_tiles, tile_size, & alt_name ) implicit none character(str_def), intent(in) :: local_mesh_name class(extrusion_type), intent(in) :: extrusion - character(str_def), intent(in), & - optional :: alt_name + logical(l_def), intent(in) :: inner_halo_tiles + integer(i_def), intent(in) :: tile_size(2) - type(local_mesh_type), pointer :: local_mesh_ptr => null() + character(str_def), intent(in), optional :: alt_name + + type(local_mesh_type), pointer :: local_mesh_ptr type(mesh_type) :: mesh integer(kind=i_def) :: mesh_id character(len=str_def) :: name - integer(kind=i_def) :: tile_size(2) - integer(kind=i_def) :: multigrid_level - integer(kind=i_def) :: max_multigrid_level - logical(kind=l_def) :: set_tile_size + nullify (local_mesh_ptr) if ( .not. present(alt_name) ) then name = local_mesh_name @@ -186,7 +178,7 @@ subroutine create_mesh_single( local_mesh_name, & end if - ! 1.0 Check if mesh_type already exists. + ! Check if mesh_type already exists. !=============================================== if ( mesh_collection%check_for(name) ) then write(log_scratch_space,'(A)') & @@ -197,7 +189,7 @@ subroutine create_mesh_single( local_mesh_name, & end if - ! 2.0 Extrude the local_mesh_object. + ! Extrude the local_mesh_object. !=============================================== local_mesh_ptr => local_mesh_collection%get_local_mesh(local_mesh_name) @@ -213,65 +205,15 @@ subroutine create_mesh_single( local_mesh_name, & end if end if - - ! 3.0 Set up tiling - !=============================================== - ! Set coarsest multigrid level that will be tiled; - ! restrict to the finest grid by default - max_multigrid_level = 1 - if ( max_tiled_multigrid_level /= imdi ) then - max_multigrid_level = max_tiled_multigrid_level - end if - - ! The tiling module uses 1x1 tiles (equivalent to colouring) by - ! default; allow user-specified tile sizes in case of 3D meshes - ! (PRIME_EXTRUSION, SHIFTED, and DOUBLE_LEVEL extrusions) and up to - ! the specified multigrid level (count levels until mesh name - ! includes the chain mesh tag). This relies on mesh name conventions - ! and a tag order from finest (level 1) to coarsest mesh (level n). - set_tile_size = .false. - if ( extrusion%get_id() == PRIME_EXTRUSION .or. & - extrusion%get_id() == SHIFTED .or. & - extrusion%get_id() == DOUBLE_LEVEL ) then - if ( allocated(chain_mesh_tags) ) then - ! Multigrid setup - use tiling if multigrid level is allowed, and - ! if mesh name includes the mesh tag at that level - do multigrid_level = 1, SIZE(chain_mesh_tags) - if ( index( trim(name), trim(chain_mesh_tags(multigrid_level)) ) > 0 & - .and. multigrid_level <= max_multigrid_level ) then - set_tile_size = .true. - exit - end if - end do - else - ! Not a multigrid setup - use tiling - set_tile_size = .true. - end if - end if - - ! Set user-specified tile size if tiling is allowed and adapt it to coarser - ! multigrid levels if requested and applicable - tile_size = 1 - if ( set_tile_size ) then - if ( tile_size_x /= imdi ) tile_size(1) = tile_size_x - if ( tile_size_y /= imdi ) tile_size(2) = tile_size_y - if ( coarsen_multigrid_tiles .and. allocated( chain_mesh_tags ) ) then - do multigrid_level = 1, SIZE(chain_mesh_tags) - if ( index( trim(name), & - trim(chain_mesh_tags(multigrid_level)) ) > 0 ) exit - tile_size = max( tile_size / 2, 1 ) - end do - end if - end if - mesh = mesh_type( local_mesh_ptr, extrusion, mesh_name=name, & - tile_size=tile_size, inner_halo_tiles=inner_halo_tiles ) + tile_size=tile_size, & + inner_halo_tiles=inner_halo_tiles ) mesh_id = mesh_collection%add_new_mesh( mesh ) call mesh%clear() - ! 4.0 Report on mesh_type creation. + ! Report on mesh_type creation. !=============================================== write(log_scratch_space,'(A,I0,A)') & ' ... "'//trim(name)//'"(id:', mesh_id,') '// & diff --git a/components/driver/source/mesh/multigrid_mod.f90 b/components/driver/source/mesh/multigrid_mod.f90 new file mode 100644 index 000000000..3c5cf2dc6 --- /dev/null +++ b/components/driver/source/mesh/multigrid_mod.f90 @@ -0,0 +1,120 @@ +!----------------------------------------------------------------------------- +! (C) Crown copyright Met Office. All rights reserved. +! The file LICENCE, distributed with this code, contains details of the terms +! under which the code may be used. +!----------------------------------------------------------------------------- + +module multigrid_mod + + use extrusion_mod, only: extrusion_type, prime_extrusion, & + shifted, double_level + use config_mod, only: config_type + use constants_mod, only: i_def, l_def, str_def, imdi + use log_mod, only: log_event, log_level_error + + implicit none + + public :: get_multigrid_tile_size + + +contains + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!> @brief Routine returns tile sizes for supplied mesh names/extrusion where +!> applicable to the multigrid configuration. +!> +!> @param[in] config Application configuration object +!> @param[in] local_mesh_names Meshes to set multigrid tile sizes +!> @param[in] extrusion Extrusion object being applied to meshes. +!> +!> @return tile_size Updated tile sizes for multigrid, if applicable. +!> Missing data indicator is returned for where the local +!> mesh tile size is not to be updated for multigrid. +! +function get_multigrid_tile_size( config, local_mesh_names, extrusion) & + result ( tile_size ) + + implicit none + + type(config_type), intent(in) :: config + character(str_def), intent(in) :: local_mesh_names(:) + class(extrusion_type), intent(in) :: extrusion + + integer(i_def), allocatable :: tile_size(:,:) + + integer(i_def) :: multigrid_level + integer(i_def) :: max_multigrid_level + logical(l_def) :: coarsen_multigrid_tiles + logical(l_def) :: set_tile_size + + character(str_def), allocatable :: chain_mesh_tags(:) + + integer(i_def) :: extrusion_id, i + character(str_def) :: name + + !========================================================================= + ! This whole section should probably be in gungho science. It allows the + ! Gungho multigrid scheme to override the tile settings in the + ! configuration. This should really be written in the gungho science, + ! though the decision to call it should be made by the application, i.e. + ! the application may wish to use it's own tileing settings. + ! + ! In partitioning namelist, should be in multigrid + ! max_tiled_multigrid_level = config%multigrid%max_tiled_multigrid_level() + ! coarsen_multigrid_tiles = config%multigrid%coarsen_multigrid_tiles() + coarsen_multigrid_tiles = config%multigrid%coarsen_multigrid_tiles() + max_multigrid_level = config%multigrid%max_tiled_multigrid_level() + chain_mesh_tags = config%multigrid%chain_mesh_tags() + + !========================================================================= + if (allocated(tile_size)) deallocate(tile_size) + allocate(tile_size(2,(size(local_mesh_names)))) + tile_size = imdi + + if (coarsen_multigrid_tiles) then + + extrusion_id = extrusion%get_id() + select case (extrusion_id) + case(prime_extrusion, shifted, double_level) + + ! Set coarsest multigrid level that will be tiled; + ! restrict to the finest grid by default + if (max_multigrid_level == imdi) then + call log_event('no max multigrid level set', log_level_error) + end if + + do i=1, size(local_mesh_names) + set_tile_size = .false. + name =local_mesh_names(i) + + ! Multigrid setup - use tiling if multigrid level is allowed, and + ! if mesh name includes the mesh tag at that level + do multigrid_level=1, size(chain_mesh_tags) + if ( index( trim(name), & + trim(chain_mesh_tags(multigrid_level)) ) > 0 & + .and. multigrid_level <= max_multigrid_level ) then + set_tile_size = .true. + exit + end if + end do + + if (set_tile_size) then + do multigrid_level=1, size(chain_mesh_tags) + if ( index( trim(name), & + trim(chain_mesh_tags(multigrid_level)) ) > 0 ) then + exit + end if + tile_size(:,i) = max( tile_size(:,i)/2, 1 ) + end do + end if ! set_tile_size + end do ! local_mesh_names + + case default + return + end select + + end if ! Coarsen multigrid_tiles + +end function get_multigrid_tile_size + +end module multigrid_mod diff --git a/components/driver/unit-test/assign_coordinate_alphabetaz_mod_test.pf b/components/driver/unit-test/assign_coordinate_alphabetaz_mod_test.pf index 2928f43ca..2fc232c3a 100644 --- a/components/driver/unit-test/assign_coordinate_alphabetaz_mod_test.pf +++ b/components/driver/unit-test/assign_coordinate_alphabetaz_mod_test.pf @@ -14,47 +14,10 @@ module assign_coordinate_alphabetaz_mod_test implicit none private - public :: set_up, tear_down, test_all - - real(kind=r_def), parameter :: radius = 104.0_r_def + public :: test_all contains - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - @before - subroutine set_up() - - use extrusion_config_mod, only : method_uniform, & - stretching_method_linear - use feign_config_mod, only : feign_extrusion_config, & - feign_planet_config - - implicit none - - call feign_extrusion_config( method=method_uniform, & - planet_radius=radius, & - domain_height=10.0_r_def, & - number_of_layers=5_i_def, & - stretching_method=stretching_method_linear, & - stretching_height=15.0_r_def, & - eta_values=(/0.5_r_def/) ) - - call feign_planet_config( scaling_factor=1.0_r_def ) - - end subroutine set_up - - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - @after - subroutine tear_down() - - use config_loader_mod, only: final_configuration - - implicit none - - call final_configuration() - - end subroutine tear_down - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! @test subroutine test_all() @@ -64,6 +27,8 @@ contains implicit none real(kind=r_def), parameter :: dh = 17.0_r_def + real(kind=r_def), parameter :: scaled_radius = 104.0_r_def + real(kind=r_def), parameter :: radius = 104.0_r_def integer(kind=i_def) :: nlayers, ndf_chi, undf_chi, ndf_pid, undf_pid integer(kind=i_def) :: map_chi(1), map_pid(1), nverts, i @@ -152,8 +117,8 @@ contains alpha_out, beta_out, height_out, & verts_XYZ, nodal_coord, verts_ref, & to_rotate, inverse_rot_matrix, & - stretch_factor, & - panel_id, ndf_pid, undf_pid, map_pid ) + stretch_factor, panel_id, scaled_radius, & + ndf_pid, undf_pid, map_pid) ! The answer should be the central point of the box ! alpha = pi/12, beta = -pi/24, h = dh / 2, panel_id = 3 @@ -176,8 +141,9 @@ contains alpha_out, beta_out, height_out, & verts_XYZ, nodal_coord, verts_ref, & to_rotate, inverse_rot_matrix, & - stretch_factor, & - panel_id, ndf_pid, undf_pid, map_pid ) + stretch_factor, panel_id, scaled_radius, & + ndf_pid, undf_pid, map_pid) + ! Alpha and beta swap from the example above ! alpha = -pi/24, beta = pi/12, h = dh / 2, panel_id = 5 @@ -220,8 +186,8 @@ contains alpha_out, beta_out, height_out, & verts_XYZ, nodal_coord, verts_ref, & to_rotate, inverse_rot_matrix, & - stretch_factor, & - panel_id, ndf_pid, undf_pid, map_pid ) + stretch_factor, panel_id, scaled_radius, & + ndf_pid, undf_pid, map_pid) ! Answer is in panel 3, with alpha and beta rotated from panel 1 ! alpha = 0.0, beta = -pi/6 (negative of the original longitude), h = dh / 2 diff --git a/components/driver/unit-test/assign_coordinate_lonlatz_mod_test.pf b/components/driver/unit-test/assign_coordinate_lonlatz_mod_test.pf index 175ba7ce2..83e72a780 100644 --- a/components/driver/unit-test/assign_coordinate_lonlatz_mod_test.pf +++ b/components/driver/unit-test/assign_coordinate_lonlatz_mod_test.pf @@ -13,47 +13,10 @@ module assign_coordinate_lonlatz_mod_test implicit none private - public :: set_up, tear_down, test_all - - real(kind=r_def), parameter :: radius = 19.0_r_def + public :: test_all contains - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - @before - subroutine set_up() - - use extrusion_config_mod, only : method_uniform, & - stretching_method_linear - use feign_config_mod, only : feign_extrusion_config, & - feign_planet_config - - implicit none - - call feign_extrusion_config( method=method_uniform, & - planet_radius=radius, & - domain_height=10.0_r_def, & - number_of_layers=5_i_def, & - stretching_method=stretching_method_linear, & - stretching_height=15.0_r_def, & - eta_values=(/0.5_r_def/) ) - - call feign_planet_config( scaling_factor=1.0_r_def ) - - end subroutine set_up - - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - @after - subroutine tear_down() - - use config_loader_mod, only : final_configuration - - implicit none - - call final_configuration() - - end subroutine tear_down - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! @test subroutine test_all() @@ -63,6 +26,8 @@ contains implicit none + real(kind=r_def), parameter :: radius = 19.0_r_def + real(kind=r_def), parameter :: scaled_radius = 19.0_r_def real(kind=r_def), parameter :: dh = 2.4_r_def real(kind=r_def), parameter :: tol = 1.0e-10_r_def ! for r_def 64-bit real(kind=r_def) :: use_tol @@ -118,17 +83,18 @@ contains ! Test with no rotation ---------------------------------------------------- to_rotate = .false. - inverse_rot_matrix = reshape([[1.0_r_def, 0.0_r_def, 0.0_r_def], & - [0.0_r_def, 1.0_r_def, 0.0_r_def], & - [0.0_r_def, 0.0_r_def, 1.0_r_def]], & + inverse_rot_matrix = reshape([[1.0_r_def, 0.0_r_def, 0.0_r_def], & + [0.0_r_def, 1.0_r_def, 0.0_r_def], & + [0.0_r_def, 0.0_r_def, 1.0_r_def]], & shape=[3,3]) - call assign_coordinate_lonlatz(nlayers, ndf_chi, nverts, & - undf_chi, map_chi, & - longitude, latitude, height, & - verts_XYZ, nodal_coord, verts_ref, & - to_rotate, inverse_rot_matrix, & - panel_id, ndf_pid, undf_pid, map_pid ) + call assign_coordinate_lonlatz(nlayers, ndf_chi, nverts, & + undf_chi, map_chi, & + longitude, latitude, height, & + verts_XYZ, nodal_coord, verts_ref, & + to_rotate, inverse_rot_matrix, & + panel_id, scaled_radius, & + ndf_pid, undf_pid, map_pid) ! The answer should be the central point of the box ! lon = 5*pi/12, lat = 5*pi/24, h = dh / 2 @@ -158,18 +124,19 @@ contains ! Rotation matrix corresponding to North pole at lon=pi/2 and lat=pi/2 to_rotate = .true. ! The corresponding rotation matrix is: - inverse_rot_matrix = reshape( & - [[ 0.0_r_def, -1.0_r_def, 0.0_r_def], & - [ 1.0_r_def, 0.0_r_def, 0.0_r_def], & - [ 0.0_r_def, 0.0_r_def, 1.0_r_def]], & + inverse_rot_matrix = reshape( & + [[ 0.0_r_def, -1.0_r_def, 0.0_r_def], & + [ 1.0_r_def, 0.0_r_def, 0.0_r_def], & + [ 0.0_r_def, 0.0_r_def, 1.0_r_def]], & shape=[3,3]) - call assign_coordinate_lonlatz(nlayers, ndf_chi, nverts, & - undf_chi, map_chi, & - longitude, latitude, height, & - verts_XYZ, nodal_coord, verts_ref, & - to_rotate, inverse_rot_matrix, & - panel_id, ndf_pid, undf_pid, map_pid ) + call assign_coordinate_lonlatz(nlayers, ndf_chi, nverts, & + undf_chi, map_chi, & + longitude, latitude, height, & + verts_XYZ, nodal_coord, verts_ref, & + to_rotate, inverse_rot_matrix, & + panel_id, scaled_radius, & + ndf_pid, undf_pid, map_pid) ! The answer should be the central point of the box ! lon = -pi/12, lat = 5*pi/24, h = dh / 2 diff --git a/components/driver/unit-test/assign_coordinate_xyz_mod_test.pf b/components/driver/unit-test/assign_coordinate_xyz_mod_test.pf index 8ad3613f9..37eaac9ef 100644 --- a/components/driver/unit-test/assign_coordinate_xyz_mod_test.pf +++ b/components/driver/unit-test/assign_coordinate_xyz_mod_test.pf @@ -8,80 +8,34 @@ module assign_coordinate_xyz_mod_test use constants_mod, only : r_def, i_def + + use base_mesh_config_mod, only: geometry_planar, topology_fully_periodic use funit implicit none private - public :: assign_coordinate_xyz_test_type, test_all - - @TestCase - type, extends(TestCase) :: assign_coordinate_xyz_test_type - private - contains - procedure setUp - procedure tearDown - procedure test_all - end type assign_coordinate_xyz_test_type + public :: test_all contains - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - subroutine setUp( this ) - - use base_mesh_config_mod, only : geometry_planar, & - topology_fully_periodic - use feign_config_mod, only : feign_base_mesh_config, & - feign_finite_element_config - use finite_element_config_mod, only : cellshape_quadrilateral, & - coord_system_xyz - - implicit none - - class(assign_coordinate_xyz_test_type), intent(inout) :: this - - call feign_base_mesh_config( file_prefix='foo', & - prime_mesh_name='unit_test', & - geometry=geometry_planar, & - prepartitioned=.false., & - topology=topology_fully_periodic, & - fplane=.false., f_lat_deg=0.0_r_def ) - - call feign_finite_element_config( cellshape=cellshape_quadrilateral, & - coord_order=0_i_def, & - coord_system=coord_system_xyz, & - element_order_h=0_i_def, & - element_order_v=0_i_def, & - rehabilitate=.true. ) - - end subroutine setUp - - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - subroutine tearDown( this ) - - use config_loader_mod, only: final_configuration - - implicit none - - class(assign_coordinate_xyz_test_type), intent(inout) :: this - - call final_configuration() - - end subroutine tearDown - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! @Test - subroutine test_all( this ) + subroutine test_all() use driver_coordinates_mod, only : assign_coordinate_xyz - implicit none - class(assign_coordinate_xyz_test_type), intent(inout) :: this + + implicit none real(kind=r_def), parameter :: tol = 1.0e-3_r_def, & one = 1.0_r_def + integer(i_def), parameter :: geometry = geometry_planar + integer(i_def), parameter :: topology = topology_fully_periodic + real(r_def), parameter :: scaled_radius = 1.0_r_def + integer(kind=i_def) :: nlayers, ndf, nverts, i, undf, ndf_pid, undf_pid integer(kind=i_def) :: map(1), map_pid(1) real(kind=r_def) :: x(1),y(1),z(1), dz(1), panel_id(1) @@ -114,8 +68,9 @@ contains call assign_coordinate_xyz( nlayers, ndf, nverts, undf, map, dz, x, y, z, & vertices_phys, nodal_coord, vertices_comp, & - 2.0_r_def, 0.0_r_def, panel_id, ndf_pid, & - undf_pid, map_pid ) + 2.0_r_def, 0.0_r_def, panel_id, & + geometry, topology, scaled_radius, & + ndf_pid, undf_pid, map_pid ) @assertEqual( one, x(1), tol ) @assertEqual( one, y(1), tol ) diff --git a/components/driver/unit-test/mesh/create_mesh_mod_test.pf b/components/driver/unit-test/mesh/create_mesh_mod_test.pf index 646c641a1..e83145d59 100644 --- a/components/driver/unit-test/mesh/create_mesh_mod_test.pf +++ b/components/driver/unit-test/mesh/create_mesh_mod_test.pf @@ -19,61 +19,46 @@ module create_mesh_mod_test use mesh_mod, only : mesh_type use extrusion_mod, only : extrusion_type, & PRIME_EXTRUSION - use extrusion_config_mod, only : METHOD_UNIFORM - use feign_config_mod, only : feign_multigrid_config, & - feign_partitioning_config use create_mesh_mod, only : create_extrusion, & create_mesh use panel_decomposition_mod, only : custom_decomposition_type + + use extrusion_config_mod, only: method_uniform + use partitioning_config_mod, only: panel_decomposition_auto, & + partitioner_planar + use pfunit implicit none private - public :: test_single_cell_tiling - - @TestCase - type, extends(TestCase), public :: create_mesh_test_type - private - contains - procedure setUp - procedure tearDown - end type create_mesh_test_type + public :: set_up, tear_down, test_single_cell_tiling contains !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - subroutine setUp( this ) + @before + subroutine set_up() implicit none - class(create_mesh_test_type), intent(inout) :: this - - mesh_collection = mesh_collection_type() + mesh_collection = mesh_collection_type() global_mesh_collection = global_mesh_collection_type() - local_mesh_collection = local_mesh_collection_type() - - call feign_multigrid_config( chain_mesh_tags = [ 'dummy' ], & - multigrid_chain_nitems = 1_i_def, & - n_coarsesmooth = 1_i_def, & - n_postsmooth = 1_i_def, & - n_presmooth = 1_i_def, & - smooth_relaxation = 0.0_r_def ) + local_mesh_collection = local_mesh_collection_type() - end subroutine setUp + end subroutine set_up !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - subroutine tearDown( this ) + @after + subroutine tear_down() implicit none - class(create_mesh_test_type), intent(inout) :: this - call mesh_collection%clear() call global_mesh_collection%clear() call local_mesh_collection%clear() - end subroutine tearDown + end subroutine tear_down !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine initialise_mesh_collections( filename, mesh_name, xyprocs, & @@ -148,14 +133,10 @@ contains ! Test 1x1 tiling only, which is equivalent to colouring - other tiling ! setups require MPI @Test - subroutine test_single_cell_tiling( this ) - - use partitioning_config_mod, only : panel_decomposition_auto, & - partitioner_planar + subroutine test_single_cell_tiling() implicit none - class(create_mesh_test_type), intent(inout) :: this class(extrusion_type), allocatable :: extrusion type(mesh_type), pointer :: mesh integer(i_def) :: xyprocs(2) @@ -176,19 +157,6 @@ contains max_stencil_depth, generate_inner_halos, & 0_i_def, 1_i_def) - call feign_partitioning_config( coarsen_multigrid_tiles = .false., & - generate_inner_halos = .true., & - inner_halo_tiles = .false., & - max_tiled_multigrid_level = 1_i_def, & - panel_decomposition = & - panel_decomposition_auto, & - panel_xproc = 1_i_def, & - panel_yproc = 1_i_def, & - partitioner = & - partitioner_planar, & - tile_size_x = 1_i_def, & - tile_size_y = 1_i_def ) - allocate( extrusion, source=create_extrusion( METHOD_UNIFORM, & 100.0_r_def, & 10000.0_r_def, & @@ -197,7 +165,12 @@ contains ! Check that exactly one mesh is added to the collection @assertEqual( mesh_collection%n_meshes(), 0 ) - call create_mesh( mesh_name, extrusion ) + + call create_mesh( mesh_name, extrusion, & + inner_halo_tiles=.false., & + tile_size=[1,1] ) + + @assertEqual( mesh_collection%n_meshes(), 1 ) ! Check that mesh name is correct diff --git a/components/lfric-xios/integration-test/lfric_xios_context_test.f90 b/components/lfric-xios/integration-test/lfric_xios_context_test.f90 index c2fbc4f59..73811cecf 100644 --- a/components/lfric-xios/integration-test/lfric_xios_context_test.f90 +++ b/components/lfric-xios/integration-test/lfric_xios_context_test.f90 @@ -32,10 +32,15 @@ program lfric_xios_context_test allocate(io_context) call io_context%initialise( "test_io_context", 1, 10 ) - call io_context%initialise_xios_context( test_db%comm, & - test_db%chi, test_db%panel_id, & + call io_context%initialise_xios_context( test_db%comm, & + test_db%chi, test_db%panel_id, & test_db%clock, test_db%calendar, & - before_close ) + before_close, & + test_db%config%base_mesh%geometry(), & + test_db%config%base_mesh%topology(), & + test_db%config%finite_element%coord_system(), & + test_db%config%planet%scaled_radius() ) + deallocate(io_context) ! ============================== Finish test ================================= diff --git a/components/lfric-xios/integration-test/lfric_xios_cyclic_temporal_test.f90 b/components/lfric-xios/integration-test/lfric_xios_cyclic_temporal_test.f90 index 03f062cbd..a67048327 100644 --- a/components/lfric-xios/integration-test/lfric_xios_cyclic_temporal_test.f90 +++ b/components/lfric-xios/integration-test/lfric_xios_cyclic_temporal_test.f90 @@ -64,10 +64,14 @@ program lfric_xios_cyclic_temporal_test fields_in_file=test_db%temporal_fields ) ) before_close => null() - call io_context%initialise_xios_context( test_db%comm, & - test_db%chi, test_db%panel_id, & + call io_context%initialise_xios_context( test_db%comm, & + test_db%chi, test_db%panel_id, & test_db%clock, test_db%calendar, & - before_close ) + before_close, & + test_db%config%base_mesh%geometry(), & + test_db%config%base_mesh%topology(), & + test_db%config%finite_element%coord_system(), & + test_db%config%planet%scaled_radius() ) context_advance => advance diff --git a/components/lfric-xios/integration-test/lfric_xios_temporal_interp_test.f90 b/components/lfric-xios/integration-test/lfric_xios_temporal_interp_test.f90 index 1d0122a23..ca1ea7b90 100644 --- a/components/lfric-xios/integration-test/lfric_xios_temporal_interp_test.f90 +++ b/components/lfric-xios/integration-test/lfric_xios_temporal_interp_test.f90 @@ -8,7 +8,7 @@ ! program lfric_xios_temporal_interp_test - use constants_mod, only: i_timestep, r_second + use constants_mod, only: i_timestep, r_second, r_def, i_def use event_mod, only: event_action use event_actor_mod, only: event_actor_type use field_mod, only: field_type, field_proxy_type @@ -38,9 +38,19 @@ program lfric_xios_temporal_interp_test type(xios_date) :: date integer(i_timestep) :: file_freq + integer(i_def) :: geometry + integer(i_def) :: topology + integer(i_def) :: coord_system + real(r_def) :: scaled_radius + call test_db%initialise() call lfric_xios_initialise( "test", test_db%comm, .false. ) + geometry = test_db%config%base_mesh%geometry() + topology = test_db%config%base_mesh%topology() + coord_system = test_db%config%finite_element%coord_system() + scaled_radius = test_db%config%planet%scaled_radius() + ! =============================== Start test ================================ allocate(io_context) @@ -66,11 +76,11 @@ program lfric_xios_temporal_interp_test fields_in_file=test_db%temporal_fields ) ) before_close => null() - call io_context%initialise_xios_context( test_db%comm, & - test_db%chi, test_db%panel_id, & - test_db%clock, test_db%calendar, & - before_close ) - + call io_context%initialise_xios_context( test_db%comm, & + test_db%chi, test_db%panel_id, & + test_db%clock, test_db%calendar, & + before_close, geometry, topology, & + coord_system, scaled_radius ) context_advance => advance context_actor => io_context diff --git a/components/lfric-xios/integration-test/lfric_xios_temporal_iodef_test.f90 b/components/lfric-xios/integration-test/lfric_xios_temporal_iodef_test.f90 index a2b8cac79..1c3e8a303 100644 --- a/components/lfric-xios/integration-test/lfric_xios_temporal_iodef_test.f90 +++ b/components/lfric-xios/integration-test/lfric_xios_temporal_iodef_test.f90 @@ -58,11 +58,14 @@ program lfric_xios_temporal_iodef_test fields_in_file=test_db%temporal_fields ) ) before_close => null() - call io_context%initialise_xios_context( test_db%comm, & - test_db%chi, test_db%panel_id, & + call io_context%initialise_xios_context( test_db%comm, & + test_db%chi, test_db%panel_id, & test_db%clock, test_db%calendar, & - before_close ) - + before_close, & + test_db%config%base_mesh%geometry(), & + test_db%config%base_mesh%topology(), & + test_db%config%finite_element%coord_system(), & + test_db%config%planet%scaled_radius() ) context_advance => advance context_actor => io_context diff --git a/components/lfric-xios/integration-test/lfric_xios_temporal_test.f90 b/components/lfric-xios/integration-test/lfric_xios_temporal_test.f90 index 60df5997e..572ce0193 100644 --- a/components/lfric-xios/integration-test/lfric_xios_temporal_test.f90 +++ b/components/lfric-xios/integration-test/lfric_xios_temporal_test.f90 @@ -62,10 +62,14 @@ program lfric_xios_temporal_test fields_in_file=test_db%temporal_fields ) ) before_close => null() - call io_context%initialise_xios_context( test_db%comm, & - test_db%chi, test_db%panel_id, & + call io_context%initialise_xios_context( test_db%comm, & + test_db%chi, test_db%panel_id, & test_db%clock, test_db%calendar, & - before_close ) + before_close, & + test_db%config%base_mesh%geometry(), & + test_db%config%base_mesh%topology(), & + test_db%config%finite_element%coord_system(), & + test_db%config%planet%scaled_radius() ) context_advance => advance diff --git a/components/lfric-xios/integration-test/lfric_xios_time_read_test.f90 b/components/lfric-xios/integration-test/lfric_xios_time_read_test.f90 index c26d1c577..2f453a381 100755 --- a/components/lfric-xios/integration-test/lfric_xios_time_read_test.f90 +++ b/components/lfric-xios/integration-test/lfric_xios_time_read_test.f90 @@ -31,10 +31,14 @@ program lfric_xios_time_read_test allocate(io_context) call io_context%initialise( "test_io_context", 1, 10 ) - call io_context%initialise_xios_context( test_db%comm, & - test_db%chi, test_db%panel_id, & + call io_context%initialise_xios_context( test_db%comm, & + test_db%chi, test_db%panel_id, & test_db%clock, test_db%calendar, & - before_close ) + before_close, & + test_db%config%base_mesh%geometry(), & + test_db%config%base_mesh%topology(), & + test_db%config%finite_element%coord_system(), & + test_db%config%planet%scaled_radius() ) allocate(check(10)) check = [ xios_date(2024, 1, 1, 15, 1, 0), & diff --git a/components/lfric-xios/integration-test/test_db_mod.f90 b/components/lfric-xios/integration-test/test_db_mod.f90 index 028dc97cb..f2f95e485 100644 --- a/components/lfric-xios/integration-test/test_db_mod.f90 +++ b/components/lfric-xios/integration-test/test_db_mod.f90 @@ -9,6 +9,7 @@ module test_db_mod use calendar_mod, only: calendar_type use cli_mod, only: parse_command_line + use config_mod, only: config_type use config_loader_mod, only: read_configuration use constants_mod, only: i_def, r_def, str_def, imdi, & r_second, i_timestep @@ -31,7 +32,6 @@ module test_db_mod finalise_logging, & log_set_level, log_event, & LOG_LEVEL_TRACE, LOG_LEVEL_ERROR - use config_mod, only: config_type use lfric_xios_read_mod, only: read_field_generic use lfric_xios_write_mod, only: write_field_generic use local_mesh_collection_mod, only: local_mesh_collection_type, & diff --git a/components/lfric-xios/source/lfric_xios_context_mod.f90 b/components/lfric-xios/source/lfric_xios_context_mod.f90 index 9f951767d..935ec50be 100644 --- a/components/lfric-xios/source/lfric_xios_context_mod.f90 +++ b/components/lfric-xios/source/lfric_xios_context_mod.f90 @@ -9,7 +9,7 @@ module lfric_xios_context_mod use calendar_mod, only : calendar_type use clock_mod, only : clock_type - use constants_mod, only : i_def, & + use constants_mod, only : i_def, r_def, & r_second, i_timestep, & l_def use field_mod, only : field_type @@ -81,18 +81,26 @@ end subroutine initialise_lfric_xios_context !> @brief Set up an XIOS context. !> - !> @param [in] communicator MPI communicator used by context. - !> @param [in] chi Array of coordinate fields - !> @param [in] panel_id Panel ID field - !> @param [in] model_clock The model clock. - !> @param [in] calendar The model calendar. - !> @param [in] before_close Routine to be called before context closes - !> @param [in] alt_coords Array of coordinate fields for alternative meshes - !> @param [in] alt_panel_ids Panel ID fields for alternative meshes - subroutine initialise_xios_context( this, communicator, & + !> @param [in] communicator MPI communicator used by context. + !> @param [in] chi Array of coordinate fields + !> @param [in] panel_id Panel ID field + !> @param [in] model_clock The model clock. + !> @param [in] calendar The model calendar. + !> @param [in] before_close Routine to be called before context closes + !> @param [in] geometry Mesh geometry enumeration value. + !> @param [in] topology Mesh topology enumeration value. + !> @param [in] coord_system Finite-element coord-system enumeration value + !> @param [in] scaled_radius Planet scaled radius + !> @param [in] alt_coords Array of coordinate fields for alternative meshes + !> @param [in] alt_panel_ids Panel ID fields for alternative meshes + subroutine initialise_xios_context( this, & + communicator, & chi, panel_id, & model_clock, calendar, & before_close, & + geometry, topology, & + coord_system, & + scaled_radius, & alt_coords, & alt_panel_ids, & start_at_zero ) @@ -100,6 +108,7 @@ subroutine initialise_xios_context( this, communicator, & implicit none class(lfric_xios_context_type), intent(inout) :: this + type(lfric_comm_type), intent(in) :: communicator type(field_type), intent(in) :: chi(:) type(field_type), intent(in) :: panel_id @@ -107,6 +116,12 @@ subroutine initialise_xios_context( this, communicator, & class(calendar_type), intent(in) :: calendar procedure(callback_clock_arg), pointer, & intent(in) :: before_close + + integer(i_def), intent(in) :: geometry + integer(i_def), intent(in) :: topology + integer(i_def), intent(in) :: coord_system + real(r_def), intent(in) :: scaled_radius + type(field_type), optional, intent(in) :: alt_coords(:,:) type(field_type), optional, intent(in) :: alt_panel_ids(:) logical, optional, intent(in) :: start_at_zero @@ -134,7 +149,9 @@ subroutine initialise_xios_context( this, communicator, & ! Run XIOS setup routines call init_xios_calendar(model_clock, calendar, zero_start, this%context_clock_step) - call init_xios_dimensions(chi, panel_id, alt_coords, alt_panel_ids) + call init_xios_dimensions( chi, panel_id, geometry, topology, & + coord_system, scaled_radius, & + alt_coords, alt_panel_ids ) if (this%filelist%get_length() > 0) call setup_xios_files(this%filelist) if (associated(before_close)) call before_close(model_clock) diff --git a/components/lfric-xios/source/lfric_xios_setup_mod.x90 b/components/lfric-xios/source/lfric_xios_setup_mod.x90 index 9c762d7bb..83ac4e602 100644 --- a/components/lfric-xios/source/lfric_xios_setup_mod.x90 +++ b/components/lfric-xios/source/lfric_xios_setup_mod.x90 @@ -155,14 +155,26 @@ contains !> !> @param[in] chi Coordinate field !> @param[in] panel_id Field with IDs of mesh panels - !> - subroutine init_xios_dimensions(chi, panel_id, alt_coords, alt_panel_ids) + !> @param[in] geometry Mesh geometry enumeration value + !> @param[in] topology Mesh topology enumeration value + !> @param[in] coord_system Finite-element coord-system enumeration value + !> @param[in] scaled_radius Planet scaled radius + subroutine init_xios_dimensions( chi, panel_id, & + geometry, topology, & + coord_system, scaled_radius, & + alt_coords, alt_panel_ids ) implicit none ! Arguments - type(field_type), intent(in) :: chi(:) - type(field_type), intent(in) :: panel_id + type(field_type), intent(in) :: chi(:) + type(field_type), intent(in) :: panel_id + + integer(i_def),intent(in) :: geometry + integer(i_def),intent(in) :: topology + integer(i_def),intent(in) :: coord_system + real(r_def), intent(in) :: scaled_radius + type(field_type), optional, intent(in) :: alt_coords(:,:) type(field_type), optional, intent(in) :: alt_panel_ids(:) @@ -170,12 +182,16 @@ contains type(mesh_type), pointer :: mesh => null() ! Initialise XIOS prime mesh - call init_xios_mesh( chi, panel_id, prime_mesh=.true. ) + call init_xios_mesh( chi, panel_id, geometry, topology, & + coord_system, scaled_radius, & + prime_mesh=.true. ) ! Initialise additional meshes if (present(alt_coords) .and. present(alt_panel_ids)) then do i = 1, size(alt_panel_ids) - call init_xios_mesh( alt_coords(i,:), alt_panel_ids(i), prime_mesh=.false. ) + call init_xios_mesh( alt_coords(i,:), alt_panel_ids(i), & + geometry, topology, coord_system, & + scaled_radius, prime_mesh=.false. ) end do end if @@ -195,16 +211,29 @@ contains !> !> @param[in] chi Coordinate field !> @param[in] panel_id Field with IDs of mesh panels + !> @param[in] geometry Mesh geometry enumeration value + !> @param[in] topology Mesh topology enumeration value + !> @param[in] coord_system Finite-element coord-system enumeration value + !> @param[in] scaled_radius Planet scaled radius !> @param[in] prime_mesh Logical flag denoting if the mesh is the primary !! I/O mesh !> - subroutine init_xios_mesh(chi, panel_id, prime_mesh) + subroutine init_xios_mesh( chi, panel_id, & + geometry, topology, & + coord_system, scaled_radius, & + prime_mesh ) implicit none ! Arguments type(field_type), intent(in) :: chi(:) type(field_type), intent(in) :: panel_id + + integer(i_def),intent(in) :: geometry + integer(i_def),intent(in) :: topology + integer(i_def),intent(in) :: coord_system + real(r_def), intent(in) :: scaled_radius + logical, optional, intent(in) :: prime_mesh ! Local variables @@ -277,7 +306,7 @@ contains r2d = radians_to_degrees else r2d = 1.0_r_def - endif + end if ! Set up fields to hold the output coordinates output_field_fs => function_space_collection%get_fs( mesh, 0, 0, W0 ) @@ -295,7 +324,7 @@ contains r2d = radians_to_degrees else r2d = 1.0_r_def - endif + end if ! Calculate the local size of a W2H fs in order to determine ! how many edge dofs for the current partition @@ -316,7 +345,9 @@ contains call sample_chi(i)%initialise( vector_space = output_field_fs ) end do ! Convert to (X,Y,Z) coordinates - call invoke(nodal_xyz_coordinates_kernel_type(sample_chi, chi, panel_id)) + call invoke(nodal_xyz_coordinates_kernel_type(sample_chi, chi, panel_id, & + geometry, topology, & + coord_system, scaled_radius)) else ! For planar geometries just re-use existing chi which are already (X,Y,Z) do i = 1,3 diff --git a/components/science/source/algorithm/sci_geometric_constants_mod.x90 b/components/science/source/algorithm/sci_geometric_constants_mod.x90 index cffdcc9d1..9bcdc1297 100644 --- a/components/science/source/algorithm/sci_geometric_constants_mod.x90 +++ b/components/science/source/algorithm/sci_geometric_constants_mod.x90 @@ -15,25 +15,25 @@ module sci_geometric_constants_mod ! Infrastructure - use constants_mod, only: i_def, r_def, l_def, str_def - use extrusion_mod, only: TWOD, PRIME_EXTRUSION - use field_mod, only: field_type - use fs_continuity_mod, only: W0, W1, W2, W2H, W3, Wtheta - use function_space_collection_mod, only: function_space_collection - use function_space_mod, only: function_space_type - use integer_field_mod, only: integer_field_type - use inventory_by_mesh_mod, only: inventory_by_mesh_type - use inventory_by_local_mesh_mod, only: inventory_by_local_mesh_type - use local_mesh_mod, only: local_mesh_type - use log_mod, only: log_event, LOG_LEVEL_ERROR - use mesh_collection_mod, only: mesh_collection - use mesh_mod, only: mesh_type - use timing_mod, only: start_timing, stop_timing, & - tik, LPROF + use constants_mod, only: i_def, r_def, l_def, str_def + use extrusion_mod, only: TWOD, PRIME_EXTRUSION + use field_mod, only: field_type + use fs_continuity_mod, only: W0, W1, W2, W2H, W3, Wtheta + use function_space_collection_mod, only: function_space_collection + use function_space_mod, only: function_space_type + use integer_field_mod, only: integer_field_type + use inventory_by_mesh_mod, only: inventory_by_mesh_type + use inventory_by_local_mesh_mod, only: inventory_by_local_mesh_type + use local_mesh_mod, only: local_mesh_type + use log_mod, only: log_event, LOG_LEVEL_ERROR + use mesh_collection_mod, only: mesh_collection + use mesh_mod, only: mesh_type + use timing_mod, only: start_timing, stop_timing, & + tik, LPROF ! Configuration - use finite_element_config_mod, only: element_order_h, & - element_order_v + use base_mesh_config_mod, only: geometry_spherical + use finite_element_config_mod, only: coord_system_native implicit none @@ -136,18 +136,29 @@ contains !> @brief Private routine for computing longitude and latitude fields !> @param[in,out] long_inventory Inventory containing longitude fields !> @param[in,out] lat_inventory Inventory containing latitude fields - !> @param[in] mesh Mesh used to determine local mesh for - !! computing the fields for - !> @param[in] fs_id Identifier for function space to compute - !! longitude and latitude fields for - !> @param[in] use_fe Flag to indicate whether to use finite - !! element or finite volume cells - subroutine compute_latlon(long_inventory, lat_inventory, mesh, fs_id, use_fe) - - use base_mesh_config_mod, only: f_lat, geometry, & - geometry_spherical - use idealised_config_mod, only: f_lon - use sci_compute_latlon_kernel_mod, only: compute_latlon_kernel_type + !> + !> @param[in] mesh Mesh used to determine local mesh for + !! computing the fields for + !> @param[in] fs_id Identifier for function space to compute + !! longitude and latitude fields for + !> @param[in] use_fe Flag to indicate whether to use finite + !! element or finite volume cells + !> @param[in] geometry Mesh geometry enumeration value + !> @param[in] topology Mesh topology enumeration value + !> @param[in] element_order_h Function space order in horizontal + !> @param[in] element_order_v Function space order in vertical + !> @param[in] coord_system Finite-Element coord-system enumeration value + !> @param[in] f_lat Latitiude of f-plane + !> @param[in] f_lon Longitude of f-plane + !> @param[in] scaled_radius Planet scaled radius + subroutine compute_latlon( long_inventory, lat_inventory, & + mesh, fs_id, use_fe, & + geometry, topology, & + element_order_h, element_order_v, & + coord_system, f_lat, f_lon, & + scaled_radius ) + + use sci_compute_latlon_kernel_mod, only: compute_latlon_kernel_type implicit none @@ -157,6 +168,12 @@ contains integer(kind=i_def), intent(in) :: fs_id logical(kind=l_def), intent(in) :: use_fe + integer(i_def), intent(in) :: geometry, topology + integer(i_def), intent(in) :: element_order_h, element_order_v + integer(i_def), intent(in) :: coord_system + real(r_def), intent(in) :: f_lat, f_lon + real(r_def), intent(in) :: scaled_radius + ! Internal variables type(mesh_type), pointer :: twod_mesh type(local_mesh_type), pointer :: local_mesh @@ -187,9 +204,11 @@ contains if ( geometry == geometry_spherical ) then chi => get_coordinates(mesh%get_id()) panel_id => get_panel_id(mesh%get_id()) - call invoke( compute_latlon_kernel_type(lat, long, chi, panel_id) ) + call invoke( compute_latlon_kernel_type(lat, long, chi, panel_id, & + geometry, topology, & + coord_system, scaled_radius) ) else - call invoke( setval_c(lat, f_lat), & + call invoke( setval_c(lat, f_lat), & setval_c(long, f_lon) ) end if @@ -302,16 +321,18 @@ contains end function get_coordinates !> @brief Returns a pointer to the extended coordinate field array - !> @param[in] mesh_id The ID of the mesh to get the object for + !> @param[in] mesh_id The ID of the mesh to get the object for + !> @param[in] coord_system Finite-Element coord-system enumeration value !> @return The coordinate field array - function get_extended_coordinates(mesh_id) result(extended_chi) + function get_extended_coordinates(mesh_id, coord_system) result(extended_chi) - use finite_element_config_mod, only: coord_system, coord_system_native use sci_extend_chi_field_kernel_mod, only: extend_chi_field_kernel_type implicit none - integer(kind=i_def), intent(in) :: mesh_id + integer(i_def), intent(in) :: coord_system + integer(i_def), intent(in) :: mesh_id + type(mesh_type), pointer :: mesh type(field_type), pointer :: extended_chi(:) logical(kind=l_def) :: constant_exists @@ -319,7 +340,8 @@ contains type(field_type), pointer :: chi(:) type(field_type), pointer :: panel_id type(function_space_type), pointer :: wchi_fs - integer(tik) :: id + + integer(tik) :: id ! Initialise inventory if this is the first time getting this constant if (.not. extended_chi_inventory%is_initialised()) then @@ -368,7 +390,7 @@ contains !> @return The dA field function get_dA_at_w2(mesh_id) result(dA_at_w2) - use sci_calc_da_at_w2_kernel_mod, only: calc_dA_at_w2_kernel_type + use sci_calc_da_at_w2_kernel_mod, only: calc_dA_at_w2_kernel_type implicit none @@ -411,24 +433,32 @@ contains end function get_dA_at_w2 !> @brief Returns the (finite element) Det(J) values at W3 dof locations - !> @param[in] mesh_id The ID of the mesh to get the object for + !> @param[in] mesh_id The ID of the mesh to get the object for + !> @param[in] element_order_h Function space order in horizontal + !> @param[in] element_order_v Function space order in vertical + !> @param[in] nqp_h_exact Number of quadrature points in horizontal + !> @param[in] nqp_v_exact Number of quadrature points in vertical !> @return The Det(J) field - function get_detj_at_w3_fe(mesh_id) result(detj_at_w3) + + function get_detj_at_w3_fe( mesh_id, element_order_h, element_order_v, & + nqp_h_exact, nqp_v_exact) & + result( detj_at_w3 ) ! @TODO #4487: update these imports ! use sci_calc_detj_at_w3_kernel_mod, only: calc_detj_at_w3_kernel_type use sci_compute_mass_matrix_kernel_w_scalar_mod, & only: compute_mass_matrix_kernel_w_scalar_type use sci_mm_diagonal_kernel_mod, only: mm_diagonal_kernel_type - use finite_element_config_mod, only: nqp_h_exact, & - nqp_v_exact use operator_mod, only: operator_type use quadrature_xyoz_mod, only: quadrature_xyoz_type use quadrature_rule_gaussian_mod, only: quadrature_rule_gaussian_type implicit none - integer(kind=i_def), intent(in) :: mesh_id + integer(i_def), intent(in) :: mesh_id + integer(i_def), intent(in) :: element_order_h, element_order_v + integer(i_def), intent(in) :: nqp_h_exact, nqp_v_exact + type(mesh_type), pointer :: mesh logical(kind=l_def) :: constant_exists type(field_type), pointer :: detj_at_w3 @@ -443,7 +473,7 @@ contains integer(tik) :: id ! If running at lowest order, use finite volume - if (element_order_h == 0 .and. element_order_v == 0) then + if ( element_order_h == 0 .and. element_order_v == 0) then detj_at_w3 => get_detj_at_w3_fv(mesh_id) return end if @@ -464,8 +494,10 @@ contains ! Create the object as it doesn't exist yet if ( LPROF ) call start_timing( id, 'runtime_constants.geometric' ) - w3_fs => function_space_collection%get_fs(mesh, element_order_h, & - element_order_v, W3) + w3_fs => function_space_collection%get_fs( mesh, & + element_order_h, & + element_order_v, & + W3 ) call detj_at_w3_inventory_fe%add_field(detj_at_w3, w3_fs, mesh) ! @TODO #4487: it is inefficient to calculate this via mass matrices @@ -564,16 +596,21 @@ contains end function get_detj_at_w3_fv !> @brief Returns the (finite element) Det(J) values at W2 dof locations - !> @param[in] mesh_id The ID of the mesh to get the object for + !> @param[in] mesh_id The ID of the mesh to get the object for + !> @param[in] element_order_h Function space order in horizontal + !> @param[in] element_order_v Function space order in vertical !> @return The Det(J) field - function get_detj_at_w2_fe(mesh_id) result(detj_at_w2) + function get_detj_at_w2_fe( mesh_id, element_order_h, element_order_v ) & + result( detj_at_w2 ) use sci_calc_detj_at_w2_kernel_mod, only: calc_detj_at_w2_kernel_type use sci_multiplicity_kernel_mod, only: multiplicity_kernel_type implicit none - integer(kind=i_def), intent(in) :: mesh_id + integer(i_def), intent(in) :: mesh_id + integer(i_def), intent(in) :: element_order_h, element_order_v + type(mesh_type), pointer :: mesh logical(kind=l_def) :: constant_exists type(field_type), pointer :: detj_at_w2 @@ -605,8 +642,10 @@ contains ! Create the object as it doesn't exist yet if ( LPROF ) call start_timing( id, 'runtime_constants.geometric' ) - w2_fs => function_space_collection%get_fs(mesh, element_order_h, & - element_order_v, W2) + w2_fs => function_space_collection%get_fs( mesh, & + element_order_h, & + element_order_v, & + W2 ) call multiplicity_w2%initialise( w2_fs ) call detj_at_w2_inventory_fe%add_field(detj_at_w2, w2_fs, mesh) @@ -681,15 +720,25 @@ contains end function get_detj_at_w2_fv !> @brief Returns a pointer to the vertical grid spacing, located at W3 DoFs - !> @param[in] mesh_id The ID of the mesh to get the object for + !> @param[in] mesh_id The ID of the mesh to get the object for + !> @param[in] geometry Mesh geometry enumeration value + !> @param[in] coord_system Finite-Element coord-system enumeration value + !> @param[in] scaled_radius Planet scaled radius !> @return The physical height difference of layers, at W3 - function get_dz_w3(mesh_id) result(dz_w3) + function get_dz_w3( mesh_id, geometry, & + coord_system, scaled_radius ) & + result( dz_w3 ) - use sci_get_dz_w3_kernel_mod, only: get_dz_w3_kernel_type + use sci_get_dz_w3_kernel_mod, only: get_dz_w3_kernel_type implicit none - integer(kind=i_def), intent(in) :: mesh_id + integer(i_def), intent(in) :: mesh_id + + integer(i_def), intent(in) :: geometry + integer(i_def), intent(in) :: coord_system + real(r_def), intent(in) :: scaled_radius + type(mesh_type), pointer :: mesh type(field_type), pointer :: dz_w3 logical(kind=l_def) :: constant_exists @@ -708,7 +757,9 @@ contains if (.not. constant_exists) then ! If this constant doesn't exist, create it ! Get height first to avoid potentially timing twice - height_w2 => get_height_fv(W2, mesh_id) + height_w2 => get_height_fv( W2, mesh_id, & + geometry, coord_system, & + scaled_radius ) if ( LPROF ) call start_timing( id, 'runtime_constants.geometric' ) @@ -815,17 +866,25 @@ contains end function get_dx_at_w2 - !> @brief Returns the 1/dz values at lowest-order Wtheta DoF locations - !> @param[in] mesh_id The ID of the mesh to get the object for + !> @param[in] mesh_id The ID of the mesh to get the object for + !> @param[in] geometry Mesh geometry enumeration value + !> @param[in] coord_system Finite-Element coord-system enumeration value + !> @param[in] scaled_radius Planet scaled radius !> @return The dz_at_wtheta field - function get_dz_at_wtheta(mesh_id) result(dz_at_wtheta) + function get_dz_at_wtheta( mesh_id, geometry, & + coord_system, scaled_radius ) & + result( dz_at_wtheta ) use sci_calc_dz_face_kernel_mod, only: calc_dz_face_kernel_type implicit none - integer(kind=i_def), intent(in) :: mesh_id + integer(i_def), intent(in) :: mesh_id + integer(i_def), intent(in) :: geometry + integer(i_def), intent(in) :: coord_system + real(r_def), intent(in) :: scaled_radius + type(mesh_type), pointer :: mesh type(field_type), pointer :: dz_at_wtheta type(function_space_type), pointer :: wtheta_k0_fs @@ -849,8 +908,10 @@ contains ! Create constant if it doesn't already exist if (.not. constant_exists) then ! NB: this assumes heights are in the lowest-order space - height_w3 => get_height_fv(W3, mesh_id) - height_wth => get_height_fv(Wtheta, mesh_id) + height_w3 => get_height_fv( W3, mesh_id, geometry, coord_system, & + scaled_radius ) + height_wth => get_height_fv( Wtheta, mesh_id,geometry, coord_system, & + scaled_radius ) if ( LPROF ) call start_timing( id, 'runtime_constants.geometric' ) @@ -871,17 +932,24 @@ contains !> @brief Returns the surface area of a cell projected to mean sea level !> i.e. ignoring the orographic effect on the area - !> @param[in] mesh_id The ID of the mesh to get the object for + !> @param[in] mesh_id The ID of the mesh to get the object for + !> @param[in] geometry Mesh geometry enumeration value + !> @param[in] planet_radius Planet radius (m) + !> @param[in] domain_height Top of atmosphere height above mean surface (m) !> @return The dA_msl_proj field - function get_dA_msl_proj(mesh_id) result(dA_msl_proj) + function get_dA_msl_proj( mesh_id, geometry, & + planet_radius, domain_height ) & + result( dA_msl_proj ) - use base_mesh_config_mod, only: geometry, geometry_spherical - use extrusion_config_mod, only: planet_radius, domain_height use sci_calc_da_msl_proj_kernel_mod, only: calc_da_msl_proj_kernel_type implicit none - integer(kind=i_def), intent(in) :: mesh_id + integer(i_def), intent(in) :: mesh_id + integer(i_def), intent(in) :: geometry + real(r_def), intent(in) :: planet_radius + real(r_def), intent(in) :: domain_height + integer(kind=i_def) :: local_mesh_id type(mesh_type), pointer :: mesh type(mesh_type), pointer :: prime_mesh @@ -928,17 +996,35 @@ contains ! ========================================================================== ! ! PHYSICAL COORDINATES OF DOFs ! ========================================================================== ! - !> @brief Returns a pointer to the longitude of finite element DoFs - !> @param[in] space_id The space for which to get the longitude of DoFs for - !> @param[in] mesh_id The ID of the mesh to get the object for + !> @param[in] space_id The space for which to get the longitude of DoFs for + !> @param[in] mesh_id The ID of the mesh to get the object for + !> @param[in] geometry Mesh geometry enumeration value + !> @param[in] topology Mesh topology enumeration value + !> @param[in] element_order_h Function space order in horizontal + !> @param[in] element_order_v Function space order in vertical + !> @param[in] coord_system Finite-Element coord-system enumeration value + !> @param[in] f_lat Latitiude of f-plane + !> @param[in] f_lon Longitude of f-plane + !> @param[in] scaled_radius Planet scaled radius !> @return The longitude field - function get_longitude_fe(space_id, mesh_id) result(long_ptr) + function get_longitude_fe( space_id, mesh_id, geometry, topology, & + element_order_h, element_order_v, & + coord_system, f_lat, f_lon,scaled_radius ) & + result( long_ptr ) implicit none - integer(kind=i_def), intent(in) :: space_id - integer(kind=i_def), intent(in) :: mesh_id + integer(i_def), intent(in) :: space_id + integer(i_def), intent(in) :: mesh_id + + integer(i_def), intent(in) :: element_order_h, element_order_v + integer(i_def), intent(in) :: geometry, topology + integer(i_def), intent(in) :: coord_system + + real(r_def), intent(in) :: f_lat, f_lon + real(r_def), intent(in) :: scaled_radius + type(mesh_type), pointer :: mesh type(local_mesh_type), pointer :: local_mesh type(inventory_by_local_mesh_type), pointer :: long_inventory @@ -947,9 +1033,14 @@ contains logical(kind=l_def) :: constant_exists character(len=str_def) :: inventory_name + logical(l_def), parameter :: use_fe = .true. + ! If running at lowest order, use finite volume if (element_order_h == 0 .and. element_order_v == 0) then - long_ptr => get_longitude_fv(space_id, mesh_id) + long_ptr => get_longitude_fv( space_id, mesh_id, geometry, topology, & + element_order_h, element_order_v, & + coord_system, f_lat, f_lon, & + scaled_radius ) return end if @@ -985,8 +1076,10 @@ contains constant_exists = long_inventory%paired_object_exists(local_mesh%get_id()) if (.not. constant_exists) then - call compute_latlon(long_inventory, lat_inventory, mesh, space_id, & - use_fe=.true.) + call compute_latlon( long_inventory, lat_inventory, & + mesh, space_id, use_fe, geometry, topology, & + element_order_h, element_order_v, coord_system, & + f_lat, f_lon, scaled_radius ) end if call long_inventory%get_field(local_mesh, long_ptr) @@ -994,15 +1087,34 @@ contains end function get_longitude_fe !> @brief Returns a pointer to the longitude of finite volume DoFs - !> @param[in] space_id The space for which to get the longitude of DoFs for - !> @param[in] mesh_id The ID of the mesh to get the object for + !> @param[in] space_id The space for which to get the longitude of DoFs for + !> @param[in] mesh_id The ID of the mesh to get the object for + !> @param[in] geometry Mesh geometry enumeration value + !> @param[in] topology Mesh topology enumeration value + !> @param[in] element_order_h Function space order in horizontal + !> @param[in] element_order_v Function space order in vertical + !> @param[in] coord_system Finite-Element coord-system enumeration value + !> @param[in] f_lat Latitiude of f-plane + !> @param[in] f_lon Longitude of f-plane + !> @param[in] scaled_radius Planet scaled radius !> @return The longitude field - function get_longitude_fv(space_id, mesh_id) result(long_ptr) + function get_longitude_fv( space_id, mesh_id, geometry, topology, & + element_order_h, element_order_v, & + coord_system, f_lat, f_lon, scaled_radius ) & + result( long_ptr ) implicit none - integer(kind=i_def), intent(in) :: space_id - integer(kind=i_def), intent(in) :: mesh_id + integer(i_def), intent(in) :: space_id + integer(i_def), intent(in) :: mesh_id + + integer(i_def), intent(in) :: element_order_h, element_order_v + integer(i_def), intent(in) :: geometry, topology + integer(i_def), intent(in) :: coord_system + + real(r_def), intent(in) :: f_lat, f_lon + real(r_def), intent(in) :: scaled_radius + type(mesh_type), pointer :: mesh type(local_mesh_type), pointer :: local_mesh type(inventory_by_local_mesh_type), pointer :: long_inventory @@ -1011,6 +1123,8 @@ contains logical(kind=l_def) :: constant_exists character(len=str_def) :: inventory_name + logical(l_def), parameter :: use_fe = .false. + ! NB: Longitude and latitude fields are computed simultaneously ! Determine inventory based on space select case (space_id) @@ -1043,24 +1157,47 @@ contains constant_exists = long_inventory%paired_object_exists(local_mesh%get_id()) if (.not. constant_exists) then - call compute_latlon(long_inventory, lat_inventory, mesh, space_id, & - use_fe=.false.) + call compute_latlon( long_inventory, lat_inventory, & + mesh, space_id, use_fe, geometry, topology, & + element_order_h, element_order_v, & + coord_system, f_lat, f_lon, scaled_radius ) end if call long_inventory%get_field(local_mesh, long_ptr) end function get_longitude_fv + !> @brief Returns a pointer to the latitude of finite element DoFs - !> @param[in] space_id The space for which to get the latitude of DoFs for - !> @param[in] mesh_id The ID of the mesh to get the object for + !> @param[in] space_id The space for which to get the latitude of DoFs for + !> @param[in] mesh_id The ID of the mesh to get the object for + !> @param[in] geometry Mesh geometry enumeration value + !> @param[in] topology Mesh topology enumeration value + !> @param[in] element_order_h Function space order in horizontal + !> @param[in] element_order_v Function space order in vertical + !> @param[in] coord_system Finite-Element coord-system enumeration value + !> @param[in] f_lat Latitiude of f-plane + !> @param[in] f_lon Longitude of f-plane + !> @param[in] scaled_radius Planet scaled radius !> @return The latitude field - function get_latitude_fe(space_id, mesh_id) result(lat_ptr) + function get_latitude_fe( space_id, mesh_id, geometry, topology, & + element_order_h, element_order_v, & + coord_system, f_lat, f_lon, scaled_radius ) & + result( lat_ptr ) + implicit none - integer(kind=i_def), intent(in) :: space_id - integer(kind=i_def), intent(in) :: mesh_id + integer(i_def), intent(in) :: space_id + integer(i_def), intent(in) :: mesh_id + + integer(i_def), intent(in) :: element_order_h, element_order_v + integer(i_def), intent(in) :: geometry, topology + integer(i_def), intent(in) :: coord_system + + real(r_def), intent(in) :: f_lat, f_lon + real(r_def), intent(in) :: scaled_radius + type(mesh_type), pointer :: mesh type(local_mesh_type), pointer :: local_mesh type(inventory_by_local_mesh_type), pointer :: long_inventory @@ -1069,9 +1206,13 @@ contains logical(kind=l_def) :: constant_exists character(len=str_def) :: inventory_name + logical(l_def), parameter :: use_fe = .true. + ! If running at lowest order, use finite volume if (element_order_h == 0 .and. element_order_v == 0) then - lat_ptr => get_latitude_fv(space_id, mesh_id) + lat_ptr => get_latitude_fv( space_id, mesh_id, geometry, topology, & + element_order_h, element_order_v, & + coord_system, f_lat, f_lon, scaled_radius ) return end if @@ -1107,24 +1248,46 @@ contains constant_exists = lat_inventory%paired_object_exists(local_mesh%get_id()) if (.not. constant_exists) then - call compute_latlon(long_inventory, lat_inventory, mesh, space_id, & - use_fe=.true.) + call compute_latlon( long_inventory, lat_inventory, & + mesh, space_id, use_fe, geometry, topology, & + element_order_h, element_order_v, & + coord_system, f_lat, f_lon, scaled_radius ) end if call lat_inventory%get_field(local_mesh, lat_ptr) end function get_latitude_fe + !> @brief Returns a pointer to the latitude of finite volume DoFs - !> @param[in] space_id The space for which to get the latitude of DoFs for - !> @param[in] mesh_id The ID of the mesh to get the object for + !> @param[in] space_id The space for which to get the latitude of DoFs for + !> @param[in] mesh_id The ID of the mesh to get the object for + !> @param[in] geometry Mesh geometry enumeration value + !> @param[in] topology Mesh topology enumeration value + !> @param[in] element_order_h Function space order in horizontal + !> @param[in] element_order_v Function space order in vertical + !> @param[in] coord_system Finite-Element coord-system enumeration value + !> @param[in] f_lat Latitiude of f-plane + !> @param[in] f_lon Longitude of f-plane + !> @param[in] scaled_radius Planet scaled radius !> @return The latitude field - function get_latitude_fv(space_id, mesh_id) result(lat_ptr) + function get_latitude_fv( space_id, mesh_id, geometry, topology, & + element_order_h, element_order_v, & + coord_system, f_lat, f_lon, scaled_radius ) & + result( lat_ptr ) implicit none - integer(kind=i_def), intent(in) :: space_id - integer(kind=i_def), intent(in) :: mesh_id + integer(i_def), intent(in) :: space_id + integer(i_def), intent(in) :: mesh_id + + integer(i_def), intent(in) :: element_order_h, element_order_v + integer(i_def), intent(in) :: geometry, topology + integer(i_def), intent(in) :: coord_system + + real(r_def), intent(in) :: f_lat, f_lon + real(r_def), intent(in) :: scaled_radius + type(mesh_type), pointer :: mesh type(local_mesh_type), pointer :: local_mesh type(inventory_by_local_mesh_type), pointer :: long_inventory @@ -1133,6 +1296,8 @@ contains logical(kind=l_def) :: constant_exists character(len=str_def) :: inventory_name + logical(l_def), parameter :: use_fe = .false. + ! NB: Longitude and latitude fields are computed simultaneously ! Determine inventory based on space select case (space_id) @@ -1165,8 +1330,10 @@ contains constant_exists = lat_inventory%paired_object_exists(local_mesh%get_id()) if (.not. constant_exists) then - call compute_latlon(long_inventory, lat_inventory, mesh, space_id, & - use_fe=.false.) + call compute_latlon( long_inventory, lat_inventory, & + mesh, space_id, use_fe, geometry, topology, & + element_order_h, element_order_v, & + coord_system, f_lat, f_lon, scaled_radius ) end if call lat_inventory%get_field(local_mesh, lat_ptr) @@ -1174,22 +1341,36 @@ contains end function get_latitude_fv !> @brief Returns a pointer to a finite element height field - !> @param[in] space The space of the desired height field - !> @param[in] mesh_id The ID of the mesh to get the object for + !> @param[in] space_id The space of the desired height field + !> @param[in] mesh_id The ID of the mesh to get the object for + !> @param[in] geometry Mesh geometry enumeration value + !> @param[in] element_order_h Function space order in horizontal + !> @param[in] element_order_v Function space order in vertical + !> @param[in] coord_system Finite-Element coord-system enumeration value + !> @param[in] scaled_radius Planet scaled radius !> @return A height field - function get_height_fe(space_id, mesh_id) result(height) + function get_height_fe( space_id, mesh_id, geometry, & + element_order_h, element_order_v, & + coord_system, scaled_radius ) & + result( height ) + use sci_height_continuous_kernel_mod, only: height_continuous_kernel_type use sci_height_discontinuous_kernel_mod, & only: height_discontinuous_kernel_type - use base_mesh_config_mod, only: geometry - use finite_element_config_mod, only: coord_system - use planet_config_mod, only: scaled_radius implicit none - integer(kind=i_def), intent(in) :: space_id - integer(kind=i_def), intent(in) :: mesh_id + integer(i_def), intent(in) :: space_id + integer(i_def), intent(in) :: mesh_id + + integer(i_def), intent(in) :: geometry + integer(i_def), intent(in) :: coord_system + integer(i_def), intent(in) :: element_order_h + integer(i_def), intent(in) :: element_order_v + + real(r_def), intent(in) :: scaled_radius + type(mesh_type), pointer :: mesh type(inventory_by_mesh_type), pointer :: inventory logical(kind=l_def) :: constant_exists @@ -1204,7 +1385,8 @@ contains ! If running at lowest order, use finite volume if (element_order_h == 0 .and. element_order_v == 0) then - height => get_height_fv(space_id, mesh_id) + height => get_height_fv( space_id, mesh_id, geometry, & + coord_system, scaled_radius ) return end if @@ -1247,9 +1429,10 @@ contains if ( LPROF ) call start_timing( id, 'runtime_constants.geometric' ) - space => function_space_collection%get_fs( & - mesh, element_order_h, element_order_v, space_id & - ) + space => function_space_collection%get_fs( mesh, & + element_order_h, & + element_order_v, & + space_id ) call inventory%add_field(height, space, mesh) select case (space_id) @@ -1289,22 +1472,30 @@ contains end function get_height_fe !> @brief Returns a pointer to a finite volume height field - !> @param[in] space The space of the desired height field - !> @param[in] mesh_id The ID of the mesh to get the object for + !> @param[in] space_id The space of the desired height field + !> @param[in] mesh_id The ID of the mesh to get the object for + !> @param[in] geometry Mesh geometry enumeration value + !> @param[in] coord_system Finite-Element coord-system enumeration value + !> @param[in] scaled_radius Planet scaled radius !> @return A height field - function get_height_fv(space_id, mesh_id) result(height) + function get_height_fv( space_id, mesh_id, geometry, & + coord_system, scaled_radius ) & + result( height ) + use sci_height_continuous_kernel_mod, only: height_continuous_kernel_type - use sci_height_discontinuous_kernel_mod, & + use sci_height_discontinuous_kernel_mod, & only: height_discontinuous_kernel_type - use base_mesh_config_mod, only: geometry - use finite_element_config_mod, only: coord_system - use planet_config_mod, only: scaled_radius implicit none - integer(kind=i_def), intent(in) :: space_id - integer(kind=i_def), intent(in) :: mesh_id + integer(i_def), intent(in) :: space_id + integer(i_def), intent(in) :: mesh_id + + integer(i_def), intent(in) :: geometry + integer(i_def), intent(in) :: coord_system + real(r_def), intent(in) :: scaled_radius + type(mesh_type), pointer :: mesh type(inventory_by_mesh_type), pointer :: inventory logical(kind=l_def) :: constant_exists diff --git a/components/science/source/algorithm/sci_mapping_constants_mod.x90 b/components/science/source/algorithm/sci_mapping_constants_mod.x90 index 594b26063..32f5d86a9 100644 --- a/components/science/source/algorithm/sci_mapping_constants_mod.x90 +++ b/components/science/source/algorithm/sci_mapping_constants_mod.x90 @@ -198,8 +198,13 @@ contains !> @brief Create the operators for projecting spherical components in !! (W3, W3, Wtheta) to a vector-valued field in W2 - !> @param[in] mesh The mesh to compute the operators for - subroutine create_spherical_components_to_w2_projection(mesh) + !> @param[in] mesh The mesh to compute the operators for + !> @param[in] geometry Mesh geometry enumeration value + !> @param[in] topology Mesh topology enumeration value + !> @param[in] coord_system Finite-Element coord-system enumeration value + !> @param[in] scaled_radius Planet scaled radius + subroutine create_spherical_components_to_w2_projection( mesh, geometry, topology, & + coord_system, scaled_radius ) use sci_compute_map_u_operators_kernel_mod, & only: compute_map_u_operators_kernel_type @@ -207,6 +212,12 @@ contains implicit none type(mesh_type), pointer, intent(in) :: mesh + + integer(i_def), intent(in) :: geometry + integer(i_def), intent(in) :: topology + integer(i_def), intent(in) :: coord_system + real(r_def), intent(in) :: scaled_radius + integer(kind=i_def) :: mesh_id type(function_space_type), pointer :: w2_fs type(function_space_type), pointer :: w3_fs @@ -242,19 +253,27 @@ contains call u_lat_map_inventory%add_operator(u_lat_map, w2_fs, w3_fs, mesh) call u_up_map_inventory%add_operator(u_up_map, w2_fs, wtheta_fs, mesh) - call invoke( name="compute_lonlatr_galerkin_operators", & - compute_map_u_operators_kernel_type(u_lon_map, & - u_lat_map, & - u_up_map, & - chi, panel_id, & + call invoke( name="compute_lonlatr_galerkin_operators", & + compute_map_u_operators_kernel_type(u_lon_map, & + u_lat_map, & + u_up_map, & + chi, panel_id, & + geometry, topology, & + coord_system, & + scaled_radius, & qr) ) end subroutine create_spherical_components_to_w2_projection !> @brief Create the operators for sampling spherical components in !! (W3, W3, Wtheta) to a vector-valued field in W2 - !> @param[in] mesh The mesh to compute the operators for - subroutine create_spherical_components_to_w2_sample(mesh) + !> @param[in] mesh The mesh to compute the operators for + !> @param[in] geometry Mesh geometry enumeration value + !> @param[in] topology Mesh topology enumeration value + !> @param[in] coord_system Finite-Element coord-system enumeration value + !> @param[in] scaled_radius Planet scaled radius + subroutine create_spherical_components_to_w2_sample( mesh, geometry, topology, & + coord_system, scaled_radius ) use sci_compute_sample_u_ops_kernel_mod, & only: compute_sample_u_ops_kernel_type @@ -262,6 +281,12 @@ contains implicit none type(mesh_type), pointer, intent(in) :: mesh + + integer(i_def), intent(in) :: geometry + integer(i_def), intent(in) :: topology + integer(i_def), intent(in) :: coord_system + real(r_def), intent(in) :: scaled_radius + integer(kind=i_def) :: mesh_id type(function_space_type), pointer :: w2_fs type(function_space_type), pointer :: w3_fs @@ -298,11 +323,14 @@ contains call u_lat_sample_inventory%add_operator(u_lat_sample, w2_fs, w3_fs, mesh) call u_up_sample_inventory%add_operator(u_up_sample, w2_fs, wtheta_fs, mesh) - call invoke( name="compute_lonlatr_sample_operators", & - compute_sample_u_ops_kernel_type(u_lon_sample, & - u_lat_sample, & - u_up_sample, & - chi, panel_id) ) + call invoke( name="compute_lonlatr_sample_operators", & + compute_sample_u_ops_kernel_type(u_lon_sample, & + u_lat_sample, & + u_up_sample, & + chi, panel_id, & + geometry, topology, & + coord_system, & + scaled_radius) ) if ( LPROF ) call stop_timing( id, 'runtime_constants.mapping' ) @@ -845,13 +873,24 @@ contains end function get_project_zdot_to_w2 !> @brief Returns a pointer to the u_lon mapping operator - !> @param[in] mesh_id The ID of the mesh to get the object for + !> @param[in] mesh_id The ID of the mesh to get the object for + !> @param[in] geometry Mesh geometry enumeration value + !> @param[in] topology Mesh topology enumeration value + !> @param[in] coord_system Finite-Element coord-system enumeration value + !> @param[in] scaled_radius Planet scaled radius !> @return The mapping operator for u_lon to W2 - function get_u_lon_map(mesh_id) result(u_lon_map_op) + function get_u_lon_map( mesh_id, geometry, topology, & + coord_system, scaled_radius ) & + result( u_lon_map_op ) implicit none - integer(kind=i_def), intent(in) :: mesh_id + integer(i_def), intent(in) :: mesh_id + integer(i_def), intent(in) :: geometry + integer(i_def), intent(in) :: topology + integer(i_def), intent(in) :: coord_system + real(r_def), intent(in) :: scaled_radius + type(mesh_type), pointer :: mesh type(operator_type), pointer :: u_lon_map_op logical(kind=l_def) :: constant_exists @@ -864,7 +903,10 @@ contains mesh => mesh_collection%get_mesh(mesh_id) constant_exists = u_lon_map_inventory%paired_object_exists(mesh_id) - if (.not. constant_exists) call create_spherical_components_to_w2_projection(mesh) + if (.not. constant_exists) then + call create_spherical_components_to_w2_projection(mesh, geometry, topology, & + coord_system, scaled_radius) + end if ! Return constant call u_lon_map_inventory%get_operator(mesh, u_lon_map_op) @@ -872,13 +914,24 @@ contains end function get_u_lon_map !> @brief Returns a pointer to the u_lat mapping operator - !> @param[in] mesh_id The ID of the mesh to get the object for + !> @param[in] mesh_id The ID of the mesh to get the object for + !> @param[in] geometry Mesh geometry enumeration value + !> @param[in] topology Mesh topology enumeration value + !> @param[in] coord_system Finite-Element coord-system enumeration value + !> @param[in] scaled_radius Planet scaled radius !> @return The mapping operator for u_lat to W2 - function get_u_lat_map(mesh_id) result(u_lat_map_op) + function get_u_lat_map( mesh_id, geometry, topology, & + coord_system, scaled_radius ) & + result( u_lat_map_op ) implicit none - integer(kind=i_def), intent(in) :: mesh_id + integer(i_def), intent(in) :: mesh_id + integer(i_def), intent(in) :: geometry + integer(i_def), intent(in) :: topology + integer(i_def), intent(in) :: coord_system + real(r_def), intent(in) :: scaled_radius + type(mesh_type), pointer :: mesh type(operator_type), pointer :: u_lat_map_op logical(kind=l_def) :: constant_exists @@ -891,7 +944,10 @@ contains mesh => mesh_collection%get_mesh(mesh_id) constant_exists = u_lat_map_inventory%paired_object_exists(mesh_id) - if (.not. constant_exists) call create_spherical_components_to_w2_projection(mesh) + if (.not. constant_exists) then + call create_spherical_components_to_w2_projection(mesh, geometry, topology, & + coord_system, scaled_radius) + end if ! Return constant call u_lat_map_inventory%get_operator(mesh, u_lat_map_op) @@ -899,13 +955,24 @@ contains end function get_u_lat_map !> @brief Returns a pointer to the u_up mapping operator - !> @param[in] mesh_id The ID of the mesh to get the object for + !> @param[in] mesh_id The ID of the mesh to get the object for + !> @param[in] geometry Mesh geometry enumeration value + !> @param[in] topology Mesh topology enumeration value + !> @param[in] coord_system Finite-Element coord-system enumeration value + !> @param[in] scaled_radius Planet scaled radius !> @return The mapping operator for u_up to W2 - function get_u_up_map(mesh_id) result(u_up_map_op) + function get_u_up_map( mesh_id, geometry, topology, & + coord_system, scaled_radius ) & + result( u_up_map_op ) implicit none - integer(kind=i_def), intent(in) :: mesh_id + integer(i_def), intent(in) :: mesh_id + integer(i_def), intent(in) :: geometry + integer(i_def), intent(in) :: topology + integer(i_def), intent(in) :: coord_system + real(r_def), intent(in) :: scaled_radius + type(mesh_type), pointer :: mesh type(operator_type), pointer :: u_up_map_op logical(kind=l_def) :: constant_exists @@ -918,7 +985,10 @@ contains mesh => mesh_collection%get_mesh(mesh_id) constant_exists = u_up_map_inventory%paired_object_exists(mesh_id) - if (.not. constant_exists) call create_spherical_components_to_w2_projection(mesh) + if (.not. constant_exists) then + call create_spherical_components_to_w2_projection(mesh, geometry, topology, & + coord_system, scaled_radius) + end if ! Return constant call u_up_map_inventory%get_operator(mesh, u_up_map_op) @@ -926,13 +996,24 @@ contains end function get_u_up_map !> @brief Returns a pointer to the u_lon sampling operator - !> @param[in] mesh_id The ID of the mesh to get the object for + !> @param[in] mesh_id The ID of the mesh to get the object for + !> @param[in] geometry Mesh geometry enumeration value + !> @param[in] topology Mesh topology enumeration value + !> @param[in] coord_system Finite-Element coord-system enumeration value + !> @param[in] scaled_radius Planet scaled radius !> @return The sampling operator for u_lon to W2 - function get_u_lon_sample(mesh_id) result(u_lon_sample_op) + function get_u_lon_sample( mesh_id, geometry, topology, & + coord_system, scaled_radius ) & + result( u_lon_sample_op ) implicit none - integer(kind=i_def), intent(in) :: mesh_id + integer(i_def), intent(in) :: mesh_id + integer(i_def), intent(in) :: geometry + integer(i_def), intent(in) :: topology + integer(i_def), intent(in) :: coord_system + real(r_def), intent(in) :: scaled_radius + type(mesh_type), pointer :: mesh type(operator_type), pointer :: u_lon_sample_op logical(kind=l_def) :: constant_exists @@ -945,7 +1026,10 @@ contains mesh => mesh_collection%get_mesh(mesh_id) constant_exists = u_lon_sample_inventory%paired_object_exists(mesh_id) - if (.not. constant_exists) call create_spherical_components_to_w2_sample(mesh) + if (.not. constant_exists) then + call create_spherical_components_to_w2_sample(mesh, geometry, topology, & + coord_system, scaled_radius) + end if ! Return constant call u_lon_sample_inventory%get_operator(mesh, u_lon_sample_op) @@ -953,13 +1037,24 @@ contains end function get_u_lon_sample !> @brief Returns a pointer to the u_lat sampling operator - !> @param[in] mesh_id The ID of the mesh to get the object for + !> @param[in] mesh_id The ID of the mesh to get the object for + !> @param[in] geometry Mesh geometry enumeration value + !> @param[in] topology Mesh topology enumeration value + !> @param[in] coord_system Finite-Element coord-system enumeration value + !> @param[in] scaled_radius Planet scaled radius !> @return The sampling operator for u_lat to W2 - function get_u_lat_sample(mesh_id) result(u_lat_sample_op) + function get_u_lat_sample( mesh_id, geometry, topology, & + coord_system, scaled_radius ) & + result( u_lat_sample_op ) implicit none - integer(kind=i_def), intent(in) :: mesh_id + integer(i_def), intent(in) :: mesh_id + integer(i_def), intent(in) :: geometry + integer(i_def), intent(in) :: topology + integer(i_def), intent(in) :: coord_system + real(r_def), intent(in) :: scaled_radius + type(mesh_type), pointer :: mesh type(operator_type), pointer :: u_lat_sample_op logical(kind=l_def) :: constant_exists @@ -972,7 +1067,10 @@ contains mesh => mesh_collection%get_mesh(mesh_id) constant_exists = u_lat_sample_inventory%paired_object_exists(mesh_id) - if (.not. constant_exists) call create_spherical_components_to_w2_sample(mesh) + if (.not. constant_exists) then + call create_spherical_components_to_w2_sample(mesh, geometry, topology, & + coord_system, scaled_radius) + end if ! Return constant call u_lat_sample_inventory%get_operator(mesh, u_lat_sample_op) @@ -980,13 +1078,24 @@ contains end function get_u_lat_sample !> @brief Returns a pointer to the u_up sampling operator - !> @param[in] mesh_id The ID of the mesh to get the object for + !> @param[in] mesh_id The ID of the mesh to get the object for + !> @param[in] geometry Mesh geometry enumeration value + !> @param[in] topology Mesh topology enumeration value + !> @param[in] coord_system Finite-Element coord-system enumeration value + !> @param[in] scaled_radius Planet scaled radius !> @return The sampling operator for u_up to W2 - function get_u_up_sample(mesh_id) result(u_up_sample_op) + function get_u_up_sample( mesh_id, geometry, topology, & + coord_system, scaled_radius) & + result( u_up_sample_op ) implicit none - integer(kind=i_def), intent(in) :: mesh_id + integer(i_def), intent(in) :: mesh_id + integer(i_def), intent(in) :: geometry + integer(i_def), intent(in) :: topology + integer(i_def), intent(in) :: coord_system + real(r_def), intent(in) :: scaled_radius + type(mesh_type), pointer :: mesh type(operator_type), pointer :: u_up_sample_op logical(kind=l_def) :: constant_exists @@ -999,7 +1108,10 @@ contains mesh => mesh_collection%get_mesh(mesh_id) constant_exists = u_up_sample_inventory%paired_object_exists(mesh_id) - if (.not. constant_exists) call create_spherical_components_to_w2_sample(mesh) + if (.not. constant_exists) then + call create_spherical_components_to_w2_sample( mesh, geometry, topology, & + coord_system, scaled_radius ) + end if ! Return constant call u_up_sample_inventory%get_operator(mesh, u_up_sample_op) @@ -1007,13 +1119,24 @@ contains end function get_u_up_sample !> @brief Returns a pointer to the operator projection from lon dot to W1 - !> @param[in] mesh_id The ID of the mesh to get the object for + !> @param[in] mesh_id The ID of the mesh to get the object for + !> @param[in] geometry Mesh geometry enumeration value + !> @param[in] topology Mesh topology enumeration value + !> @param[in] coord_system Finite-Element coord-system enumeration value + !> @param[in] scaled_radius Planet scaled radius !> @return The projection operator - function get_project_lon_dot_to_w1(mesh_id) result(proj_op) + function get_project_lon_dot_to_w1( mesh_id, geometry, topology, & + coord_system, scaled_radius ) & + result( proj_op ) implicit none - integer(kind=i_def), intent(in) :: mesh_id + integer(i_def), intent(in) :: mesh_id + integer(i_def), intent(in) :: geometry + integer(i_def), intent(in) :: topology + integer(i_def), intent(in) :: coord_system + real(r_def), intent(in) :: scaled_radius + type(mesh_type), pointer :: mesh type(operator_type), pointer :: proj_op logical(kind=l_def) :: constant_exists @@ -1050,10 +1173,14 @@ contains proj_op, w1_fs, w3_fs, mesh & ) - call invoke( name='proj_lon_dot_to_w1_op', & - project_ws_to_w1_operator_kernel_type(proj_op, & - chi, panel_id, & - xdirection, qr) ) + call invoke( name='proj_lon_dot_to_w1_op', & + project_ws_to_w1_operator_kernel_type(proj_op, & + chi, panel_id, & + xdirection, & + geometry, topology, & + coord_system, & + scaled_radius, & + qr) ) if ( LPROF ) call stop_timing( id, 'runtime_constants.mapping' ) end if @@ -1064,13 +1191,24 @@ contains end function get_project_lon_dot_to_w1 !> @brief Returns a pointer to the operator projection from lat dot to W1 - !> @param[in] mesh_id The ID of the mesh to get the object for + !> @param[in] mesh_id The ID of the mesh to get the object for + !> @param[in] geometry Mesh geometry enumeration value + !> @param[in] topology Mesh topology enumeration value + !> @param[in] coord_system Finite-Element coord-system enumeration value + !> @param[in] scaled_radius Planet scaled radius !> @return The projection operator - function get_project_lat_dot_to_w1(mesh_id) result(proj_op) + function get_project_lat_dot_to_w1( mesh_id, geometry, topology, & + coord_system, scaled_radius ) & + result( proj_op ) implicit none - integer(kind=i_def), intent(in) :: mesh_id + integer(i_def), intent(in) :: mesh_id + integer(i_def), intent(in) :: geometry + integer(i_def), intent(in) :: topology + integer(i_def), intent(in) :: coord_system + real(r_def), intent(in) :: scaled_radius + type(mesh_type), pointer :: mesh type(operator_type), pointer :: proj_op logical(kind=l_def) :: constant_exists @@ -1107,10 +1245,14 @@ contains proj_op, w1_fs, w3_fs, mesh & ) - call invoke( name='proj_lat_dot_to_w1_op', & - project_ws_to_w1_operator_kernel_type(proj_op, & - chi, panel_id, & - ydirection, qr) ) + call invoke( name='proj_lat_dot_to_w1_op', & + project_ws_to_w1_operator_kernel_type(proj_op, & + chi, panel_id, & + ydirection, & + geometry, topology, & + coord_system, & + scaled_radius, & + qr) ) if ( LPROF ) call stop_timing( id, 'runtime_constants.mapping' ) end if @@ -1121,14 +1263,25 @@ contains end function get_project_lat_dot_to_w1 !> @brief Returns a pointer to the operator projection from r dot to W1 - !> @param[in] mesh_id The ID of the mesh to get the object for + !> @param[in] mesh_id The ID of the mesh to get the object for + !> @param[in] geometry Mesh geometry enumeration value + !> @param[in] topology Mesh topology enumeration value + !> @param[in] coord_system Finite-Element coord-system enumeration value + !> @param[in] scaled_radius Planet scaled radius !> @return The projection operator - function get_project_r_dot_to_w1(mesh_id) result(proj_op) + function get_project_r_dot_to_w1( mesh_id, geometry, topology, & + coord_system, scaled_radius ) & + result( proj_op ) implicit none - integer(kind=i_def), intent(in) :: mesh_id - type(mesh_type), pointer :: mesh + integer(i_def), intent(in) :: mesh_id + integer(i_def), intent(in) :: geometry + integer(i_def), intent(in) :: topology + integer(i_def), intent(in) :: coord_system + real(r_def), intent(in) :: scaled_radius + + type(mesh_type), pointer :: mesh type(operator_type), pointer :: proj_op logical(kind=l_def) :: constant_exists type(function_space_type), pointer :: w1_fs @@ -1162,10 +1315,14 @@ contains proj_op, w1_fs, w3_fs, mesh & ) - call invoke( name='proj_r_dot_to_w1_op', & - project_ws_to_w1_operator_kernel_type(proj_op, & - chi, panel_id, & - zdirection, qr) ) + call invoke( name='proj_r_dot_to_w1_op', & + project_ws_to_w1_operator_kernel_type(proj_op, & + chi, panel_id, & + zdirection, & + geometry, topology, & + coord_system, & + scaled_radius, & + qr) ) if ( LPROF ) call stop_timing( id, 'runtime_constants.mapping' ) end if @@ -1176,15 +1333,26 @@ contains end function get_project_r_dot_to_w1 !> @brief Returns the displacement when averaging from W3 to W2 - !> @param[in] mesh_id The ID of the mesh to get the object for + !> @param[in] mesh_id The ID of the mesh to get the object for + !> @param[in] geometry Mesh geometry enumeration value + !> @param[in] topology Mesh topology enumeration value + !> @param[in] coord_system Finite-Element coord-system enumeration value + !> @param[in] scaled_radius Planet scaled radius !> @return The displacement field used for correcting mappings from W3 to W2 - function get_w3_to_w2_displacement(mesh_id) result(w3_to_w2_displacement) + function get_w3_to_w2_displacement( mesh_id, geometry, topology, & + coord_system, scaled_radius ) & + result( w3_to_w2_displacement ) use sci_w3_to_w2_displacement_kernel_mod, & only: w3_to_w2_displacement_kernel_type implicit none - integer(kind=i_def), intent(in) :: mesh_id + integer(i_def), intent(in) :: mesh_id + integer(i_def), intent(in) :: geometry + integer(i_def), intent(in) :: topology + integer(i_def), intent(in) :: coord_system + real(r_def), intent(in) :: scaled_radius + type(mesh_type), pointer :: mesh type(local_mesh_type), pointer :: local_mesh type(field_type), pointer :: w3_to_w2_displacement @@ -1226,9 +1394,12 @@ contains w2h_k0_fs, local_mesh) call dummy_w3%initialise( w3_k0_fs ) - call invoke( setval_c(w3_to_w2_displacement, 0.0_r_def), & - w3_to_w2_displacement_kernel_type(w3_to_w2_displacement, & - chi, panel_id, dummy_w3) ) + call invoke( setval_c(w3_to_w2_displacement, 0.0_r_def), & + w3_to_w2_displacement_kernel_type(w3_to_w2_displacement, & + chi, panel_id, dummy_w3, & + geometry, topology, & + coord_system, & + scaled_radius) ) if ( LPROF ) call stop_timing( id, 'runtime_constants.mapping' ) end if diff --git a/components/science/source/kernel/fem/sci_gp_vector_rhs_kernel_mod.F90 b/components/science/source/kernel/fem/sci_gp_vector_rhs_kernel_mod.F90 index 0db2fbe49..38a827c21 100644 --- a/components/science/source/kernel/fem/sci_gp_vector_rhs_kernel_mod.F90 +++ b/components/science/source/kernel/fem/sci_gp_vector_rhs_kernel_mod.F90 @@ -220,7 +220,8 @@ subroutine gp_vector_rhs_code(nlayers, & end do ! Obtain (X,Y,Z) coordinates for converting components of u - call chi2xyz(coords(1), coords(2), coords(3), ipanel, & + call chi2xyz(coords(1), coords(2), coords(3), ipanel, & + geometry, topology, coord_system, scaled_radius, & x_at_quad(1), x_at_quad(2), x_at_quad(3)) u_physical(:) = cart2sphere_vector(x_at_quad, u_at_quad) diff --git a/components/science/source/kernel/geometry/sci_chi_transform_mod.F90 b/components/science/source/kernel/geometry/sci_chi_transform_mod.F90 index 39ccfe0fb..d939d8c6f 100644 --- a/components/science/source/kernel/geometry/sci_chi_transform_mod.F90 +++ b/components/science/source/kernel/geometry/sci_chi_transform_mod.F90 @@ -30,15 +30,12 @@ module sci_chi_transform_mod LOG_LEVEL_WARNING use matrix_invert_mod, only : matrix_invert_3x3 -use base_mesh_config_mod, only : geometry, & - geometry_spherical, & - geometry_planar, & - topology, & - topology_fully_periodic -use finite_element_config_mod, only : coord_system, & - coord_system_xyz, & - coord_system_native -use planet_config_mod, only : scaled_radius +! Configuration modules +use base_mesh_config_mod, only: geometry_spherical, & + geometry_planar, & + topology_fully_periodic +use finite_element_config_mod, only: coord_system_xyz, & + coord_system_native implicit none @@ -77,6 +74,8 @@ module sci_chi_transform_mod !------------------------------------------------------------------------------ !> @brief Initialise the coordinate transform information !! +!> @param[in] geometry Mesh geometry enumeration value +!> @param[in] topology Mesh topology enumeration value !> @param[in] mesh_collection Optional: a collection of meshes, which contain !! metadata used to determine the rotation matrix !! and stretching factors. @@ -92,11 +91,12 @@ module sci_chi_transform_mod !------------------------------------------------------------------------------ subroutine init_chi_transforms( geometry, topology, & mesh_collection, & - north_pole_arg, equator_lat_arg ) + north_pole_arg, & + equator_lat_arg ) - use local_mesh_mod, only : local_mesh_type - use mesh_collection_mod, only : mesh_collection_type - use mesh_mod, only : mesh_type + use local_mesh_mod, only: local_mesh_type + use mesh_collection_mod, only: mesh_collection_type + use mesh_mod, only: mesh_type implicit none @@ -187,7 +187,7 @@ subroutine init_chi_transforms( geometry, topology, & LOG_LEVEL_WARNING & ) end if - end if + end if ! present(mesh_collection) if (present(north_pole_arg)) north_pole = north_pole_arg if (present(equator_lat_arg)) equatorial_latitude = equator_lat_arg @@ -243,15 +243,21 @@ end subroutine final_chi_transforms !> will be added to the height to give the radius before the coordinates !> are transformed to (X,Y,Z) coordinates. !! -!! @param[in] chi_1 The first coordinate field in -!! @param[in] chi_2 The second coordinate field in -!! @param[in] chi_3 The third coordinate field in -!! @param[in] panel_id The mesh panel ID -!! @param[out] x The first coordinate field out (global Cartesian X) -!! @param[out] y The second coordinate field out (global Cartesian Y) -!! @param[out] z The third coordinate field out (global Cartesian Z) +!! @param[in] chi_1 The first coordinate field in +!! @param[in] chi_2 The second coordinate field in +!! @param[in] chi_3 The third coordinate field in +!! @param[in] panel_id The mesh panel ID +!> @param[in] geometry Mesh geometry enumeration value +!> @param[in] topology Mesh topology enumeration value +!> @param[in] coord_system Finite-Element coord-system enumeration value +!> @param[in] scaled_radius Planet scaled radius +!! @param[out] x The first coordinate field out (global Cartesian X) +!! @param[out] y The second coordinate field out (global Cartesian Y) +!! @param[out] z The third coordinate field out (global Cartesian Z) !------------------------------------------------------------------------------- -subroutine chi2xyz(chi_1, chi_2, chi_3, panel_id, x, y, z) +subroutine chi2xyz( chi_1, chi_2, chi_3, panel_id, & + geometry, topology, coord_system, scaled_radius, & + x, y, z ) implicit none @@ -261,6 +267,11 @@ subroutine chi2xyz(chi_1, chi_2, chi_3, panel_id, x, y, z) real(kind=r_def) :: xyz(3) + integer(i_def), intent(in) :: geometry + integer(i_def), intent(in) :: topology + integer(i_def), intent(in) :: coord_system + real(r_def), intent(in) :: scaled_radius + if (geometry == geometry_planar .or. coord_system == coord_system_xyz) then ! chi already uses (geocentric) Cartesian coordinates x = chi_1 @@ -325,21 +336,30 @@ end subroutine chi2xyz !> function from chi2xyz above). Therefore this will not add the !> scaled_radius to transform. !! -!! @param[in] chi_1 The first coordinate field in -!! @param[in] chi_2 The second coordinate field in -!! @param[in] chi_3 The third coordinate field in -!! @param[in] panel_id The mesh panel ID -!! @param[out] x The first coordinate field out (global Cartesian X) -!! @param[out] y The second coordinate field out (global Cartesian Y) -!! @param[out] z The third coordinate field out (global Cartesian Z) +!! @param[in] chi_1 The first coordinate field in +!! @param[in] chi_2 The second coordinate field in +!! @param[in] chi_3 The third coordinate field in +!! @param[in] panel_id The mesh panel ID +!> @param[in] geometry Mesh geometry enumeration value +!> @param[in] topology Mesh topology enumeration value +!> @param[in] coord_system Finite-Element coord-system enumeration value +!! @param[out] x The first coordinate field out (global Cartesian X) +!! @param[out] y The second coordinate field out (global Cartesian Y) +!! @param[out] z The third coordinate field out (global Cartesian Z) !------------------------------------------------------------------------------- -subroutine chir2xyz(chi_1, chi_2, chi_3, panel_id, x, y, z) +subroutine chir2xyz( chi_1, chi_2, chi_3, panel_id, & + geometry, topology, coord_system, & + x, y, z ) implicit none - integer(kind=i_def), intent(in) :: panel_id - real(kind=r_def), intent(in) :: chi_1, chi_2, chi_3 - real(kind=r_def), intent(out) :: x, y, z + integer(kind=i_def), intent(in) :: panel_id + real(kind=r_def), intent(in) :: chi_1, chi_2, chi_3 + integer(kind=i_def), intent(in) :: geometry + integer(kind=i_def), intent(in) :: topology + integer(kind=i_def), intent(in) :: coord_system + + real(kind=r_def), intent(out) :: x, y, z real(kind=r_def) :: xyz(3) @@ -404,21 +424,32 @@ end subroutine chir2xyz !> @brief Transforms a coordinate field chi from any system into spherical polar !> (longitude, latitude, radius) coordinates !! -!! @param[in] chi_1 The first coordinate field in -!! @param[in] chi_2 The second coordinate field in -!! @param[in] chi_3 The third coordinate field in -!! @param[in] panel_id The mesh panel ID -!! @param[out] longitude The first coordinate field out (longitude) -!! @param[out] latitude The second coordinate field out (latitude) -!! @param[out] radius The third coordinate field out (radius) +!! @param[in] chi_1 The first coordinate field in +!! @param[in] chi_2 The second coordinate field in +!! @param[in] chi_3 The third coordinate field in +!! @param[in] panel_id The mesh panel ID +!> @param[in] geometry Mesh geometry enumeration value +!> @param[in] topology Mesh topology enumeration value +!> @param[in] coord_system Finite-Element coord-system enumeration value +!> @param[in] scaled_radius Planet scaled radius +!! @param[out] longitude The first coordinate field out (longitude) +!! @param[out] latitude The second coordinate field out (latitude) +!! @param[out] radius The third coordinate field out (radius) !------------------------------------------------------------------------------- -subroutine chi2llr(chi_1, chi_2, chi_3, panel_id, lon, lat, radius) +subroutine chi2llr( chi_1, chi_2, chi_3, panel_id, & + geometry, topology, coord_system, scaled_radius, & + lon, lat, radius ) implicit none - integer(kind=i_def), intent(in) :: panel_id - real(kind=r_def), intent(in) :: chi_1, chi_2, chi_3 - real(kind=r_def), intent(out) :: lon, lat, radius + integer(kind=i_def), intent(in) :: panel_id + real(kind=r_def), intent(in) :: chi_1, chi_2, chi_3 + integer(kind=i_def), intent(in) :: geometry + integer(kind=i_def), intent(in) :: topology + integer(kind=i_def), intent(in) :: coord_system + real(kind=r_def), intent(in) :: scaled_radius + + real(kind=r_def), intent(out) :: lon, lat, radius real(kind=r_def) :: xyz(3) @@ -476,28 +507,39 @@ end subroutine chi2llr !> @brief Transforms a coordinate field chi from any system into *native* !! equiangular cubed sphere (alpha,beta,radius) coordinates !! -!! @param[in] chi_1 The first coordinate field in -!! @param[in] chi_2 The second coordinate field in -!! @param[in] chi_3 The third coordinate field in -!! @param[in] panel_id The mesh panel ID -!! @param[out] alpha The first coordinate field out (alpha) -!! @param[out] beta The second coordinate field out (beta) -!! @param[out] radius The third coordinate field out (radius) +!! @param[in] chi_1 The first coordinate field in +!! @param[in] chi_2 The second coordinate field in +!! @param[in] chi_3 The third coordinate field in +!! @param[in] panel_id The mesh panel ID +!> @param[in] geometry Mesh geometry enumeration value +!> @param[in] topology Mesh topology enumeration value +!> @param[in] coord_system Finite-Element coord-system enumeration value +!> @param[in] scaled_radius Planet scaled radius +!! @param[out] alpha The first coordinate field out (alpha) +!! @param[out] beta The second coordinate field out (beta) +!! @param[out] radius The third coordinate field out (radius) !------------------------------------------------------------------------------- -subroutine chi2abr(chi_1, chi_2, chi_3, panel_id, alpha, beta, radius) +subroutine chi2abr( chi_1, chi_2, chi_3, panel_id, & + geometry, topology, coord_system, scaled_radius, & + alpha, beta, radius ) implicit none - integer(kind=i_def), intent(in) :: panel_id - real(kind=r_def), intent(in) :: chi_1, chi_2, chi_3 - real(kind=r_def), intent(out) :: alpha, beta, radius + integer(kind=i_def), intent(in) :: panel_id + real(kind=r_def), intent(in) :: chi_1, chi_2, chi_3 + integer(kind=i_def), intent(in) :: geometry + integer(kind=i_def), intent(in) :: topology + integer(kind=i_def), intent(in) :: coord_system + real(kind=r_def), intent(in) :: scaled_radius + + real(kind=r_def), intent(out) :: alpha, beta, radius real(kind=r_def) :: xyz(3) if (topology /= topology_fully_periodic .or. geometry /= geometry_spherical) then - call log_event( & - 'chi2abr can only be used on cubed-sphere meshes', LOG_LEVEL_ERROR & - ) + + call log_event( 'chi2abr can only be used on cubed-sphere meshes', & + LOG_LEVEL_ERROR ) else if (coord_system == coord_system_native) then alpha = chi_1 @@ -531,6 +573,7 @@ end subroutine chi2abr !! native Cartesian coordinates to the physical Cartesian coordinates !------------------------------------------------------------------------------- function get_mesh_rotation_matrix() result(rot_mat) + implicit none real(kind=r_def) :: rot_mat(3,3) @@ -543,6 +586,7 @@ end function get_mesh_rotation_matrix !! physical Cartesian coordinates to native Cartesian coordinates !------------------------------------------------------------------------------- function get_inverse_mesh_rotation_matrix() result(rot_mat) + implicit none real(kind=r_def) :: rot_mat(3,3) @@ -554,6 +598,7 @@ end function get_inverse_mesh_rotation_matrix !> @brief Returns the Schmidt transform stretch factor !------------------------------------------------------------------------------- function get_stretch_factor() result(stretch_factor_out) + implicit none real(kind=r_def) :: stretch_factor_out @@ -565,6 +610,7 @@ end function get_stretch_factor !> @brief Returns whether coordinates are rotated !------------------------------------------------------------------------------- function get_to_rotate() result(to_rotate_out) + implicit none logical(kind=l_def) :: to_rotate_out @@ -576,6 +622,7 @@ end function get_to_rotate !> @brief Returns whether coordinates are stretched !------------------------------------------------------------------------------- function get_to_stretch() result(to_stretch_out) + implicit none logical(kind=l_def) :: to_stretch_out diff --git a/components/science/source/kernel/geometry/sci_compute_latlon_kernel_mod.F90 b/components/science/source/kernel/geometry/sci_compute_latlon_kernel_mod.F90 index 1ca3f2776..35da135a3 100644 --- a/components/science/source/kernel/geometry/sci_compute_latlon_kernel_mod.F90 +++ b/components/science/source/kernel/geometry/sci_compute_latlon_kernel_mod.F90 @@ -8,14 +8,16 @@ module sci_compute_latlon_kernel_mod use argument_mod, only: arg_type, func_type, & - GH_FIELD, GH_REAL, & + GH_FIELD, GH_SCALAR, & + GH_INTEGER, GH_REAL, & GH_WRITE, GH_READ, & - ANY_SPACE_1, & + ANY_SPACE_1, & ANY_DISCONTINUOUS_SPACE_3, & ANY_SPACE_9, GH_BASIS, & CELL_COLUMN, GH_EVALUATOR use constants_mod, only: r_def, i_def use kernel_mod, only: kernel_type + use sci_chi_transform_mod, only: chi2llr implicit none @@ -29,12 +31,17 @@ module sci_compute_latlon_kernel_mod !> type, public, extends(kernel_type) :: compute_latlon_kernel_type private - type(arg_type) :: meta_args(4) = (/ & - arg_type(GH_FIELD, GH_REAL, GH_WRITE, ANY_SPACE_1), & - arg_type(GH_FIELD, GH_REAL, GH_WRITE, ANY_SPACE_1), & - arg_type(GH_FIELD*3, GH_REAL, GH_READ, ANY_SPACE_9), & - arg_type(GH_FIELD, GH_REAL, GH_READ, ANY_DISCONTINUOUS_SPACE_3) & + type(arg_type) :: meta_args(8) = (/ & + arg_type(GH_FIELD, GH_REAL, GH_WRITE, ANY_SPACE_1), &! latitude + arg_type(GH_FIELD, GH_REAL, GH_WRITE, ANY_SPACE_1), &! longitude + arg_type(GH_FIELD*3, GH_REAL, GH_READ, ANY_SPACE_9), &! chi_1, chi_2, chi_3 + arg_type(GH_FIELD, GH_REAL, GH_READ, ANY_DISCONTINUOUS_SPACE_3), &! panel_id + arg_type(GH_SCALAR, GH_INTEGER, GH_READ), &! geometry + arg_type(GH_SCALAR, GH_INTEGER, GH_READ), &! topology + arg_type(GH_SCALAR, GH_INTEGER, GH_READ), &! coord_system + arg_type(GH_SCALAR, GH_REAL, GH_READ) &! scaled_radius /) + type(func_type) :: meta_funcs(1) = (/ & func_type(ANY_SPACE_9, GH_BASIS) & /) @@ -61,6 +68,10 @@ module sci_compute_latlon_kernel_mod !> @param[in] chi_2 Second component of the coordinate field !> @param[in] chi_3 Third component of the coordinate field !> @param[in] panel_id A field giving the ID for mesh panels +!> @param[in] geometry Mesh geometry enumeration value +!> @param[in] topology Mesh topology enumeration value +!> @param[in] coord_system Finite-Element coord-system enumeration value +!> @param[in] scaled_radius Planet scaled radius !> @param[in] ndf_x Number of degrees of freedom per cell for height !> @param[in] undf_x Number of unique degrees of freedom for height !> @param[in] map_x Dofmap for the cell at the base of the column for height @@ -71,15 +82,16 @@ module sci_compute_latlon_kernel_mod !> @param[in] ndf_pid Number of degrees of freedom per cell for panel_id !> @param[in] undf_pid Number of unique degrees of freedom for panel_id !> @param[in] map_pid Dofmap for the cell at the base of the column for panel_id -subroutine compute_latlon_code(nlayers, & - latitude, longitude, & - chi_1, chi_2, chi_3, & - panel_id, & - ndf_x, undf_x, map_x, & - ndf_chi, undf_chi, map_chi, & - basis_chi, & - ndf_pid, undf_pid, map_pid & - ) +subroutine compute_latlon_code( nlayers, & + latitude, longitude, & + chi_1, chi_2, chi_3, & + panel_id, & + geometry, topology, & + coord_system, scaled_radius, & + ndf_x, undf_x, map_x, & + ndf_chi, undf_chi, map_chi, & + basis_chi, & + ndf_pid, undf_pid, map_pid ) implicit none @@ -93,6 +105,11 @@ subroutine compute_latlon_code(nlayers, & real(kind=r_def), dimension(undf_chi), intent(in) :: chi_1, chi_2, chi_3 real(kind=r_def), dimension(undf_pid), intent(in) :: panel_id + integer(kind=i_def), intent(in) :: geometry + integer(kind=i_def), intent(in) :: topology + integer(kind=i_def), intent(in) :: coord_system + real(kind=r_def), intent(in) :: scaled_radius + integer(kind=i_def), dimension(ndf_x), intent(in) :: map_x integer(kind=i_def), dimension(ndf_chi), intent(in) :: map_chi integer(kind=i_def), dimension(ndf_pid), intent(in) :: map_pid @@ -112,7 +129,9 @@ subroutine compute_latlon_code(nlayers, & coords(2) = coords(2) + chi_2(map_chi(df_chi)+k)*basis_chi(1,df_chi,df_x) coords(3) = coords(3) + chi_3(map_chi(df_chi)+k)*basis_chi(1,df_chi,df_x) end do - call chi2llr(coords(1), coords(2), coords(3), ipanel, lon, lat, radius) + call chi2llr(coords(1), coords(2), coords(3), ipanel, & + geometry, topology, coord_system, scaled_radius, & + lon, lat, radius) latitude(map_x(df_x) + k) = lat longitude(map_x(df_x) + k) = lon end do diff --git a/components/science/source/kernel/geometry/sci_nodal_xyz_coordinates_kernel_mod.F90 b/components/science/source/kernel/geometry/sci_nodal_xyz_coordinates_kernel_mod.F90 index c3558e2d9..1d5a1b5f9 100644 --- a/components/science/source/kernel/geometry/sci_nodal_xyz_coordinates_kernel_mod.F90 +++ b/components/science/source/kernel/geometry/sci_nodal_xyz_coordinates_kernel_mod.F90 @@ -11,8 +11,9 @@ module sci_nodal_xyz_coordinates_kernel_mod use kernel_mod, only : kernel_type use argument_mod, only : arg_type, func_type, & - GH_FIELD, GH_WRITE, & - GH_READ, GH_REAL, & + GH_SCALAR, GH_FIELD, & + GH_WRITE, GH_READ, & + GH_REAL, GH_INTEGER, & ANY_SPACE_9, ANY_SPACE_1, & ANY_DISCONTINUOUS_SPACE_3, & GH_BASIS, CELL_COLUMN, & @@ -28,10 +29,14 @@ module sci_nodal_xyz_coordinates_kernel_mod ! The type declaration for the kernel. Contains the metadata needed by the Psy layer type, public, extends(kernel_type) :: nodal_xyz_coordinates_kernel_type private - type(arg_type) :: meta_args(3) = (/ & - arg_type(GH_FIELD*3, GH_REAL, GH_WRITE, ANY_SPACE_1), & - arg_type(GH_FIELD*3, GH_REAL, GH_READ, ANY_SPACE_9), & - arg_type(GH_FIELD, GH_REAL, GH_READ, ANY_DISCONTINUOUS_SPACE_3) & + type(arg_type) :: meta_args(7) = (/ & + arg_type(GH_FIELD*3, GH_REAL, GH_WRITE, ANY_SPACE_1), &! nodal_x, nodal_y, nodal_z + arg_type(GH_FIELD*3, GH_REAL, GH_READ, ANY_SPACE_9), &! chi_1, chi_2, chi_3 + arg_type(GH_FIELD, GH_REAL, GH_READ, ANY_DISCONTINUOUS_SPACE_3), &! panel_id + arg_type(GH_SCALAR, GH_INTEGER, GH_READ), &! geometry + arg_type(GH_SCALAR, GH_INTEGER, GH_READ), &! topology + arg_type(GH_SCALAR, GH_INTEGER, GH_READ), &! coord_system + arg_type(GH_SCALAR, GH_REAL, GH_READ) &! scaled_radius /) type(func_type) :: meta_funcs(1) = (/ & func_type(ANY_SPACE_9, GH_BASIS) & @@ -58,6 +63,10 @@ module sci_nodal_xyz_coordinates_kernel_mod !> @param[in] chi2 Coordinates in the second direction !> @param[in] chi3 Coordinates in the third direction !> @param[in] panel_id A field giving the ID for mesh panels. +!! @param[in] geometry +!! @param[in] topology +!! @param[in] coord_system +!! @param[in] scaled_radius !> @param[in] ndf_x Number of degrees of freedom per cell for the output field !> @param[in] undf_x Number of unique degrees of freedom for the output field !> @param[in] map_x Dofmap for the output field @@ -69,15 +78,16 @@ module sci_nodal_xyz_coordinates_kernel_mod !> @param[in] ndf_pid Number of degrees of freedom per cell for panel_id !> @param[in] undf_pid Number of unique degrees of freedom for panel_id !> @param[in] map_pid Dofmap for the panel_id field -subroutine nodal_xyz_coordinates_code(nlayers, & - nodal_x, nodal_y, nodal_z, & - chi1, chi2, chi3, & - panel_id, & - ndf_x, undf_x, map_x, & - ndf_chi, undf_chi, map_chi, & - basis_chi, & - ndf_pid, undf_pid, map_pid & - ) +subroutine nodal_xyz_coordinates_code( nlayers, & + nodal_x, nodal_y, nodal_z, & + chi1, chi2, chi3, & + panel_id, & + geometry, topology, & + coord_system, scaled_radius, & + ndf_x, undf_x, map_x, & + ndf_chi, undf_chi, map_chi, & + basis_chi, & + ndf_pid, undf_pid, map_pid ) implicit none @@ -92,6 +102,11 @@ subroutine nodal_xyz_coordinates_code(nlayers, & real(kind=r_def), dimension(undf_pid), intent(in) :: panel_id real(kind=r_def), dimension(1,ndf_chi,ndf_x), intent(in) :: basis_chi + integer(i_def), intent(in) :: geometry + integer(i_def), intent(in) :: topology + integer(i_def), intent(in) :: coord_system + real(r_def), intent(in) :: scaled_radius + ! Internal variables integer(kind=i_def) :: df_x, df_chi, k, ipanel real(kind=r_def) :: coords(3) @@ -107,11 +122,13 @@ subroutine nodal_xyz_coordinates_code(nlayers, & coords(3) = coords(3) + chi3(map_chi(df_chi)+k)*basis_chi(1,df_chi,df_x) end do - call chi2xyz(coords(1), coords(2), & - coords(3), ipanel, & - nodal_x(map_x(df_x)+k), & - nodal_y(map_x(df_x)+k), & - nodal_z(map_x(df_x)+k) ) + call chi2xyz( coords(1), coords(2), & + coords(3), ipanel, & + geometry, topology, & + coord_system, scaled_radius, & + nodal_x(map_x(df_x)+k), & + nodal_y(map_x(df_x)+k), & + nodal_z(map_x(df_x)+k) ) end do end do diff --git a/components/science/source/kernel/inter_function_space/sci_compute_map_u_operators_kernel_mod.F90 b/components/science/source/kernel/inter_function_space/sci_compute_map_u_operators_kernel_mod.F90 index a78f203d4..a1fcfdfd6 100644 --- a/components/science/source/kernel/inter_function_space/sci_compute_map_u_operators_kernel_mod.F90 +++ b/components/science/source/kernel/inter_function_space/sci_compute_map_u_operators_kernel_mod.F90 @@ -16,7 +16,8 @@ module sci_compute_map_u_operators_kernel_mod use argument_mod, only : arg_type, func_type, & - GH_FIELD, GH_REAL, & + GH_FIELD, GH_SCALAR, & + GH_REAL, GH_INTEGER, & GH_OPERATOR, & GH_INC, GH_READ, GH_WRITE, & ANY_SPACE_9, & @@ -28,6 +29,8 @@ module sci_compute_map_u_operators_kernel_mod use kernel_mod, only : kernel_type use log_mod, only : log_event, LOG_LEVEL_ERROR, LOG_LEVEL_INFO + use base_mesh_config_mod, only: geometry_spherical, geometry_planar + implicit none private @@ -40,13 +43,17 @@ module sci_compute_map_u_operators_kernel_mod !> type, public, extends(kernel_type) :: compute_map_u_operators_kernel_type private - type(arg_type) :: meta_args(5) = (/ & - arg_type(GH_OPERATOR, GH_REAL, GH_WRITE, W2, W3), & - arg_type(GH_OPERATOR, GH_REAL, GH_WRITE, W2, W3), & - arg_type(GH_OPERATOR, GH_REAL, GH_WRITE, W2, WTHETA), & - arg_type(GH_FIELD*3, GH_REAL, GH_READ, ANY_SPACE_9), & - arg_type(GH_FIELD, GH_REAL, GH_READ, ANY_DISCONTINUOUS_SPACE_3) & - /) + type(arg_type) :: meta_args(9) = (/ & + arg_type(GH_OPERATOR, GH_REAL, GH_WRITE, W2, W3), &! u_lon_op + arg_type(GH_OPERATOR, GH_REAL, GH_WRITE, W2, W3), &! u_lat_op + arg_type(GH_OPERATOR, GH_REAL, GH_WRITE, W2, WTHETA), &! u_up_op + arg_type(GH_FIELD*3, GH_REAL, GH_READ, ANY_SPACE_9), &! chi_sph_1, chi_sph_2, chi_sph_3 + arg_type(GH_FIELD, GH_REAL, GH_READ, ANY_DISCONTINUOUS_SPACE_3), &! panel_id + arg_type(GH_SCALAR, GH_INTEGER, GH_READ), &! geometry + arg_type(GH_SCALAR, GH_INTEGER, GH_READ), &! topology + arg_type(GH_SCALAR, GH_INTEGER, GH_READ), &! coord_system + arg_type(GH_SCALAR, GH_REAL, GH_READ) &! scaled_radius + /) type(func_type) :: meta_funcs(4) = (/ & func_type(W2, GH_BASIS), & func_type(W3, GH_BASIS), & @@ -80,6 +87,10 @@ module sci_compute_map_u_operators_kernel_mod !! @param[in] chi_sph_2 2nd coordinate in spherical Wchi !! @param[in] chi_sph_3 3rd coordinate in spherical Wchi !! @param[in] panel_id Field giving the ID for mesh panels +!! @param[in] geometry +!! @param[in] topology +!! @param[in] coord_system +!! @param[in] scaled_radius !! @param[in] ndf_w2 Number of degrees of freedom per cell for w2 !! @param[in] basis_w2 W2 basis functions evaluated at quadrature points !! @param[in] ndf_w3 Number of degrees of freedom per cell for w3 @@ -100,29 +111,24 @@ module sci_compute_map_u_operators_kernel_mod !! @param[in] nqp_v Number of quadrature points in the vertical !! @param[in] wqp_h Horizontal quadrature weights !! @param[in] wqp_v Vertical quadrature weights -subroutine compute_map_u_operators_code(cell, nlayers, ncell_3d_1, & - u_lon_op, ncell_3d_2, u_lat_op, & - ncell_3d_3, u_up_op, & - chi_sph_1, chi_sph_2, chi_sph_3, panel_id, & - ndf_w2, basis_w2, & - ndf_w3, basis_w3, & - ndf_wt, basis_wt, & - ndf_chi_sph, undf_chi_sph, map_chi_sph, & - chi_sph_basis, chi_sph_diff_basis, & - ndf_pid, undf_pid, map_pid, & - nqp_h, nqp_v, wqp_h, wqp_v & - ) +subroutine compute_map_u_operators_code( cell, nlayers, ncell_3d_1, & + u_lon_op, ncell_3d_2, u_lat_op, & + ncell_3d_3, u_up_op, & + chi_sph_1, chi_sph_2, chi_sph_3, panel_id, & + geometry, topology, & + coord_system, scaled_radius, & + ndf_w2, basis_w2, & + ndf_w3, basis_w3, & + ndf_wt, basis_wt, & + ndf_chi_sph, undf_chi_sph, map_chi_sph, & + chi_sph_basis, chi_sph_diff_basis, & + ndf_pid, undf_pid, map_pid, & + nqp_h, nqp_v, wqp_h, wqp_v ) use sci_chi_transform_mod, only : chi2llr use sci_coordinate_jacobian_mod, only : coordinate_jacobian use coord_transform_mod, only : sphere2cart_vector - use finite_element_config_mod, only: coord_system - use base_mesh_config_mod, only: geometry, topology, & - geometry_spherical, & - geometry_planar - use planet_config_mod, only: scaled_radius - implicit none ! Arguments @@ -151,6 +157,11 @@ subroutine compute_map_u_operators_code(cell, nlayers, ncell_3d_1, & real(kind=r_def), dimension(nqp_h), intent(in) :: wqp_h real(kind=r_def), dimension(nqp_v), intent(in) :: wqp_v + integer(i_def), intent(in) :: geometry + integer(i_def), intent(in) :: topology + integer(i_def), intent(in) :: coord_system + real(r_def), intent(in) :: scaled_radius + ! Internal variables integer(kind=i_def) :: df, df2, df3, dft integer(kind=i_def) :: k, ik, qp1, qp2, ipanel @@ -208,8 +219,10 @@ subroutine compute_map_u_operators_code(cell, nlayers, ncell_3d_1, & llr(:) = 0.0_r_def - call chi2llr(coords(1), coords(2), coords(3), & - ipanel, llr(1), llr(2), llr(3)) + call chi2llr( coords(1), coords(2), coords(3), & + ipanel, geometry, topology, & + coord_system, scaled_radius, & + llr(1), llr(2), llr(3) ) end if diff --git a/components/science/source/kernel/inter_function_space/sci_compute_sample_u_ops_kernel_mod.F90 b/components/science/source/kernel/inter_function_space/sci_compute_sample_u_ops_kernel_mod.F90 index 4bfbb79e1..4c9992d82 100644 --- a/components/science/source/kernel/inter_function_space/sci_compute_sample_u_ops_kernel_mod.F90 +++ b/components/science/source/kernel/inter_function_space/sci_compute_sample_u_ops_kernel_mod.F90 @@ -15,7 +15,8 @@ module sci_compute_sample_u_ops_kernel_mod use argument_mod, only : arg_type, func_type, & - GH_FIELD, GH_REAL, & + GH_FIELD, GH_SCALAR, & + GH_REAL, GH_INTEGER, & GH_OPERATOR, & GH_INC, GH_READ, GH_WRITE, & ANY_DISCONTINUOUS_SPACE_3, & @@ -32,11 +33,7 @@ module sci_compute_sample_u_ops_kernel_mod use coord_transform_mod, only : sphere2cart_vector use reference_element_mod, only : W, S, N, E, T, B - use finite_element_config_mod, only: coord_system - use base_mesh_config_mod, only: geometry, topology, & - geometry_spherical, & - geometry_planar - use planet_config_mod, only: scaled_radius + use base_mesh_config_mod, only: geometry_spherical, geometry_planar implicit none @@ -50,18 +47,22 @@ module sci_compute_sample_u_ops_kernel_mod !> type, public, extends(kernel_type) :: compute_sample_u_ops_kernel_type private - type(arg_type) :: meta_args(5) = (/ & - arg_type(GH_OPERATOR, GH_REAL, GH_WRITE, W2broken, W3), & - arg_type(GH_OPERATOR, GH_REAL, GH_WRITE, W2broken, W3), & - arg_type(GH_OPERATOR, GH_REAL, GH_WRITE, W2broken, WTHETA), & - arg_type(GH_FIELD*3, GH_REAL, GH_READ, Wchi), & - arg_type(GH_FIELD, GH_REAL, GH_READ, ANY_DISCONTINUOUS_SPACE_3) & + type(arg_type) :: meta_args(9) = (/ & + arg_type(GH_OPERATOR, GH_REAL, GH_WRITE, W2broken, W3), &! u_lon_op + arg_type(GH_OPERATOR, GH_REAL, GH_WRITE, W2broken, W3), &! u_lat_op + arg_type(GH_OPERATOR, GH_REAL, GH_WRITE, W2broken, WTHETA), &! u_rad_op + arg_type(GH_FIELD*3, GH_REAL, GH_READ, Wchi), &! chi_1, chi_2, chi_3 + arg_type(GH_FIELD, GH_REAL, GH_READ, ANY_DISCONTINUOUS_SPACE_3), &! panel_id + arg_type(GH_SCALAR, GH_INTEGER, GH_READ), &! geometry + arg_type(GH_SCALAR, GH_INTEGER, GH_READ), &! topology + arg_type(GH_SCALAR, GH_INTEGER, GH_READ), &! coord_system + arg_type(GH_SCALAR, GH_REAL, GH_READ) &! scaled_radius /) - type(func_type) :: meta_funcs(1) = (/ & - func_type(Wchi, GH_BASIS, GH_DIFF_BASIS) & + type(func_type) :: meta_funcs(1) = (/ & + func_type(Wchi, GH_BASIS, GH_DIFF_BASIS) & /) - type(reference_element_data_type), dimension(1) :: & - meta_reference_element = & + type(reference_element_data_type), dimension(1) :: & + meta_reference_element = & (/ reference_element_data_type(normals_to_faces) /) integer :: operates_on = CELL_COLUMN integer :: gh_shape = GH_EVALUATOR @@ -89,6 +90,10 @@ module sci_compute_sample_u_ops_kernel_mod !> @param[in] chi2 Coordinates in the second direction !> @param[in] chi3 Coordinates in the third direction !> @param[in] panel_id A field giving the ID for mesh panels +!> @param[in] geometry +!> @param[in] topology +!> @param[in] coord_system +!> @param[in] scaled_radius !> @param[in] ndf_w2b Number of DoFs per cell for broken W2 !> @param[in] ndf_w3 Number of DoFs per cell for W3 !> @param[in] ndf_wt Number of DoFs per cell for Wtheta @@ -104,18 +109,18 @@ module sci_compute_sample_u_ops_kernel_mod !> @param[in] map_pid DoF map for the column's base cell for panel ID field !> @param[in] nfaces Number of cell faces !> @param[in] face_normals The normal vectors to each face -subroutine compute_sample_u_ops_code( col, nlayers, & - ncell_3d_1, u_lon_op, & - ncell_3d_2, u_lat_op, & - ncell_3d_3, u_rad_op, & - chi1, chi2, chi3, & - panel_id, & - ndf_w2b, ndf_w3, ndf_wt, & - ndf_chi, undf_chi, map_chi, & - chi_basis, chi_diff_basis, & - ndf_pid, undf_pid, map_pid, & - nfaces, face_normals & - ) +subroutine compute_sample_u_ops_code( col, nlayers, & + ncell_3d_1, u_lon_op, & + ncell_3d_2, u_lat_op, & + ncell_3d_3, u_rad_op, & + chi1, chi2, chi3, & + panel_id, geometry, topology, & + coord_system, scaled_radius, & + ndf_w2b, ndf_w3, ndf_wt, & + ndf_chi, undf_chi, map_chi, & + chi_basis, chi_diff_basis, & + ndf_pid, undf_pid, map_pid, & + nfaces, face_normals ) implicit none @@ -143,6 +148,11 @@ subroutine compute_sample_u_ops_code( col, nlayers, & real(kind=r_def), dimension(ncell_3d_2,ndf_w2b,ndf_w3), intent(inout) :: u_lat_op real(kind=r_def), dimension(ncell_3d_3,ndf_w2b,ndf_wt), intent(inout) :: u_rad_op + integer(kind=i_def), intent(in) :: geometry + integer(kind=i_def), intent(in) :: topology + integer(kind=i_def), intent(in) :: coord_system + real(kind=r_def), intent(in) :: scaled_radius + ! Internal variables integer(kind=i_def) :: df_w2, df_wt, df_chi, k, ipanel, cell_3d real(kind=r_def), dimension(3,3,ndf_w2b) :: jacobian, jac_inv @@ -249,8 +259,9 @@ subroutine compute_sample_u_ops_code( col, nlayers, & end do ! Calculate (lon,lat,r) coordinates for W2 points in this cell - call chi2llr(chi1_w2(df_w2), chi2_w2(df_w2), chi3_w2(df_w2), ipanel, & - llr(1), llr(2), llr(3)) + call chi2llr( chi1_w2(df_w2), chi2_w2(df_w2), chi3_w2(df_w2), ipanel, & + geometry, topology, coord_system, scaled_radius, & + llr(1), llr(2), llr(3) ) ! Rotate (lon,lat,r) unit vectors to (X,Y,Z) coordinates lon_vector_xyz(df_w2,:) = sphere2cart_vector(lon_vector_llr, llr) diff --git a/components/science/source/kernel/inter_function_space/sci_convert_phys_to_hdiv_kernel_mod.F90 b/components/science/source/kernel/inter_function_space/sci_convert_phys_to_hdiv_kernel_mod.F90 index fece67189..bd2b30bfe 100644 --- a/components/science/source/kernel/inter_function_space/sci_convert_phys_to_hdiv_kernel_mod.F90 +++ b/components/science/source/kernel/inter_function_space/sci_convert_phys_to_hdiv_kernel_mod.F90 @@ -10,9 +10,10 @@ module sci_convert_phys_to_hdiv_kernel_mod use argument_mod, only : arg_type, func_type, & - GH_FIELD, GH_REAL, & + GH_FIELD, GH_SCALAR, & + GH_REAL, GH_INTEGER, & GH_WRITE, GH_READ, & - ANY_SPACE_9, GH_SCALAR, & + ANY_SPACE_9, & ANY_DISCONTINUOUS_SPACE_3, & GH_BASIS, GH_DIFF_BASIS, & CELL_COLUMN, GH_EVALUATOR, & @@ -21,6 +22,8 @@ module sci_convert_phys_to_hdiv_kernel_mod use fs_continuity_mod, only : W2 use kernel_mod, only : kernel_type + use base_mesh_config_mod, only: geometry_spherical + implicit none private @@ -33,14 +36,17 @@ module sci_convert_phys_to_hdiv_kernel_mod !> type, public, extends(kernel_type) :: convert_phys_to_hdiv_kernel_type private - type(arg_type) :: meta_args(7) = (/ & - arg_type(GH_FIELD, GH_REAL, GH_WRITE, W2), & - arg_type(GH_FIELD, GH_REAL, GH_READ, W2), & - arg_type(GH_FIELD, GH_REAL, GH_READ, W2), & - arg_type(GH_FIELD, GH_REAL, GH_READ, W2), & - arg_type(GH_FIELD*3, GH_REAL, GH_READ, ANY_SPACE_9), & - arg_type(GH_FIELD, GH_REAL, GH_READ, ANY_DISCONTINUOUS_SPACE_3), & - arg_type(GH_SCALAR, GH_INTEGER, GH_READ) & + type(arg_type) :: meta_args(10) = (/ & + arg_type(GH_FIELD, GH_REAL, GH_WRITE, W2), &! u_hdiv + arg_type(GH_FIELD, GH_REAL, GH_READ, W2), &! u_lon + arg_type(GH_FIELD, GH_REAL, GH_READ, W2), &! u_lat + arg_type(GH_FIELD, GH_REAL, GH_READ, W2), &! u_up + arg_type(GH_FIELD*3, GH_REAL, GH_READ, ANY_SPACE_9), &! chi_1, chi_2, chi_3 + arg_type(GH_FIELD, GH_REAL, GH_READ, ANY_DISCONTINUOUS_SPACE_3), &! panel_id + arg_type(GH_SCALAR, GH_INTEGER, GH_READ), &! geometry + arg_type(GH_SCALAR, GH_INTEGER, GH_READ), &! topology + arg_type(GH_SCALAR, GH_INTEGER, GH_READ), &! coord_system + arg_type(GH_SCALAR, GH_REAL, GH_READ) &! scaled_radius /) type(func_type) :: meta_funcs(2) = (/ & func_type(W2, GH_BASIS), & @@ -70,6 +76,9 @@ module sci_convert_phys_to_hdiv_kernel_mod !> @param[in] chi_3 3rd coordinate field !> @param[in] panel_id Field giving the ID for mesh panels !> @param[in] geometry Integer indicating the domain geometry +!> @param[in] topology +!> @param[in] coord_system +!> @param[in] scaled_radius !> @param[in] ndf_w2 Number of DoFs per cell for W2 !> @param[in] undf_w2 Number of DoFs for W2 for this partition !> @param[in] map_w2 Map of DoFs for lowest-layer cells for W2 @@ -93,6 +102,9 @@ subroutine convert_phys_to_hdiv_code( nlayers, & chi_3, & panel_id, & geometry, & + topology, & + coord_system, & + scaled_radius, & ndf_w2, & undf_w2, & map_w2, & @@ -111,18 +123,12 @@ subroutine convert_phys_to_hdiv_code( nlayers, & pointwise_coordinate_jacobian_inverse use coord_transform_mod, only : sphere2cart_vector - use base_mesh_config_mod, only: topology, & - geometry_spherical - use finite_element_config_mod, only: coord_system - use planet_config_mod, only: scaled_radius - implicit none ! Arguments integer(kind=i_def), intent(in) :: nlayers integer(kind=i_def), intent(in) :: ndf_w2, ndf_pid, ndf_chi integer(kind=i_def), intent(in) :: undf_w2, undf_pid, undf_chi - integer(kind=i_def), intent(in) :: geometry integer(kind=i_def), intent(in) :: map_w2(ndf_w2) integer(kind=i_def), intent(in) :: map_chi(ndf_chi) @@ -141,6 +147,11 @@ subroutine convert_phys_to_hdiv_code( nlayers, & real(kind=r_def), intent(in) :: chi_2(undf_chi) real(kind=r_def), intent(in) :: chi_3(undf_chi) + integer(i_def), intent(in) :: geometry + integer(i_def), intent(in) :: topology + integer(i_def), intent(in) :: coord_system + real(r_def), intent(in) :: scaled_radius + ! Internal variables integer(kind=i_def) :: df_w2, df_chi, k, ipanel real(kind=r_def) :: detj @@ -195,7 +206,9 @@ subroutine convert_phys_to_hdiv_code( nlayers, & ! Convert coordinates from whatever coordinate system the model uses ! into spherical-polar coordinates call chi2llr(coords(1), coords(2), coords(3), & - ipanel, llr(1), llr(2), llr(3)) + ipanel, geometry, topology, & + coord_system, scaled_radius, & + llr(1), llr(2), llr(3)) u_spherical(1) = u_lon(map_w2(df_w2) + k) u_spherical(2) = u_lat(map_w2(df_w2) + k) diff --git a/components/science/source/kernel/inter_function_space/sci_project_ws_to_w1_operator_kernel_mod.F90 b/components/science/source/kernel/inter_function_space/sci_project_ws_to_w1_operator_kernel_mod.F90 index 607c2782b..36ad67b83 100644 --- a/components/science/source/kernel/inter_function_space/sci_project_ws_to_w1_operator_kernel_mod.F90 +++ b/components/science/source/kernel/inter_function_space/sci_project_ws_to_w1_operator_kernel_mod.F90 @@ -26,6 +26,8 @@ module sci_project_ws_to_w1_operator_kernel_mod use fs_continuity_mod, only : W1, Wchi use log_mod, only : log_event, LOG_LEVEL_ERROR +use base_mesh_config_mod, only: geometry_spherical, geometry_planar + implicit none private @@ -35,12 +37,16 @@ module sci_project_ws_to_w1_operator_kernel_mod !> The type declaration for the kernel. Contains the metadata needed by the PSy layer type, public, extends(kernel_type) :: project_ws_to_w1_operator_kernel_type private - type(arg_type) :: meta_args(4) = (/ & - arg_type(GH_OPERATOR, GH_REAL, GH_WRITE, W1, ANY_DISCONTINUOUS_SPACE_1), & - arg_type(GH_FIELD*3, GH_REAL, GH_READ, Wchi), & - arg_type(GH_FIELD, GH_REAL, GH_READ, ANY_DISCONTINUOUS_SPACE_3), & - arg_type(GH_SCALAR, GH_INTEGER, GH_READ) & - /) + type(arg_type) :: meta_args(8) = (/ & + arg_type(GH_OPERATOR, GH_REAL, GH_WRITE, W1, ANY_DISCONTINUOUS_SPACE_1), &! projection_operator + arg_type(GH_FIELD*3, GH_REAL, GH_READ, Wchi), &! chi_1, chi_2, chi_3 + arg_type(GH_FIELD, GH_REAL, GH_READ, ANY_DISCONTINUOUS_SPACE_3), &! panel_id + arg_type(GH_SCALAR, GH_INTEGER, GH_READ), &! direction + arg_type(GH_SCALAR, GH_INTEGER, GH_READ), &! geometry + arg_type(GH_SCALAR, GH_INTEGER, GH_READ), &! topology + arg_type(GH_SCALAR, GH_INTEGER, GH_READ), &! coord_system + arg_type(GH_SCALAR, GH_REAL, GH_READ) &! scaled_radius + /) type(func_type) :: meta_funcs(3) = (/ & func_type(W1, GH_BASIS), & func_type(ANY_DISCONTINUOUS_SPACE_1, GH_BASIS), & @@ -73,6 +79,10 @@ module sci_project_ws_to_w1_operator_kernel_mod !> @param[in] chi3 3rd coordinate field in Wchi !> @param[in] panel_id Field giving the ID for mesh panels. !> @param[in] direction Index of the vector component (1,2 or 3) to project +!> @param[in] geometry +!> @param[in] topology +!> @param[in] coord_system +!> @param[in] scaled_radius !> @param[in] ndf_w1 Number of degrees of freedom per cell for vector space !> @param[in] basis_w1 Basis functions for the vector space at quadrature points !> @param[in] ndf_ws Number of degrees of freedom per cell for scalar space @@ -95,6 +105,10 @@ subroutine project_ws_to_w1_operator_code( cell, nlayers, & chi1, chi2, chi3, & panel_id, & direction, & + geometry, & + topology, & + coord_system, & + scaled_radius, & ndf_w1, basis_w1, & ndf_ws, basis_ws, & ndf_wx, undf_wx, map_wx, & @@ -107,12 +121,6 @@ subroutine project_ws_to_w1_operator_code( cell, nlayers, & use sci_chi_transform_mod, only: chi2llr use coord_transform_mod, only: sphere2cart_vector - use base_mesh_config_mod, only: geometry, topology, & - geometry_spherical, & - geometry_planar - use finite_element_config_mod, only: coord_system - use planet_config_mod, only: scaled_radius - implicit none ! Arguments @@ -136,6 +144,11 @@ subroutine project_ws_to_w1_operator_code( cell, nlayers, & real(kind=r_def), intent(in) :: wqp_h(nqp_h) real(kind=r_def), intent(in) :: wqp_v(nqp_v) + integer(i_def), intent(in) :: geometry + integer(i_def), intent(in) :: topology + integer(i_def), intent(in) :: coord_system + real(r_def), intent(in) :: scaled_radius + ! Internal variables integer(kind=i_def) :: df_x, df_s, df_1, ik, k, qp_h, qp_v real(kind=r_def), dimension(ndf_wx) :: chi1_e, chi2_e, chi3_e @@ -172,8 +185,9 @@ subroutine project_ws_to_w1_operator_code( cell, nlayers, & llr(:) = 0.0_r_def - call chi2llr(coords(1), coords(2), coords(3), & - ipanel, llr(1), llr(2), llr(3)) + call chi2llr(coords(1), coords(2), coords(3), ipanel, & + geometry, topology, coord_system, scaled_radius, & + llr(1), llr(2), llr(3)) end if call pointwise_coordinate_jacobian(coord_system, geometry, & diff --git a/components/science/source/kernel/inter_function_space/sci_w3_to_w2_displacement_kernel_mod.F90 b/components/science/source/kernel/inter_function_space/sci_w3_to_w2_displacement_kernel_mod.F90 index 64cc4ad24..7d6390122 100644 --- a/components/science/source/kernel/inter_function_space/sci_w3_to_w2_displacement_kernel_mod.F90 +++ b/components/science/source/kernel/inter_function_space/sci_w3_to_w2_displacement_kernel_mod.F90 @@ -12,7 +12,8 @@ module sci_w3_to_w2_displacement_kernel_mod use argument_mod, only : arg_type, func_type, & - GH_FIELD, GH_REAL, & + GH_FIELD, GH_SCALAR, & + GH_REAL, GH_INTEGER, & GH_READ, GH_INC, & ANY_DISCONTINUOUS_SPACE_3, & GH_BASIS, GH_EVALUATOR, & @@ -33,14 +34,18 @@ module sci_w3_to_w2_displacement_kernel_mod !> The type declaration for the kernel. Contains the metadata needed by the PSy layer type, public, extends(kernel_type) :: w3_to_w2_displacement_kernel_type private - type(arg_type) :: meta_args(4) = (/ & - arg_type(GH_FIELD, GH_REAL, GH_INC, W2H), & - arg_type(GH_FIELD*3, GH_REAL, GH_READ, Wchi), & - arg_type(GH_FIELD, GH_REAL, GH_READ, ANY_DISCONTINUOUS_SPACE_3), & - arg_type(GH_FIELD, GH_REAL, GH_READ, W3) & + type(arg_type) :: meta_args(8) = (/ & + arg_type(GH_FIELD, GH_REAL, GH_INC, W2H), &! displacement + arg_type(GH_FIELD*3, GH_REAL, GH_READ, Wchi), &! chi_1, chi_2, chi_3 + arg_type(GH_FIELD, GH_REAL, GH_READ, ANY_DISCONTINUOUS_SPACE_3), &! panel_id + arg_type(GH_FIELD, GH_REAL, GH_READ, W3), &! dummy_w3 + arg_type(GH_SCALAR, GH_INTEGER, GH_READ), &! geometry + arg_type(GH_SCALAR, GH_INTEGER, GH_READ), &! topology + arg_type(GH_SCALAR, GH_INTEGER, GH_READ), &! coord_system + arg_type(GH_SCALAR, GH_REAL, GH_READ) &! scaled_radius /) - type(func_type) :: meta_funcs(1) = (/ & - func_type(Wchi, GH_BASIS) & + type(func_type) :: meta_funcs(1) = (/ & + func_type(Wchi, GH_BASIS) & /) integer :: operates_on = CELL_COLUMN integer :: gh_shape = GH_EVALUATOR @@ -67,6 +72,10 @@ module sci_w3_to_w2_displacement_kernel_mod !> @param[in] chi_3 The third coordinate field !> @param[in] panel_id ID for panels of the underlying mesh !> @param[in] dummy_w3 An unused dummy field in W3 + !> @param[in] geometry + !> @param[in] topology + !> @param[in] coord_system + !> @param[in] scaled_radius !> @param[in] ndf_w2h Number of DoFs for W2H per cell !> @param[in] undf_w2h Number of unique DoFs for W2H per partition !> @param[in] map_w2h The DoF map for bottom layer cells for W2H @@ -88,6 +97,10 @@ subroutine w3_to_w2_displacement_code( nlayers, & chi_3, & panel_id, & dummy_w3, & + geometry, & + topology, & + coord_system, & + scaled_radius, & ndf_w2h, & undf_w2h, & map_w2h, & @@ -127,6 +140,11 @@ subroutine w3_to_w2_displacement_code( nlayers, & real(kind=r_def), intent(in) :: basis_chi_w2h(1,ndf_chi,ndf_w2h) real(kind=r_def), intent(in) :: basis_chi_w3(1,ndf_chi,ndf_w3) + integer(kind=i_def), intent(in) :: geometry + integer(kind=i_def), intent(in) :: topology + integer(kind=i_def), intent(in) :: coord_system + real(kind=r_def), intent(in) :: scaled_radius + ! Vertical cell index integer(kind=i_def) :: df_w2h, df_w3, df_chi integer(kind=i_def) :: ipanel @@ -156,7 +174,8 @@ subroutine w3_to_w2_displacement_code( nlayers, & chi3_at_dof = chi3_at_dof + & basis_chi_w3(1,df_chi,df_w3) * chi_3(map_chi(df_chi)) end do - call chi2abr(chi1_at_dof, chi2_at_dof, chi3_at_dof, ipanel, & + call chi2abr(chi1_at_dof, chi2_at_dof, chi3_at_dof, ipanel, & + geometry, topology, coord_system, scaled_radius, & alpha_w3, beta_w3, dummy_r) ! W2H points --------------------------------------------------------------- @@ -174,7 +193,8 @@ subroutine w3_to_w2_displacement_code( nlayers, & chi3_at_dof = chi3_at_dof + & basis_chi_w2h(1,df_chi,df_w2h) * chi_3(map_chi(df_chi)) end do - call chi2abr(chi1_at_dof, chi2_at_dof, chi3_at_dof, ipanel, & + call chi2abr(chi1_at_dof, chi2_at_dof, chi3_at_dof, ipanel, & + geometry, topology, coord_system, scaled_radius, & alpha_w2h(df_w2h), beta_w2h(df_w2h), dummy_r) end do diff --git a/components/science/unit-test/kernel/geometry/chi_transform_mod_test.pf b/components/science/unit-test/kernel/geometry/chi_transform_mod_test.pf index fd1f60bbd..6ce2dcf94 100644 --- a/components/science/unit-test/kernel/geometry/chi_transform_mod_test.pf +++ b/components/science/unit-test/kernel/geometry/chi_transform_mod_test.pf @@ -8,7 +8,15 @@ module chi_transform_mod_test use, intrinsic :: iso_fortran_env, only : real64 - use constants_mod, only : i_def, r_def, str_long, PI + use constants_mod, only : i_def, r_def, str_long, PI, rmdi + + use base_mesh_config_mod, only: geometry_spherical, & + geometry_planar, & + topology_fully_periodic, & + topology_non_periodic + + use finite_element_config_mod, only: coord_system_native, & + coord_system_xyz use funit @@ -44,6 +52,12 @@ module chi_transform_mod_test real(r_def) :: target_chi_1 real(r_def) :: target_chi_2 real(r_def) :: target_chi_3 + + integer(i_def) :: src_coord_system + integer(i_def) :: topology + integer(i_def) :: geometry + real(r_def) :: scaled_radius + contains procedure setUp procedure tearDown @@ -51,7 +65,6 @@ module chi_transform_mod_test end type chi_transform_mod_test_type ! Add my own parameters for the different coordinate system cases - ! This is done here because it is before the feign_config is initialised integer(i_def), parameter :: ABH = 1 integer(i_def), parameter :: LLH = 2 integer(i_def), parameter :: XYZ = 3 @@ -69,7 +82,7 @@ contains implicit none type(chi_parameters_type), intent(in) :: test_parameter - type(chi_transform_mod_test_type) :: new_test + type(chi_transform_mod_test_type) :: new_test new_test%source_coord_system = test_parameter%source_coord_system new_test%target_coord_system = test_parameter%target_coord_system @@ -90,30 +103,30 @@ contains class( chi_parameters_type ), intent( in ) :: this character(:), allocatable :: output_string - character(str_long) :: source_string, target_string + character(:), allocatable :: source_string, target_string select case ( this%source_coord_system ) case ( XYZ ) - write( source_string, '(A)') 'XYZ' + source_string = 'XYZ' case ( LLH ) - write( source_string, '(A)') 'LLH' + source_string = 'LLH' case ( ABH ) - write( source_string, '(A)') 'ABH' + source_string = 'ABH' case ( LLH_rot ) - write( source_string, '(A)') 'LLH rot' + source_string = 'LLH rot' case ( ABH_stretch_rot ) - write( source_string, '(A)') 'ABH stretch+rot' + source_string = 'ABH stretch+rot' end select select case ( this%target_coord_system ) case ( XYZ ) - write( target_string, '(A)') 'XYZ' + target_string = 'XYZ' case ( LLH ) - write( target_string, '(A)') 'LLH' + target_string = 'LLH' case ( ABH ) - write( target_string, '(A)') 'ABH' + target_string = 'ABH' case ( R2XYZ ) - write( target_string, '(A)') 'R2XYZ' + target_string = 'R2XYZ' end select output_string = trim( source_string // '2' // target_string ) @@ -221,82 +234,47 @@ contains !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine setUp( this ) - use base_mesh_config_mod, only : geometry_spherical, & - geometry_planar, & - topology_fully_periodic, & - topology_non_periodic - use extrusion_config_mod, only : method_uniform, & - stretching_method_linear - use finite_element_config_mod, only : cellshape_quadrilateral, & - coord_system_native, & - coord_system_xyz - use feign_config_mod, only : feign_base_mesh_config, & - feign_extrusion_config, & - feign_finite_element_config, & - feign_planet_config - use sci_chi_transform_mod, only : init_chi_transforms + use sci_chi_transform_mod, only: init_chi_transforms implicit none class(chi_transform_mod_test_type), intent(inout) :: this - integer(i_def) :: coord_system, topology - real(r_def) :: north_pole(2), equatorial_latitude - - call feign_extrusion_config( method=method_uniform, & - planet_radius=planet_radius, & - domain_height=10.0_r_def, & - number_of_layers=5_i_def, & - stretching_method=stretching_method_linear, & - stretching_height=15.0_r_def, & - eta_values=(/0.5_r_def/) ) + real(r_def) :: north_pole(2), equatorial_latitude select case ( this%source_coord_system ) case ( XYZ ) - coord_system = coord_system_xyz - topology = topology_fully_periodic + this%src_coord_system = coord_system_xyz + this%topology = topology_fully_periodic + case ( LLH, LLH_rot ) - coord_system = coord_system_native - topology = topology_non_periodic + this%src_coord_system = coord_system_native + this%topology = topology_non_periodic + case ( ABH, ABH_stretch_rot ) - coord_system = coord_system_native - topology = topology_fully_periodic + this%src_coord_system = coord_system_native + this%topology = topology_fully_periodic end select - call feign_base_mesh_config( file_prefix='foo', & - prime_mesh_name='unit_test', & - geometry=geometry_spherical, & - prepartitioned=.false., & - topology=topology, & - fplane=.false., f_lat_deg=0.0_r_def ) - - call feign_finite_element_config( & - cellshape=cellshape_quadrilateral, & - element_order_h=0_i_def, & - element_order_v=0_i_def, & - rehabilitate=.true., & - coord_order=0_i_def, & - coord_system=coord_system ) - - call feign_planet_config( scaling_factor=scaling ) + this%geometry = geometry_spherical + this%scaled_radius = planet_radius*scaling if ( this%source_coord_system == LLH_rot ) then north_pole(1) = PI/2.0_r_def north_pole(2) = 0.0_r_def - call init_chi_transforms(geometry_spherical, & - topology, & + call init_chi_transforms(this%geometry, this%topology, & north_pole_arg=north_pole) else if ( this%source_coord_system == ABH_stretch_rot ) then north_pole(1) = -PI/2.0_r_def north_pole(2) = 0.0_r_def equatorial_latitude = PI/6.0_r_def - call init_chi_transforms(geometry_spherical, & - topology, & + call init_chi_transforms(this%geometry, this%topology, & north_pole_arg=north_pole, & equator_lat_arg=equatorial_latitude) else ! Non-rotated or stretched case - call init_chi_transforms(geometry_spherical, topology) + call init_chi_transforms(this%geometry, this%topology) + end if end subroutine setUp @@ -320,8 +298,7 @@ contains @Test subroutine test_all( this ) - use sci_chi_transform_mod, only : chi2abr, chi2llr, chi2xyz, chir2xyz - use finite_element_config_mod, only : coord_system + use sci_chi_transform_mod, only: chi2abr, chi2llr, chi2xyz, chir2xyz implicit none @@ -332,16 +309,24 @@ contains select case ( this%target_coord_system ) case ( ABH ) call chi2abr(this%source_chi_1, this%source_chi_2, this%source_chi_3, & - this%panel_id, new_coord_1, new_coord_2, new_coord_3 ) + this%panel_id, this%geometry, this%topology, & + this%src_coord_system, this%scaled_radius, & + new_coord_1, new_coord_2, new_coord_3) case ( LLH ) call chi2llr(this%source_chi_1, this%source_chi_2, this%source_chi_3, & - this%panel_id, new_coord_1, new_coord_2, new_coord_3 ) + this%panel_id, this%geometry, this%topology, & + this%src_coord_system, this%scaled_radius, & + new_coord_1, new_coord_2, new_coord_3 ) case ( XYZ ) call chi2xyz(this%source_chi_1, this%source_chi_2, this%source_chi_3, & - this%panel_id, new_coord_1, new_coord_2, new_coord_3 ) + this%panel_id, this%geometry, this%topology, & + this%src_coord_system, this%scaled_radius, & + new_coord_1, new_coord_2, new_coord_3 ) case ( R2XYZ ) call chir2xyz(this%source_chi_1, this%source_chi_2, this%source_chi_3, & - this%panel_id, new_coord_1, new_coord_2, new_coord_3 ) + this%panel_id, this%geometry, this%topology, & + this%src_coord_system, & + new_coord_1, new_coord_2, new_coord_3 ) end select ! Check if answers are correct diff --git a/components/science/unit-test/kernel/geometry/compute_latlon_kernel_mod_test.pf b/components/science/unit-test/kernel/geometry/compute_latlon_kernel_mod_test.pf index ee65c0293..4b7899cbf 100644 --- a/components/science/unit-test/kernel/geometry/compute_latlon_kernel_mod_test.pf +++ b/components/science/unit-test/kernel/geometry/compute_latlon_kernel_mod_test.pf @@ -6,7 +6,7 @@ module compute_latlon_kernel_mod_test - use constants_mod, only : i_def, r_def, pi, imdi + use constants_mod, only : i_def, r_def, pi, imdi, rmdi use get_unit_test_m3x3_dofmap_mod, & only : get_w3_m3x3_dofmap, get_wchi_m3x3_dofmap use get_unit_test_m3x3_q3x3x3_sizes_mod, & @@ -16,77 +16,48 @@ module compute_latlon_kernel_mod_test use funit + use finite_element_config_mod, only: coord_system_xyz + implicit none private - public test_all + public :: set_up, tear_down, test_all - @TestCase - type, extends(TestCase), public :: compute_latlon_kernel_test_type - private - contains - procedure setUp - procedure tearDown - procedure test_all - end type compute_latlon_kernel_test_type contains !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - subroutine setUp( this ) - - use sci_chi_transform_mod, only : init_chi_transforms - use feign_config_mod, only : feign_finite_element_config - use finite_element_config_mod, only : cellshape_quadrilateral, & - coord_system_xyz - + @before + subroutine set_up() - use extrusion_config_mod, only: method_uniform + use sci_chi_transform_mod, only: init_chi_transforms implicit none - class(compute_latlon_kernel_test_type), intent(inout) :: this - - integer(kind=i_def) :: nlayers - - call feign_finite_element_config( & - cellshape=cellshape_quadrilateral, & - coord_order=0_i_def, & - coord_system=coord_system_xyz, & - element_order_h=0_i_def, & - element_order_v=0_i_def, & - rehabilitate=.true. ) - call init_chi_transforms(imdi, imdi) - end subroutine setUp + end subroutine set_up !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - subroutine tearDown( this ) + @after + subroutine tear_down() use sci_chi_transform_mod, only: final_chi_transforms - use config_loader_mod, only: final_configuration implicit none - class(compute_latlon_kernel_test_type), intent(inout) :: this - - ! Finalise namelists - call final_configuration() call final_chi_transforms() - end subroutine tearDown + end subroutine tear_down !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! @Test - subroutine test_all( this ) + subroutine test_all() use sci_compute_latlon_kernel_mod, only: compute_latlon_code implicit none - class(compute_latlon_kernel_test_type), intent(inout) :: this - real(r_def), parameter :: tol = 1.0e-12_r_def integer(i_def), parameter :: nlayers = 1 integer(i_def) :: k, df_w3 @@ -100,6 +71,11 @@ contains real(r_def), allocatable :: latitude(:), longitude(:) real(r_def), allocatable :: lat_answer(:), lon_answer(:) + integer(i_def), parameter :: geometry = imdi + integer(i_def), parameter :: topology = imdi + integer(i_def), parameter :: coord_system = coord_system_xyz + real(r_def), parameter :: scaled_radius = rmdi + call get_w3_m3x3_q3x3x3_size( ndf_w3, undf_w3, unused, & unused, unused, unused, & unused, nlayers=nlayers) @@ -131,6 +107,8 @@ contains latitude, longitude, & chi_1, chi_2, chi_3, & panel_id, & + geometry, topology, & + coord_system, scaled_radius, & ndf_w3, undf_w3, map_w3(:,1), & ndf_chi, undf_chi, map_chi(:,1), & basis_chi(:,:,1,:), & diff --git a/components/science/unit-test/kernel/geometry/coordinate_jacobian_alphabetaz_mod_test.pf b/components/science/unit-test/kernel/geometry/coordinate_jacobian_alphabetaz_mod_test.pf index 756a4c84b..ae1aa557a 100644 --- a/components/science/unit-test/kernel/geometry/coordinate_jacobian_alphabetaz_mod_test.pf +++ b/components/science/unit-test/kernel/geometry/coordinate_jacobian_alphabetaz_mod_test.pf @@ -8,11 +8,11 @@ module coordinate_jacobian_alphabetaz_mod_test use funit - use constants_mod, only : r_def, i_def + use constants_mod, only: r_def, i_def - use base_mesh_config_mod, only : geometry_spherical, & - topology_fully_periodic - use finite_element_config_mod, only : coord_system_native + use base_mesh_config_mod, only: geometry_spherical, & + topology_fully_periodic + use finite_element_config_mod, only: coord_system_native implicit none @@ -38,12 +38,10 @@ contains @after subroutine tear_down() - use config_loader_mod, only: final_configuration use sci_chi_transform_mod, only: final_chi_transforms implicit none - call final_configuration() call final_chi_transforms() end subroutine tear_down diff --git a/components/science/unit-test/kernel/geometry/coordinate_jacobian_lonlatz_mod_test.pf b/components/science/unit-test/kernel/geometry/coordinate_jacobian_lonlatz_mod_test.pf index 3d717c8b0..e4a2c021f 100644 --- a/components/science/unit-test/kernel/geometry/coordinate_jacobian_lonlatz_mod_test.pf +++ b/components/science/unit-test/kernel/geometry/coordinate_jacobian_lonlatz_mod_test.pf @@ -10,6 +10,10 @@ module coordinate_jacobian_lonlatz_mod_test use funit use constants_mod, only : r_def, i_def, PI + use base_mesh_config_mod, only: geometry_spherical, & + topology_non_periodic + use finite_element_config_mod, only: coord_system_native + implicit none public :: set_up, tear_down, test_all @@ -30,12 +34,10 @@ contains @after subroutine tear_down() - use config_loader_mod, only: final_configuration use sci_chi_transform_mod, only: final_chi_transforms implicit none - call final_configuration() call final_chi_transforms() end subroutine tear_down @@ -52,10 +54,6 @@ contains pointwise_coordinate_jacobian, & pointwise_coordinate_jacobian_inverse - use base_mesh_config_mod, only: geometry_spherical, & - topology_non_periodic - use finite_element_config_mod, only: coord_system_native - implicit none real(kind=r_def), parameter :: tol = 1.0e-12_r_def ! r_def 64bit diff --git a/components/science/unit-test/kernel/geometry/native_jacobian_alphabetaz_mod_test.pf b/components/science/unit-test/kernel/geometry/native_jacobian_alphabetaz_mod_test.pf index 53aed1296..162bc1642 100644 --- a/components/science/unit-test/kernel/geometry/native_jacobian_alphabetaz_mod_test.pf +++ b/components/science/unit-test/kernel/geometry/native_jacobian_alphabetaz_mod_test.pf @@ -37,12 +37,10 @@ contains @after subroutine tear_down() - use config_loader_mod, only: final_configuration use sci_chi_transform_mod, only: final_chi_transforms implicit none - call final_configuration() call final_chi_transforms() end subroutine tear_down diff --git a/components/science/unit-test/kernel/geometry/native_jacobian_lonlatz_mod_test.pf b/components/science/unit-test/kernel/geometry/native_jacobian_lonlatz_mod_test.pf index e63daeeab..53ed27607 100644 --- a/components/science/unit-test/kernel/geometry/native_jacobian_lonlatz_mod_test.pf +++ b/components/science/unit-test/kernel/geometry/native_jacobian_lonlatz_mod_test.pf @@ -38,12 +38,10 @@ contains @after subroutine tear_down() - use config_loader_mod, only: final_configuration use sci_chi_transform_mod, only: final_chi_transforms implicit none - call final_configuration() call final_chi_transforms() end subroutine tear_down diff --git a/components/science/unit-test/kernel/geometry/nodal_xyz_coordinates_kernel_mod_test.pf b/components/science/unit-test/kernel/geometry/nodal_xyz_coordinates_kernel_mod_test.pf index 01bccfc8e..ed263deff 100644 --- a/components/science/unit-test/kernel/geometry/nodal_xyz_coordinates_kernel_mod_test.pf +++ b/components/science/unit-test/kernel/geometry/nodal_xyz_coordinates_kernel_mod_test.pf @@ -6,7 +6,7 @@ module nodal_xyz_coordinates_kernel_mod_test - use constants_mod, only : i_def, r_def, imdi + use constants_mod, only : i_def, r_def, imdi, rmdi use get_unit_test_m3x3_q3x3x3_sizes_mod, only : get_w0_m3x3_q3x3x3_size, & get_w2_m3x3_q3x3x3_size, & @@ -20,82 +20,60 @@ module nodal_xyz_coordinates_kernel_mod_test use get_unit_test_3x3x3_chi_mod, only : get_w0_3x3x3_field + use finite_element_config_mod, only: coord_system_xyz + use funit implicit none private - public :: test_all - - @TestCase - type, extends(TestCase), public :: nodal_xyz_coordinates_test_type - private - contains - procedure setUp - procedure tearDown - procedure test_all - end type nodal_xyz_coordinates_test_type + public :: set_up, tear_down, test_all - integer(i_def), parameter :: element_order_h = 0 - integer(i_def), parameter :: element_order_v = 0 contains !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - subroutine setUp( this ) + @before + subroutine set_up() - use sci_chi_transform_mod, only : init_chi_transforms - use feign_config_mod, only : feign_finite_element_config - use finite_element_config_mod, only : cellshape_quadrilateral, & - coord_system_xyz + use sci_chi_transform_mod, only: init_chi_transforms implicit none - class(nodal_xyz_coordinates_test_type), intent(inout) :: this - - call feign_finite_element_config( & - cellshape=cellshape_quadrilateral, & - coord_order=0_i_def, & - coord_system=coord_system_xyz, & - element_order_h=0_i_def, & - element_order_v=0_i_def, & - rehabilitate=.true. ) - call init_chi_transforms(imdi, imdi) - end subroutine setUp + end subroutine set_up !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - subroutine tearDown( this ) + @after + subroutine tear_down() use sci_chi_transform_mod, only: final_chi_transforms - use config_loader_mod, only: final_configuration implicit none - class(nodal_xyz_coordinates_test_type), intent(inout) :: this - - ! Finalise namelists - call final_configuration() call final_chi_transforms() - end subroutine tearDown + end subroutine tear_down !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! @Test - subroutine test_all( this ) + subroutine test_all() use sci_nodal_xyz_coordinates_kernel_mod, only: nodal_xyz_coordinates_code implicit none - class(nodal_xyz_coordinates_test_type), intent(inout) :: this - real(r_def), parameter :: dx = 6000.0_r_def real(r_def), parameter :: dy = 1000.0_r_def real(r_def), parameter :: dz = 2000.0_r_def real(r_def), parameter :: tol = 1.0e-6_r_def + integer(i_def), parameter :: geometry = imdi + integer(i_def), parameter :: topology = imdi + integer(i_def), parameter :: coord_system = coord_system_xyz + real(r_def), parameter :: scaled_radius = rmdi + ! Fields real(r_def), allocatable :: nodal1(:), nodal2(:), nodal3(:) real(r_def), allocatable :: chi1(:), chi2(:), chi3(:), panel_id(:) @@ -171,6 +149,10 @@ contains chi2, & chi3, & panel_id, & + geometry, & + topology, & + coord_system, & + scaled_radius, & ndf_w2, undf_w2, & map_w2(:,cell), & ndf_w0, undf_w0, & diff --git a/components/science/unit-test/kernel/inter_function_space/compute_map_u_operators_kernel_mod_test.pf b/components/science/unit-test/kernel/inter_function_space/compute_map_u_operators_kernel_mod_test.pf index 2a5c4dd49..5686d794b 100644 --- a/components/science/unit-test/kernel/inter_function_space/compute_map_u_operators_kernel_mod_test.pf +++ b/components/science/unit-test/kernel/inter_function_space/compute_map_u_operators_kernel_mod_test.pf @@ -7,60 +7,30 @@ module compute_map_u_operators_kernel_mod_test use constants_mod, only : i_def, pi, r_def + + use base_mesh_config_mod, only: geometry_spherical, & + topology_non_periodic + use finite_element_config_mod, only: coord_system_native + use funit + + implicit none private public :: set_up, tear_down, test_all - real(r_def), parameter :: planet_radius = 6000000.0_r_def - contains !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! @before subroutine set_up() - use base_mesh_config_mod, only : geometry_spherical, & - topology_non_periodic - use extrusion_config_mod, only : method_uniform, & - stretching_method_linear - use sci_chi_transform_mod, only : init_chi_transforms - use finite_element_config_mod, only : cellshape_quadrilateral, & - coord_system_native - use feign_config_mod, only : feign_extrusion_config, & - feign_finite_element_config, & - feign_base_mesh_config, & - feign_planet_config + use sci_chi_transform_mod, only: init_chi_transforms implicit none - call feign_base_mesh_config( file_prefix='foo', & - prime_mesh_name='unit_test', & - geometry=geometry_spherical, & - prepartitioned=.false., & - topology=topology_non_periodic, & - fplane=.false., f_lat_deg=0.0_r_def ) - - call feign_extrusion_config( method=method_uniform, & - planet_radius=planet_radius, & - domain_height=10.0_r_def, & - number_of_layers=4_i_def, & - stretching_height=0.5_r_def, & - stretching_method=stretching_method_linear, & - eta_values=(/0.5_r_def/) ) - - call feign_finite_element_config( & - cellshape=cellshape_quadrilateral, & - element_order_h=0_i_def, & - element_order_v=0_i_def, & - rehabilitate=.true., & - coord_order = 0_i_def, & - coord_system=coord_system_native ) - - call feign_planet_config( scaling_factor=1.0_r_def ) - call init_chi_transforms(geometry_spherical,topology_non_periodic) end subroutine set_up @@ -69,12 +39,10 @@ contains @after subroutine tear_down() - use config_loader_mod, only: final_configuration use sci_chi_transform_mod, only: final_chi_transforms implicit none - call final_configuration() call final_chi_transforms() end subroutine tear_down @@ -114,13 +82,21 @@ contains only : get_gaussian_q3x3x3_quadrature_weights_xy, & get_gaussian_q3x3x3_quadrature_weights_z - implicit none + implicit none + + real(r_def), parameter :: planet_radius = 6000000.0_r_def + real(r_def), parameter :: scaling = 1.0_r_def real(r_def), parameter :: tol = 1.0e-3_r_def real(r_def), parameter :: dlon = 1.0_r_def/planet_radius, & dlat = 1.0_r_def/planet_radius, & dz = 1.0_r_def + integer(i_def), parameter :: geometry = geometry_spherical + integer(i_def), parameter :: topology = topology_non_periodic + integer(i_def), parameter :: coord_system = coord_system_native + real(r_def), parameter :: scaled_radius = scaling*planet_radius + real(r_def) :: answer integer(i_def) :: cell, ncells, ncell_3d @@ -220,13 +196,17 @@ contains u_up_data(:) = 0.0_r_def panel_id_data(:) = 1.0_r_def - call compute_map_u_operators_code( cell, nlayers, ncell_3d, & + call compute_map_u_operators_code( cell, nlayers, ncell_3d, & u_lon_map, ncell_3d, u_lat_map, & ncell_3d, u_up_map, & chi_data(:,1), & chi_data(:,2), & chi_data(:,3), & panel_id_data, & + geometry, & + topology, & + coord_system, & + scaled_radius, & ndf_w2, basis_w2, & ndf_w3, basis_w3, & ndf_wt, basis_wt, & diff --git a/components/science/unit-test/kernel/inter_function_space/compute_sample_u_ops_kernel_mod_test.pf b/components/science/unit-test/kernel/inter_function_space/compute_sample_u_ops_kernel_mod_test.pf index f5fda906e..3778ec60b 100644 --- a/components/science/unit-test/kernel/inter_function_space/compute_sample_u_ops_kernel_mod_test.pf +++ b/components/science/unit-test/kernel/inter_function_space/compute_sample_u_ops_kernel_mod_test.pf @@ -8,6 +8,11 @@ module compute_sample_u_ops_kernel_mod_test use constants_mod, only : i_def, r_def use reference_element_mod, only : W, S, E, N, B, T + + use base_mesh_config_mod, only: geometry_spherical, & + topology_non_periodic + use finite_element_config_mod, only: coord_system_native + use funit implicit none @@ -15,54 +20,16 @@ module compute_sample_u_ops_kernel_mod_test private public :: set_up, tear_down, test_all - real(r_def), parameter :: radius = 6000000.0_r_def - real(r_def), parameter :: scaling = 1.0_r_def - contains !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! @before subroutine set_up() - use base_mesh_config_mod, only : geometry_spherical, & - topology_non_periodic - use extrusion_config_mod, only : method_uniform, & - stretching_method_linear - use finite_element_config_mod, only : cellshape_quadrilateral, & - coord_system_native - use feign_config_mod, only : feign_base_mesh_config, & - feign_extrusion_config, & - feign_finite_element_config, & - feign_planet_config - use sci_chi_transform_mod, only : init_chi_transforms + use sci_chi_transform_mod, only: init_chi_transforms implicit none - call feign_base_mesh_config( file_prefix='foo', & - prime_mesh_name='unit_test', & - geometry=geometry_spherical, & - prepartitioned=.false., & - topology=topology_non_periodic, & - fplane=.false., f_lat_deg=0.0_r_def ) - - call feign_extrusion_config( method=method_uniform, & - planet_radius=radius, & - domain_height=10.0_r_def, & - number_of_layers=5_i_def, & - stretching_method=stretching_method_linear, & - stretching_height=15.0_r_def, & - eta_values=(/0.5_r_def/) ) - - call feign_finite_element_config( & - cellshape=cellshape_quadrilateral, & - element_order_h=0_i_def, & - element_order_v=0_i_def, & - rehabilitate=.true., & - coord_order = 0_i_def, & - coord_system=coord_system_native ) - - call feign_planet_config( scaling_factor=scaling ) - call init_chi_transforms(geometry_spherical, topology_non_periodic) end subroutine set_up @@ -71,12 +38,10 @@ contains @after subroutine tear_down() - use config_loader_mod, only: final_configuration use sci_chi_transform_mod, only: final_chi_transforms implicit none - call final_configuration() call final_chi_transforms() end subroutine tear_down @@ -102,6 +67,7 @@ contains implicit none + real(r_def), parameter :: radius = 6000000.0_r_def real(r_def), parameter :: tol = 1.0e-2_r_def real(r_def), parameter :: dlon = 2.0_r_def/radius real(r_def), parameter :: dlat = 3.0_r_def/radius @@ -117,6 +83,14 @@ contains real(r_def), parameter :: du_lat0_dz = -0.5_r_def real(r_def), parameter :: du_up0_dz = 0.001_r_def + + real(r_def), parameter :: scaling = 1.0_r_def + + integer(i_def),parameter :: geometry = geometry_spherical + integer(i_def),parameter :: topology = topology_non_periodic + integer(i_def),parameter :: coord_system = coord_system_native + real(r_def), parameter :: scaled_radius = scaling*radius + integer(i_def) :: cell, ncols_x, ncols_y, ncols_2d, ncells_3d integer(i_def) :: nlayers, k, nfaces @@ -275,6 +249,10 @@ contains chi_data(:,2), & chi_data(:,3), & panel_id_data, & + geometry, & + topology, & + coord_system, & + scaled_radius, & ndf_w2b, ndf_w3, ndf_wt, & ndf_wchi, undf_wchi, map_wchi(:,cell), & basis_wchi, diff_basis_wchi, & diff --git a/components/science/unit-test/kernel/inter_function_space/convert_phys_to_hdiv_kernel_mod_test.pf b/components/science/unit-test/kernel/inter_function_space/convert_phys_to_hdiv_kernel_mod_test.pf index 247e482e9..571181d93 100644 --- a/components/science/unit-test/kernel/inter_function_space/convert_phys_to_hdiv_kernel_mod_test.pf +++ b/components/science/unit-test/kernel/inter_function_space/convert_phys_to_hdiv_kernel_mod_test.pf @@ -7,77 +7,48 @@ !! W2 wind field module convert_phys_to_hdiv_kernel_mod_test - use constants_mod, only : i_def, r_def + use constants_mod, only: i_def, r_def, rmdi + + use base_mesh_config_mod, only: geometry_planar, & + topology_fully_periodic + use finite_element_config_mod, only: coord_system_xyz + use funit implicit none private - public :: convert_phys_to_hdiv_test_type, test_all - - @TestCase - type, extends(TestCase) :: convert_phys_to_hdiv_test_type - private - contains - procedure setUp - procedure tearDown - procedure test_all - end type convert_phys_to_hdiv_test_type + public :: set_up, tear_down, test_all contains !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - subroutine setUp( this ) + @before + subroutine set_up() - use base_mesh_config_mod, only : geometry_planar, & - topology_fully_periodic - use sci_chi_transform_mod, only : init_chi_transforms - use finite_element_config_mod, only : cellshape_quadrilateral, & - coord_system_xyz - use feign_config_mod, only : feign_finite_element_config, & - feign_base_mesh_config + use sci_chi_transform_mod, only: init_chi_transforms implicit none - class(convert_phys_to_hdiv_test_type), intent(inout) :: this - - call feign_base_mesh_config( file_prefix='foo', & - prime_mesh_name='unit_test', & - geometry=geometry_planar, & - prepartitioned=.false., & - topology=topology_fully_periodic, & - fplane=.false., f_lat_deg=0.0_r_def ) - - call feign_finite_element_config( & - cellshape=cellshape_quadrilateral, & - coord_order=0_i_def, & - coord_system=coord_system_xyz, & - element_order_h=0_i_def, & - element_order_v=0_i_def, & - rehabilitate=.true. ) - call init_chi_transforms(geometry_planar, topology_fully_periodic) - end subroutine setUp + end subroutine set_up !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - subroutine tearDown( this ) + @after + subroutine tear_down() - use config_loader_mod, only: final_configuration use sci_chi_transform_mod, only: final_chi_transforms implicit none - class(convert_phys_to_hdiv_test_type), intent(inout) :: this - - call final_configuration() call final_chi_transforms() - end subroutine tearDown + end subroutine tear_down !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! @Test - subroutine test_all( this ) + subroutine test_all() use sci_convert_phys_to_hdiv_kernel_mod, only : convert_phys_to_hdiv_code @@ -94,11 +65,8 @@ contains get_w2_w2nodal_basis use get_unit_test_3x3x3_chi_mod, only: get_w0_3x3x3_field - use base_mesh_config_mod, only: geometry_planar - - implicit none - class(convert_phys_to_hdiv_test_type), intent(inout) :: this + implicit none real(r_def), parameter :: tol = 1.0e-3_r_def real(r_def), parameter :: dx = 6000.0_r_def @@ -108,6 +76,11 @@ contains real(r_def), parameter :: u_meridional = 3.0_r_def real(r_def), parameter :: u_radial = 0.4_r_def + integer(i_def), parameter :: geometry = geometry_planar + integer(i_def), parameter :: topology = topology_fully_periodic + integer(i_def), parameter :: coord_system = coord_system_xyz + real(r_def), parameter :: scaled_radius = rmdi + real(r_def) :: answer integer(i_def) :: cell @@ -192,7 +165,10 @@ contains chi_data(:,2), & chi_data(:,3), & panel_id_data, & - geometry_planar, & + geometry, & + topology, & + coord_system, & + scaled_radius, & ndf_w2, undf_w2, & map_w2(:,cell), basis_w2, & ndf_w0, undf_w0, & diff --git a/components/science/unit-test/kernel/inter_function_space/project_ws_to_w1_operator_kernel_mod_test.pf b/components/science/unit-test/kernel/inter_function_space/project_ws_to_w1_operator_kernel_mod_test.pf index 3c8c17e91..5f1d0f284 100644 --- a/components/science/unit-test/kernel/inter_function_space/project_ws_to_w1_operator_kernel_mod_test.pf +++ b/components/science/unit-test/kernel/inter_function_space/project_ws_to_w1_operator_kernel_mod_test.pf @@ -7,7 +7,7 @@ !> Test the projection from a scalar space to W1 module project_ws_to_w1_operator_kernel_mod_test - use constants_mod, only : i_def, r_def + use constants_mod, only : i_def, r_def, rmdi use funit use get_unit_test_q3x3x3_quadrature_mod, only : get_gaussian_q3x3x3_quadrature_weights_xy, & get_gaussian_q3x3x3_quadrature_weights_z @@ -23,83 +23,56 @@ module project_ws_to_w1_operator_kernel_mod_test get_w3_m3x3_dofmap use get_unit_test_3x3x3_chi_mod, only : get_w0_3x3x3_field + use base_mesh_config_mod, only: geometry_planar, & + topology_fully_periodic + use finite_element_config_mod, only: coord_system_xyz + implicit none private - public :: test_all - - @TestCase - type, extends(TestCase), public :: project_ws_to_w1_operator_test_type - private - contains - procedure setUp - procedure tearDown - procedure test_all - end type project_ws_to_w1_operator_test_type + public :: set_up, tear_down, test_all contains !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - subroutine setUp( this ) + @before + subroutine set_up() - use finite_element_config_mod, only : cellshape_quadrilateral, & - coord_system_xyz - use feign_config_mod, only : feign_finite_element_config - use base_mesh_config_mod, only : geometry_planar, & - topology_fully_periodic - use feign_config_mod, only : feign_base_mesh_config - use sci_chi_transform_mod, only : init_chi_transforms + use sci_chi_transform_mod, only: init_chi_transforms implicit none - class(project_ws_to_w1_operator_test_type), intent(inout) :: this - - call feign_base_mesh_config( file_prefix='foo', & - prime_mesh_name='unit_test', & - geometry=geometry_planar, & - prepartitioned=.false., & - topology=topology_fully_periodic, & - fplane=.false., f_lat_deg=0.0_r_def ) - - call feign_finite_element_config( cellshape=cellshape_quadrilateral, & - coord_order=0_i_def, & - element_order_h=0_i_def, & - element_order_v=0_i_def, & - rehabilitate=.true., & - coord_system=coord_system_xyz ) - call init_chi_transforms(geometry_planar,topology_fully_periodic) - end subroutine setUp + end subroutine set_up !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - subroutine tearDown( this ) + @after + subroutine tear_down() - use config_loader_mod, only: final_configuration use sci_chi_transform_mod, only: final_chi_transforms implicit none - class(project_ws_to_w1_operator_test_type), intent(inout) :: this - - call final_configuration() call final_chi_transforms() - end subroutine tearDown + end subroutine tear_down !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - @Test - subroutine test_all( this ) + subroutine test_all() use sci_project_ws_to_w1_operator_kernel_mod, only : project_ws_to_w1_operator_code implicit none - class(project_ws_to_w1_operator_test_type), intent(inout) :: this - real(kind=r_def), parameter :: tol = 1.0e-6_r_def + integer(i_def), parameter :: geometry = geometry_planar + integer(i_def), parameter :: topology = topology_fully_periodic + integer(i_def), parameter :: coord_system = coord_system_xyz + real(r_def), parameter :: scaled_radius = rmdi + ! Mesh integer(i_def) :: nlayers, ncells, ncell_3d, cell, icell integer(i_def) :: ndf_w0, undf_w0, ndf_w1, undf_w1 @@ -210,6 +183,10 @@ contains chi1, chi2, chi3, & panel_id, & direction, & + geometry, & + topology, & + coord_system, & + scaled_radius, & ndf_w1, & basis_w1, & ndf_w3, & diff --git a/components/science/unit-test/kernel/inter_function_space/w3_to_w2_displacement_kernel_mod_test.pf b/components/science/unit-test/kernel/inter_function_space/w3_to_w2_displacement_kernel_mod_test.pf index c24a3b329..6d01ea321 100644 --- a/components/science/unit-test/kernel/inter_function_space/w3_to_w2_displacement_kernel_mod_test.pf +++ b/components/science/unit-test/kernel/inter_function_space/w3_to_w2_displacement_kernel_mod_test.pf @@ -7,9 +7,15 @@ !> Test the kernel to compute errors in W3 to W2 averaging module w3_to_w2_displacement_kernel_mod_test - use constants_mod, only: i_def, r_def, PI, l_def - use reference_element_mod, only: S, E, N, W - use funit + use constants_mod, only: i_def, r_def, PI, l_def + use reference_element_mod, only: S, E, N, W + + use base_mesh_config_mod, only: geometry_spherical, & + topology_fully_periodic + use finite_element_config_mod, only: coord_system_native + + use funit + implicit none @@ -22,45 +28,10 @@ contains @before subroutine set_up() - use base_mesh_config_mod, only : geometry_spherical, & - topology_fully_periodic - use extrusion_config_mod, only : method_uniform, & - stretching_method_linear - use finite_element_config_mod, only : cellshape_quadrilateral, & - coord_system_native - use feign_config_mod, only : feign_base_mesh_config, & - feign_extrusion_config, & - feign_finite_element_config, & - feign_planet_config - use sci_chi_transform_mod, only : init_chi_transforms + use sci_chi_transform_mod, only: init_chi_transforms implicit none - call feign_base_mesh_config( file_prefix='foo', & - prime_mesh_name='unit_test', & - geometry=geometry_spherical, & - prepartitioned=.false., & - topology=topology_fully_periodic, & - fplane=.false., f_lat_deg=0.0_r_def ) - - call feign_extrusion_config( method=method_uniform, & - planet_radius=1900000.0_r_def, & - domain_height=10.0_r_def, & - number_of_layers=5_i_def, & - stretching_method=stretching_method_linear, & - stretching_height=15.0_r_def, & - eta_values=(/0.5_r_def/) ) - - call feign_finite_element_config( & - cellshape=cellshape_quadrilateral, & - element_order_h=0_i_def, & - element_order_v=0_i_def, & - rehabilitate=.true., & - coord_order = 0_i_def, & - coord_system=coord_system_native ) - - call feign_planet_config( scaling_factor=1.0_r_def ) - call init_chi_transforms(geometry_spherical, topology_fully_periodic) end subroutine set_up @@ -69,12 +40,10 @@ contains @after subroutine tear_down() - use config_loader_mod, only: final_configuration use sci_chi_transform_mod, only: final_chi_transforms implicit none - call final_configuration() call final_chi_transforms() end subroutine tear_down @@ -84,8 +53,8 @@ contains subroutine test_all() use sci_w3_to_w2_displacement_kernel_mod, only: w3_to_w2_displacement_code - use get_unit_test_w2hnodal_basis_mod, only: get_wchi_w2hnodal_basis - use get_unit_test_w3nodal_basis_mod, only: get_wchi_w3nodal_basis + use get_unit_test_w2hnodal_basis_mod, only: get_wchi_w2hnodal_basis + use get_unit_test_w3nodal_basis_mod, only: get_wchi_w3nodal_basis implicit none @@ -96,6 +65,11 @@ contains real(r_def), parameter :: beta0 = PI/6.0_r_def real(r_def), parameter :: dz = 2.0_r_def + integer(i_def), parameter :: geometry = geometry_spherical + integer(i_def), parameter :: topology = topology_fully_periodic + integer(i_def), parameter :: coord_system = coord_system_native + real(r_def), parameter :: scaled_radius = 1900000.0_r_def + ! Non-periodic 3x3x3 domain integer(i_def), parameter :: ncells = 2 integer(i_def), parameter :: nlayers = 1 @@ -120,7 +94,7 @@ contains real(r_def) :: dummy_w3(undf_w3) real(r_def) :: answer, phi - integer(i_def) :: df_w2h, cell + integer(i_def) :: cell ! ------------------------------------------------------------------------ ! ! Make DoF maps @@ -177,6 +151,8 @@ contains chi_1, chi_2, chi_3, & panel_id, & dummy_w3, & + geometry, topology, & + coord_system, scaled_radius, & ndf_w2h, undf_w2h, map_w2h(:,cell), & ndf_chi, undf_chi, map_chi(:,cell), & basis_chi_w2h, basis_chi_w3, & diff --git a/rose-meta/lfric-lbc_demo b/rose-meta/lfric-lbc_demo new file mode 120000 index 000000000..c4b8f3c23 --- /dev/null +++ b/rose-meta/lfric-lbc_demo @@ -0,0 +1 @@ +../applications/lbc_demo/rose-meta/lfric-lbc_demo/ \ No newline at end of file diff --git a/rose-stem/app/coupled/rose-app.conf b/rose-stem/app/coupled/rose-app.conf index 38b1839b6..e9531bc67 100644 --- a/rose-stem/app/coupled/rose-app.conf +++ b/rose-stem/app/coupled/rose-app.conf @@ -70,10 +70,6 @@ write_diag=.false. log_to_rank_zero_only=.false. run_log_level='info' -[namelist:multigrid] -chain_mesh_tags='' -multigrid_chain_nitems=1 - [namelist:partitioning] generate_inner_halos=.true. panel_decomposition='auto' diff --git a/rose-stem/app/io_demo/rose-app.conf b/rose-stem/app/io_demo/rose-app.conf index d594e1c86..3199a7631 100644 --- a/rose-stem/app/io_demo/rose-app.conf +++ b/rose-stem/app/io_demo/rose-app.conf @@ -1,4 +1,4 @@ -meta=lfric-io_demo/HEAD +meta=lfric-io_demo/vn3.1 [command] default=$CORE_ROOT_DIR/bin/tweak_iodef ; \ @@ -84,10 +84,6 @@ temporal_reading=.false. log_to_rank_zero_only=.false. run_log_level='info' -[namelist:multigrid] -chain_mesh_tags='' -multigrid_chain_nitems=1 - [namelist:partitioning] generate_inner_halos=.false. panel_decomposition='auto' diff --git a/rose-stem/app/lbc_demo/opt/rose-app-mesh_lbc_demo.conf b/rose-stem/app/lbc_demo/opt/rose-app-mesh_lbc_demo.conf index 5225b2288..009f0da92 100644 --- a/rose-stem/app/lbc_demo/opt/rose-app-mesh_lbc_demo.conf +++ b/rose-stem/app/lbc_demo/opt/rose-app-mesh_lbc_demo.conf @@ -5,8 +5,8 @@ source=$MESH_DIR/mesh_lbc_demo.nc [namelist:base_mesh] file_prefix='mesh_lbc_demo' geometry='spherical' -prepartitioned=.true., -prime_mesh_name='primary', -topology='non_periodic', +prepartitioned=.true. +prime_mesh_name='primary' +topology='non_periodic' [!!namelist:partitioning] diff --git a/rose-stem/app/lbc_demo/rose-app.conf b/rose-stem/app/lbc_demo/rose-app.conf index d979a47c6..b81c4d29f 100644 --- a/rose-stem/app/lbc_demo/rose-app.conf +++ b/rose-stem/app/lbc_demo/rose-app.conf @@ -1,4 +1,4 @@ -meta=lfric-lbc_demo/HEAD +meta=lfric-lbc_demo/vn3.1 [command] default=$CORE_ROOT_DIR/bin/tweak_iodef ; \ diff --git a/rose-stem/app/simple_diffusion/rose-app.conf b/rose-stem/app/simple_diffusion/rose-app.conf index f710dfe34..09600a626 100644 --- a/rose-stem/app/simple_diffusion/rose-app.conf +++ b/rose-stem/app/simple_diffusion/rose-app.conf @@ -71,10 +71,6 @@ write_diag=.false. log_to_rank_zero_only=.false. run_log_level='info' -[namelist:multigrid] -chain_mesh_tags='' -multigrid_chain_nitems=1 - [namelist:partitioning] generate_inner_halos=.true. panel_decomposition='auto' diff --git a/rose-stem/app/skeleton/rose-app.conf b/rose-stem/app/skeleton/rose-app.conf index 7e966d0b7..11307c9d3 100644 --- a/rose-stem/app/skeleton/rose-app.conf +++ b/rose-stem/app/skeleton/rose-app.conf @@ -68,10 +68,6 @@ write_diag=.false. log_to_rank_zero_only=.false. run_log_level='info' -[namelist:multigrid] -chain_mesh_tags='' -multigrid_chain_nitems=1 - [namelist:partitioning] generate_inner_halos=.true. panel_decomposition='auto'