diff --git a/applications/io_demo/source/driver/io_demo_checkpoint_mod.f90 b/applications/io_demo/source/driver/io_demo_checkpoint_mod.f90 index 0fa29fcb9..bc6fabcee 100644 --- a/applications/io_demo/source/driver/io_demo_checkpoint_mod.f90 +++ b/applications/io_demo/source/driver/io_demo_checkpoint_mod.f90 @@ -25,11 +25,6 @@ module io_demo_checkpoint_mod LOG_LEVEL_DEBUG, LOG_LEVEL_ERROR use mesh_mod, only: mesh_type - use base_mesh_config_mod, only: geometry_spherical, & - geometry_planar, & - topology_fully_periodic, & - topology_non_periodic - implicit none private @@ -71,21 +66,9 @@ subroutine setup_checkpoint_io(modeldb, chi, panel_id) call log_event( 'io_demo: Setting up checkpoint I/O', LOG_LEVEL_DEBUG ) mesh => chi(1)%get_mesh() - if (mesh%is_geometry_spherical()) then - geometry = geometry_spherical - else - geometry = geometry_planar - end if - - if (mesh%is_topology_periodic()) then - topology = topology_fully_periodic - else if (mesh%is_topology_non_periodic()) then - topology = topology_non_periodic - else - call log_event( 'Unsupported mesh topology', & - log_level_error ) - end if + geometry = mesh%geometry() + topology = mesh%topology() coord_system = modeldb%config%finite_element%coord_system() scaled_radius = modeldb%config%planet%scaled_radius() diff --git a/applications/lbc_demo/source/kernel/set_lbc_int_kernel_mod.F90 b/applications/lbc_demo/source/kernel/set_lbc_int_kernel_mod.F90 index 0a7c37ca7..4cb31acf2 100644 --- a/applications/lbc_demo/source/kernel/set_lbc_int_kernel_mod.F90 +++ b/applications/lbc_demo/source/kernel/set_lbc_int_kernel_mod.F90 @@ -14,7 +14,7 @@ module set_lbc_int_kernel_mod use constants_mod, only: r_def, i_def, l_def, radians_to_degrees use kernel_mod, only: kernel_type - use base_mesh_config_mod, only: geometry_spherical + use mesh_mod, only: geometry_spherical implicit none diff --git a/applications/lbc_demo/source/kernel/set_lbc_real_kernel_mod.F90 b/applications/lbc_demo/source/kernel/set_lbc_real_kernel_mod.F90 index 02d8d983e..3c0765829 100644 --- a/applications/lbc_demo/source/kernel/set_lbc_real_kernel_mod.F90 +++ b/applications/lbc_demo/source/kernel/set_lbc_real_kernel_mod.F90 @@ -14,7 +14,7 @@ module set_lbc_real_kernel_mod use constants_mod, only: r_def, i_def, l_def, radians_to_degrees use kernel_mod, only: kernel_type - use base_mesh_config_mod, only: geometry_spherical + use mesh_mod, only: geometry_spherical implicit none diff --git a/applications/skeleton/source/algorithm/skeleton_constants_mod.x90 b/applications/skeleton/source/algorithm/skeleton_constants_mod.x90 index 4717f9b7c..9157664d9 100644 --- a/applications/skeleton/source/algorithm/skeleton_constants_mod.x90 +++ b/applications/skeleton/source/algorithm/skeleton_constants_mod.x90 @@ -15,7 +15,7 @@ module skeleton_constants_mod ! Infrastructure - use constants_mod, only: str_def, i_def, l_def + use constants_mod, only: str_def, r_def, i_def, l_def use driver_modeldb_mod, only: modeldb_type use field_mod, only: field_type use fs_continuity_mod, only: W0, W1, W2, W2broken, & @@ -80,7 +80,7 @@ contains ! Arguments type(modeldb_type), intent(in) :: modeldb - type(mesh_type), pointer, intent(in) :: mesh + type(mesh_type), intent(in) :: mesh type(field_type), target, intent(in) :: chi(:) type(field_type), target, intent(in) :: panel_id @@ -109,6 +109,10 @@ contains integer(i_def) :: order_h, order_v integer(i_def) :: nqp_h_exact, nqp_v_exact + integer(i_def) :: geometry + integer(i_def) :: topology + integer(i_def) :: coord_system + real(r_def) :: scaled_radius integer(tik) :: id @@ -116,10 +120,15 @@ contains call log_event( "Skeleton: creating runtime constants", LOG_LEVEL_INFO ) - order_h = modeldb%config%finite_element%element_order_h() - order_v = modeldb%config%finite_element%element_order_v() - nqp_h_exact = modeldb%config%finite_element%nqp_h_exact() - nqp_v_exact = modeldb%config%finite_element%nqp_v_exact() + geometry = mesh%geometry() + topology = mesh%topology() + + order_h = modeldb%config%finite_element%element_order_h() + order_v = modeldb%config%finite_element%element_order_v() + nqp_h_exact = modeldb%config%finite_element%nqp_h_exact() + nqp_v_exact = modeldb%config%finite_element%nqp_v_exact() + coord_system = modeldb%config%finite_element%coord_system() + scaled_radius = modeldb%config%planet%scaled_radius() !=========== Create function spaces required for setup ==================! @@ -167,13 +176,16 @@ contains call curl_inventory%add_operator( curl, w2_fs, w1_fs, mesh ) call broken_div_inventory%add_operator( broken_div, w3_fs, w2b_fs, mesh ) - call invoke( name = "create_de_rham_matrices", & - compute_derham_matrices_kernel_type(mm_w0, mm_w1, & - mm_w2, mm_w2b, & - mm_w3, mm_wtheta, & - grad, curl, & - div, broken_div, & - chi, panel_id, qr) ) + call invoke( name = "create_de_rham_matrices", & + compute_derham_matrices_kernel_type(mm_w0, mm_w1, & + mm_w2, mm_w2b, & + mm_w3, mm_wtheta, & + grad, curl, & + div, broken_div, & + chi, panel_id, & + geometry, topology, & + coord_system, & + scaled_radius, qr) ) !===================== Add operators to inventories =======================! diff --git a/components/driver/source/driver_coordinates_mod.F90 b/components/driver/source/driver_coordinates_mod.F90 index d5c5d98a4..b6595073b 100644 --- a/components/driver/source/driver_coordinates_mod.F90 +++ b/components/driver/source/driver_coordinates_mod.F90 @@ -18,11 +18,12 @@ module driver_coordinates_mod schmidt_transform_xyz, & inverse_schmidt_transform_xyz + use mesh_mod, only: geometry_planar, & + geometry_spherical, & + topology_periodic, & + topology_non_periodic + ! 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 @@ -64,7 +65,7 @@ subroutine assign_coordinate_field(config, mesh, chi, panel_id) implicit none type(config_type), intent(in) :: config - type(mesh_type), intent(in), pointer :: mesh + type(mesh_type), intent(in) :: mesh type(field_type), intent(inout) :: chi(3) type(field_type), intent(inout) :: panel_id @@ -108,18 +109,8 @@ subroutine assign_coordinate_field(config, mesh, chi, panel_id) nullify( map, map_pid, dof_coords, reference_element ) - if (mesh%is_geometry_spherical()) then - geometry = geometry_spherical - else - geometry = geometry_planar - end if - - if (mesh%is_topology_periodic()) then - topology = topology_fully_periodic - else - topology = topology_non_periodic - end if - + geometry = mesh%geometry() + topology = mesh%topology() coord_system = config%finite_element%coord_system() scaled_radius = config%planet%scaled_radius() @@ -178,7 +169,7 @@ subroutine assign_coordinate_field(config, mesh, chi, panel_id) ! Throw an error if stretching factor is not 1 and not on cubed-sphere if ( abs(stretch_factor - 1.0_r_def) > eps .and. .not. & (geometry == geometry_spherical .and. & - topology == topology_fully_periodic) ) then + topology == topology_periodic) ) then call log_event( & 'driver_coordinates: Cannot determine coordinates if Schmidt ' // & 'stretching factor is not 1 and mesh is not cubed-sphere', & @@ -228,7 +219,7 @@ subroutine assign_coordinate_field(config, mesh, chi, panel_id) end do else if ( geometry == geometry_spherical .and. & - topology /= topology_fully_periodic ) then + topology /= topology_periodic ) then do cell = 1,chi_proxy(1)%vspace%get_ncell() @@ -263,7 +254,7 @@ subroutine assign_coordinate_field(config, mesh, chi, panel_id) end do else if ( geometry == geometry_spherical .and. & - topology == topology_fully_periodic ) then + topology == topology_periodic ) then do cell = 1,chi_proxy(1)%vspace%get_ncell() @@ -352,7 +343,7 @@ subroutine calc_panel_id( nlayers, & integer(kind=i_def) :: vert, k if ( geometry == geometry_spherical .and. & - topology == topology_fully_periodic ) then + topology == topology_periodic ) then ! The following code assumes that the mesh generator has ordered the ! global cell ids panel-by-panel. If this is ever not the case, the @@ -451,7 +442,7 @@ subroutine assign_coordinate_xyz( nlayers, & end if ! Domain does not have N-S boundaries only if topology completely periodic if ( column_coords(2,SWB,k+1) > column_coords(2,NWB,k+1) .and. & - topology == topology_fully_periodic ) then + topology == topology_periodic ) then ! On y boundary vertex_local_coords(2,SWB) = domain_y vertex_local_coords(2,SEB) = domain_y diff --git a/components/driver/source/driver_fem_mod.f90 b/components/driver/source/driver_fem_mod.f90 index 8d2bfe5e4..88f9803ea 100644 --- a/components/driver/source/driver_fem_mod.f90 +++ b/components/driver/source/driver_fem_mod.f90 @@ -26,6 +26,10 @@ module driver_fem_mod use sci_chi_transform_mod, only: init_chi_transforms, & final_chi_transforms + use mesh_mod, only: geometry_spherical, & + geometry_planar, & + topology_non_periodic + ! Object types use config_mod, only: config_type use field_mod, only: field_type @@ -33,10 +37,6 @@ module driver_fem_mod use inventory_by_mesh_mod, only: inventory_by_mesh_type ! Configuration modules - use base_mesh_config_mod, only: geometry_spherical, & - geometry_planar, & - topology_non_periodic, & - topology_fully_periodic use finite_element_config_mod, only: coord_system_xyz, & coord_space_W0, & coord_space_Wchi, & @@ -121,17 +121,8 @@ subroutine init_fem(config, chi_inventory, panel_id_inventory) mesh => mesh_collection%get_mesh(all_mesh_names(i)) mesh_name = mesh%get_mesh_name() - if (mesh%is_geometry_spherical()) then - geometry = geometry_spherical - else - geometry = geometry_planar - end if - - if (mesh%is_topology_periodic()) then - topology = topology_fully_periodic - else - topology = topology_non_periodic - end if + geometry = mesh%geometry() + topology = mesh%topology() ! Initialise coordinate transformations call init_chi_transforms( geometry, topology, & diff --git a/components/driver/source/driver_mesh_mod.f90 b/components/driver/source/driver_mesh_mod.f90 index 40157459a..13c0823eb 100644 --- a/components/driver/source/driver_mesh_mod.f90 +++ b/components/driver/source/driver_mesh_mod.f90 @@ -51,11 +51,11 @@ module driver_mesh_mod use global_mesh_collection_mod, only: global_mesh_collection use local_mesh_collection_mod, only: local_mesh_collection + use mesh_mod, only: geometry_spherical, & + topology_periodic ! Configuration modules use finite_element_config_mod, only: cellshape_quadrilateral - use base_mesh_config_mod, only: geometry_spherical, & - topology_fully_periodic implicit none @@ -287,7 +287,7 @@ subroutine init_mesh( config, & generate_inner_halos = config%partitioning%generate_inner_halos() if ( geometry == geometry_spherical .and. & - topology == topology_fully_periodic ) then + topology == topology_periodic ) then mesh_selection = mesh_cubedsphere call log_event( "Setting up cubed-sphere partition mesh(es)", & log_level_debug ) 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 2fc232c3a..e7445c4ab 100644 --- a/components/driver/unit-test/assign_coordinate_alphabetaz_mod_test.pf +++ b/components/driver/unit-test/assign_coordinate_alphabetaz_mod_test.pf @@ -31,7 +31,7 @@ contains 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 + integer(kind=i_def) :: map_chi(1), map_pid(1), nverts real(kind=r_def) :: alpha(2), beta(2), panel_rho(4) real(kind=r_def) :: alpha_out(1), beta_out(1), height_out(1), panel_id(1) real(kind=r_def) :: verts_XYZ(3,8,1), verts_ref(8,3), nodal_coord(3,1) 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 37eaac9ef..4434c9128 100644 --- a/components/driver/unit-test/assign_coordinate_xyz_mod_test.pf +++ b/components/driver/unit-test/assign_coordinate_xyz_mod_test.pf @@ -9,7 +9,7 @@ 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 mesh_mod, only: geometry_planar, topology_periodic use funit implicit none @@ -33,7 +33,7 @@ contains one = 1.0_r_def integer(i_def), parameter :: geometry = geometry_planar - integer(i_def), parameter :: topology = topology_fully_periodic + integer(i_def), parameter :: topology = topology_periodic real(r_def), parameter :: scaled_radius = 1.0_r_def integer(kind=i_def) :: nlayers, ndf, nverts, i, undf, ndf_pid, undf_pid diff --git a/components/science/Makefile b/components/science/Makefile index ac6023dd5..52826618a 100644 --- a/components/science/Makefile +++ b/components/science/Makefile @@ -71,9 +71,7 @@ document-api: api-documentation # unit-tests/%: export BIN_DIR ?= $(PROJECT_DIR)/test unit-tests/%: export CXX_LINK = TRUE -#unit-tests/%: export EXTERNAL_STATIC_LIBRARIES += pfunit unit-tests/%: export IMPORT_PARTS = $(CORE_ROOT_DIR)/infrastructure \ -# $(CORE_ROOT_DIR)/components/lfric-xios \ $(CORE_ROOT_DIR)/components/inventory \ $(CORE_ROOT_DIR)/components/driver unit-tests/%: export META_FILE_DIR = rose-meta/lfric-science/HEAD diff --git a/components/science/source/algorithm/sci_fem_constants_mod.x90 b/components/science/source/algorithm/sci_fem_constants_mod.x90 index 28c31d995..edc0fc1ec 100644 --- a/components/science/source/algorithm/sci_fem_constants_mod.x90 +++ b/components/science/source/algorithm/sci_fem_constants_mod.x90 @@ -15,6 +15,7 @@ module sci_fem_constants_mod ! Infrastructure + use config_mod, only: config_type use constants_mod, only: i_def, r_def, r_second, & str_def, l_def use field_mod, only: field_type @@ -25,21 +26,15 @@ module sci_fem_constants_mod use inventory_by_mesh_mod, only: inventory_by_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 mesh_mod, only: mesh_type, topology_periodic use operator_mod, only: operator_type use quadrature_xyoz_mod, only: quadrature_xyoz_type use quadrature_rule_gaussian_mod, only: quadrature_rule_gaussian_type use timing_mod, only: start_timing, stop_timing, & tik, LPROF - ! Configuration - use finite_element_config_mod, only: element_order_h, & - element_order_v, & - nqp_h_exact, & - nqp_v_exact - ! Other algorithms - use sci_geometric_constants_mod, only: get_coordinates, & + use sci_geometric_constants_mod, only: get_coordinates, & get_panel_id ! Kernels use sci_compute_div_operator_kernel_mod, & @@ -48,6 +43,12 @@ module sci_fem_constants_mod use dg_inc_matrix_vector_kernel_mod, only: dg_inc_matrix_vector_kernel_type use sci_multiplicity_kernel_mod, only: multiplicity_kernel_type + ! Configuration modules + use finite_element_config_mod, only: element_order_h, & + element_order_v, & + nqp_h_exact, & + nqp_v_exact + implicit none private @@ -197,10 +198,11 @@ contains end function get_qr_fv !> @brief Returns a pointer to a finite element mass matrix + !> @param[in] config Application configuration object + !> @param[in] mesh Mesh to get the object for !> @param[in] space The space of the desired mass matrix - !> @param[in] mesh_id The ID of the mesh to get the object for !> @return The mass matrix operator - function get_mass_matrix_fe(space, mesh_id) result(mm_op) + function get_mass_matrix_fe(config, mesh, space) result(mm_op) use sci_compute_mass_matrix_kernel_w_scalar_mod, & only: compute_mass_matrix_kernel_w_scalar_type @@ -210,13 +212,13 @@ contains only: compute_mass_matrix_kernel_w1_type use sci_edge_lump_w2_mass_matrix_kernel_mod, & only: edge_lump_w2_mass_matrix_kernel_type - use base_mesh_config_mod, only: topology, topology_fully_periodic implicit none - integer(kind=i_def), intent(in) :: space - integer(kind=i_def), intent(in) :: mesh_id - type(mesh_type), pointer :: mesh + type(config_type), intent(in) :: config + type(mesh_type), intent(in) :: mesh + integer(kind=i_def), intent(in) :: space + type(operator_type), pointer :: mm_op logical(kind=l_def) :: constant_exists type(field_type), pointer :: chi(:) @@ -230,9 +232,20 @@ contains integer(kind=i_def), parameter :: stencil_depth = 1_i_def integer(tik) :: id + integer(i_def) :: geometry + integer(i_def) :: topology + integer(i_def) :: coord_system + real(r_def) :: scaled_radius + + geometry = mesh%geometry() + topology = mesh%topology() + + coord_system = config%finite_element%coord_system() + scaled_radius = config%planet%scaled_radius() + ! If running at lowest order, use finite volume if (element_order_h == 0 .and. element_order_v == 0) then - mm_op => get_mass_matrix_fv(space, mesh_id) + mm_op => get_mass_matrix_fv(config, mesh, space) return end if @@ -266,8 +279,7 @@ contains call inventory%initialise(name=inventory_name) end if - mesh => mesh_collection%get_mesh(mesh_id) - constant_exists = inventory%paired_object_exists(mesh_id) + constant_exists = inventory%paired_object_exists(mesh%get_id()) if (.not. constant_exists) then ! Create constant if it doesn't already exist @@ -285,21 +297,33 @@ contains case (W1) call invoke( name="compute_w1_mass_matrix_fe", & compute_mass_matrix_kernel_w1_type(mm_op, chi, & - panel_id, qr_ptr) ) + panel_id, & + geometry, topology, & + coord_system, & + scaled_radius, & + qr_ptr) ) case (W2, W2broken, W2H, W2V) call invoke( name="compute_w2_mass_matrix_fe", & compute_mass_matrix_kernel_w2_type(mm_op, chi, & - panel_id, qr_ptr) ) + panel_id, & + geometry, topology, & + coord_system, & + scaled_radius, & + qr_ptr) ) case (W3, Wtheta) call invoke( name= "compute_scalar_mass_matrix_fe", & compute_mass_matrix_kernel_w_scalar_type(mm_op, chi, & panel_id, & extend_mesh, & + geometry, & + topology, & + coord_system, & + scaled_radius, & qr_ptr) ) end select ! Lump W2 mass matrix along the edge of the domain to avoid boundary ! errors propagating into the domain with each iteration - if ( (space == W2 .or. space == W2h) .and. topology /= topology_fully_periodic ) then + if ( (space == W2 .or. space == W2h) .and. topology /= topology_periodic ) then call dummy_field%initialise( fs, halo_depth=2 ) call invoke( edge_lump_w2_mass_matrix_kernel_type(mm_op, & dummy_field, & @@ -315,10 +339,11 @@ contains end function get_mass_matrix_fe !> @brief Returns a pointer to a finite volume mass matrix + !> @param[in] config Application configuration object + !> @param[in] mesh Mesh to get the object for !> @param[in] space The space of the desired mass matrix - !> @param[in] mesh_id The ID of the mesh to get the object for !> @return The mass matrix operator - function get_mass_matrix_fv(space, mesh_id) result(mm_op) + function get_mass_matrix_fv(config, mesh, space) result(mm_op) use sci_compute_mass_matrix_kernel_w_scalar_mod, & only: compute_mass_matrix_kernel_w_scalar_type @@ -328,13 +353,13 @@ contains only: compute_mass_matrix_kernel_w1_type use sci_edge_lump_w2_mass_matrix_kernel_mod, & only: edge_lump_w2_mass_matrix_kernel_type - use base_mesh_config_mod, only: topology, topology_fully_periodic implicit none - integer(kind=i_def), intent(in) :: space - integer(kind=i_def), intent(in) :: mesh_id - type(mesh_type), pointer :: mesh + type(config_type), intent(in) :: config + type(mesh_type), intent(in) :: mesh + integer(kind=i_def), intent(in) :: space + type(operator_type), pointer :: mm_op logical(kind=l_def) :: constant_exists type(field_type), pointer :: chi(:) @@ -348,6 +373,17 @@ contains integer(kind=i_def), parameter :: stencil_depth = 1_i_def integer(tik) :: id + integer(i_def) :: geometry + integer(i_def) :: topology + integer(i_def) :: coord_system + real(r_def) :: scaled_radius + + geometry = mesh%geometry() + topology = mesh%topology() + + coord_system = config%finite_element%coord_system() + scaled_radius = config%planet%scaled_radius() + ! Point to appropriate inventory for this space select case (space) case (W1) @@ -378,8 +414,7 @@ contains call inventory%initialise(name=inventory_name) end if - mesh => mesh_collection%get_mesh(mesh_id) - constant_exists = inventory%paired_object_exists(mesh_id) + constant_exists = inventory%paired_object_exists(mesh%get_id()) if (.not. constant_exists) then ! Create constant if it doesn't already exist @@ -396,21 +431,33 @@ contains case (W1) call invoke( name="compute_w1_mass_matrix_fv", & compute_mass_matrix_kernel_w1_type(mm_op, chi, & - panel_id, qr_ptr) ) + panel_id, & + geometry, topology, & + coord_system, & + scaled_radius, & + qr_ptr) ) case (W2, W2broken, W2H, W2V) call invoke( name="compute_w2_mass_matrix_fv", & compute_mass_matrix_kernel_w2_type(mm_op, chi, & - panel_id, qr_ptr) ) + panel_id, & + geometry, topology, & + coord_system, & + scaled_radius, & + qr_ptr) ) case (W3, Wtheta) call invoke( name= "compute_scalar_mass_matrix_fv", & compute_mass_matrix_kernel_w_scalar_type(mm_op, chi, & panel_id, & extend_mesh, & + geometry, & + topology, & + coord_system, & + scaled_radius, & qr_ptr) ) end select ! Lump W2 mass matrix along the edge of the domain to avoid boundary ! errors propagating into the domain with each iteration - if ( (space == W2 .or. space == W2h) .and. topology /= topology_fully_periodic ) then + if ( (space == W2 .or. space == W2h) .and. topology /= topology_periodic ) then call dummy_field%initialise( fs, halo_depth=2 ) call invoke( edge_lump_w2_mass_matrix_kernel_type(mm_op, & dummy_field, & @@ -426,18 +473,21 @@ contains end function get_mass_matrix_fv !> @brief Returns a pointer to a diagonal finite element mass matrix + !> @param[in] config Application configuration object + !> @param[in] mesh Mesh to get the object for !> @param[in] space The space of the desired diagonal mass matrix - !> @param[in] mesh_id The ID of the mesh to get the object for + !> @return The mass matrix diagonal field - function get_mass_matrix_diagonal_fe(space, mesh_id) result(diagonal_mm) + function get_mass_matrix_diagonal_fe(config, mesh, space) result(diagonal_mm) - use sci_mm_diagonal_kernel_mod, only: mm_diagonal_kernel_type + use sci_mm_diagonal_kernel_mod, only: mm_diagonal_kernel_type implicit none - integer(kind=i_def), intent(in) :: space - integer(kind=i_def), intent(in) :: mesh_id - type(mesh_type), pointer :: mesh + type(config_type), intent(in) :: config + type(mesh_type), intent(in) :: mesh + integer(kind=i_def), intent(in) :: space + type(field_type), pointer :: diagonal_mm type(operator_type), pointer :: mass_matrix logical(kind=l_def) :: constant_exists @@ -448,7 +498,7 @@ contains ! If running at lowest order, use finite volume if (element_order_h == 0 .and. element_order_v == 0) then - diagonal_mm => get_mass_matrix_diagonal_fv(space, mesh_id) + diagonal_mm => get_mass_matrix_diagonal_fv(config, mesh, space) return end if @@ -482,19 +532,18 @@ contains call inventory%initialise(name=inventory_name) end if - mesh => mesh_collection%get_mesh(mesh_id) - constant_exists = inventory%paired_object_exists(mesh_id) + constant_exists = inventory%paired_object_exists(mesh%get_id()) if (.not. constant_exists) then ! Create constant if it doesn't already exist - mass_matrix => get_mass_matrix_fe(space, mesh_id) + mass_matrix => get_mass_matrix_fe(config, mesh, space) if ( LPROF ) call start_timing( id, 'runtime_constants.fem' ) fs => function_space_collection%get_fs(mesh, element_order_h, & element_order_v, space) call inventory%add_field(diagonal_mm, fs, mesh) - call invoke ( setval_c(diagonal_mm, 0.0_r_def), & + call invoke ( setval_c(diagonal_mm, 0.0_r_def), & mm_diagonal_kernel_type(diagonal_mm, mass_matrix) ) if ( LPROF ) call stop_timing( id, 'runtime_constants.fem' ) @@ -506,18 +555,20 @@ contains end function get_mass_matrix_diagonal_fe !> @brief Returns a pointer to a diagonal finite volume mass matrix + !> @param[in] config Application configuration object + !> @param[in] mesh Mesh to get the object for !> @param[in] space The space of the desired diagonal mass matrix - !> @param[in] mesh_id The ID of the mesh to get the object for !> @return The mass matrix diagonal field - function get_mass_matrix_diagonal_fv(space, mesh_id) result(diagonal_mm) + function get_mass_matrix_diagonal_fv(config, mesh, space) result(diagonal_mm) use sci_mm_diagonal_kernel_mod, only: mm_diagonal_kernel_type implicit none - integer(kind=i_def), intent(in) :: space - integer(kind=i_def), intent(in) :: mesh_id - type(mesh_type), pointer :: mesh + type(config_type), intent(in) :: config + type(mesh_type), intent(in) :: mesh + integer(kind=i_def), intent(in) :: space + type(field_type), pointer :: diagonal_mm type(operator_type), pointer :: mass_matrix logical(kind=l_def) :: constant_exists @@ -556,18 +607,17 @@ contains call inventory%initialise(name=inventory_name) end if - mesh => mesh_collection%get_mesh(mesh_id) - constant_exists = inventory%paired_object_exists(mesh_id) + constant_exists = inventory%paired_object_exists(mesh%get_id()) if (.not. constant_exists) then ! Create constant if it doesn't already exist - mass_matrix => get_mass_matrix_fv(space, mesh_id) + mass_matrix => get_mass_matrix_fv(config, mesh, space) if ( LPROF ) call start_timing( id, 'runtime_constants.fem' ) fs => function_space_collection%get_fs(mesh, 0, 0, space) call inventory%add_field(diagonal_mm, fs, mesh) - call invoke ( setval_c(diagonal_mm, 0.0_r_def), & + call invoke ( setval_c(diagonal_mm, 0.0_r_def), & mm_diagonal_kernel_type(diagonal_mm, mass_matrix) ) if ( LPROF ) call stop_timing( id, 'runtime_constants.fem' ) @@ -580,18 +630,20 @@ contains !> @brief Returns a pointer to a field for a finite element inverse lumped !! mass matrix + !> @param[in] config Application configuration object + !> @param[in] mesh Mesh to get the object for !> @param[in] space The space of the desired lumped inverse mass matrix - !> @param[in] mesh_id The ID of the mesh to get the object for !> @return The lumped inverse mass matrix - function get_inverse_lumped_mass_matrix_fe(space, mesh_id) result(mm_linv) + function get_inverse_lumped_mass_matrix_fe(config, mesh, space) result(mm_linv) use dg_matrix_vector_kernel_mod, only: dg_matrix_vector_kernel_type implicit none - integer(kind=i_def), intent(in) :: space - integer(kind=i_def), intent(in) :: mesh_id - type(mesh_type), pointer :: mesh + type(config_type), intent(in) :: config + type(mesh_type), intent(in) :: mesh + integer(kind=i_def), intent(in) :: space + type(field_type) :: ones type(field_type), pointer :: mm_linv type(operator_type), pointer :: mass_matrix @@ -604,7 +656,7 @@ contains ! If running at lowest order, use finite volume if (element_order_h == 0 .and. element_order_v == 0) then - mm_linv => get_inverse_lumped_mass_matrix_fv(space, mesh_id) + mm_linv => get_inverse_lumped_mass_matrix_fv(config, mesh, space) return end if @@ -626,12 +678,11 @@ contains call inventory%initialise(name=inventory_name) end if - mesh => mesh_collection%get_mesh(mesh_id) - constant_exists = inventory%paired_object_exists(mesh_id) + constant_exists = inventory%paired_object_exists(mesh%get_id()) if (.not. constant_exists) then ! Create constant if it doesn't already exist - mass_matrix => get_mass_matrix_fe(space, mesh_id) + mass_matrix => get_mass_matrix_fe(config, mesh, space) if ( LPROF ) call start_timing( id, 'runtime_constants.fem' ) @@ -664,18 +715,20 @@ contains !> @brief Returns a pointer to a field for a finite volume inverse lumped !! mass matrix + !> @param[in] config Application configuration object + !> @param[in] mesh Mesh to get the object for !> @param[in] space The space of the desired lumped inverse mass matrix - !> @param[in] mesh_id The ID of the mesh to get the object for !> @return The lumped inverse mass matrix - function get_inverse_lumped_mass_matrix_fv(space, mesh_id) result(mm_linv) + function get_inverse_lumped_mass_matrix_fv(config, mesh, space) result(mm_linv) use dg_matrix_vector_kernel_mod, only: dg_matrix_vector_kernel_type implicit none - integer(kind=i_def), intent(in) :: space - integer(kind=i_def), intent(in) :: mesh_id - type(mesh_type), pointer :: mesh + type(config_type), intent(in) :: config + type(mesh_type), intent(in) :: mesh + integer(kind=i_def), intent(in) :: space + type(field_type) :: ones type(field_type), pointer :: mm_linv type(operator_type), pointer :: mass_matrix @@ -704,12 +757,11 @@ contains call inventory%initialise(name=inventory_name) end if - mesh => mesh_collection%get_mesh(mesh_id) - constant_exists = inventory%paired_object_exists(mesh_id) + constant_exists = inventory%paired_object_exists(mesh%get_id()) if (.not. constant_exists) then ! Create constant if it doesn't already exist - mass_matrix => get_mass_matrix_fv(space, mesh_id) + mass_matrix => get_mass_matrix_fv(config, mesh, space) if ( LPROF ) call start_timing( id, 'runtime_constants.fem' ) @@ -741,17 +793,20 @@ contains !> @brief Returns a pointer to the finite element inverse W3 mass matrix !! operator - !> @param[in] mesh_id The ID of the mesh to get the object for + !> @param[in] config Application configuration object + !> @param[in] mesh Mesh to get the object for + !> @param[in] space The space of the desired inverse mass matrix !> @return The inverse W3 mass matrix operator - function get_inverse_mass_matrix_fe(space, mesh_id) result(mm_inv) + function get_inverse_mass_matrix_fe(config, mesh, space) result(mm_inv) use sci_invert_local_operator_kernel_mod, only: invert_local_operator_kernel_type implicit none - integer(kind=i_def), intent(in) :: space - integer(kind=i_def), intent(in) :: mesh_id - type(mesh_type), pointer :: mesh + type(config_type), intent(in) :: config + type(mesh_type), intent(in) :: mesh + integer(kind=i_def), intent(in) :: space + type(operator_type), pointer :: mass_matrix type(operator_type), pointer :: mm_inv logical(kind=l_def) :: constant_exists @@ -762,7 +817,7 @@ contains ! If running at lowest order, use finite volume if (element_order_h == 0 .and. element_order_v == 0) then - mm_inv => get_inverse_mass_matrix_fv(space, mesh_id) + mm_inv => get_inverse_mass_matrix_fv(config, mesh, space) return end if @@ -784,12 +839,11 @@ contains call inventory%initialise(name=inventory_name) end if - mesh => mesh_collection%get_mesh(mesh_id) - constant_exists = inventory%paired_object_exists(mesh_id) + constant_exists = inventory%paired_object_exists(mesh%get_id()) if (.not. constant_exists) then ! Create constant if it doesn't already exist - mass_matrix => get_mass_matrix_fe(space, mesh_id) + mass_matrix => get_mass_matrix_fe(config, mesh, space) if ( LPROF ) call start_timing( id, 'runtime_constants.fem' ) @@ -810,17 +864,20 @@ contains !> @brief Returns a pointer to the finite volume inverse W3 mass matrix !! operator - !> @param[in] mesh_id The ID of the mesh to get the object for + !> @param[in] config Application configuration object + !> @param[in] mesh Mesh to get the object for + !> @param[in] space The space of the desired operator !> @return The inverse W3 mass matrix operator - function get_inverse_mass_matrix_fv(space, mesh_id) result(mm_inv) + function get_inverse_mass_matrix_fv(config, mesh, space) result(mm_inv) use sci_invert_local_operator_kernel_mod, only: invert_local_operator_kernel_type implicit none - integer(kind=i_def), intent(in) :: space - integer(kind=i_def), intent(in) :: mesh_id - type(mesh_type), pointer :: mesh + type(config_type), intent(in) :: config + type(mesh_type), intent(in) :: mesh + integer(kind=i_def), intent(in) :: space + type(operator_type), pointer :: mass_matrix type(operator_type), pointer :: mm_inv logical(kind=l_def) :: constant_exists @@ -847,12 +904,11 @@ contains call inventory%initialise(name=inventory_name) end if - mesh => mesh_collection%get_mesh(mesh_id) - constant_exists = inventory%paired_object_exists(mesh_id) + constant_exists = inventory%paired_object_exists(mesh%get_id()) if (.not. constant_exists) then ! Create constant if it doesn't already exist - mass_matrix => get_mass_matrix_fv(space, mesh_id) + mass_matrix => get_mass_matrix_fv(config, mesh, space) if ( LPROF ) call start_timing( id, 'runtime_constants.fem' ) @@ -871,16 +927,18 @@ contains end function get_inverse_mass_matrix_fv !> @brief Returns a pointer to the finite element curl operator - !> @param[in] mesh_id The ID of the mesh to get the object for + !> @param[in] config Application configuration object + !> @param[in] mesh Mesh to get the object for !> @return The curl operator - function get_curl_fe(mesh_id) result(curl) + function get_curl_fe(config, mesh) result(curl) use sci_compute_curl_operator_kernel_mod, only: compute_curl_operator_kernel_type implicit none - integer(kind=i_def), intent(in) :: mesh_id - type(mesh_type), pointer :: mesh + type(config_type), intent(in) :: config + type(mesh_type), intent(in) :: mesh + type(operator_type), pointer :: curl logical(kind=l_def) :: constant_exists type(field_type), pointer :: chi(:) @@ -890,9 +948,24 @@ contains type(function_space_type), pointer :: w1_fs integer(tik) :: id + integer(i_def) :: geometry + integer(i_def) :: topology + integer(i_def) :: coord_system + integer(i_def) :: order_h + integer(i_def) :: order_v + real(r_def) :: scaled_radius + + geometry = mesh%geometry() + topology = mesh%topology() + + coord_system = config%finite_element%coord_system() + scaled_radius = config%planet%scaled_radius() + order_h = config%finite_element%element_order_h() + order_v = config%finite_element%element_order_v() + ! If running at lowest order, use finite volume - if (element_order_h == 0 .and. element_order_v == 0) then - curl => get_curl_fv(mesh_id) + if (order_h == 0 .and. order_v == 0) then + curl => get_curl_fv(config, mesh) return end if @@ -901,8 +974,7 @@ contains call curl_inventory_fe%initialise(name="curl_fe") end if - mesh => mesh_collection%get_mesh(mesh_id) - constant_exists = curl_inventory_fe%paired_object_exists(mesh_id) + constant_exists = curl_inventory_fe%paired_object_exists(mesh%get_id()) ! Create constant if it doesn't already exist if (.not. constant_exists) then @@ -910,16 +982,16 @@ contains panel_id => get_panel_id(mesh) qr_ptr => get_qr_fe() - w2_fs => function_space_collection%get_fs( mesh, element_order_h, & - element_order_v, W2 ) - w1_fs => function_space_collection%get_fs( mesh, element_order_h, & - element_order_v, W1 ) + w2_fs => function_space_collection%get_fs( mesh, order_h, order_v, W2 ) + w1_fs => function_space_collection%get_fs( mesh, order_h, order_v, W1 ) if ( LPROF ) call start_timing( id, 'runtime_constants.fem' ) call curl_inventory_fe%add_operator(curl, w2_fs, w1_fs, mesh) - call invoke( name='calculate_curl_fe', & - compute_curl_operator_kernel_type(curl, chi, & - panel_id, qr_ptr) ) + call invoke( name='calculate_curl_fe', & + compute_curl_operator_kernel_type(curl, chi, panel_id, & + geometry, topology, & + coord_system, & + scaled_radius, qr_ptr) ) if ( LPROF ) call stop_timing( id, 'runtime_constants.fem' ) end if @@ -928,16 +1000,18 @@ contains end function get_curl_fe !> @brief Returns a pointer to the finite volume curl operator - !> @param[in] mesh_id The ID of the mesh to get the object for + !> @param[in] config Application configuration object + !> @param[in] mesh Mesh to get the object for !> @return The curl operator - function get_curl_fv(mesh_id) result(curl) + function get_curl_fv(config, mesh) result(curl) use sci_compute_curl_operator_kernel_mod, only: compute_curl_operator_kernel_type implicit none - integer(kind=i_def), intent(in) :: mesh_id - type(mesh_type), pointer :: mesh + type(config_type), intent(in) :: config + type(mesh_type), intent(in) :: mesh + type(operator_type), pointer :: curl logical(kind=l_def) :: constant_exists type(field_type), pointer :: chi(:) @@ -947,13 +1021,23 @@ contains type(function_space_type), pointer :: w1_fs integer(tik) :: id + integer(i_def) :: geometry + integer(i_def) :: topology + integer(i_def) :: coord_system + real(r_def) :: scaled_radius + + geometry = mesh%geometry() + topology = mesh%topology() + + coord_system = config%finite_element%coord_system() + scaled_radius = config%planet%scaled_radius() + ! Initialise inventory if it hasn't been done so already if (.not. curl_inventory_fv%is_initialised()) then call curl_inventory_fv%initialise(name="curl_fv") end if - mesh => mesh_collection%get_mesh(mesh_id) - constant_exists = curl_inventory_fv%paired_object_exists(mesh_id) + constant_exists = curl_inventory_fv%paired_object_exists(mesh%get_id()) ! Create constant if it doesn't already exist if (.not. constant_exists) then @@ -966,9 +1050,11 @@ contains if ( LPROF ) call start_timing( id, 'runtime_constants.fem' ) call curl_inventory_fv%add_operator(curl, w2_fs, w1_fs, mesh) - call invoke( name='calculate_curl_fv', & - compute_curl_operator_kernel_type(curl, chi, & - panel_id, qr_ptr) ) + call invoke( name='calculate_curl_fv', & + compute_curl_operator_kernel_type(curl, chi, panel_id, & + geometry, topology, & + coord_system, & + scaled_radius, qr_ptr) ) if ( LPROF ) call stop_timing( id, 'runtime_constants.fem' ) end if @@ -977,14 +1063,16 @@ contains end function get_curl_fv !> @brief Returns a pointer to the div operator - !> @param[in] mesh_id The ID of the mesh to get the object for + !> @param[in] config Application configuration object + !> @param[in] mesh Mesh to get the object for !> @return The div operator - function get_div(mesh_id) result(div) + function get_div(config, mesh) result(div) implicit none - integer(kind=i_def), intent(in) :: mesh_id - type(mesh_type), pointer :: mesh + type(config_type), intent(in) :: config + type(mesh_type), intent(in) :: mesh + type(operator_type), pointer :: div logical(kind=l_def) :: constant_exists type(field_type), pointer :: chi(:) @@ -994,13 +1082,29 @@ contains type(function_space_type), pointer :: w2_fs integer(tik) :: id + integer(i_def) :: geometry + integer(i_def) :: topology + integer(i_def) :: coord_system + integer(i_def) :: order_h + integer(i_def) :: order_v + real(r_def) :: scaled_radius + logical(l_def) :: rehabilitate + + geometry = mesh%geometry() + topology = mesh%topology() + + coord_system = config%finite_element%coord_system() + scaled_radius = config%planet%scaled_radius() + rehabilitate = config%finite_element%rehabilitate() + order_h = config%finite_element%element_order_h() + order_v = config%finite_element%element_order_v() + ! Initialise inventory if it hasn't been done so already if (.not. div_inventory%is_initialised()) then call div_inventory%initialise(name="div") end if - mesh => mesh_collection%get_mesh(mesh_id) - constant_exists = div_inventory%paired_object_exists(mesh_id) + constant_exists = div_inventory%paired_object_exists(mesh%get_id()) ! Create constant if it doesn't already exist if (.not. constant_exists) then @@ -1008,16 +1112,19 @@ contains panel_id => get_panel_id(mesh) qr_ptr => get_qr_fe() - w2_fs => function_space_collection%get_fs( mesh, element_order_h, & - element_order_v, W2 ) - w3_fs => function_space_collection%get_fs( mesh, element_order_h, & - element_order_v, W3 ) + w2_fs => function_space_collection%get_fs( mesh, order_h, order_v, W2 ) + w3_fs => function_space_collection%get_fs( mesh, order_h, order_v, W3 ) if ( LPROF ) call start_timing( id, 'runtime_constants.fem' ) call div_inventory%add_operator(div, w3_fs, w2_fs, mesh) - call invoke( name='calculate_div', & - compute_div_operator_kernel_type(div, chi, & - panel_id, qr_ptr) ) + call invoke( name='calculate_div', & + compute_div_operator_kernel_type(div, chi, & + panel_id, & + geometry, topology, & + coord_system, & + scaled_radius, & + rehabilitate, & + qr_ptr) ) if ( LPROF ) call stop_timing( id, 'runtime_constants.fem' ) end if @@ -1026,14 +1133,16 @@ contains end function get_div !> @brief Returns a pointer to the div operator - !> @param[in] mesh_id The ID of the mesh to get the object for + !> @param[in] config Application configuration object + !> @param[in] mesh Mesh to get the object for !> @return The div operator for W2H - function get_div_h(mesh_id) result(div_h) + function get_div_h(config, mesh) result(div_h) implicit none - integer(kind=i_def), intent(in) :: mesh_id - type(mesh_type), pointer :: mesh + type(config_type), intent(in) :: config + type(mesh_type), intent(in) :: mesh + type(operator_type), pointer :: div_h logical(kind=l_def) :: constant_exists type(field_type), pointer :: chi(:) @@ -1043,13 +1152,29 @@ contains type(function_space_type), pointer :: w2h_fs integer(tik) :: id + integer(i_def) :: geometry + integer(i_def) :: topology + integer(i_def) :: coord_system + integer(i_def) :: order_h + integer(i_def) :: order_v + real(r_def) :: scaled_radius + logical(l_def) :: rehabilitate + + geometry = mesh%geometry() + topology = mesh%topology() + + coord_system = config%finite_element%coord_system() + scaled_radius = config%planet%scaled_radius() + rehabilitate = config%finite_element%rehabilitate() + order_h = config%finite_element%element_order_h() + order_v = config%finite_element%element_order_v() + ! Initialise inventory if it hasn't been done so already if (.not. div_h_inventory%is_initialised()) then call div_h_inventory%initialise(name="div_h") end if - mesh => mesh_collection%get_mesh(mesh_id) - constant_exists = div_h_inventory%paired_object_exists(mesh_id) + constant_exists = div_h_inventory%paired_object_exists(mesh%get_id()) ! Create constant if it doesn't already exist if (.not. constant_exists) then @@ -1057,16 +1182,19 @@ contains panel_id => get_panel_id(mesh) qr_ptr => get_qr_fe() - w2h_fs => function_space_collection%get_fs( mesh, element_order_h, & - element_order_v, W2H ) - w3_fs => function_space_collection%get_fs( mesh, element_order_h, & - element_order_v, W3 ) + w2h_fs => function_space_collection%get_fs( mesh, order_h, order_v, W2H ) + w3_fs => function_space_collection%get_fs( mesh, order_h, order_v, W3 ) if ( LPROF ) call start_timing( id, 'runtime_constants.fem' ) call div_h_inventory%add_operator(div_h, w3_fs, w2h_fs, mesh) - call invoke( name='calculate_div_h', & - compute_div_operator_kernel_type(div_h, chi, & - panel_id, qr_ptr) ) + call invoke( name='calculate_div_h', & + compute_div_operator_kernel_type(div_h, chi, & + panel_id, & + geometry, topology, & + coord_system, & + scaled_radius, & + rehabilitate, & + qr_ptr) ) if ( LPROF ) call stop_timing( id, 'runtime_constants.fem' ) end if @@ -1075,16 +1203,18 @@ contains end function get_div_h !> @brief Returns a pointer to the finite element inv_m3 * div operator - !> @param[in] mesh_id The ID of the mesh to get the object for + !> @param[in] config Application configuration object + !> @param[in] mesh Mesh to get the object for !> @return The div operator - function get_im3_div_fe(mesh_id) result(im3_div) + function get_im3_div_fe(config, mesh) result(im3_div) use sci_operator_x_times_y_kernel_mod, only: operator_x_times_y_kernel_type implicit none - integer(kind=i_def), intent(in) :: mesh_id - type(mesh_type), pointer :: mesh + type(config_type), intent(in) :: config + type(mesh_type), intent(in) :: mesh + type(operator_type), pointer :: im3_div type(operator_type), pointer :: mm_w3_inv type(operator_type) :: div @@ -1096,9 +1226,26 @@ contains type(function_space_type), pointer :: w2_fs integer(tik) :: id + integer(i_def) :: geometry + integer(i_def) :: topology + integer(i_def) :: coord_system + integer(i_def) :: order_h + integer(i_def) :: order_v + real(r_def) :: scaled_radius + logical(l_def) :: rehabilitate + + geometry = mesh%geometry() + topology = mesh%topology() + + coord_system = config%finite_element%coord_system() + scaled_radius = config%planet%scaled_radius() + rehabilitate = config%finite_element%rehabilitate() + order_h = config%finite_element%element_order_h() + order_v = config%finite_element%element_order_v() + ! If running at lowest order, use finite volume - if (element_order_h == 0 .and. element_order_v == 0) then - im3_div => get_im3_div_fv(mesh_id) + if (order_h == 0 .and. order_v == 0) then + im3_div => get_im3_div_fv(config, mesh) return end if @@ -1107,26 +1254,28 @@ contains call im3_div_inventory_fe%initialise(name="im3_div_fe") end if - mesh => mesh_collection%get_mesh(mesh_id) - constant_exists = im3_div_inventory_fe%paired_object_exists(mesh_id) + constant_exists = im3_div_inventory_fe%paired_object_exists(mesh%get_id()) ! Create constant if it doesn't already exist if (.not. constant_exists) then chi => get_coordinates(mesh) panel_id => get_panel_id(mesh) qr_ptr => get_qr_fe() - mm_w3_inv => get_inverse_mass_matrix_fe(W3, mesh_id) + mm_w3_inv => get_inverse_mass_matrix_fe(config, mesh, W3) - w2_fs => function_space_collection%get_fs( mesh, element_order_h, & - element_order_v, W2 ) - w3_fs => function_space_collection%get_fs( mesh, element_order_h, & - element_order_v, W3 ) + w2_fs => function_space_collection%get_fs( mesh, order_h, order_v, W2 ) + w3_fs => function_space_collection%get_fs( mesh, order_h, order_v, W3 ) if ( LPROF ) call start_timing( id, 'runtime_constants.fem' ) call div%initialise( w3_fs, w2_fs ) call im3_div_inventory_fe%add_operator(im3_div, w3_fs, w2_fs, mesh) - call invoke( compute_div_operator_kernel_type(div, chi, & - panel_id, qr_ptr), & + call invoke( compute_div_operator_kernel_type(div, chi, & + panel_id, & + geometry, topology, & + coord_system, & + scaled_radius, & + rehabilitate, & + qr_ptr), & operator_x_times_y_kernel_type(im3_div, mm_w3_inv, div) ) if ( LPROF ) call stop_timing( id, 'runtime_constants.fem' ) end if @@ -1136,16 +1285,18 @@ contains end function get_im3_div_fe !> @brief Returns a pointer to the finite volume inv_m3 * div operator - !> @param[in] mesh_id The ID of the mesh to get the object for + !> @param[in] config Application configuration object + !> @param[in] mesh Mesh to get the object for !> @return The div operator - function get_im3_div_fv(mesh_id) result(im3_div) + function get_im3_div_fv(config, mesh) result(im3_div) use sci_operator_x_times_y_kernel_mod, only: operator_x_times_y_kernel_type implicit none - integer(kind=i_def), intent(in) :: mesh_id - type(mesh_type), pointer :: mesh + type(config_type), intent(in) :: config + type(mesh_type), intent(in) :: mesh + type(operator_type), pointer :: im3_div type(operator_type), pointer :: mm_w3_inv type(operator_type) :: div @@ -1157,20 +1308,32 @@ contains type(function_space_type), pointer :: w2_fs integer(tik) :: id + integer(i_def) :: geometry + integer(i_def) :: topology + integer(i_def) :: coord_system + real(r_def) :: scaled_radius + logical(l_def) :: rehabilitate + + geometry = mesh%geometry() + topology = mesh%topology() + + coord_system = config%finite_element%coord_system() + scaled_radius = config%planet%scaled_radius() + rehabilitate = config%finite_element%rehabilitate() + ! Initialise inventory if it hasn't been done so already if (.not. im3_div_inventory_fv%is_initialised()) then call im3_div_inventory_fv%initialise(name="im3_div_fv") end if - mesh => mesh_collection%get_mesh(mesh_id) - constant_exists = im3_div_inventory_fv%paired_object_exists(mesh_id) + constant_exists = im3_div_inventory_fv%paired_object_exists(mesh%get_id()) ! Create constant if it doesn't already exist if (.not. constant_exists) then chi => get_coordinates(mesh) panel_id => get_panel_id(mesh) qr_ptr => get_qr_fv() - mm_w3_inv => get_inverse_mass_matrix_fv(W3, mesh_id) + mm_w3_inv => get_inverse_mass_matrix_fv(config, mesh, W3) w2_fs => function_space_collection%get_fs( mesh, 0, 0, W2 ) w3_fs => function_space_collection%get_fs( mesh, 0, 0, W3 ) @@ -1178,8 +1341,13 @@ contains if ( LPROF ) call start_timing( id, 'runtime_constants.fem' ) call div%initialise( w3_fs, w2_fs ) call im3_div_inventory_fv%add_operator(im3_div, w3_fs, w2_fs, mesh) - call invoke( compute_div_operator_kernel_type(div, chi, & - panel_id, qr_ptr), & + call invoke( compute_div_operator_kernel_type(div, chi, & + panel_id, & + geometry, topology, & + coord_system, & + scaled_radius, & + rehabilitate, & + qr_ptr), & operator_x_times_y_kernel_type(im3_div, mm_w3_inv, div) ) if ( LPROF ) call stop_timing( id, 'runtime_constants.fem' ) end if diff --git a/components/science/source/algorithm/sci_galerkin_projection_alg_mod.x90 b/components/science/source/algorithm/sci_galerkin_projection_alg_mod.x90 index c726e3cf3..5c5fe4da1 100644 --- a/components/science/source/algorithm/sci_galerkin_projection_alg_mod.x90 +++ b/components/science/source/algorithm/sci_galerkin_projection_alg_mod.x90 @@ -9,6 +9,7 @@ !> module sci_galerkin_projection_alg_mod + use config_mod, only: config_type use constants_mod, only: r_def, i_def use function_space_collection_mod, only: function_space_collection use field_mod, only: field_type @@ -17,6 +18,7 @@ module sci_galerkin_projection_alg_mod use sci_gp_vector_rhs_kernel_mod, only: gp_vector_rhs_kernel_type use log_mod, only: log_event, log_scratch_space, & log_level_debug, log_level_error + use mesh_mod, only: mesh_type use operator_mod, only: operator_type use quadrature_xyoz_mod, only: quadrature_xyoz_type use sci_mass_matrix_solver_alg_mod, & @@ -43,6 +45,7 @@ contains !> into orthogonal components and the galerkin projection of each !> component is computed. !> + !> @param[in] config Application configuration object !> @param[out] f_out Field to project onto. !> @param[in] f_in Field to project. !> @param[in] chi Field entity co-ordinates. @@ -52,7 +55,8 @@ contains !> @todo Ideally this function would not take chi and panel_id, these would !> be available through the input field. That is some way off though. !> - subroutine galerkin_projection_algorithm( f_out, & + subroutine galerkin_projection_algorithm( config, & + f_out, & f_in, & chi, & panel_id, & @@ -60,6 +64,8 @@ contains implicit none + type(config_type), intent(in) :: config + ! Field to output type(field_type), intent(inout) :: f_out(:) @@ -76,13 +82,29 @@ contains type(field_type) :: w2_field integer(i_def) :: out_fs integer(i_def) :: idx + + type(mesh_type), pointer :: mesh + integer(i_def) :: element_order_h integer(i_def) :: element_order_v + integer(i_def) :: geometry + integer(i_def) :: topology + integer(i_def) :: coord_system + real(r_def) :: scaled_radius ! Create continuous fields to project data into out_fs = f_out(1)%which_function_space() + + mesh => f_in%get_mesh() element_order_h = f_in%get_element_order_h() element_order_v = f_in%get_element_order_v() + + geometry = mesh%geometry() + topology = mesh%topology() + + coord_system = config%finite_element%coord_system() + scaled_radius = config%planet%scaled_radius() + do idx = 1, size(f_out) ! A Galerkin projection doesn't make sense between meshes. ! This check could usefully be an assertion if Fortran gains support @@ -108,7 +130,9 @@ contains if ( size(f_out) == 1 ) then write( log_scratch_space, '(A)' ) ' scalar field ... ' call log_event( log_scratch_space, log_level_debug ) - call invoke( gp_rhs_kernel_type(rhs(1), f_in, chi, panel_id, qr) ) + call invoke( gp_rhs_kernel_type(rhs(1), f_in, chi, panel_id, & + geometry, topology, coord_system, & + scaled_radius, qr) ) else write( log_scratch_space, '(A)' ) ' vector field ... ' call log_event( log_scratch_space, log_level_debug ) @@ -121,14 +145,13 @@ contains ) & ) - call invoke( & - gp_vector_rhs_kernel_type( & - rhs, f_in, chi, panel_id, w2_field, qr & - ) & - ) + call invoke( & + gp_vector_rhs_kernel_type( rhs, f_in, chi, panel_id, w2_field, & + geometry, topology, coord_system, & + scaled_radius, qr ) ) end if do idx = 1, size(f_out) - call mass_matrix_solver_alg( f_out(idx), rhs(idx) ) + call mass_matrix_solver_alg( config, f_out(idx), rhs(idx) ) end do end if diff --git a/components/science/source/algorithm/sci_geometric_constants_mod.x90 b/components/science/source/algorithm/sci_geometric_constants_mod.x90 index decdc4332..68bb8a73c 100644 --- a/components/science/source/algorithm/sci_geometric_constants_mod.x90 +++ b/components/science/source/algorithm/sci_geometric_constants_mod.x90 @@ -21,6 +21,7 @@ module sci_geometric_constants_mod use function_space_collection_mod, only: function_space_collection use log_mod, only: log_event, LOG_LEVEL_ERROR use mesh_collection_mod, only: mesh_collection + use mesh_mod, only: geometry_spherical use timing_mod, only: start_timing, stop_timing, & tik, LPROF @@ -35,10 +36,6 @@ module sci_geometric_constants_mod use local_mesh_mod, only: local_mesh_type ! Configuration - 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 implicit none @@ -175,6 +172,7 @@ contains integer(kind=i_def) :: k_h, k_v integer(tik) :: id + integer(kind=i_def) :: geometry integer(kind=i_def) :: topology integer(kind=i_def) :: element_order_h integer(kind=i_def) :: element_order_v @@ -190,6 +188,9 @@ contains element_order_h = config%finite_element%element_order_h() element_order_v = config%finite_element%element_order_v() + geometry = mesh%geometry() + topology = mesh%topology() + f_lat = config%base_mesh%f_lat() f_lon = config%idealised%f_lon() @@ -214,18 +215,9 @@ contains chi => get_coordinates(mesh) panel_id => get_panel_id(mesh) - if (mesh%is_topology_periodic()) then - topology = topology_fully_periodic - else if (mesh%is_topology_non_periodic()) then - topology = topology_non_periodic - else - call log_event('compute_latlon: Unsupported topology selection',log_level_error) - end if - call invoke( compute_latlon_kernel_type(lat, long, chi, panel_id, & - geometry_spherical, & - topology, coord_system, & - scaled_radius) ) + geometry, topology, & + coord_system, scaled_radius) ) else call invoke( setval_c(lat, f_lat), & setval_c(long, f_lon) ) @@ -519,17 +511,28 @@ contains type(quadrature_rule_gaussian_type) :: quadrature_rule integer(tik) :: id - integer(kind=i_def) :: element_order_h, element_order_v - integer(kind=i_def) :: nqp_h_exact, nqp_v_exact - - element_order_h = config%finite_element%element_order_h() - element_order_v = config%finite_element%element_order_v() - nqp_h_exact = config%finite_element%nqp_h_exact() - nqp_v_exact = config%finite_element%nqp_v_exact() + integer(i_def) :: order_h + integer(i_def) :: order_v + integer(i_def) :: nqp_h_exact + integer(i_def) :: nqp_v_exact + integer(i_def) :: geometry + integer(i_def) :: topology + integer(i_def) :: coord_system + real(r_def) :: scaled_radius + + geometry = mesh%geometry() + topology = mesh%topology() + + order_h = config%finite_element%element_order_h() + order_v = config%finite_element%element_order_v() + nqp_h_exact = config%finite_element%nqp_h_exact() + nqp_v_exact = config%finite_element%nqp_v_exact() + coord_system = config%finite_element%coord_system() + scaled_radius = config%planet%scaled_radius() ! If running at lowest order, use finite volume - if (element_order_h == 0 .and. element_order_v == 0) then - detj_at_w3 => get_detj_at_w3_fv(mesh) + if (order_h == 0 .and. order_v == 0) then + detj_at_w3 => get_detj_at_w3_fv(config, mesh) return end if @@ -548,10 +551,7 @@ 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, order_h, 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 @@ -565,6 +565,10 @@ contains chi, & panel_id, & extended_mesh, & + geometry, & + topology, & + coord_system, & + scaled_radius, & qr), & setval_c(detj_at_w3, 0.0_r_def), & mm_diagonal_kernel_type(detj_at_w3, mm_w3) ) @@ -578,9 +582,10 @@ contains end function get_detj_at_w3_fe !> @brief Returns the (finite volume) Det(J) values at W3 dof locations - !> @param[in] mesh Mesh to get the object for + !> @param[in] config Application namelist configuration object + !> @param[in] mesh Mesh to get the object for !> @return The Det(J) field - function get_detj_at_w3_fv(mesh) result(detj_at_w3) + function get_detj_at_w3_fv(config, mesh) result(detj_at_w3) ! @TODO #4487: update these imports ! use sci_calc_detj_at_w3_kernel_mod, & @@ -594,7 +599,8 @@ contains implicit none - type(mesh_type), intent(in) :: mesh + type(config_type), intent(in) :: config + type(mesh_type), intent(in) :: mesh logical(kind=l_def) :: constant_exists type(field_type), pointer :: detj_at_w3 @@ -608,6 +614,17 @@ contains type(quadrature_rule_gaussian_type) :: quadrature_rule integer(tik) :: id + integer(i_def) :: geometry + integer(i_def) :: topology + integer(i_def) :: coord_system + real(r_def) :: scaled_radius + + geometry = mesh%geometry() + topology = mesh%topology() + + coord_system = config%finite_element%coord_system() + scaled_radius = config%planet%scaled_radius() + ! Check inventory is initialised if (.not. detj_at_w3_inventory_fv%is_initialised()) then ! Initialise all inventories together @@ -636,6 +653,10 @@ contains chi, & panel_id, & extended_mesh, & + geometry, & + topology, & + coord_system, & + scaled_radius, & qr), & setval_c(detj_at_w3, 0.0_r_def), & mm_diagonal_kernel_type(detj_at_w3, mm_w3) ) @@ -782,7 +803,7 @@ contains implicit none - type(config_type), intent(in) :: config + type(config_type), intent(in) :: config type(mesh_type), intent(in) :: mesh type(field_type), pointer :: dz_w3 @@ -972,8 +993,8 @@ contains implicit none - type(config_type), intent(in) :: config - type(mesh_type), intent(in), pointer :: mesh + type(config_type), intent(in) :: config + type(mesh_type), intent(in) :: mesh integer(kind=i_def) :: local_mesh_id @@ -995,11 +1016,7 @@ contains call dA_msl_proj_inventory%initialise(name="dA_msl_proj") end if - if (mesh%is_geometry_spherical()) then - geometry = geometry_spherical - else - geometry = geometry_planar - end if + geometry = mesh%geometry() planet_radius = config%extrusion%planet_radius() domain_height = config%extrusion%domain_height() @@ -1340,11 +1357,7 @@ contains real(kind=r_def) :: scaled_radius - if (mesh%is_geometry_spherical()) then - geometry= geometry_spherical - else - geometry= geometry_planar - end if + geometry = mesh%geometry() element_order_h = config%finite_element%element_order_h() element_order_v = config%finite_element%element_order_v() @@ -1470,11 +1483,7 @@ contains integer(kind=i_def) :: coord_system real(kind=r_def) :: scaled_radius - if (mesh%is_geometry_spherical()) then - geometry = geometry_spherical - else - geometry = geometry_planar - end if + geometry = mesh%geometry() coord_system = config%finite_element%coord_system() scaled_radius = config%planet%scaled_radius() diff --git a/components/science/source/algorithm/sci_mapping_constants_mod.x90 b/components/science/source/algorithm/sci_mapping_constants_mod.x90 index e295d9b75..5e2c67c0c 100644 --- a/components/science/source/algorithm/sci_mapping_constants_mod.x90 +++ b/components/science/source/algorithm/sci_mapping_constants_mod.x90 @@ -15,6 +15,7 @@ module sci_mapping_constants_mod ! Infrastructure + use config_mod, only: config_type use constants_mod, only: i_def, r_def, l_def, str_def use extrusion_mod, only: PRIME_EXTRUSION, & DOUBLE_LEVEL, & @@ -40,16 +41,6 @@ module sci_mapping_constants_mod use timing_mod, only: start_timing, stop_timing, & tik, LPROF - ! Object types - use config_mod, only: config_type - - ! Configuration modules - use base_mesh_config_mod, only: geometry_spherical, & - geometry_planar, & - topology_fully_periodic, & - topology_non_periodic - - ! Other algorithms use sci_geometric_constants_mod, only: get_coordinates, & get_panel_id @@ -230,18 +221,8 @@ contains integer(kind=i_def) :: coord_system real(kind=r_def) :: scaled_radius - if (mesh%is_geometry_spherical()) then - geometry = geometry_spherical - else - geometry = geometry_planar - end if - - if (mesh%is_topology_periodic()) then - topology = topology_fully_periodic - else - topology = topology_non_periodic - end if - + geometry = mesh%geometry() + topology = mesh%topology() coord_system = config%finite_element%coord_system() scaled_radius = config%planet%scaled_radius() @@ -309,18 +290,8 @@ contains integer(kind=i_def) :: coord_system real(kind=r_def) :: scaled_radius - if (mesh%is_geometry_spherical()) then - geometry = geometry_spherical - else - geometry = geometry_planar - end if - - if (mesh%is_topology_periodic()) then - topology = topology_fully_periodic - else - topology = topology_non_periodic - end if - + geometry = mesh%geometry() + topology = mesh%topology() coord_system = config%finite_element%coord_system() scaled_radius = config%planet%scaled_radius() @@ -696,16 +667,20 @@ contains end function get_intermesh_weights_w2 !> @brief Returns a pointer to the weights for conservative W3 mapping + !> @param[in] config Application configuration object !> @param[in] fine_mesh The fine mesh for the transform !> @param[in] coarse_mesh The coarse mesh for the transform !> @return The field containing weights for conservative W3 mapping - function get_intermesh_weights_w3_rdef(fine_mesh, coarse_mesh) result(weights_rdef) + function get_intermesh_weights_w3_rdef(config, fine_mesh, coarse_mesh) & + result(weights_rdef) use sci_fem_constants_mod, only: get_mass_matrix_diagonal_fv use sci_weights_intermesh_w3_kernel_mod, only: weights_intermesh_w3_kernel_type implicit none + type(config_type), intent(in) :: config + type(mesh_type), pointer, intent(in) :: coarse_mesh type(mesh_type), pointer, intent(in) :: fine_mesh integer(kind=i_def) :: intermesh_id @@ -741,8 +716,8 @@ contains end if ! @TODO #4487: this would be better as just the cell volumes - mm_w3_fine => get_mass_matrix_diagonal_fv(W3, fine_mesh%get_id()) - mm_w3_coarse => get_mass_matrix_diagonal_fv(W3, coarse_mesh%get_id()) + mm_w3_fine => get_mass_matrix_diagonal_fv(config, fine_mesh, W3) + mm_w3_coarse => get_mass_matrix_diagonal_fv(config, coarse_mesh, W3) if ( LPROF ) call start_timing( id, 'runtime_constants.mapping' ) @@ -768,13 +743,17 @@ contains end function get_intermesh_weights_w3_rdef !> @brief Returns a pointer to the weights for conservative W3 mapping + !> @param[in] config Application configuration object !> @param[in] fine_mesh The fine mesh for the transform !> @param[in] coarse_mesh The coarse mesh for the transform !> @return The field containing weights for conservative W3 mapping - function get_intermesh_weights_w3_rtran(fine_mesh, coarse_mesh) result(weights_rtran) + function get_intermesh_weights_w3_rtran(config, fine_mesh, coarse_mesh) & + result(weights_rtran) implicit none + type(config_type), intent(in) :: config + type(mesh_type), pointer, intent(in) :: coarse_mesh type(mesh_type), pointer, intent(in) :: fine_mesh integer(kind=i_def) :: intermesh_id @@ -800,7 +779,9 @@ contains if (.not. constant_exists) then ! Create the object as it doesn't exist yet - weights_rdef => get_intermesh_weights_w3_rdef(fine_mesh, coarse_mesh) + weights_rdef => get_intermesh_weights_w3_rdef(config, & + fine_mesh, & + coarse_mesh) if ( LPROF ) call start_timing( id, 'runtime_constants.mapping' ) @@ -1109,18 +1090,8 @@ contains integer(kind=i_def), parameter :: xdirection = 1_i_def integer(tik) :: id - if (mesh%is_geometry_spherical()) then - geometry = geometry_spherical - else - geometry = geometry_planar - end if - - if (mesh%is_topology_periodic()) then - topology = topology_fully_periodic - else - topology = topology_non_periodic - end if - + geometry = mesh%geometry() + topology = mesh%topology() coord_system = config%finite_element%coord_system() scaled_radius = config%planet%scaled_radius() @@ -1191,18 +1162,8 @@ contains integer(kind=i_def), parameter :: ydirection = 2_i_def integer(tik) :: id - if (mesh%is_geometry_spherical()) then - geometry = geometry_spherical - else - geometry = geometry_planar - end if - - if (mesh%is_topology_periodic()) then - topology = topology_fully_periodic - else - topology = topology_non_periodic - end if - + geometry = mesh%geometry() + topology = mesh%topology() coord_system = config%finite_element%coord_system() scaled_radius = config%planet%scaled_radius() @@ -1274,18 +1235,8 @@ contains integer(kind=i_def), parameter :: zdirection = 3_i_def integer(tik) :: id - if (mesh%is_geometry_spherical()) then - geometry = geometry_spherical - else - geometry = geometry_planar - end if - - if (mesh%is_topology_periodic()) then - topology = topology_fully_periodic - else - topology = topology_non_periodic - end if - + geometry = mesh%geometry() + topology = mesh%topology() coord_system = config%finite_element%coord_system() scaled_radius = config%planet%scaled_radius() @@ -1358,18 +1309,8 @@ contains integer(kind=i_def) :: coord_system real(kind=r_def) :: scaled_radius - if (mesh%is_geometry_spherical()) then - geometry = geometry_spherical - else - geometry = geometry_planar - end if - - if (mesh%is_topology_periodic()) then - topology = topology_fully_periodic - else - topology = topology_non_periodic - end if - + geometry = mesh%geometry() + topology = mesh%topology() coord_system = config%finite_element%coord_system() scaled_radius = config%planet%scaled_radius() diff --git a/components/science/source/algorithm/sci_project_output_mod.F90 b/components/science/source/algorithm/sci_project_output_mod.F90 index 8074c742d..6783b6409 100644 --- a/components/science/source/algorithm/sci_project_output_mod.F90 +++ b/components/science/source/algorithm/sci_project_output_mod.F90 @@ -20,16 +20,18 @@ module sci_project_output_mod !> @details This procedure uses the galerkin projection and a precomputed !> mass matrix to project a field !> + !> @param[in] config Application configuration object !> @param[in] field To be projected. !> @param[inout] projected_field Receives projection. !> @param[in] chi Field entity co-ordinates. !> @param[in] panel_id Cell orientation map. !> @param[in] output_fs Desired output function space. !> - subroutine project_output( field, projected_field, & + subroutine project_output( config, field, projected_field, & chi, panel_id, & output_fs ) + use config_mod, only: config_type use constants_mod, only: r_def, str_max_filename, i_def use field_mod, only: field_type use field_parent_mod, only: write_interface @@ -42,13 +44,18 @@ subroutine project_output( field, projected_field, & implicit none + type(config_type), intent(in) :: config + ! Input field to project from type(field_type), intent(in) :: field + ! Output field to project to type(field_type), intent(inout) :: projected_field(:) + ! Co-ordinate system type(field_type), intent(in) :: chi(:) type(field_type), intent(in) :: panel_id + ! Output function space integer(i_def), intent(in) :: output_fs @@ -84,9 +91,8 @@ subroutine project_output( field, projected_field, & end do ! do the projection - call galerkin_projection_algorithm( & - projected_field, field, chi, panel_id, qr & - ) + call galerkin_projection_algorithm(config, projected_field, field, & + chi, panel_id, qr) end subroutine project_output diff --git a/components/science/source/algorithm/solver/sci_dense_operator_alg_mod.x90 b/components/science/source/algorithm/solver/sci_dense_operator_alg_mod.x90 index de9e14b18..bfeaaea70 100644 --- a/components/science/source/algorithm/solver/sci_dense_operator_alg_mod.x90 +++ b/components/science/source/algorithm/solver/sci_dense_operator_alg_mod.x90 @@ -21,7 +21,8 @@ module sci_dense_operator_alg_mod only : compute_mass_matrix_kernel_w_scalar_type use sci_compute_mass_matrix_kernel_w3_mod, & only : compute_mass_matrix_kernel_w3_type - use constants_mod, only : r_def, i_def, l_def + use config_mod, only : config_type + use constants_mod, only : r_def, i_def, l_def, imdi use field_mod, only : field_type use function_space_mod, only : function_space_type use function_space_collection_mod, only : function_space_collection @@ -71,15 +72,18 @@ module sci_dense_operator_alg_mod contains !> @brief Construct a dense_operator object. !> + !> @param[in] config Application configuration object !> @param[in] fv A field_vector of size 2 !> @param[in] chi The coordinate field !> @param[in] panel_id The field with mesh panel IDs !> @param[in] qr The quadrature object !> @return self The constructed dense operator - function dense_operator_constructor( fv, qr, chi, panel_id ) result(self) + function dense_operator_constructor( config, fv, qr, chi, panel_id ) result(self) implicit none + type(config_type), intent(in) :: config + class(abstract_vector_type), intent(in) :: fv type( field_type ), dimension(3), intent(in) :: chi type( field_type ), intent(in) :: panel_id @@ -87,13 +91,23 @@ contains type(dense_operator_type) :: self + type(mesh_type), pointer :: mesh + type(function_space_type), pointer :: fs_wa => null() type(function_space_type), pointer :: fs_wb => null() - type(mesh_type), pointer :: mesh => null() + integer(kind=i_def) :: element_order_h, element_order_v, fs_label logical(kind=l_def) :: extended_mesh + logical(l_def) :: rehabilitate + integer(i_def) :: geometry, topology, coord_system + real(r_def) :: scaled_radius + + rehabilitate = config%finite_element%rehabilitate() + coord_system = config%finite_element%coord_system() + scaled_radius = config%planet%scaled_radius() + select type (fv) type is (field_vector_type) @@ -105,6 +119,8 @@ contains end if mesh => fv%vector(1)%get_mesh() + geometry = mesh%geometry() + topology = mesh%topology() ! Make the mass_w0 operator element_order_h = fv%vector(1)%get_element_order_h() @@ -125,8 +141,12 @@ contains extended_mesh = .false. call invoke( & compute_mass_matrix_kernel_w_scalar_type(self%mass_w0, chi, panel_id, & - extended_mesh, qr), & - compute_mass_matrix_kernel_w3_type(self%mass_w3, chi, panel_id, qr) ) + extended_mesh, geometry, & + topology, coord_system, & + scaled_radius, qr), & + compute_mass_matrix_kernel_w3_type(self%mass_w3, chi, panel_id, & + geometry, topology, coord_system, & + scaled_radius, rehabilitate, qr) ) class default @@ -141,15 +161,19 @@ contains !> @brief Applies the LMA mass matrices to the vector, y = M x. !> - !> @param[in,out] self The dense operator - !> @param[in] x Field vector to be read, size 2 - !> @param[in,out] y Field vector to be written, size 2 - subroutine apply_dense_op(self, x, y) + !> @param[in] config Application configuration object + !> @param[in] x Field vector to be read, size 2 + !> @param[in,out] y Field vector to be written, size 2 + subroutine apply_dense_op(self, config, x, y) implicit none - class(dense_operator_type), intent(inout) :: self - class(abstract_vector_type), intent(in) :: x - class(abstract_vector_type), intent(inout) :: y + + class(dense_operator_type), intent(inout) :: self + + type(config_type), intent(in) :: config + class(abstract_vector_type), intent(in) :: x + class(abstract_vector_type), intent(inout) :: y + type(field_type), pointer :: x_field1 => null(), & y_field1 => null(), & x_field2 => null(), & diff --git a/components/science/source/algorithm/solver/sci_diagonal_preconditioner_alg_mod.x90 b/components/science/source/algorithm/solver/sci_diagonal_preconditioner_alg_mod.x90 index a04d49a95..fcaad004e 100644 --- a/components/science/source/algorithm/solver/sci_diagonal_preconditioner_alg_mod.x90 +++ b/components/science/source/algorithm/solver/sci_diagonal_preconditioner_alg_mod.x90 @@ -7,6 +7,7 @@ !> @brief Application of the diagonal of an operator. module sci_diagonal_preconditioner_alg_mod + use config_mod, only: config_type use constants_mod, only: i_def, r_def use field_mod, only: field_type use log_mod, only: log_event, & @@ -65,16 +66,18 @@ contains !> @brief Apply diagonal preconditioner to a field to obtain \f$y=D^-1x\f$. !> - !> @param[in] self Instance of diagonal_preconditioner_type - !> @param[in] x Field \f$x\f$ to apply preconditioner to - !> @param[in,out] y Resulting field \f$y=D^-1x\f$ - subroutine apply_diagonal_preconditioner(self, x, y) + !> @param[in] config Application configuration object + !> @param[in] x Field \f$x\f$ to apply preconditioner to + !> @param[in,out] y Resulting field \f$y=D^-1x\f$ + subroutine apply_diagonal_preconditioner(self, config, x, y) implicit none class(diagonal_preconditioner_type), intent(inout) :: self - class(abstract_vector_type), intent(in) :: x - class(abstract_vector_type), intent(inout) :: y + + type(config_type), intent(in) :: config + class(abstract_vector_type), intent(in) :: x + class(abstract_vector_type), intent(inout) :: y integer(kind=i_def) :: i, nfields ! Workaround for PSyclone to get pointers of the correct type for x and y diff --git a/components/science/source/algorithm/solver/sci_hori_mass_matrix_solver_alg_mod.x90 b/components/science/source/algorithm/solver/sci_hori_mass_matrix_solver_alg_mod.x90 index 366d22313..a5a4f31b8 100644 --- a/components/science/source/algorithm/solver/sci_hori_mass_matrix_solver_alg_mod.x90 +++ b/components/science/source/algorithm/solver/sci_hori_mass_matrix_solver_alg_mod.x90 @@ -12,6 +12,7 @@ module sci_hori_mass_matrix_solver_alg_mod ! Constants and types + use config_mod, only: config_type use constants_mod, only: i_def use integer_field_mod, only: integer_field_type use field_mod, only: field_type @@ -35,6 +36,7 @@ contains ! ======================================================= !> @brief Mass matrix solver for only the horizontal component of a W2 field + !> @param[in] config Application configuration object !> @param[in,out] field_solve The W2 solved field !> @param[in] field The W2 field to apply the horizontal mass !! matrix solve @@ -42,10 +44,12 @@ contains !! not inputted then the vertical solved !! component is set to zero - subroutine hori_mass_matrix_solver_alg( field_solve, field, vert_field_solve ) + subroutine hori_mass_matrix_solver_alg( config, field_solve, field, vert_field_solve ) implicit none + type(config_type), intent(in) :: config + ! Arguments type(field_type), intent(inout) :: field_solve type(field_type), intent(in) :: field @@ -76,7 +80,7 @@ contains ! Take horizontal part of field call split_w2_field_alg( hori_field, vert_field, field ) ! Mass matrix solve of horizontal part - call mass_matrix_solver_alg(hori_field_solve, hori_field) + call mass_matrix_solver_alg(config, hori_field_solve, hori_field) ! Convert to W2 field if (present(vert_field_solve)) then @@ -100,4 +104,4 @@ contains end subroutine hori_mass_matrix_solver_alg -end module sci_hori_mass_matrix_solver_alg_mod \ No newline at end of file +end module sci_hori_mass_matrix_solver_alg_mod diff --git a/components/science/source/algorithm/solver/sci_mass_matrix_operator_alg_mod.x90 b/components/science/source/algorithm/solver/sci_mass_matrix_operator_alg_mod.x90 index 8a56ff1e0..928ea0854 100644 --- a/components/science/source/algorithm/solver/sci_mass_matrix_operator_alg_mod.x90 +++ b/components/science/source/algorithm/solver/sci_mass_matrix_operator_alg_mod.x90 @@ -9,6 +9,7 @@ !> @details Applies the appropriate mass matrix operator to a field. module sci_mass_matrix_operator_alg_mod + use config_mod, only: config_type use constants_mod, only: i_def, r_def, l_def use field_mod, only: field_type use log_mod, only: log_event, & @@ -66,10 +67,10 @@ contains !> @brief Apply mass matrix operator to a pressure field to obtain \f$y=Mx\f$. !> - !> @param[in,out] self Instance of mass_matrix_operator_type - !> @param[in] x Field \f$x\f$ to apply operator to - !> @param[in,out] y Resulting field \f$y=Mx\f$ - subroutine apply_mass_matrix_operator(self, x, y) + !> @param[in] config Application configuration object + !> @param[in] x Field \f$x\f$ to apply operator to + !> @param[in,out] y Resulting field \f$y=Mx\f$ + subroutine apply_mass_matrix_operator(self, config, x, y) use sci_enforce_bc_kernel_mod, only: enforce_bc_kernel_type use matrix_vector_kernel_mod, only: matrix_vector_kernel_type @@ -77,8 +78,10 @@ contains implicit none class(mass_matrix_operator_type), intent(inout) :: self - class(abstract_vector_type), intent(in) :: x - class(abstract_vector_type), intent(inout) :: y + + type(config_type), intent(in) :: config + class(abstract_vector_type), intent(in) :: x + class(abstract_vector_type), intent(inout) :: y ! Workaround for PSyclone to get pointers of the correct type for x and y type(field_type), pointer :: x_vec => null(), & diff --git a/components/science/source/algorithm/solver/sci_mass_matrix_solver_alg_mod.x90 b/components/science/source/algorithm/solver/sci_mass_matrix_solver_alg_mod.x90 index 20c1246b8..1c1ce9a4d 100644 --- a/components/science/source/algorithm/solver/sci_mass_matrix_solver_alg_mod.x90 +++ b/components/science/source/algorithm/solver/sci_mass_matrix_solver_alg_mod.x90 @@ -14,12 +14,14 @@ !> created when needed module sci_mass_matrix_solver_alg_mod + use config_mod, only: config_type use constants_mod, only: i_def, r_def, l_def use log_mod, only: log_event, & log_scratch_space, & LOG_LEVEL_DEBUG, & LOG_LEVEL_ERROR, & LOG_LEVEL_TRACE + use mesh_mod, only: mesh_type ! Derived Types use field_mod, only: field_type @@ -54,11 +56,13 @@ module sci_mass_matrix_solver_alg_mod contains !=============================================================================! !> @details Combined initialisation, run and finalise procedure for the mass matrix solver to solve My = x + !> @param[in] config Application configuration object !> @param[inout] y result field !> @param[in] x input field !> @param[in] bc_flag optional flag to overide the default application of !> boundary conditions - subroutine mass_matrix_solver_alg(y, x, bc_flag) + subroutine mass_matrix_solver_alg(config, y, x, bc_flag) + use solver_config_mod, only: maximum_iterations, & tolerance, & method, & @@ -91,9 +95,12 @@ contains only: operator_tri_solve_kernel_type implicit none + type(config_type), intent(in) :: config + ! Prognostic fields - type( field_type), intent(inout) :: y - type( field_type), intent(in) :: x + type( field_type), intent(inout) :: y + type( field_type), intent(in) :: x + ! Optional flag to overide default boundary condition application logical(kind=l_def), optional, intent(in) :: bc_flag @@ -111,7 +118,9 @@ contains integer(kind=i_def) :: fs logical(kind=l_def) :: apply_bc logical(kind=l_def) :: lowest_order - integer(kind=i_def) :: mesh_id + + type(mesh_type), pointer :: mesh + type(operator_type), pointer :: mass_matrix => null() type(field_type), pointer :: mass_matrix_diagonal => null() type(field_vector_type) :: vec_mm_diagonal @@ -121,7 +130,7 @@ contains if ( LPROF ) call start_timing( id, 'mass_matrix_solver_alg' ) - mesh_id = y%get_mesh_id() + mesh => y%get_mesh() lowest_order = x%get_element_order_h() == 0 .and. & x%get_element_order_v() == 0 @@ -131,14 +140,14 @@ contains ! Direct solve for fully discontinuous spaces !> todo this needs to be extended to all DG spaces if ( lowest_order ) then - mass_matrix => get_inverse_mass_matrix_fv(fs, mesh_id) + mass_matrix => get_inverse_mass_matrix_fv(config, mesh, fs) else - mass_matrix => get_inverse_mass_matrix_fe(fs, mesh_id) + mass_matrix => get_inverse_mass_matrix_fe(config, mesh, fs) end if call invoke( dg_matrix_vector_kernel_type(y, x, mass_matrix) ) else if ( fs == Wtheta .and. lowest_order ) then - mass_matrix => get_mass_matrix_fv(fs, mesh_id) + mass_matrix => get_mass_matrix_fv(config, mesh, fs) call invoke( operator_tri_solve_kernel_type(y, x, mass_matrix) ) else ! Iterative solve for (semi-)continuous spaces @@ -192,12 +201,12 @@ contains if ( lowest_order ) then ! Use mass matrix operator for FV spaces - mass_matrix => get_mass_matrix_fv(fs, mesh_id) - mass_matrix_diagonal => get_mass_matrix_diagonal_fv(fs, mesh_id) + mass_matrix => get_mass_matrix_fv(config, mesh, fs) + mass_matrix_diagonal => get_mass_matrix_diagonal_fv(config, mesh, fs) else ! Use mass matrix operator for higher order spaces - mass_matrix => get_mass_matrix_fe(fs, mesh_id) - mass_matrix_diagonal => get_mass_matrix_diagonal_fe(fs, mesh_id) + mass_matrix => get_mass_matrix_fe(config, mesh, fs) + mass_matrix_diagonal => get_mass_matrix_diagonal_fe(config, mesh, fs) end if select case ( fs ) @@ -307,7 +316,7 @@ contains vec_y = field_vector_type(1) call vec_x%import_field(x, 1) call vec_y%import_field(y, 1) - call mass_matrix_solver%apply(vec_y, vec_x) + call mass_matrix_solver%apply(config, vec_y, vec_x) call vec_y%export_field(y, 1) end if ! direct or iterative solve diff --git a/components/science/source/algorithm/solver/sci_null_preconditioner_alg_mod.x90 b/components/science/source/algorithm/solver/sci_null_preconditioner_alg_mod.x90 index c229fff55..e2e9b6321 100644 --- a/components/science/source/algorithm/solver/sci_null_preconditioner_alg_mod.x90 +++ b/components/science/source/algorithm/solver/sci_null_preconditioner_alg_mod.x90 @@ -11,6 +11,7 @@ module sci_null_preconditioner_alg_mod + use config_mod, only : config_type use constants_mod, only : i_def, r_def use log_mod, only : log_event, & LOG_LEVEL_INFO, & @@ -67,15 +68,18 @@ contains !> !>@details Apply the trivial (null) preconditioner by setting \f$y=x\f$ !> - !>@param[inout] self instance of type null_preconditioner_type - !>@param[inout] x field-vector containing the right hand side of the pressure - !>@param[inout] y field-vector containing the solution - subroutine apply_null_preconditioner(self, x, y) + !>@param[in] config Application configuration object + !>@param[inout] x field-vector containing the right hand side of the pressure + !>@param[inout] y field-vector containing the solution + subroutine apply_null_preconditioner(self, config, x, y) implicit none + class(null_preconditioner_type), intent(inout) :: self - class(abstract_vector_type), intent(in) :: x - class(abstract_vector_type), intent(inout) :: y + + type(config_type), intent(in) :: config + class(abstract_vector_type), intent(in) :: x + class(abstract_vector_type), intent(inout) :: y select type(x) type is(field_vector_type) diff --git a/components/science/source/kernel/fem/sci_compute_broken_div_operator_kernel_mod.F90 b/components/science/source/kernel/fem/sci_compute_broken_div_operator_kernel_mod.F90 index 3d00ff143..f09402ef7 100644 --- a/components/science/source/kernel/fem/sci_compute_broken_div_operator_kernel_mod.F90 +++ b/components/science/source/kernel/fem/sci_compute_broken_div_operator_kernel_mod.F90 @@ -7,21 +7,19 @@ module sci_compute_broken_div_operator_kernel_mod use argument_mod, only: arg_type, func_type, & - GH_OPERATOR, GH_FIELD, & + GH_OPERATOR, & + GH_FIELD, GH_SCALAR, & GH_READ, GH_WRITE, & - GH_REAL, ANY_SPACE_1, & + GH_REAL, GH_INTEGER, & + GH_LOGICAL, ANY_SPACE_1, & ANY_DISCONTINUOUS_SPACE_3, & GH_BASIS, GH_DIFF_BASIS, & CELL_COLUMN, GH_QUADRATURE_XYoZ - use constants_mod, only: r_def, i_def + use constants_mod, only: r_def, i_def, l_def use sci_coordinate_jacobian_mod, only: coordinate_jacobian use fs_continuity_mod, only: W2broken, W3 use kernel_mod, only: kernel_type - use base_mesh_config_mod, only: geometry, topology - use finite_element_config_mod, only: coord_system, rehabilitate - use planet_config_mod, only: scaled_radius - implicit none private @@ -32,10 +30,15 @@ module sci_compute_broken_div_operator_kernel_mod type, public, extends(kernel_type) :: compute_broken_div_operator_kernel_type private - type(arg_type) :: meta_args(3) = (/ & - arg_type(GH_OPERATOR, GH_REAL, GH_WRITE, W3, W2broken), & - arg_type(GH_FIELD*3, GH_REAL, GH_READ, ANY_SPACE_1), & - arg_type(GH_FIELD, GH_REAL, GH_READ, ANY_DISCONTINUOUS_SPACE_3) & + type(arg_type) :: meta_args(8) = (/ & + arg_type(GH_OPERATOR, GH_REAL, GH_WRITE, W3, W2broken), & ! broken_div + arg_type(GH_FIELD*3, GH_REAL, GH_READ, ANY_SPACE_1), & ! chi1, chi2, chi3 + 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 + arg_type(GH_SCALAR, GH_LOGICAL, GH_READ) & ! rehabilitate /) type(func_type) :: meta_funcs(3) = (/ & func_type(W3, GH_BASIS), & @@ -64,6 +67,11 @@ module sci_compute_broken_div_operator_kernel_mod !! @param[in] chi2 2nd coordinate field in Wchi !! @param[in] chi3 3rd coordinate field in Wchi !! @param[in] panel_id Field giving the ID for mesh panels. + !! @param[in] geometry Mesh geometry enumeration + !! @param[in] topology Mesh topology enumeration + !! @param[in] coord_system Finite-element coordinate system enumeration + !! @param[in] scaled_radius Scaled planet radius + !! @param[in] rehabilitate Apply rehabilitation !! @param[in] ndf_w3 Number of degrees of freedom per cell for W3 space. !! @param[in] basis_w3 Scalar basis functions evaluated at quadrature points !! for W3 space. @@ -89,6 +97,9 @@ module sci_compute_broken_div_operator_kernel_mod subroutine compute_broken_div_operator_code(cell, nlayers, ncell_3d, & broken_div, & chi1, chi2, chi3, panel_id, & + geometry, topology, & + coord_system, scaled_radius, & + rehabilitate, & ndf_w3, basis_w3, & ndf_w2b, diff_basis_w2b, & ndf_chi, undf_chi, map_chi, & @@ -122,6 +133,12 @@ subroutine compute_broken_div_operator_code(cell, nlayers, ncell_3d, & real(kind=r_def), dimension(nqp_h), intent(in) :: wqp_h real(kind=r_def), dimension(nqp_v), intent(in) :: wqp_v + 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 + logical(kind=l_def), intent(in) :: rehabilitate + ! Internal variables integer(kind=i_def) :: df, df2, df3, k, ik, ipanel integer(kind=i_def) :: qp1, qp2 diff --git a/components/science/source/kernel/fem/sci_compute_curl_operator_kernel_mod.F90 b/components/science/source/kernel/fem/sci_compute_curl_operator_kernel_mod.F90 index ddcb67d05..b504249fd 100644 --- a/components/science/source/kernel/fem/sci_compute_curl_operator_kernel_mod.F90 +++ b/components/science/source/kernel/fem/sci_compute_curl_operator_kernel_mod.F90 @@ -6,21 +6,20 @@ module sci_compute_curl_operator_kernel_mod use argument_mod, only: arg_type, func_type, & - GH_OPERATOR, GH_FIELD, & + GH_OPERATOR, & + GH_FIELD, GH_SCALAR, & GH_READ, GH_WRITE, & - GH_REAL, ANY_SPACE_1, & + GH_REAL, GH_INTEGER, & + ANY_SPACE_1, & ANY_DISCONTINUOUS_SPACE_3, & GH_BASIS, GH_DIFF_BASIS, & CELL_COLUMN, GH_QUADRATURE_XYoZ + use constants_mod, only: r_def, i_def use sci_coordinate_jacobian_mod, only: coordinate_jacobian use fs_continuity_mod, only: W1, W2 use kernel_mod, only: kernel_type - use base_mesh_config_mod, only: geometry, topology - use finite_element_config_mod, only: coord_system - use planet_config_mod, only: scaled_radius - implicit none private @@ -31,10 +30,14 @@ module sci_compute_curl_operator_kernel_mod type, public, extends(kernel_type) :: compute_curl_operator_kernel_type private - type(arg_type) :: meta_args(3) = (/ & - arg_type(GH_OPERATOR, GH_REAL, GH_WRITE, W2, W1), & - arg_type(GH_FIELD*3, GH_REAL, GH_READ, ANY_SPACE_1), & - arg_type(GH_FIELD, GH_REAL, GH_READ, ANY_DISCONTINUOUS_SPACE_3) & + type(arg_type) :: meta_args(7) = (/ & + arg_type(GH_OPERATOR, GH_REAL, GH_WRITE, W2, W1), & ! curl + arg_type(GH_FIELD*3, GH_REAL, GH_READ, ANY_SPACE_1), & ! chi1, chi2, chi3 + 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(3) = (/ & func_type(W2, GH_BASIS), & @@ -63,6 +66,10 @@ module sci_compute_curl_operator_kernel_mod !! @param[in] chi2 2nd component of coordinate field !! @param[in] chi3 3rd component of coordinate field !! @param[in] panel_id A field giving the ID for mesh panels. +!! @param[in] geometry Mesh geometry enumeration +!! @param[in] topology Mesh topology enumeration +!! @param[in] coord_system Finite-element coordinate system enumeration +!! @param[in] scaled_radius Scaled planet radius !! @param[in] ndf_w2 Number of degrees of freedom per cell for W2. !! @param[in] basis_w2 W2 vector basis functions evaluated at quadrature points. !! @param[in] ndf_w1 Number of degrees of freedom per cell for W1. @@ -83,6 +90,8 @@ module sci_compute_curl_operator_kernel_mod subroutine compute_curl_operator_code(cell, nlayers, ncell_3d, & curl, & chi1, chi2, chi3, panel_id, & + geometry, topology, & + coord_system, scaled_radius, & ndf_w2, basis_w2, & ndf_w1, diff_basis_w1, & ndf_chi, undf_chi, map_chi, & @@ -115,6 +124,11 @@ subroutine compute_curl_operator_code(cell, nlayers, ncell_3d, & real(kind=r_def), dimension(nqp_h), intent(in) :: wqp_h real(kind=r_def), dimension(nqp_v), intent(in) :: wqp_v + 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, df1, df2, k, ik, ipanel integer(kind=i_def) :: qp1, qp2 diff --git a/components/science/source/kernel/fem/sci_compute_derham_matrices_kernel_mod.F90 b/components/science/source/kernel/fem/sci_compute_derham_matrices_kernel_mod.F90 index dcbf2607b..45ada61c4 100644 --- a/components/science/source/kernel/fem/sci_compute_derham_matrices_kernel_mod.F90 +++ b/components/science/source/kernel/fem/sci_compute_derham_matrices_kernel_mod.F90 @@ -16,9 +16,11 @@ module sci_compute_derham_matrices_kernel_mod use argument_mod, only: arg_type, func_type, & - GH_OPERATOR, GH_FIELD, & - GH_READ, GH_WRITE, & - GH_REAL, ANY_SPACE_9, & + GH_OPERATOR, & + GH_FIELD, GH_SCALAR, & + GH_READ, GH_WRITE, & + GH_REAL, GH_INTEGER, & + ANY_SPACE_9, & ANY_DISCONTINUOUS_SPACE_3, & GH_BASIS, GH_DIFF_BASIS, & CELL_COLUMN, GH_QUADRATURE_XYoZ @@ -28,10 +30,6 @@ module sci_compute_derham_matrices_kernel_mod use fs_continuity_mod, only: W0, W1, W2, W2broken, W3, Wtheta use kernel_mod, only: kernel_type - use base_mesh_config_mod, only: geometry, topology - use finite_element_config_mod, only: coord_system - use planet_config_mod, only: scaled_radius - implicit none private @@ -42,19 +40,23 @@ module sci_compute_derham_matrices_kernel_mod type, public, extends(kernel_type) :: compute_derham_matrices_kernel_type private - type(arg_type) :: meta_args(12) = (/ & - arg_type(GH_OPERATOR, GH_REAL, GH_WRITE, W0, W0), & - arg_type(GH_OPERATOR, GH_REAL, GH_WRITE, W1, W1), & - arg_type(GH_OPERATOR, GH_REAL, GH_WRITE, W2, W2), & - arg_type(GH_OPERATOR, GH_REAL, GH_WRITE, W2broken, W2broken), & - arg_type(GH_OPERATOR, GH_REAL, GH_WRITE, W3, W3), & - arg_type(GH_OPERATOR, GH_REAL, GH_WRITE, Wtheta, Wtheta), & - arg_type(GH_OPERATOR, GH_REAL, GH_WRITE, W1, W0), & - arg_type(GH_OPERATOR, GH_REAL, GH_WRITE, W2, W1), & - arg_type(GH_OPERATOR, GH_REAL, GH_WRITE, W3, W2), & - arg_type(GH_OPERATOR, GH_REAL, GH_WRITE, W3, W2broken), & - 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(16) = (/ & + arg_type(GH_OPERATOR, GH_REAL, GH_WRITE, W0, W0), & ! mm0 + arg_type(GH_OPERATOR, GH_REAL, GH_WRITE, W1, W1), & ! mm1 + arg_type(GH_OPERATOR, GH_REAL, GH_WRITE, W2, W2), & ! mm2 + arg_type(GH_OPERATOR, GH_REAL, GH_WRITE, W2broken, W2broken), & ! mm2b + arg_type(GH_OPERATOR, GH_REAL, GH_WRITE, W3, W3), & ! mm3 + arg_type(GH_OPERATOR, GH_REAL, GH_WRITE, Wtheta, Wtheta), & ! mmt + arg_type(GH_OPERATOR, GH_REAL, GH_WRITE, W1, W0), & ! grad + arg_type(GH_OPERATOR, GH_REAL, GH_WRITE, W2, W1), & ! curl + arg_type(GH_OPERATOR, GH_REAL, GH_WRITE, W3, W2), & ! div + arg_type(GH_OPERATOR, GH_REAL, GH_WRITE, W3, W2broken), & ! broken_div + arg_type(GH_FIELD*3, GH_REAL, GH_READ, ANY_SPACE_9), & ! chi1, chi2, chi3 + 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(7) = (/ & func_type(W0, GH_BASIS, GH_DIFF_BASIS), & @@ -105,6 +107,11 @@ module sci_compute_derham_matrices_kernel_mod !! @param[in] chi2 Physical coordinates in the 2nd dir. !! @param[in] chi3 Physical coordinates in the 3rd dir. !! @param[in] panel_id Field giving the ID for mesh panels. +!! @param[in] geometry Mesh geometry enumeration +!! @param[in] topology Mesh topology enumeration +!! @param[in] coord_system Finite-element coordinate system enumeration +!! @param[in] scaled_radius Scaled planet radius +!! @param[in] rehabilitate Apply rehabilitation !! @param[in] ndf_w0 Number of degrees of freedom per cell for W0 space. !! @param[in] basis_w0 Basis functions evaluated at quadrature points for W0 space. !! @param[in] diff_basis_w0 Differential of basis functions evaluated at quadrature points for W0 space. @@ -148,6 +155,8 @@ subroutine compute_derham_matrices_code(cell, nlayers, & ncell_3d7, broken_div, & chi1, chi2, chi3, & panel_id, & + geometry, topology, & + coord_system, scaled_radius, & ndf_w0, basis_w0, diff_basis_w0, & ndf_w1, basis_w1, diff_basis_w1, & ndf_w2, basis_w2, diff_basis_w2, & @@ -205,6 +214,11 @@ subroutine compute_derham_matrices_code(cell, nlayers, & real(kind=r_def), intent(in) :: wqp_h(nqp_h) real(kind=r_def), intent(in) :: wqp_v(nqp_v) + 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, df2, k, ik integer(kind=i_def) :: qp1, qp2 diff --git a/components/science/source/kernel/fem/sci_compute_div_operator_kernel_mod.F90 b/components/science/source/kernel/fem/sci_compute_div_operator_kernel_mod.F90 index 05b84b344..93eec65a7 100644 --- a/components/science/source/kernel/fem/sci_compute_div_operator_kernel_mod.F90 +++ b/components/science/source/kernel/fem/sci_compute_div_operator_kernel_mod.F90 @@ -6,21 +6,19 @@ module sci_compute_div_operator_kernel_mod use argument_mod, only: arg_type, func_type, & - GH_OPERATOR, GH_FIELD, & + GH_OPERATOR, & + GH_FIELD, GH_SCALAR, & GH_READ, GH_WRITE, & - GH_REAL, ANY_SPACE_1, & + GH_REAL, GH_INTEGER, & + GH_LOGICAL, ANY_SPACE_1, & GH_BASIS, GH_DIFF_BASIS, & ANY_DISCONTINUOUS_SPACE_3, & GH_QUADRATURE_XYoZ, CELL_COLUMN - use constants_mod, only: r_def, i_def + use constants_mod, only: r_def, i_def, l_def use sci_coordinate_jacobian_mod, only: coordinate_jacobian use fs_continuity_mod, only: W2, W3 use kernel_mod, only: kernel_type - use base_mesh_config_mod, only: geometry, topology - use finite_element_config_mod, only: coord_system, rehabilitate - use planet_config_mod, only: scaled_radius - implicit none private @@ -31,10 +29,16 @@ module sci_compute_div_operator_kernel_mod type, public, extends(kernel_type) :: compute_div_operator_kernel_type private - type(arg_type) :: meta_args(3) = (/ & - arg_type(GH_OPERATOR, GH_REAL, GH_WRITE, W3, W2), & - arg_type(GH_FIELD*3, GH_REAL, GH_READ, ANY_SPACE_1), & - arg_type(GH_FIELD, GH_REAL, GH_READ, ANY_DISCONTINUOUS_SPACE_3) & + type(arg_type) :: meta_args(8) = (/ & + arg_type(GH_OPERATOR, GH_REAL, GH_WRITE, W3, W2), & ! div + arg_type(GH_FIELD*3, GH_REAL, GH_READ, ANY_SPACE_1), & ! chi1, chi2, chi3 + 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 + arg_type(GH_SCALAR, GH_LOGICAL, GH_READ) & ! rehabilitate + /) type(func_type) :: meta_funcs(3) = (/ & func_type(W3, GH_BASIS), & @@ -63,6 +67,11 @@ module sci_compute_div_operator_kernel_mod !! @param[in] chi_2 2nd coordinate field in Wchi !! @param[in] chi_3 3rd coordinate field in Wchi !! @param[in] panel_id Field giving the ID for mesh panels. +!! @param[in] geometry Mesh geometry enumeration +!! @param[in] topology Mesh topology enumeration +!! @param[in] coord_system Finite-element coordinate system enumeration +!! @param[in] scaled_radius Scaled planet radius +!! @param[in] rehabilitate Apply rehabilitation !! @param[in] ndf_w3 Number of degrees of freedom per cell. !! @param[in] basis_w3 Scalar basis functions !! evaluated at quadrature points. @@ -87,6 +96,9 @@ subroutine compute_div_operator_code(cell, nlayers, ncell_3d, & div, & chi1, chi2, chi3, & panel_id, & + geometry, topology, & + coord_system, scaled_radius, & + rehabilitate, & ndf_w3, & basis_w3, & ndf_w2, diff_basis_w2, & @@ -124,6 +136,12 @@ subroutine compute_div_operator_code(cell, nlayers, ncell_3d, & real(kind=r_def), dimension(nqp_h), intent(in) :: wqp_h real(kind=r_def), dimension(nqp_v), intent(in) :: wqp_v + 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 + logical(kind=l_def), intent(in) :: rehabilitate + ! Internal variables integer(kind=i_def) :: df, df2, df3, k, ik integer(kind=i_def) :: qp1, qp2, ipanel diff --git a/components/science/source/kernel/fem/sci_compute_grad_operator_kernel_mod.F90 b/components/science/source/kernel/fem/sci_compute_grad_operator_kernel_mod.F90 index 9c06d929f..2006ebc4a 100644 --- a/components/science/source/kernel/fem/sci_compute_grad_operator_kernel_mod.F90 +++ b/components/science/source/kernel/fem/sci_compute_grad_operator_kernel_mod.F90 @@ -6,9 +6,11 @@ module sci_compute_grad_operator_kernel_mod use argument_mod, only: arg_type, func_type, & - GH_OPERATOR, GH_FIELD, & + GH_OPERATOR, & + GH_FIELD, GH_SCALAR, & GH_READ, GH_WRITE, & - GH_REAL, ANY_SPACE_1, & + GH_REAL, GH_INTEGER, & + ANY_SPACE_1, & ANY_DISCONTINUOUS_SPACE_3, & GH_BASIS, GH_DIFF_BASIS, & CELL_COLUMN, GH_QUADRATURE_XYoZ @@ -19,10 +21,6 @@ module sci_compute_grad_operator_kernel_mod use fs_continuity_mod, only: W0, W1 use kernel_mod, only: kernel_type - use base_mesh_config_mod, only: geometry, topology - use finite_element_config_mod, only: coord_system - use planet_config_mod, only: scaled_radius - implicit none private @@ -33,10 +31,14 @@ module sci_compute_grad_operator_kernel_mod type, public, extends(kernel_type) :: compute_grad_operator_kernel_type private - type(arg_type) :: meta_args(3) = (/ & - arg_type(GH_OPERATOR, GH_REAL, GH_WRITE, W1, W0), & - arg_type(GH_FIELD*3, GH_REAL, GH_READ, ANY_SPACE_1), & - arg_type(GH_FIELD, GH_REAL, GH_READ, ANY_DISCONTINUOUS_SPACE_3) & + type(arg_type) :: meta_args(7) = (/ & + arg_type(GH_OPERATOR, GH_REAL, GH_WRITE, W1, W0), & ! grad + arg_type(GH_FIELD*3, GH_REAL, GH_READ, ANY_SPACE_1), & ! chi1, chi2, chi3 + 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(3) = (/ & func_type(W1, GH_BASIS), & @@ -65,6 +67,10 @@ module sci_compute_grad_operator_kernel_mod !! @param[in] chi2 2nd coordinate field in Wchi !! @param[in] chi3 3rd coordinate field in Wchi !! @param[in] panel_id Field giving the ID for mesh panels +!! @param[in] geometry Mesh geometry enumeration +!! @param[in] topology Mesh topology enumeration +!! @param[in] coord_system Finite-element coordinate system enumeration +!! @param[in] scaled_radius Scaled planet radius !! @param[in] ndf_w1 Number of degrees of freedom per cell !! @param[in] basis_w1 Vector basis functions !! evaluated at quadrature points @@ -85,14 +91,16 @@ module sci_compute_grad_operator_kernel_mod !! @param[in] nqp_v Number of vertical quadrature points !! @param[in] wqp_h Horizontal quadrature weights !! @param[in] wqp_v Vertical quadrature weights -subroutine compute_grad_operator_code(cell, nlayers, ncell_3d, & - grad, & - chi1, chi2, chi3, panel_id, & - ndf_w1, basis_w1, & - ndf_w0, diff_basis_w0, & - ndf_chi, undf_chi, map_chi, & - basis_chi, diff_basis_chi, & - ndf_pid, undf_pid, map_pid, & +subroutine compute_grad_operator_code(cell, nlayers, ncell_3d, & + grad, & + chi1, chi2, chi3, panel_id, & + geometry, topology, & + coord_system, scaled_radius, & + ndf_w1, basis_w1, & + ndf_w0, diff_basis_w0, & + ndf_chi, undf_chi, map_chi, & + basis_chi, diff_basis_chi, & + ndf_pid, undf_pid, map_pid, & nqp_h, nqp_v, wqp_h, wqp_v ) implicit none @@ -120,6 +128,11 @@ subroutine compute_grad_operator_code(cell, nlayers, ncell_3d, & real(kind=r_def), dimension(nqp_h), intent(in) :: wqp_h real(kind=r_def), dimension(nqp_v), intent(in) :: wqp_v + 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, df0, df1, k, ik integer(kind=i_def) :: qp1, qp2, ipanel diff --git a/components/science/source/kernel/fem/sci_compute_mass_matrix_kernel_w1_mod.F90 b/components/science/source/kernel/fem/sci_compute_mass_matrix_kernel_w1_mod.F90 index b0beabbbb..045c228fd 100644 --- a/components/science/source/kernel/fem/sci_compute_mass_matrix_kernel_w1_mod.F90 +++ b/components/science/source/kernel/fem/sci_compute_mass_matrix_kernel_w1_mod.F90 @@ -11,22 +11,21 @@ module sci_compute_mass_matrix_kernel_w1_mod use argument_mod, only: arg_type, func_type, & - GH_OPERATOR, GH_FIELD, & + GH_OPERATOR, & + GH_FIELD, GH_SCALAR, & GH_READ, GH_WRITE, & - GH_REAL, ANY_SPACE_9, & + GH_REAL, GH_INTEGER, & + ANY_SPACE_9, & ANY_DISCONTINUOUS_SPACE_3, & GH_BASIS, GH_DIFF_BASIS, & CELL_COLUMN, GH_QUADRATURE_XYoZ + use sci_coordinate_jacobian_mod, only: coordinate_jacobian, & coordinate_jacobian_inverse use constants_mod, only: r_def, i_def use fs_continuity_mod, only: W1 use kernel_mod, only: kernel_type - use base_mesh_config_mod, only: geometry, topology - use finite_element_config_mod, only: coord_system - use planet_config_mod, only: scaled_radius - implicit none private @@ -37,10 +36,14 @@ module sci_compute_mass_matrix_kernel_w1_mod type, public, extends(kernel_type) :: compute_mass_matrix_kernel_w1_type private - type(arg_type) :: meta_args(3) = (/ & - arg_type(GH_OPERATOR, GH_REAL, GH_WRITE, W1, W1), & - 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_OPERATOR, GH_REAL, GH_WRITE, W1, W1), & ! mm + arg_type(GH_FIELD*3, GH_REAL, GH_READ, ANY_SPACE_9), & ! chi1, chi2, chi3 + 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(ANY_SPACE_9, GH_BASIS, GH_DIFF_BASIS), & @@ -68,6 +71,10 @@ module sci_compute_mass_matrix_kernel_w1_mod !! @param[in] chi2 2nd coordinate field in Wchi !! @param[in] chi3 3rd coordinate field in Wchi !! @param[in] panel_id Field giving the ID for mesh panels +!! @param[in] geometry Mesh geometry enumeration +!! @param[in] topology Mesh topology enumeration +!! @param[in] coord_system Finite-element coordinate system enumeration +!! @param[in] scaled_radius Scaled planet radius !! @param[in] ndf_w1 Number of degrees of freedom per cell !! @param[in] basis_w1 Vector basis functions evaluated at quadrature points !! @param[in] ndf_chi Number of degrees of freedom per cell for chi field @@ -88,6 +95,8 @@ module sci_compute_mass_matrix_kernel_w1_mod subroutine compute_mass_matrix_w1_code(cell, nlayers, ncell_3d, & mm, & chi1, chi2, chi3, panel_id, & + geometry, topology, & + coord_system, scaled_radius, & ndf_w1, basis_w1, & ndf_chi, undf_chi, map_chi, & basis_chi, diff_basis_chi, & @@ -116,6 +125,11 @@ subroutine compute_mass_matrix_w1_code(cell, nlayers, ncell_3d, & real(kind=r_def), intent(in) :: wqp_h(nqp_h) real(kind=r_def), intent(in) :: wqp_v(nqp_v) + 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, df2, k, ik, ipanel integer(kind=i_def) :: qp1, qp2 diff --git a/components/science/source/kernel/fem/sci_compute_mass_matrix_kernel_w2_mod.F90 b/components/science/source/kernel/fem/sci_compute_mass_matrix_kernel_w2_mod.F90 index f1ebbcaae..a97369886 100644 --- a/components/science/source/kernel/fem/sci_compute_mass_matrix_kernel_w2_mod.F90 +++ b/components/science/source/kernel/fem/sci_compute_mass_matrix_kernel_w2_mod.F90 @@ -10,22 +10,19 @@ !> module sci_compute_mass_matrix_kernel_w2_mod - use argument_mod, only: arg_type, func_type, & - GH_OPERATOR, GH_FIELD, & - GH_READ, GH_WRITE, & - GH_REAL, ANY_W2, & - ANY_DISCONTINUOUS_SPACE_3, & - ANY_SPACE_9, & - GH_BASIS, GH_DIFF_BASIS, & + use argument_mod, only: arg_type, func_type, & + GH_OPERATOR, & + GH_FIELD, GH_SCALAR, & + GH_READ, GH_WRITE, & + GH_REAL, GH_INTEGER, ANY_W2, & + ANY_DISCONTINUOUS_SPACE_3, & + ANY_SPACE_9, & + GH_BASIS, GH_DIFF_BASIS, & CELL_COLUMN, GH_QUADRATURE_XYoZ use constants_mod, only: i_def, r_def use sci_coordinate_jacobian_mod, only: coordinate_jacobian use kernel_mod, only: kernel_type - use base_mesh_config_mod, only: geometry, topology - use finite_element_config_mod, only: coord_system - use planet_config_mod, only: scaled_radius - implicit none private @@ -36,10 +33,14 @@ module sci_compute_mass_matrix_kernel_w2_mod type, public, extends(kernel_type) :: compute_mass_matrix_kernel_w2_type private - type(arg_type) :: meta_args(3) = (/ & - arg_type(GH_OPERATOR, GH_REAL, GH_WRITE, ANY_W2, ANY_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) & + type(arg_type) :: meta_args(7) = (/ & + arg_type(GH_OPERATOR, GH_REAL, GH_WRITE, ANY_W2, ANY_W2), & ! mm + arg_type(GH_FIELD*3, GH_REAL, GH_READ, ANY_SPACE_9), & ! chi1, chi2, chi3 + 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(ANY_W2, GH_BASIS), & @@ -68,6 +69,10 @@ module sci_compute_mass_matrix_kernel_w2_mod !! @param[in] chi2 2nd coordinate field in Wchi !! @param[in] chi3 3rd coordinate field in Wchi !! @param[in] panel_id Field giving the ID for mesh panels +!! @param[in] geometry Mesh geometry enumeration +!! @param[in] topology Mesh topology enumeration +!! @param[in] coord_system Finite-element coordinate system enumeration +!! @param[in] scaled_radius Scaled planet radius !! @param[in] ndf_w2 Degrees of freedom per cell !! @param[in] basis_w2 Vector basis functions evaluated at quadrature points !! @param[in] ndf_chi Degrees of freedom per cell for chi field @@ -88,6 +93,8 @@ subroutine compute_mass_matrix_w2_code(cell, nlayers, ncell_3d, & mm, & chi1, chi2, chi3, & panel_id, & + geometry, topology, & + coord_system, scaled_radius, & ndf_w2, basis_w2, & ndf_chi, undf_chi, map_chi, & basis_chi, & @@ -120,6 +127,11 @@ subroutine compute_mass_matrix_w2_code(cell, nlayers, ncell_3d, & real(kind=r_def), intent(in) :: wqp_h(nqp_h) real(kind=r_def), intent(in) :: wqp_v(nqp_v) + 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, df2, k, ik, ipanel integer(kind=i_def) :: qp1, qp2 diff --git a/components/science/source/kernel/fem/sci_compute_mass_matrix_kernel_w3_mod.F90 b/components/science/source/kernel/fem/sci_compute_mass_matrix_kernel_w3_mod.F90 index 5423fbe34..cf48acde9 100644 --- a/components/science/source/kernel/fem/sci_compute_mass_matrix_kernel_w3_mod.F90 +++ b/components/science/source/kernel/fem/sci_compute_mass_matrix_kernel_w3_mod.F90 @@ -11,21 +11,20 @@ module sci_compute_mass_matrix_kernel_w3_mod use, intrinsic :: iso_fortran_env, only: real32, real64 - use argument_mod, only: arg_type, func_type, & - GH_OPERATOR, GH_FIELD, & - GH_READ, GH_WRITE, & - GH_REAL, ANY_SPACE_9, & - ANY_DISCONTINUOUS_SPACE_3, & - GH_BASIS, GH_DIFF_BASIS, & - CELL_COLUMN, GH_QUADRATURE_XYoZ - use sci_coordinate_jacobian_mod, only: coordinate_jacobian - use constants_mod, only: i_def - use fs_continuity_mod, only: W3 - use kernel_mod, only: kernel_type + use argument_mod, only: arg_type, func_type, & + GH_OPERATOR, GH_FIELD, GH_SCALAR, & + GH_READ, GH_WRITE, & + GH_REAL, GH_INTEGER, GH_LOGICAL, & + ANY_SPACE_9, & + ANY_DISCONTINUOUS_SPACE_3, & + GH_BASIS, GH_DIFF_BASIS, & + CELL_COLUMN, GH_QUADRATURE_XYoZ + + use constants_mod, only: i_def, r_def, l_def + use fs_continuity_mod, only: W3 + use kernel_mod, only: kernel_type - use base_mesh_config_mod, only: geometry, topology - use finite_element_config_mod, only: coord_system, rehabilitate - use planet_config_mod, only: scaled_radius + use sci_coordinate_jacobian_mod, only: coordinate_jacobian implicit none @@ -36,10 +35,15 @@ module sci_compute_mass_matrix_kernel_w3_mod !--------------------------------------------------------------------------- type, public, extends(kernel_type) :: compute_mass_matrix_kernel_w3_type private - type(arg_type) :: meta_args(3) = (/ & - arg_type(GH_OPERATOR, GH_REAL, GH_WRITE, W3, W3), & - 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_OPERATOR, GH_REAL, GH_WRITE, W3, W3), & ! mm + arg_type(GH_FIELD*3, GH_REAL, GH_READ, ANY_SPACE_9), & ! chi1, chi2, chi3 + 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 + arg_type(GH_SCALAR, GH_LOGICAL, GH_READ) & ! rehabilitate /) type(func_type) :: meta_funcs(2) = (/ & func_type(W3, GH_BASIS), & @@ -72,6 +76,11 @@ module sci_compute_mass_matrix_kernel_w3_mod !! @param[in] chi2 2nd coordinate field in Wchi !! @param[in] chi3 3rd coordinate field in Wchi !! @param[in] panel_id Field giving the ID for mesh panels + !! @param[in] geometry Mesh geometry enumeration + !! @param[in] topology Mesh topology enumeration + !! @param[in] coord_system Finite-element coordinate system enumeration + !! @param[in] scaled_radius Scaled planet radius + !! @param[in] rehabilitate Apply rehabilitation !! @param[in] ndf_w3 Number of degrees of freedom per cell for the operator space !! @param[in] basis_w3 Scalar basis functions evaluated at quadrature points !! @param[in] ndf_chi Number of degrees of freedom per cell for the coordinate field @@ -90,9 +99,12 @@ module sci_compute_mass_matrix_kernel_w3_mod ! REAL32 PRECISION ! ================== - subroutine compute_mass_matrix_w3_code_real32( & + subroutine compute_mass_matrix_w3_code_real32( & cell, nlayers, ncell_3d, mm, & chi1, chi2, chi3, panel_id, & + geometry, topology, & + coord_system, scaled_radius, & + rehabilitate, & ndf_w3, basis_w3, & ndf_chi, undf_chi, map_chi, & basis_chi, diff_basis_chi, & @@ -110,17 +122,22 @@ subroutine compute_mass_matrix_w3_code_real32( & integer(kind=i_def), dimension(ndf_chi), intent(in) :: map_chi integer(kind=i_def), dimension(ndf_pid), intent(in) :: map_pid - real(kind=real32), dimension(ncell_3d,ndf_w3,ndf_w3), intent(inout) :: mm + real(kind=real32), dimension(ncell_3d,ndf_w3,ndf_w3), intent(inout) :: mm + real(kind=real32), dimension(1,ndf_chi,nqp_h,nqp_v), intent(in) :: basis_chi + real(kind=real32), dimension(3,ndf_chi,nqp_h,nqp_v), intent(in) :: diff_basis_chi + real(kind=real32), dimension(1,ndf_w3,nqp_h,nqp_v), intent(in) :: basis_w3 - real(kind=real32), dimension(1,ndf_chi,nqp_h,nqp_v), intent(in) :: basis_chi - real(kind=real32), dimension(3,ndf_chi,nqp_h,nqp_v), intent(in) :: diff_basis_chi - real(kind=real32), dimension(1,ndf_w3,nqp_h,nqp_v), intent(in) :: basis_w3 + real(kind=real32), dimension(undf_chi), intent(in) :: chi1, chi2, chi3 + real(kind=real32), dimension(undf_pid), intent(in) :: panel_id - real(kind=real32), dimension(undf_chi), intent(in) :: chi1, chi2, chi3 - real(kind=real32), dimension(undf_pid), intent(in) :: panel_id + real(kind=real32), dimension(nqp_h), intent(in) :: wqp_h + real(kind=real32), dimension(nqp_v), intent(in) :: wqp_v - real(kind=real32), dimension(nqp_h), intent(in) :: wqp_h - real(kind=real32), dimension(nqp_v), intent(in) :: wqp_v + 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 + logical(kind=l_def), intent(in) :: rehabilitate !Internal variables integer(kind=i_def) :: df, df2, k, ik, ipanel @@ -186,6 +203,9 @@ end subroutine compute_mass_matrix_w3_code_real32 subroutine compute_mass_matrix_w3_code_mixed_precision( & cell, nlayers, ncell_3d, mm, & chi1, chi2, chi3, panel_id, & + geometry, topology, & + coord_system, scaled_radius, & + rehabilitate, & ndf_w3, basis_w3, & ndf_chi, undf_chi, map_chi, & basis_chi, diff_basis_chi, & @@ -215,6 +235,12 @@ subroutine compute_mass_matrix_w3_code_mixed_precision( & real(kind=real64), dimension(nqp_h), intent(in) :: wqp_h real(kind=real64), dimension(nqp_v), intent(in) :: wqp_v + 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 + logical, intent(in) :: rehabilitate + !Internal variables integer(kind=i_def) :: df, df2, k, ik, ipanel integer(kind=i_def) :: qp1, qp2 @@ -279,6 +305,9 @@ end subroutine compute_mass_matrix_w3_code_mixed_precision subroutine compute_mass_matrix_w3_code_real64( & cell, nlayers, ncell_3d, mm, & chi1, chi2, chi3, panel_id, & + geometry, topology, & + coord_system, scaled_radius, & + rehabilitate, & ndf_w3, basis_w3, & ndf_chi, undf_chi, map_chi, & basis_chi, diff_basis_chi, & @@ -308,6 +337,12 @@ subroutine compute_mass_matrix_w3_code_real64( & real(kind=real64), dimension(nqp_h), intent(in) :: wqp_h real(kind=real64), dimension(nqp_v), intent(in) :: wqp_v + 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 + logical, intent(in) :: rehabilitate + !Internal variables integer(kind=i_def) :: df, df2, k, ik, ipanel integer(kind=i_def) :: qp1, qp2 diff --git a/components/science/source/kernel/fem/sci_compute_mass_matrix_kernel_w_scalar_mod.F90 b/components/science/source/kernel/fem/sci_compute_mass_matrix_kernel_w_scalar_mod.F90 index b54a6ba43..c9876d4f1 100644 --- a/components/science/source/kernel/fem/sci_compute_mass_matrix_kernel_w_scalar_mod.F90 +++ b/components/science/source/kernel/fem/sci_compute_mass_matrix_kernel_w_scalar_mod.F90 @@ -17,19 +17,17 @@ module sci_compute_mass_matrix_kernel_w_scalar_mod GH_OPERATOR, GH_FIELD, & GH_LOGICAL, GH_SCALAR, & GH_READ, GH_WRITE, & - GH_REAL, ANY_SPACE_2, & - ANY_SPACE_9, & + GH_REAL, GH_INTEGER, & + ANY_SPACE_2, ANY_SPACE_9, & ANY_DISCONTINUOUS_SPACE_3, & GH_BASIS, GH_DIFF_BASIS, & CELL_COLUMN, GH_QUADRATURE_XYoZ - use constants_mod, only: i_def, l_def - use sci_coordinate_jacobian_mod, only: coordinate_jacobian + + use constants_mod, only: i_def, r_def, l_def use fs_continuity_mod, only: W0, Wtheta use kernel_mod, only: kernel_type - use base_mesh_config_mod, only: geometry, topology - use finite_element_config_mod, only: coord_system - use planet_config_mod, only: scaled_radius + use sci_coordinate_jacobian_mod, only: coordinate_jacobian implicit none @@ -40,11 +38,15 @@ module sci_compute_mass_matrix_kernel_w_scalar_mod !--------------------------------------------------------------------------- type, public, extends(kernel_type) :: compute_mass_matrix_kernel_w_scalar_type private - type(arg_type) :: meta_args(4) = (/ & - arg_type(GH_OPERATOR, GH_REAL, GH_WRITE, ANY_SPACE_2, ANY_SPACE_2), & - 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_LOGICAL, GH_READ) & + type(arg_type) :: meta_args(8) = (/ & + arg_type(GH_OPERATOR, GH_REAL, GH_WRITE, ANY_SPACE_2, ANY_SPACE_2), & ! mm + arg_type(GH_FIELD*3, GH_REAL, GH_READ, ANY_SPACE_9), & ! chi1, chi2, chi3 + arg_type(GH_FIELD, GH_REAL, GH_READ, ANY_DISCONTINUOUS_SPACE_3), & ! panel_id + arg_type(GH_SCALAR, GH_LOGICAL, GH_READ), & ! extended_mesh + 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(ANY_SPACE_2, GH_BASIS), & @@ -78,6 +80,10 @@ module sci_compute_mass_matrix_kernel_w_scalar_mod !! @param[in] chi3 3rd coordinate field in Wchi !! @param[in] panel_id Field giving the ID for mesh panels. !! @param[in] extended_mesh Compute on an extended mesh + !! @param[in] geometry Mesh geometry enumeration + !! @param[in] topology Mesh topology enumeration + !! @param[in] coord_system Finite-element coordinate system enumeration + !! @param[in] scaled_radius Scaled planet radius !! @param[in] ndf_w_scalar The number of degrees of freedom per cell for w_scalar. !! @param[in] ndf_chi The number of degrees of freedom per cell for chi. !! @param[in] basis_w_scalar 4-dim array holding SCALAR basis functions evaluated at @@ -103,6 +109,8 @@ subroutine compute_mass_matrix_w_scalar_code_r32( & cell, nlayers, ncell_3d, mm, & chi1, chi2, chi3, panel_id, & extended_mesh, & + geometry, topology, & + coord_system, scaled_radius, & ndf_w_scalar, basis_w_scalar, & ndf_chi, undf_chi, map_chi, & basis_chi, diff_basis_chi, & @@ -123,6 +131,11 @@ subroutine compute_mass_matrix_w_scalar_code_r32( & logical(kind=l_def), intent(in) :: extended_mesh + 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=real32), dimension(ncell_3d,ndf_w_scalar,ndf_w_scalar), & intent(inout) :: mm @@ -203,6 +216,8 @@ subroutine compute_mass_matrix_w_scalar_code_r32r64( & cell, nlayers, ncell_3d, mm, & chi1, chi2, chi3, panel_id, & extended_mesh, & + geometry, topology, & + coord_system, scaled_radius, & ndf_w_scalar, basis_w_scalar, & ndf_chi, undf_chi, map_chi, & basis_chi, diff_basis_chi, & @@ -223,6 +238,11 @@ subroutine compute_mass_matrix_w_scalar_code_r32r64( & logical(kind=l_def), intent(in) :: extended_mesh + 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=real32), dimension(ncell_3d,ndf_w_scalar,ndf_w_scalar), & intent(inout) :: mm @@ -302,6 +322,8 @@ subroutine compute_mass_matrix_w_scalar_code_r64( & cell, nlayers, ncell_3d, mm, & chi1, chi2, chi3, panel_id, & extended_mesh, & + geometry, topology, & + coord_system, scaled_radius, & ndf_w_scalar, basis_w_scalar, & ndf_chi, undf_chi, map_chi, & basis_chi, diff_basis_chi, & @@ -322,6 +344,11 @@ subroutine compute_mass_matrix_w_scalar_code_r64( & logical(kind=l_def), intent(in) :: extended_mesh + 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=real64), dimension(ncell_3d,ndf_w_scalar,ndf_w_scalar), & intent(inout) :: mm diff --git a/components/science/source/kernel/fem/sci_gp_rhs_kernel_mod.F90 b/components/science/source/kernel/fem/sci_gp_rhs_kernel_mod.F90 index ca6d73d1e..9323af616 100644 --- a/components/science/source/kernel/fem/sci_gp_rhs_kernel_mod.F90 +++ b/components/science/source/kernel/fem/sci_gp_rhs_kernel_mod.F90 @@ -7,14 +7,16 @@ !> @brief Kernel which projects a field into into a given space module sci_gp_rhs_kernel_mod + use kernel_mod, only : kernel_type use constants_mod, only : r_def, i_def -use argument_mod, only : arg_type, func_type, & - GH_FIELD, GH_REAL, GH_INC, & - GH_READ, ANY_SPACE_9, & - ANY_SPACE_1, ANY_SPACE_2, & - ANY_DISCONTINUOUS_SPACE_3, & - GH_BASIS, GH_DIFF_BASIS, & +use argument_mod, only : arg_type, func_type, & + GH_FIELD, GH_SCALAR, & + GH_REAL, GH_INTEGER, & + GH_INC, GH_READ, ANY_SPACE_9, & + ANY_SPACE_1, ANY_SPACE_2, & + ANY_DISCONTINUOUS_SPACE_3, & + GH_BASIS, GH_DIFF_BASIS, & CELL_COLUMN, GH_QUADRATURE_XYoZ implicit none @@ -27,11 +29,15 @@ module sci_gp_rhs_kernel_mod !> The type declaration for the kernel. Contains the metadata needed by the Psy layer type, public, extends(kernel_type) :: gp_rhs_kernel_type private - type(arg_type) :: meta_args(4) = (/ & - arg_type(GH_FIELD, GH_REAL, GH_INC, ANY_SPACE_1), & - arg_type(GH_FIELD, GH_REAL, GH_READ, ANY_SPACE_2), & - 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_INC, ANY_SPACE_1), & ! rhs + arg_type(GH_FIELD, GH_REAL, GH_READ, ANY_SPACE_2), & ! field + arg_type(GH_FIELD*3, GH_REAL, GH_READ, ANY_SPACE_9), & ! chi1, chi2, chi3 + 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(3) = (/ & func_type(ANY_SPACE_1, GH_BASIS), & @@ -63,6 +69,10 @@ module sci_gp_rhs_kernel_mod !! @param[in] chi_2 2nd coordinate field in Wchi !! @param[in] chi_3 3rd coordinate field in Wchi !! @param[in] panel_id Field giving the ID for mesh panels +!! @param[in] geometry Mesh geometry enumeration +!! @param[in] topology Mesh topology enumeration +!! @param[in] coord_system Finite-element coordinate system enumeration +!! @param[in] scaled_radius Scaled planet radius !! @param[in] ndf Number of degrees of freedom per cell !! @param[in] undf Number of (local) unique degrees of freedom of the field rhs !! @param[in] map Dofmap for the cell at the base of the column @@ -87,6 +97,8 @@ module sci_gp_rhs_kernel_mod subroutine gp_rhs_code(nlayers, & rhs, field, & chi_1, chi_2, chi_3, panel_id, & + geometry, topology, & + coord_system, scaled_radius, & ndf, undf, map, basis, & ndf_f, undf_f, map_f, f_basis, & ndf_chi, undf_chi, map_chi, & @@ -96,10 +108,6 @@ subroutine gp_rhs_code(nlayers, & use sci_coordinate_jacobian_mod, only: coordinate_jacobian - use base_mesh_config_mod, only: geometry, topology - use finite_element_config_mod, only: coord_system - use planet_config_mod, only: scaled_radius - implicit none ! Arguments @@ -124,6 +132,11 @@ subroutine gp_rhs_code(nlayers, & real(kind=r_def), dimension(nqp_h), intent(in) :: wqp_h real(kind=r_def), dimension(nqp_v), intent(in) :: wqp_v + 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, df2, k, qp1, qp2, ipanel real(kind=r_def), dimension(nqp_h,nqp_v) :: dj 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 38a827c21..9d872b33e 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 @@ -8,8 +8,10 @@ module sci_gp_vector_rhs_kernel_mod use argument_mod, only : arg_type, func_type, & - GH_FIELD, GH_REAL, GH_INC, & - GH_READ, ANY_SPACE_1, & + GH_FIELD, GH_SCALAR, & + GH_REAL, GH_INTEGER, & + GH_INC, GH_READ, & + ANY_SPACE_1, & ANY_SPACE_2, ANY_SPACE_9, & ANY_DISCONTINUOUS_SPACE_3, & GH_BASIS, GH_DIFF_BASIS, & @@ -21,11 +23,7 @@ module sci_gp_vector_rhs_kernel_mod use coord_transform_mod, only : cart2sphere_vector use fs_continuity_mod, only : W0, W2 use kernel_mod, only : kernel_type - - use base_mesh_config_mod, only: geometry, topology, & - geometry_spherical - use finite_element_config_mod, only: coord_system - use planet_config_mod, only: scaled_radius + use mesh_mod, only : geometry_spherical implicit none @@ -39,12 +37,16 @@ module sci_gp_vector_rhs_kernel_mod !> type, public, extends(kernel_type) :: gp_vector_rhs_kernel_type private - type(arg_type) :: meta_args(5) = (/ & - arg_type(GH_FIELD*3, GH_REAL, GH_INC, ANY_SPACE_1), & - arg_type(GH_FIELD, GH_REAL, GH_READ, ANY_SPACE_2), & - 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_FIELD, GH_REAL, GH_READ, W2) & + type(arg_type) :: meta_args(9) = (/ & + arg_type(GH_FIELD*3, GH_REAL, GH_INC, ANY_SPACE_1), & ! rhs1, rhs2, rhs3 + arg_type(GH_FIELD, GH_REAL, GH_READ, ANY_SPACE_2), & ! field + arg_type(GH_FIELD*3, GH_REAL, GH_READ, ANY_SPACE_9), & ! chi1, chi2, chi3 + arg_type(GH_FIELD, GH_REAL, GH_READ, ANY_DISCONTINUOUS_SPACE_3), & ! panel_id + arg_type(GH_FIELD, GH_REAL, GH_READ, W2), & ! w2_field + 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(ANY_SPACE_1, GH_BASIS), & @@ -82,6 +84,10 @@ module sci_gp_vector_rhs_kernel_mod !! @param[in] chi_3 3rd coordinate field !! @param[in] panel_id Field giving the ID for mesh panels. !! @param[in] w2_field W2_field needed to get function space components +!! @param[in] geometry Mesh geometry enumeration +!! @param[in] topology Mesh topology enumeration +!! @param[in] coord_system Finite-element coordinate system enumeration +!! @param[in] scaled_radius Scaled planet radius !! @param[in] ndf Number of degrees of freedom per cell !! @param[in] undf Number of degrees of freedom !! @param[in] map Dofmap for the cell at the base of the column @@ -111,6 +117,8 @@ subroutine gp_vector_rhs_code(nlayers, & rhs1, rhs2, rhs3, field, & chi_1, chi_2, chi_3, & panel_id, w2_field, & + geometry, topology, coord_system, & + scaled_radius, & ndf, undf, map, basis, & ndf_f, undf_f, map_f, f_basis, & ndf_chi, undf_chi, & @@ -149,6 +157,11 @@ subroutine gp_vector_rhs_code(nlayers, & real(kind=r_def), dimension(nqp_h), intent(in) :: wqp_h real(kind=r_def), dimension(nqp_v), intent(in) :: wqp_v + 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, df2, k, qp1, qp2 real(kind=r_def), dimension(nqp_h,nqp_v) :: dj diff --git a/components/science/source/kernel/fem/sci_mg_derham_mat_kernel_mod.F90 b/components/science/source/kernel/fem/sci_mg_derham_mat_kernel_mod.F90 index 04329bd5b..3143289f2 100644 --- a/components/science/source/kernel/fem/sci_mg_derham_mat_kernel_mod.F90 +++ b/components/science/source/kernel/fem/sci_mg_derham_mat_kernel_mod.F90 @@ -15,8 +15,10 @@ module sci_mg_derham_mat_kernel_mod use argument_mod, only: arg_type, func_type, & - GH_OPERATOR, GH_FIELD, & - GH_REAL, GH_READ, GH_WRITE, & + GH_OPERATOR, & + GH_FIELD, GH_SCALAR, & + GH_REAL, GH_INTEGER, & + GH_READ, GH_WRITE, & ANY_SPACE_1, ANY_SPACE_9, & ANY_DISCONTINUOUS_SPACE_3, & GH_BASIS, GH_DIFF_BASIS, & @@ -27,10 +29,6 @@ module sci_mg_derham_mat_kernel_mod use fs_continuity_mod, only: W2, W3, wtheta use kernel_mod, only: kernel_type - use base_mesh_config_mod, only: geometry, topology - use finite_element_config_mod, only: coord_system - use planet_config_mod, only: scaled_radius - implicit none private @@ -41,13 +39,17 @@ module sci_mg_derham_mat_kernel_mod type, public, extends(kernel_type) :: mg_derham_mat_kernel_type private - type(arg_type) :: meta_args(6) = (/ & - arg_type(GH_OPERATOR, GH_REAL, GH_WRITE, W2, W2), & - arg_type(GH_OPERATOR, GH_REAL, GH_WRITE, W3, W3), & - arg_type(GH_OPERATOR, GH_REAL, GH_WRITE, Wtheta, Wtheta), & - arg_type(GH_OPERATOR, GH_REAL, GH_WRITE, W3, 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) & + type(arg_type) :: meta_args(10) = (/ & + arg_type(GH_OPERATOR, GH_REAL, GH_WRITE, W2, W2), & ! mm2 + arg_type(GH_OPERATOR, GH_REAL, GH_WRITE, W3, W3), & ! mm3 + arg_type(GH_OPERATOR, GH_REAL, GH_WRITE, Wtheta, Wtheta), & ! mmt + arg_type(GH_OPERATOR, GH_REAL, GH_WRITE, W3, W2), & ! div + arg_type(GH_FIELD*3, GH_REAL, GH_READ, ANY_SPACE_9), & ! chi1, chi2, chi3 + 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, GH_DIFF_BASIS), & @@ -83,6 +85,10 @@ module sci_mg_derham_mat_kernel_mod !! @param[in] chi2 Physical coordinates in the 2nd dir !! @param[in] chi3 Physical coordinates in the 3rd dir !! @param[in] panel_id Field giving the ID for mesh panels +!! @param[in] geometry Mesh geometry enumeration +!! @param[in] topology Mesh topology enumeration +!! @param[in] coord_system Finite-element coordinate system enumeration +!! @param[in] scaled_radius Scaled planet radius !! @param[in] ndf_w2 Number of degrees of freedom per cell for W2 space !! @param[in] basis_w2 Basis functions evaluated at quadrature points for W2 space !! @param[in] diff_basis_w2 Differential of basis functions evaluated at @@ -112,6 +118,8 @@ subroutine mg_derham_mat_code(cell, nlayers, & ncell_3d6, div, & chi1, chi2, chi3, & panel_id, & + geometry, topology, & + coord_system, scaled_radius, & ndf_w2, basis_w2, diff_basis_w2, & ndf_w3, basis_w3, & ndf_wt, basis_wt, & @@ -149,6 +157,11 @@ subroutine mg_derham_mat_code(cell, nlayers, & real(kind=r_def), intent(in) :: wqp_h(nqp_h) real(kind=r_def), intent(in) :: wqp_v(nqp_v) + 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, df2, k, ik integer(kind=i_def) :: qp1, qp2 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 22e5bd04b..0523f7556 100644 --- a/components/science/source/kernel/geometry/sci_chi_transform_mod.F90 +++ b/components/science/source/kernel/geometry/sci_chi_transform_mod.F90 @@ -29,11 +29,11 @@ module sci_chi_transform_mod LOG_LEVEL_DEBUG, & LOG_LEVEL_WARNING use matrix_invert_mod, only : matrix_invert_3x3 +use mesh_mod, only : geometry_spherical, & + geometry_planar, & + topology_periodic ! 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 @@ -180,7 +180,7 @@ subroutine init_chi_transforms( geometry, topology, & ) end if if ( abs(equatorial_latitude - rmdi) < EPS .or. & - geometry == geometry_planar .or. topology /= topology_fully_periodic ) then + geometry == geometry_planar .or. topology /= topology_periodic ) then equatorial_latitude = 0.0_r_def call log_event( & 'Equatorial latitude for mesh not set, so using 0.0 as default', & @@ -278,7 +278,7 @@ subroutine chi2xyz( chi_1, chi_2, chi_3, panel_id, & y = chi_2 z = chi_3 - else if (topology /= topology_fully_periodic) then + else if (topology /= topology_periodic) then ! domain is a spherical LAM, using (lon,lat,z) coordinates call llr2xyz(chi_1, chi_2, chi_3+scaled_radius, x, y, z) @@ -369,7 +369,7 @@ subroutine chir2xyz( chi_1, chi_2, chi_3, panel_id, & y = chi_2 z = chi_3 - else if (topology /= topology_fully_periodic) then + else if (topology /= topology_periodic) then ! domain is a spherical LAM, using (lon,lat,z) coordinates call llr2xyz(chi_1, chi_2, chi_3, x, y, z) @@ -457,7 +457,7 @@ subroutine chi2llr( chi_1, chi_2, chi_3, panel_id, & ! chi uses (geocentric) Cartesian coordinates call xyz2llr(chi_1, chi_2, chi_3, lon, lat, radius) - else if (topology /= topology_fully_periodic) then + else if (topology /= topology_periodic) then ! domain is a spherical LAM, already using (lon,lat,z) coordinates ! may need to rotate these to the physical (lon,lat) coordinates @@ -536,7 +536,7 @@ subroutine chi2abr( chi_1, chi_2, chi_3, panel_id, & real(kind=r_def) :: xyz(3) - if (topology /= topology_fully_periodic .or. geometry /= geometry_spherical) then + if (topology /= topology_periodic .or. geometry /= geometry_spherical) then call log_event( & 'chi2abr can only be used on cubed-sphere meshes', LOG_LEVEL_ERROR & ) diff --git a/components/science/source/kernel/geometry/sci_coordinate_jacobian_mod.F90 b/components/science/source/kernel/geometry/sci_coordinate_jacobian_mod.F90 index 852afbfa6..d1207e33e 100644 --- a/components/science/source/kernel/geometry/sci_coordinate_jacobian_mod.F90 +++ b/components/science/source/kernel/geometry/sci_coordinate_jacobian_mod.F90 @@ -19,15 +19,14 @@ module sci_coordinate_jacobian_mod xyz2ll, & llr2xyz, & schmidt_transform_lat - + use mesh_mod, only: geometry_planar, & + topology_periodic use sci_chi_transform_mod, only: get_mesh_rotation_matrix, & get_to_stretch, & get_to_rotate, & get_stretch_factor ! Configuration modules - use base_mesh_config_mod, only: geometry_planar, & - topology_fully_periodic use finite_element_config_mod, only: coord_system_xyz, & coord_system_native @@ -209,7 +208,7 @@ subroutine coordinate_jacobian_quadrature_real32( & ! Using (X,Y,Z) coordinates or on a plane jac = jac_ref2sph - else if (topology == topology_fully_periodic) then + else if (topology == topology_periodic) then ! Native coordinates for a cubed-sphere mesh to_rotate = get_to_rotate() to_stretch = get_to_stretch() @@ -372,7 +371,7 @@ subroutine coordinate_jacobian_quadrature_real64( & ! Using (X,Y,Z) coordinates or on a plane jac = jac_ref2sph - else if (topology == topology_fully_periodic) then + else if (topology == topology_periodic) then ! Native coordinates for a cubed-sphere mesh to_rotate = get_to_rotate() to_stretch = get_to_stretch() @@ -549,7 +548,7 @@ subroutine coordinate_jacobian_evaluator_real32( & ! Using (X,Y,Z) coordinates or on a plane jac = jac_ref2sph - else if (topology == topology_fully_periodic) then + else if (topology == topology_periodic) then ! Native coordinates for a cubed-sphere mesh to_rotate = get_to_rotate() to_stretch = get_to_stretch() @@ -690,7 +689,7 @@ subroutine coordinate_jacobian_evaluator_real64( & ! Using (X,Y,Z) coordinates or on a plane jac = jac_ref2sph - else if (topology == topology_fully_periodic) then + else if (topology == topology_periodic) then ! Native coordinates for a cubed-sphere mesh to_rotate = get_to_rotate() to_stretch = get_to_stretch() @@ -981,7 +980,7 @@ subroutine pointwise_coordinate_jacobian_real32( & ! Using (X,Y,Z) coordinates or on a plane jac = jac_ref2sph - else if (topology == topology_fully_periodic) then + else if (topology == topology_periodic) then ! Native coordinates for a cubed-sphere mesh alpha = 0.0_real32 beta = 0.0_real32 @@ -1100,7 +1099,7 @@ subroutine pointwise_coordinate_jacobian_real64( & ! Using (X,Y,Z) coordinates or on a plane jac = jac_ref2sph - else if (topology == topology_fully_periodic) then + else if (topology == topology_periodic) then ! Native coordinates for a cubed-sphere mesh alpha = 0.0_real64 beta = 0.0_real64 diff --git a/components/science/source/kernel/geometry/sci_height_continuous_kernel_mod.F90 b/components/science/source/kernel/geometry/sci_height_continuous_kernel_mod.F90 index 4964cb6aa..326cf565c 100644 --- a/components/science/source/kernel/geometry/sci_height_continuous_kernel_mod.F90 +++ b/components/science/source/kernel/geometry/sci_height_continuous_kernel_mod.F90 @@ -15,10 +15,10 @@ module sci_height_continuous_kernel_mod GH_READ, GH_INC, & ANY_SPACE_1, ANY_SPACE_9, & CELL_COLUMN, GH_BASIS, GH_EVALUATOR - use base_mesh_config_mod, only: geometry_spherical use constants_mod, only: r_def, i_def, l_def use finite_element_config_mod, only: coord_system_xyz use kernel_mod, only: kernel_type + use mesh_mod, only: geometry_spherical implicit none private diff --git a/components/science/source/kernel/geometry/sci_height_discontinuous_kernel_mod.F90 b/components/science/source/kernel/geometry/sci_height_discontinuous_kernel_mod.F90 index 101e8aa4b..af7f9d1b2 100644 --- a/components/science/source/kernel/geometry/sci_height_discontinuous_kernel_mod.F90 +++ b/components/science/source/kernel/geometry/sci_height_discontinuous_kernel_mod.F90 @@ -16,10 +16,10 @@ module sci_height_discontinuous_kernel_mod GH_READ, GH_WRITE, & ANY_DISCONTINUOUS_SPACE_1, ANY_SPACE_9, & CELL_COLUMN, GH_BASIS, GH_EVALUATOR - use base_mesh_config_mod, only: geometry_spherical use constants_mod, only: r_def, i_def, l_def use finite_element_config_mod, only: coord_system_xyz use kernel_mod, only: kernel_type + use mesh_mod, only: geometry_spherical implicit none private diff --git a/components/science/source/kernel/geometry/sci_native_jacobian_mod.F90 b/components/science/source/kernel/geometry/sci_native_jacobian_mod.F90 index 10116b3ee..866bea469 100644 --- a/components/science/source/kernel/geometry/sci_native_jacobian_mod.F90 +++ b/components/science/source/kernel/geometry/sci_native_jacobian_mod.F90 @@ -26,11 +26,11 @@ module sci_native_jacobian_mod get_to_stretch, & get_to_rotate, & get_stretch_factor + use mesh_mod, only: geometry_planar, & + topology_periodic use finite_element_config_mod, only: coord_system_xyz, & coord_system_native - use base_mesh_config_mod, only: geometry_planar, & - topology_fully_periodic implicit none @@ -131,7 +131,7 @@ subroutine native_jacobian(coord_system, geometry, topology, scaled_radius, & ! Using (X,Y,Z) coordinates or on a plane jac = jac_ref2sph - else if (topology == topology_fully_periodic) then + else if (topology == topology_periodic) then radius = real(scaled_radius, kind=r_def) jac_sph2XYZ = jacobian_abr2XYZ(nlayers, chi_1_df, chi_2_df, chi_3_df+radius, panel_id) 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 f05d600f8..f8c44e203 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 @@ -28,8 +28,7 @@ module sci_compute_map_u_operators_kernel_mod use fs_continuity_mod, only : W2, W3, Wtheta 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 + use mesh_mod, only : geometry_spherical, geometry_planar implicit none 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 eb8514b3c..127ae2900 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 @@ -28,13 +28,14 @@ module sci_compute_sample_u_ops_kernel_mod use constants_mod, only : r_def, i_def use fs_continuity_mod, only : W2broken, W3, Wtheta use kernel_mod, only : kernel_type + use sci_chi_transform_mod, only : chi2llr use sci_coordinate_jacobian_mod, only : coordinate_jacobian, & coordinate_jacobian_inverse use coord_transform_mod, only : sphere2cart_vector use reference_element_mod, only : W, S, N, E, T, B + use mesh_mod, only : geometry_spherical, geometry_planar - use base_mesh_config_mod, only: geometry_spherical, geometry_planar implicit none 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 2888e3577..9afedf8bc 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 @@ -21,8 +21,7 @@ module sci_convert_phys_to_hdiv_kernel_mod use constants_mod, only : r_def, i_def use fs_continuity_mod, only : W2 use kernel_mod, only : kernel_type - - use base_mesh_config_mod, only: geometry_spherical + use mesh_mod, only : geometry_spherical implicit none 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 d00964b45..51cf41bb0 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,8 +26,7 @@ module sci_project_ws_to_w1_operator_kernel_mod use constants_mod, only : r_def, i_def use fs_continuity_mod, only : W1 use log_mod, only : log_event, LOG_LEVEL_ERROR - -use base_mesh_config_mod, only: geometry_spherical, geometry_planar +use mesh_mod, only : geometry_spherical, geometry_planar implicit none diff --git a/components/science/source/solver/sci_hierarchical_linear_operator_mod.f90 b/components/science/source/solver/sci_hierarchical_linear_operator_mod.f90 index 69cb77928..39dd957bc 100644 --- a/components/science/source/solver/sci_hierarchical_linear_operator_mod.f90 +++ b/components/science/source/solver/sci_hierarchical_linear_operator_mod.f90 @@ -12,6 +12,7 @@ module sci_hierarchical_linear_operator_mod use sci_linear_operator_mod, only : abstract_linear_operator_type + use config_mod, only: config_type implicit none @@ -28,12 +29,15 @@ module sci_hierarchical_linear_operator_mod !> Abstract interface defined for the coarsening linear operator to the !> next level of the multigrid hierarchy. !> - !>@param[in] self a hierarchical linear operator - !>@param[in] fs_coarse coarse level function space - !>@param[inout] coarse_operator coarsened version on the next multigrid level - subroutine coarsen_interface(self,coarse_operator) + !>@param[in] config Application configuration object + !>@param[inout] coarse_operator Coarsened version on the next multigrid level + subroutine coarsen_interface(self, config, coarse_operator) import :: abstract_hierarchical_linear_operator_type + import :: config_type + class(abstract_hierarchical_linear_operator_type), intent(in) :: self + + type(config_type), intent(in) :: config class(abstract_hierarchical_linear_operator_type), allocatable, intent(inout) :: coarse_operator end subroutine coarsen_interface end interface diff --git a/components/science/source/solver/sci_hierarchical_preconditioner_mod.f90 b/components/science/source/solver/sci_hierarchical_preconditioner_mod.f90 index e63d8f42f..040be6cf7 100644 --- a/components/science/source/solver/sci_hierarchical_preconditioner_mod.f90 +++ b/components/science/source/solver/sci_hierarchical_preconditioner_mod.f90 @@ -14,6 +14,7 @@ module sci_hierarchical_preconditioner_mod use sci_preconditioner_mod, only: abstract_preconditioner_type + use config_mod, only: config_type implicit none private @@ -29,11 +30,15 @@ module sci_hierarchical_preconditioner_mod !> Abstract interface defined for the coarsening preconditioner to the !> next level of the multigrid hierarchy. !> - !>@param[inout] self a hierarchical preconditioner - !>@param[inout] other coarsened version on the next multigrid level - subroutine coarsen_interface(self, other) + !>@param[in] config Application configuration object + !>@param[inout] other Coarsened version on the next multigrid level + subroutine coarsen_interface(self, config, other) import :: abstract_hierarchical_preconditioner_type + import :: config_type + class(abstract_hierarchical_preconditioner_type), intent(inout) :: self + + type(config_type), intent(in) :: config class(abstract_hierarchical_preconditioner_type), allocatable, intent(inout) :: other end subroutine coarsen_interface end interface diff --git a/components/science/source/solver/sci_iterative_solver_mod.f90 b/components/science/source/solver/sci_iterative_solver_mod.f90 index f85ef074f..b955242f2 100644 --- a/components/science/source/solver/sci_iterative_solver_mod.f90 +++ b/components/science/source/solver/sci_iterative_solver_mod.f90 @@ -14,6 +14,7 @@ module sci_iterative_solver_mod + use config_mod, only : config_type use constants_mod, only : r_def, i_def, EPS, l_def use vector_mod, only : abstract_vector_type use sci_linear_operator_mod, & @@ -63,15 +64,22 @@ module sci_iterative_solver_mod !> @details Apply the iterative solver for a given right hand side !> \f$b\f$. !> - !> @param[inout] x Resulting solution \f$x\f$ - !> @param[in] b Right hand side vector \f$b\f$ + !> @param[in] config Application configuration object + !> @param[inout] x Resulting solution \f$x\f$ + !> @param[in] b Right hand side vector \f$b\f$ !> - subroutine apply_interface(self, x, b) + subroutine apply_interface(self, config, x, b) + import :: abstract_vector_type import :: abstract_iterative_solver_type + import :: config_type + class(abstract_iterative_solver_type), intent(inout) :: self + + type(config_type), intent(in) :: config class(abstract_vector_type), intent(inout) :: x class(abstract_vector_type), intent(inout) :: b + end subroutine apply_interface end interface @@ -113,8 +121,9 @@ module function cg_constructor( lin_op, prec, r_tol, a_tol, max_iter, & end interface interface - module subroutine cg_solve(self, x, b) + module subroutine cg_solve(self, config, x, b) class(conjugate_gradient_type), intent(inout) :: self + type(config_type), intent(in) :: config class(abstract_vector_type), intent(inout) :: x class(abstract_vector_type), intent(inout) :: b end subroutine @@ -149,8 +158,9 @@ module function bicgstab_constructor( lin_op, prec, r_tol, a_tol, max_iter, & end interface interface - module subroutine bicgstab_solve(self, x, b) + module subroutine bicgstab_solve(self, config, x, b) class(bicgstab_type), intent(inout) :: self + type(config_type), intent(in) :: config class(abstract_vector_type), intent(inout) :: x class(abstract_vector_type), intent(inout) :: b end subroutine @@ -188,8 +198,9 @@ end function gmres_constructor end interface interface - module subroutine gmres_solve(self, x, b) + module subroutine gmres_solve(self, config, x, b) class(gmres_type), intent(inout) :: self + type(config_type), intent(in) :: config class(abstract_vector_type), intent(inout) :: x class(abstract_vector_type), intent(inout) :: b end subroutine gmres_solve @@ -226,8 +237,9 @@ end function fgmres_constructor end interface interface - module subroutine fgmres_solve(self, x, b) - class(fgmres_type), intent(inout) :: self + module subroutine fgmres_solve(self, config, x, b) + class(fgmres_type), intent(inout) :: self + type(config_type), intent(in) :: config class(abstract_vector_type), intent(inout) :: x class(abstract_vector_type), intent(inout) :: b end subroutine fgmres_solve @@ -264,8 +276,9 @@ end function gcr_constructor end interface interface - module subroutine gcr_solve(self, x, b) + module subroutine gcr_solve(self, config, x, b) class(gcr_type), intent(inout) :: self + type(config_type), intent(in) :: config class(abstract_vector_type), intent(inout) :: x class(abstract_vector_type), intent(inout) :: b end subroutine gcr_solve @@ -304,8 +317,9 @@ end function block_gcr_constructor end interface interface - module subroutine block_gcr_solve(self, x, b) + module subroutine block_gcr_solve(self, config, x, b) class(block_gcr_type), intent(inout) :: self + type(config_type), intent(in) :: config class(abstract_vector_type), intent(inout) :: x class(abstract_vector_type), intent(inout) :: b end subroutine block_gcr_solve @@ -339,8 +353,9 @@ module function precondition_only_constructor( lin_op, prec, & end interface interface - module subroutine precondition_only_solve(self, x, b) + module subroutine precondition_only_solve(self, config, x, b) class(precondition_only_type), intent(inout) :: self + type(config_type), intent(in) :: config class(abstract_vector_type), intent(inout) :: x class(abstract_vector_type), intent(inout) :: b end subroutine @@ -380,10 +395,11 @@ module function jacobi_constructor( lin_op, prec, r_tol, a_tol, max_iter, end interface interface - module subroutine jacobi_solve(self, x, b) - class(jacobi_type), intent(inout) :: self - class(abstract_vector_type), intent(inout) :: x - class(abstract_vector_type), intent(inout) :: b + module subroutine jacobi_solve(self, config, x, b) + class(jacobi_type), intent(inout) :: self + type(config_type), intent(in) :: config + class(abstract_vector_type), intent(inout) :: x + class(abstract_vector_type), intent(inout) :: b end subroutine end interface @@ -421,8 +437,9 @@ module function chebyshev_constructor( lin_op, prec, r_tol, a_tol, max_iter, & end interface interface - module subroutine chebyshev_solve(self, x, b) + module subroutine chebyshev_solve(self, config, x, b) class(chebyshev_type), intent(inout) :: self + type(config_type), intent(in) :: config class(abstract_vector_type), intent(inout) :: x class(abstract_vector_type), intent(inout) :: b end subroutine @@ -490,12 +507,16 @@ module function cg_constructor(lin_op, prec, r_tol, a_tol, max_iter, & !> !> Over-rides the abstract interface to do the actual solve. !> - !> @param[inout] b "RHS" or boundary conditions. - !> @param[inout] x Solution. + !> @param[in] config Application configuration object + !> @param[inout] b "RHS" or boundary conditions. + !> @param[inout] x Solution. !> - module subroutine cg_solve(self, x, b) + module subroutine cg_solve(self, config, x, b) + implicit none + class(conjugate_gradient_type), intent(inout) :: self + type(config_type), intent(in) :: config class(abstract_vector_type), intent(inout) :: x class(abstract_vector_type), intent(inout) :: b @@ -519,7 +540,7 @@ module subroutine cg_solve(self, x, b) converged=.false. !set up the algorithm - call self%lin_op%apply(x,r) ! r = A.x + call self%lin_op%apply(config, x, r) ! r = A.x call r%scale(-1.0_r_def) ! r = -A.x call r%axpy(1.0_r_def, b) ! r = b - A.x r_nrm_0 = r%norm() ! r_0 = ||r||_2 @@ -539,8 +560,8 @@ module subroutine cg_solve(self, x, b) end if call z%set_scalar(0.0_r_def) - call self%prec%apply(r,z) ! z = P^{-1}.r - rz = r%dot(z) ! rz = + call self%prec%apply(config, r, z) ! z = P^{-1}.r + rz = r%dot(z) ! rz = r_nrm_old = r_nrm_0 call p%copy(z) @@ -550,10 +571,10 @@ module subroutine cg_solve(self, x, b) end if ! iterate until maximal number of iterations is reached do iter=1, self%max_iter - call self%lin_op%apply(p,z) ! z = A.p - alpha = rz / p%dot(z) ! alpha = / - call x%axpy(alpha,p) ! x -> x + alpha*p - call r%axpy(-alpha,z) ! r -> r - alpha*A.p + call self%lin_op%apply(config, p, z) ! z = A.p + alpha = rz / p%dot(z) ! alpha = / + call x%axpy(alpha,p) ! x -> x + alpha*p + call r%axpy(-alpha,z) ! r -> r - alpha*A.p if ( self%monitor_convergence ) then r_nrm = r%norm() ! r = ||r||_2 @@ -567,10 +588,10 @@ module subroutine cg_solve(self, x, b) exit end if end if - call self%prec%apply(r,z) ! z = P^{-1}.r - rz_new = r%dot(z) ! rz_new = - beta = rz_new/rz ! beta = / - call p%aypx(beta,z) ! p -> z + beta*p + call self%prec%apply(config, r, z) ! z = P^{-1}.r + rz_new = r%dot(z) ! rz_new = + beta = rz_new/rz ! beta = / + call p%aypx(beta,z) ! p -> z + beta*p rz = rz_new r_nrm_old = r_nrm end do @@ -640,15 +661,19 @@ module function bicgstab_constructor( lin_op, prec, r_tol, a_tol, max_iter, & end function bicgstab_constructor !> bicgstab solve. Over-rides the abstract interface to do the actual solve. - !> @param[inout] b an abstract vector which will be an actual vector of unkown extended type - !> This the "RHS" or boundary conditions, - !> @param[inout] x an abstract vector which is the solution + !> @param[in] config Application configuration object + !> @param[inout] b An abstract vector which will be an actual vector of + !> unkown extended type. This the "RHS" or boundary conditions, + !> @param[inout] x An abstract vector which is the solution !> @param[self] The solver which has pointers to the lin_op and preconditioner - module subroutine bicgstab_solve(self, x, b) + module subroutine bicgstab_solve(self, config, x, b) + implicit none - class(bicgstab_type), intent(inout) :: self - class(abstract_vector_type), intent(inout) :: x - class(abstract_vector_type), intent(inout) :: b + + class(bicgstab_type), intent(inout) :: self + type(config_type), intent(in) :: config + class(abstract_vector_type), intent(inout) :: x + class(abstract_vector_type), intent(inout) :: b ! tempory vectors class(abstract_vector_type), allocatable :: r @@ -674,7 +699,7 @@ module subroutine bicgstab_solve(self, x, b) ! v = Ax call x%duplicate(v) call v%set_scalar(0.0_r_def) - call self%lin_op%apply(x,v) + call self%lin_op%apply(config, x, v) ! r = b - Ax call r%axpy(-1.0_r_def,v) ! store initial residual @@ -721,20 +746,20 @@ module subroutine bicgstab_solve(self, x, b) call t%copy(r) call t%axpy(-beta*omega, v) ! stage 2 post-condition - call self%prec%apply(t,y) + call self%prec%apply(config, t, y) ! now add on beta P call p%aypx(beta, y) ! apply the matrix - call self%lin_op%apply(p,v) + call self%lin_op%apply(config, p, v) alpha = rho/r0%dot(v) ! s = r - alpha * v call s%copy(r) call s%axpy(-alpha,v) ! apply the preconditioner - call self%prec%apply(s,z) + call self%prec%apply(config, s, z) ! apply the operator - call self%lin_op%apply(z,t) + call self%lin_op%apply(config, z, t) ! final scalars tt = t%dot(t) @@ -826,15 +851,19 @@ end function gmres_constructor !> gmres_solve. Over-rides the abstract interface to do the actual solve. !> @detail The solver implements left-preconditioning, i.e. is solving M{-1}.A.x = M{-1}.b + !> @param[in] config Application configuration object !> @param[inout] b an abstract vector which will be an actual vector of unkown extended type !> This the "RHS" or boundary conditions, !> @param[inout] x an abstract vector which is the solution !> @param[self] The solver which has pointers to the lin_op and preconditioner - module subroutine gmres_solve(self, x, b) + module subroutine gmres_solve(self, config, x, b) + implicit none - class(gmres_type), intent(inout) :: self - class(abstract_vector_type), intent(inout) :: x - class(abstract_vector_type), intent(inout) :: b + + class(gmres_type), intent(inout) :: self + type(config_type), intent(in) :: config + class(abstract_vector_type), intent(inout) :: x + class(abstract_vector_type), intent(inout) :: b ! temporary vectors class(abstract_vector_type), allocatable :: s @@ -859,7 +888,7 @@ module subroutine gmres_solve(self, x, b) call b%duplicate(res) ! compute res = b -Ax ... in stages - call self%lin_op%apply(x,Ax) + call self%lin_op%apply(config, x, Ax) call res%copy(b) call res%axpy(-1.0_r_def,Ax) @@ -899,7 +928,7 @@ module subroutine gmres_solve(self, x, b) ! initialisation complete, lets go to work. do iter = 1, self%max_iter call s%set_scalar(0.0_r_def) - call self%prec%apply(res,s) + call self%prec%apply(config, res, s) beta = s%norm() call v(1)%vt%copy(s) call v(1)%vt%scale(1.0_r_def/beta) @@ -911,9 +940,9 @@ module subroutine gmres_solve(self, x, b) call w%copy(v(iv)%vt) ! apply the operator - call self%lin_op%apply( w, s ) + call self%lin_op%apply(config, w, s) ! apply the preconditioner - call self%prec%apply( s, w ) + call self%prec%apply(config, s, w) ! compute the h values do ivj = 1, iv @@ -963,7 +992,7 @@ module subroutine gmres_solve(self, x, b) end do call Ax%set_scalar(0.0_r_def) - call self%lin_op%apply(x, Ax) + call self%lin_op%apply(config, x, Ax) call res%copy(Ax) call res%aypx(-1.0_r_def, b) @@ -1052,15 +1081,17 @@ end function fgmres_constructor !> @detail The solver implements flexible right-preconditioning, i.e. is solving A.M{-1}.M.x = b !> as 1) A.M{-1}.y = b !> 2) M{-1}y = x + !> @param[in] config Application configuration object !> @param[inout] b an abstract vector which will be an actual vector of unkown extended type !> This the "RHS" or boundary conditions, !> @param[inout] x an abstract vector which is the solution !> @param[self] The solver which has pointers to the lin_op and preconditioner - module subroutine fgmres_solve(self, x, b) + module subroutine fgmres_solve(self, config, x, b) implicit none - class(fgmres_type), intent(inout) :: self - class(abstract_vector_type), intent(inout) :: x - class(abstract_vector_type), intent(inout) :: b + class(fgmres_type), intent(inout) :: self + type(config_type), intent(in) :: config + class(abstract_vector_type), intent(inout) :: x + class(abstract_vector_type), intent(inout) :: b ! temporary vectors class(abstract_vector_type), allocatable :: s @@ -1138,9 +1169,9 @@ module subroutine fgmres_solve(self, x, b) do iter = 1, self%max_iter do iv = 1, self%gcrk - call self%prec%apply(v(iv)%vt, Pv(iv)%vt) + call self%prec%apply(config, v(iv)%vt, Pv(iv)%vt) ! apply the operator - call self%lin_op%apply( Pv(iv)%vt, s ) + call self%lin_op%apply(config, Pv(iv)%vt, s) call w%copy(s) ! compute the h values do ivj = 1, iv @@ -1189,7 +1220,7 @@ module subroutine fgmres_solve(self, x, b) end do ! check for convergence - call self%lin_op%apply(dx, Ax) + call self%lin_op%apply(config, dx, Ax) call res%copy(Ax) call res%aypx(-1.0_r_def, b) @@ -1286,15 +1317,19 @@ end function gcr_constructor !> @detail The solver implements right-preconditioning, i.e. is solving A.M{-1}.M.x = b !> as 1) A.M{-1}.y = b !> 2) M{-1}y = x + !> @param[in] config Application configuration object !> @param[inout] b an abstract vector which will be an actual vector of unkown extended type !> This the "RHS" or boundary conditions, !> @param[inout] x an abstract vector which is the solution !> @param[self] The solver which has pointers to the lin_op and preconditioner - module subroutine gcr_solve(self, x, b) + module subroutine gcr_solve(self, config, x, b) + implicit none - class(gcr_type), intent(inout) :: self - class(abstract_vector_type), intent(inout) :: x - class(abstract_vector_type), intent(inout) :: b + + class(gcr_type), intent(inout) :: self + type(config_type), intent(in) :: config + class(abstract_vector_type), intent(inout) :: x + class(abstract_vector_type), intent(inout) :: b ! temporary vectors class(abstract_vector_type), allocatable :: dx @@ -1317,8 +1352,8 @@ module subroutine gcr_solve(self, x, b) ! initial guess call dx%set_scalar(0.0_r_def) - call self%prec%apply( b, dx ) - call self%lin_op%apply( dx, Ax ) + call self%prec%apply(config, b, dx) + call self%lin_op%apply(config, dx, Ax) !res = b - Ax call res%copy(b) call res%axpy(-1.0_r_def, Ax) @@ -1356,9 +1391,9 @@ module subroutine gcr_solve(self, x, b) do iv = 1, self%gcrk ! apply the preconditioner - call self%prec%apply( res, Pv(iv)%vt ) + call self%prec%apply(config, res, Pv(iv)%vt) ! apply the operator - call self%lin_op%apply( Pv(iv)%vt, v(iv)%vt ) + call self%lin_op%apply(config, Pv(iv)%vt, v(iv)%vt) do ivj = 1, iv-1 alpha = v(iv)%vt%dot(v(ivj)%vt) @@ -1460,15 +1495,17 @@ module function block_gcr_constructor( lin_op, prec, gcrk, r_tol, a_tol, max_ite !> @detail The solver implements right-preconditioning, i.e. is solving A.M{-1}.M.x = b !> as 1) A.M{-1}.y = b !> 2) M{-1}y = x + !> @param[in] config Application configuration object !> @param[inout] b an abstract vector which will be an actual vector of unkown extended type !> This the "RHS" or boundary conditions, !> @param[inout] x an abstract vector which is the solution !> @param[self] The solver which has pointers to the lin_op and preconditioner - module subroutine block_gcr_solve(self, x, b) + module subroutine block_gcr_solve(self, config, x, b) implicit none - class(block_gcr_type), intent(inout) :: self - class(abstract_vector_type), intent(inout) :: x - class(abstract_vector_type), intent(inout) :: b + class(block_gcr_type), intent(inout) :: self + type(config_type), intent(in) :: config + class(abstract_vector_type), intent(inout) :: x + class(abstract_vector_type), intent(inout) :: b ! temporary vectors class(abstract_vector_type), allocatable :: Ax @@ -1501,8 +1538,8 @@ module subroutine block_gcr_solve(self, x, b) ! initial guess call x%set_scalar(0.0_r_def) - call self%prec%apply( b, x ) - call self%lin_op%apply( x, Ax ) + call self%prec%apply(config, b, x) + call self%lin_op%apply(config, x, Ax) call res%copy(b) call res%axpy(-1.0_r_def, Ax) @@ -1536,9 +1573,9 @@ module subroutine block_gcr_solve(self, x, b) do iv = 1, self%gcrk ! apply the preconditioner - call self%prec%apply( res, Pv(iv)%vt ) + call self%prec%apply(config, res, Pv(iv)%vt) ! apply the operator - call self%lin_op%apply( Pv(iv)%vt, v(iv)%vt ) + call self%lin_op%apply(config, Pv(iv)%vt, v(iv)%vt) do ivj = 1, iv-1 alpha = v(iv)%vt%dot(v(ivj)%vt) @@ -1661,13 +1698,17 @@ module function precondition_only_constructor(lin_op, prec, monitor_convergence) end function !> Precondition only solve. Over-rides the abstract interface to do the actual solve. + !> @param[in] config Application configuration object !> @param[inout] b an abstract vector which will be an actual vector of unkown extended type !> This the "RHS" or boundary conditions, !> @param[inout] x an abstract vector which is the solution !> @param[self] The solver which has pointers to the lin_op and preconditioner - module subroutine precondition_only_solve(self, x, b) + module subroutine precondition_only_solve(self, config, x, b) + implicit none + class(precondition_only_type), intent(inout) :: self + type(config_type), intent(in) :: config class(abstract_vector_type), intent(inout) :: x class(abstract_vector_type), intent(inout) :: b @@ -1677,7 +1718,7 @@ module subroutine precondition_only_solve(self, x, b) call log_event("Precondition only starting", LOG_LEVEL_DEBUG) - call self%prec%apply(b,x) ! x = P^{-1}.b + call self%prec%apply(config, b, x) ! x = P^{-1}.b if( self%monitor_convergence ) then ! Compute initial and final error @@ -1685,7 +1726,7 @@ module subroutine precondition_only_solve(self, x, b) call res%copy(b) e0 = res%norm() call x%duplicate(Ax) - call self%lin_op%apply(x, Ax) + call self%lin_op%apply(config, x, Ax) call res%axpy(-1.0_r_def, Ax) e = res%norm() write(log_scratch_space,'(A,3E15.8)') & @@ -1743,13 +1784,17 @@ module function jacobi_constructor(lin_op, prec, r_tol, a_tol, max_iter, & end function !> Jacobi solve. Over-rides the abstract interface to do the actual solve. + !> @param[in] config Application configuration object !> @param[inout] b an abstract vector which will be an actual vector of unkown extended type !> This the "RHS" or boundary conditions, !> @param[inout] x an abstract vector which is the solution !> @param[self] The solver which has pointers to the lin_op and preconditioner - module subroutine jacobi_solve(self, x, b) + module subroutine jacobi_solve(self, config, x, b) + implicit none + class(jacobi_type), intent(inout) :: self + type(config_type), intent(in) :: config class(abstract_vector_type), intent(inout) :: x class(abstract_vector_type), intent(inout) :: b @@ -1774,7 +1819,7 @@ module subroutine jacobi_solve(self, x, b) ! Solve Ax = b do iter=1, self%max_iter - call self%lin_op%apply(x, r) ! r = Ax + call self%lin_op%apply(config, x, r) ! r = Ax call r%axpy(-1.0_r_def, b) ! r = r - b = Ax - b ! Check for convergence @@ -1792,7 +1837,7 @@ module subroutine jacobi_solve(self, x, b) end if end if - call self%prec%apply(r, z) ! z = P^{-1}.r = P^{-1}.(Ax - b) + call self%prec%apply(config, r, z) ! z = P^{-1}.r = P^{-1}.(Ax - b) call x%axpy(const, z) ! x = x + const*z = x + const*P^{-1}(Ax - b) end do @@ -1867,13 +1912,17 @@ module function chebyshev_constructor(lin_op, prec, r_tol, a_tol, max_iter, end function !> chebyshev solve. Over-rides the abstract interface to do the actual solve. + !> @param[in] config Application configuration object !> @param[inout] b an abstract vector which will be an actual vector of unkown extended type !> This the "RHS" or boundary conditions, !> @param[inout] x an abstract vector which is the solution !> @param[self] The solver which has pointers to the lin_op and preconditioner - module subroutine chebyshev_solve(self, x, b) + module subroutine chebyshev_solve(self, config, x, b) + implicit none + class(chebyshev_type), intent(inout) :: self + type(config_type), intent(in) :: config class(abstract_vector_type), intent(inout) :: x class(abstract_vector_type), intent(inout) :: b @@ -1908,11 +1957,11 @@ module subroutine chebyshev_solve(self, x, b) do iter = 1, self%max_iter ! r = b-M*xo - call self%lin_op%apply(xo,z) + call self%lin_op%apply(config, xo, z) call r%axpby(1.0_r_def, b, -1.0_r_def, z) ! z = D^{-1}.r - call self%prec%apply(r,z) + call self%prec%apply(config, r, z) ! x = w*(a/b*z+xo) + (1-w)*xp w = 1.0_r_def/(1.0_r_def - w/(4.0_r_def*a2**2)) @@ -1931,7 +1980,7 @@ module subroutine chebyshev_solve(self, x, b) end do ! residiual = norm(b - M*x) if ( self%monitor_convergence ) then - call self%lin_op%apply(x,z) ! z = M.x + call self%lin_op%apply(config, x, z) ! z = M.x call z%axpy(-1.0_r_def, b) ! z = M.x-b final_norm = z%norm() write(log_scratch_space, & diff --git a/components/science/source/solver/sci_linear_operator_mod.f90 b/components/science/source/solver/sci_linear_operator_mod.f90 index 87ef31dab..78c8bef41 100644 --- a/components/science/source/solver/sci_linear_operator_mod.f90 +++ b/components/science/source/solver/sci_linear_operator_mod.f90 @@ -9,8 +9,11 @@ !! defines an interface for the linear operator application y = A.x module sci_linear_operator_mod - use vector_mod, only : abstract_vector_type - use function_space_mod, only : function_space_type + + use config_mod, only: config_type + use function_space_mod, only: function_space_type + use vector_mod, only: abstract_vector_type + implicit none private @@ -22,16 +25,21 @@ module sci_linear_operator_mod end type abstract_linear_operator_type abstract interface - !> Abstract interface defined for the apply method. - !> param[in] self a linear operator - !> param[in] x a vector the linear operator is applied to - !> param[inout] y a vector, the result. - subroutine apply_interface(self, x, y) + !> @brief Abstract interface defined for the apply method. + !> @param[in] config Application configuration object + !> @param[in] x a vector the linear operator is applied to + !> @param[inout] y a vector, the result. + subroutine apply_interface(self, config, x, y) import :: abstract_linear_operator_type import :: abstract_vector_type + import :: config_type + class(abstract_linear_operator_type), intent(inout) :: self - class(abstract_vector_type), intent(in) :: x - class(abstract_vector_type), intent(inout) :: y + + type(config_type), intent(in) :: config + class(abstract_vector_type), intent(in) :: x + class(abstract_vector_type), intent(inout) :: y + end subroutine apply_interface end interface diff --git a/components/science/source/solver/sci_preconditioner_mod.f90 b/components/science/source/solver/sci_preconditioner_mod.f90 index 2b88c9a6a..a602a2961 100644 --- a/components/science/source/solver/sci_preconditioner_mod.f90 +++ b/components/science/source/solver/sci_preconditioner_mod.f90 @@ -10,7 +10,9 @@ !> defines an interface for the preconditioner application y = P^{-1}.x module sci_preconditioner_mod - use vector_mod, only : abstract_vector_type + + use config_mod, only: config_type + use vector_mod, only: abstract_vector_type implicit none private @@ -25,15 +27,20 @@ module sci_preconditioner_mod abstract interface !> abstract interface defined for the apply procedure of a preconditioner !! y = P.x - !> @param[in] self a preconditioner - !> @param[in] x a vector that the preconditioner is applied to. - !> @param[inout] y a vector, the result. - subroutine apply_interface(self, x, y) + !> @param[in] config Application configuration object + !> @param[in] x A vector that the preconditioner is applied to. + !> @param[inout] y A vector, the result. + subroutine apply_interface(self, config, x, y) + import :: abstract_vector_type import :: abstract_preconditioner_type + import :: config_type + class(abstract_preconditioner_type), intent(inout) :: self - class(abstract_vector_type), intent(in) :: x - class(abstract_vector_type), intent(inout) :: y + + type(config_type), intent(in) :: config + class(abstract_vector_type), intent(in) :: x + class(abstract_vector_type), intent(inout) :: y end subroutine apply_interface end interface diff --git a/components/science/source/utilities/sci_query_mod.f90 b/components/science/source/utilities/sci_query_mod.f90 index 7613aabcc..1db2491f2 100644 --- a/components/science/source/utilities/sci_query_mod.f90 +++ b/components/science/source/utilities/sci_query_mod.f90 @@ -16,9 +16,9 @@ module sci_query_mod implicit none private - public :: valid_for_global_model, & - check_lbc, & - is_lbc + public :: valid_for_global_model + public :: check_lbc + public :: is_lbc interface is_lbc module procedure is_global_lbc @@ -109,7 +109,6 @@ function check_lbc( mesh_name ) result ( answer ) character(4), parameter :: lbc_tag = '-lbc' integer(i_def) :: mesh_name_length -! character(4) :: tmp_str ! Check to see if the loaded mesh is an LBC mesh ! where the lbc_suffix would be appended. diff --git a/components/science/unit-test/kernel/fem/compute_broken_div_operator_kernel_mod_test.pf b/components/science/unit-test/kernel/fem/compute_broken_div_operator_kernel_mod_test.pf index ea0f476ef..e22e9ce2d 100644 --- a/components/science/unit-test/kernel/fem/compute_broken_div_operator_kernel_mod_test.pf +++ b/components/science/unit-test/kernel/fem/compute_broken_div_operator_kernel_mod_test.pf @@ -7,7 +7,9 @@ !> module compute_broken_div_operator_kernel_mod_test - use constants_mod, only : i_def, r_def, imdi + use constants_mod, only : l_def, i_def, r_def, imdi, rmdi + use finite_element_config_mod, only: coord_system_xyz + use funit use get_unit_test_m3x3_dofmap_mod, & only : get_w0_m3x3_dofmap, & get_w3_m3x3_dofmap @@ -25,76 +27,47 @@ module compute_broken_div_operator_kernel_mod_test get_gaussian_q3x3x3_quadrature_weights_z use get_unit_test_3x3x3_chi_mod, & only : get_w0_3x3x3_field - use funit implicit none private - public :: test_all - - @TestCase - type, extends(TestCase), public :: compute_broken_div_operator_test_type - private - contains - procedure setUp - procedure tearDown - procedure test_all - end type compute_broken_div_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, coord_space_wchi - use feign_config_mod, only : feign_finite_element_config - use sci_chi_transform_mod, only : init_chi_transforms + use sci_chi_transform_mod, only: init_chi_transforms implicit none - class(compute_broken_div_operator_test_type), intent(inout) :: this - - call feign_finite_element_config( & - cellshape=cellshape_quadrilateral, & - coord_order=0_i_def, & - coord_order_nonprime=1_i_def, & - coord_space=coord_space_wchi, & - coord_system=coord_system_xyz, & - element_order_h=1_i_def, & - element_order_v=1_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 config_loader_mod, only: final_configuration use sci_chi_transform_mod, only: final_chi_transforms implicit none - class(compute_broken_div_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_compute_broken_div_operator_kernel_mod, & only : compute_broken_div_operator_code implicit none - class(compute_broken_div_operator_test_type), intent(inout) :: this - real(r_def), parameter :: tol = 1.0e-6_r_def real(r_def), parameter :: dx = 6000.0_r_def real(r_def), parameter :: dy = 1000.0_r_def @@ -102,7 +75,12 @@ contains integer(i_def), parameter :: nlayers = 3 - integer(i_def) :: i, j, k, icell + 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 + logical(l_def), parameter :: rehabilitate = .true. + integer(i_def) :: ndf_w0, undf_w0, ndf_w2b integer(i_def) :: ndf_w3, ndf_pid, undf_pid integer(i_def) :: nqp_h, nqp_v @@ -165,6 +143,11 @@ contains chi2, & chi3, & panel_id, & + geometry, & + topology, & + coord_system, & + scaled_radius, & + rehabilitate, & ndf_w3, & basis_w3, & ndf_w2b, & diff --git a/components/science/unit-test/kernel/fem/compute_curl_operator_kernel_mod_test.pf b/components/science/unit-test/kernel/fem/compute_curl_operator_kernel_mod_test.pf index 12ae5c4c5..06b536476 100644 --- a/components/science/unit-test/kernel/fem/compute_curl_operator_kernel_mod_test.pf +++ b/components/science/unit-test/kernel/fem/compute_curl_operator_kernel_mod_test.pf @@ -8,7 +8,7 @@ !> module compute_curl_operator_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_dofmap_mod, only : get_w0_m3x3_dofmap, & get_w3_m3x3_dofmap use get_unit_test_m3x3_q3x3x3_sizes_mod, only : get_w0_m3x3_q3x3x3_size, & @@ -25,87 +25,62 @@ module compute_curl_operator_kernel_mod_test use get_unit_test_3x3x3_chi_mod, only : get_w0_3x3x3_field use funit + use finite_element_config_mod, only: coord_system_xyz + implicit none private - public :: test_all - @TestCase - type, extends(TestCase), public :: compute_curl_operator_test_type - private - contains - procedure setUp - procedure tearDown - procedure test_all - end type compute_curl_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, coord_space_wchi - use feign_config_mod, only : feign_finite_element_config - use sci_chi_transform_mod, only : init_chi_transforms + use sci_chi_transform_mod, only: init_chi_transforms implicit none - class(compute_curl_operator_test_type), intent(inout) :: this - - integer(i_def) :: nlayers, dim, diff_dim, dim_w2, diff_dim_w1 - integer(i_def) :: ndf_w0, ndf_w2, ndf_w1 - integer(i_def) :: nqp_h, nqp_v - integer(i_def) :: i - - call feign_finite_element_config( & - cellshape=cellshape_quadrilateral, & - coord_order=0_i_def, & - coord_order_nonprime=1_i_def, & - coord_space=coord_space_wchi, & - 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 config_loader_mod, only: final_configuration use sci_chi_transform_mod, only: final_chi_transforms implicit none - class(compute_curl_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_compute_curl_operator_kernel_mod, & only : compute_curl_operator_code implicit none - class(compute_curl_operator_test_type), intent(inout) :: this - real(kind=r_def), parameter :: tol = 1.0e-6_r_def - real(kind=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 - integer(i_def), parameter :: nlayers = 3 + real(r_def), parameter :: tol = 1.0e-6_r_def + 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 + integer(i_def), parameter :: nlayers = 3 + + 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 - integer :: i, j, k, cell - integer(i_def) :: err + integer :: cell integer(i_def) :: ndf_w0, undf_w0, ndf_w1, ndf_w2, ndf_pid, undf_pid integer(i_def) :: ncells, ncell_3d, nqp_h, nqp_v @@ -164,6 +139,10 @@ contains chi2, & chi3, & panel_id, & + geometry, & + topology, & + coord_system, & + scaled_radius, & ndf_w2, & basis_w2, & ndf_w1, & diff --git a/components/science/unit-test/kernel/fem/compute_derham_matrices_kernel_mod_test.pf b/components/science/unit-test/kernel/fem/compute_derham_matrices_kernel_mod_test.pf index 1674c37e8..31e7f52af 100644 --- a/components/science/unit-test/kernel/fem/compute_derham_matrices_kernel_mod_test.pf +++ b/components/science/unit-test/kernel/fem/compute_derham_matrices_kernel_mod_test.pf @@ -8,77 +8,50 @@ !> module compute_derham_matrices_kernel_mod_test - use constants_mod, only: i_def, r_def, imdi + use constants_mod, only: i_def, r_def, imdi, rmdi + use finite_element_config_mod, only: coord_system_xyz use funit implicit none private - public :: test_all - - @TestCase - type, extends(TestCase), public :: compute_derham_matrices_test_type - private - contains - procedure setUp - procedure tearDown - procedure test_all - end type compute_derham_matrices_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, coord_space_wchi - use feign_config_mod, only : feign_finite_element_config - use sci_chi_transform_mod, only : init_chi_transforms + use sci_chi_transform_mod, only: init_chi_transforms implicit none - class(compute_derham_matrices_test_type), intent(inout) :: this - - call feign_finite_element_config( & - cellshape=cellshape_quadrilateral, & - coord_order=0_i_def, & - coord_order_nonprime=1_i_def, & - coord_space=coord_space_wchi, & - 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 config_loader_mod, only: final_configuration use sci_chi_transform_mod, only: final_chi_transforms implicit none - class(compute_derham_matrices_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_compute_derham_matrices_kernel_mod, only : compute_derham_matrices_code + use sci_compute_derham_matrices_kernel_mod, only: compute_derham_matrices_code implicit none - class(compute_derham_matrices_test_type), intent(inout) :: this - real(r_def), parameter :: tol = 1.0e-12_r_def ! Operators to test @@ -107,6 +80,11 @@ contains integer(i_def) :: i real(kind=r_def) :: 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 + ! Set cell coordinates dx = 2.0_r_def dy = 3.0_r_def @@ -157,6 +135,9 @@ contains 1_i_def, broken_div, & x, y, z, & panel_id, & + geometry, topology, & + coord_system, & + scaled_radius, & ndf, s, ds, & ndf, v, ds, & ndf, v, dv, & diff --git a/components/science/unit-test/kernel/fem/compute_div_operator_kernel_mod_test.pf b/components/science/unit-test/kernel/fem/compute_div_operator_kernel_mod_test.pf index 95e85be0a..9f48a1452 100644 --- a/components/science/unit-test/kernel/fem/compute_div_operator_kernel_mod_test.pf +++ b/components/science/unit-test/kernel/fem/compute_div_operator_kernel_mod_test.pf @@ -8,7 +8,9 @@ !> module compute_div_operator_kernel_mod_test - use constants_mod, only : i_def, r_def, imdi + use constants_mod, only : l_def, i_def, r_def, imdi, rmdi + use finite_element_config_mod, only: coord_system_xyz + use funit use get_unit_test_m3x3_dofmap_mod, & only : get_w0_m3x3_dofmap, & get_w3_m3x3_dofmap @@ -26,82 +28,59 @@ module compute_div_operator_kernel_mod_test get_gaussian_q3x3x3_quadrature_weights_z use get_unit_test_3x3x3_chi_mod, & only : get_w0_3x3x3_field - use funit implicit none private - public :: test_all - - @TestCase - type, extends(TestCase), public :: compute_div_operator_test_type - private - contains - procedure setUp - procedure tearDown - procedure test_all - end type compute_div_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, coord_space_wchi - use feign_config_mod, only : feign_finite_element_config - use sci_chi_transform_mod, only : init_chi_transforms + use sci_chi_transform_mod, only: init_chi_transforms implicit none - class(compute_div_operator_test_type), intent(inout) :: this - - call feign_finite_element_config( & - cellshape=cellshape_quadrilateral, & - coord_order=0_i_def, & - coord_order_nonprime=1_i_def, & - coord_space=coord_space_wchi, & - coord_system=coord_system_xyz, & - element_order_h=1_i_def, & - element_order_v=1_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 config_loader_mod, only: final_configuration use sci_chi_transform_mod, only: final_chi_transforms implicit none - class(compute_div_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_compute_div_operator_kernel_mod, only : compute_div_operator_code implicit none - class(compute_div_operator_test_type), intent(inout) :: this - real(kind=r_def), parameter :: tol = 1.0e-6_r_def real(kind=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 integer(i_def), parameter :: nlayers = 3 - integer :: i, j, k, cell + 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 + logical(l_def), parameter :: rehabilitate = .true. + + integer :: cell integer(i_def) :: ndf_w0, undf_w0, ndf_w2, ndf_w3, ndf_pid, undf_pid integer(i_def) :: nqp_h, nqp_v @@ -165,6 +144,11 @@ contains chi2, & chi3, & panel_id, & + geometry, & + topology, & + coord_system, & + scaled_radius, & + rehabilitate, & ndf_w3, & basis_w3, & ndf_w2, & diff --git a/components/science/unit-test/kernel/fem/compute_grad_operator_kernel_mod_test.pf b/components/science/unit-test/kernel/fem/compute_grad_operator_kernel_mod_test.pf index 3ee9177c2..57c7b4ecc 100644 --- a/components/science/unit-test/kernel/fem/compute_grad_operator_kernel_mod_test.pf +++ b/components/science/unit-test/kernel/fem/compute_grad_operator_kernel_mod_test.pf @@ -8,7 +8,9 @@ !> module compute_grad_operator_kernel_mod_test - use constants_mod, only : i_def, r_def, imdi + use constants_mod, only : i_def, r_def, imdi, rmdi + use finite_element_config_mod, only : coord_system_xyz + use funit use get_unit_test_m3x3_q3x3x3_sizes_mod, only : get_w0_m3x3_q3x3x3_size, & get_w1_m3x3_q3x3x3_size, & get_w3_m3x3_q3x3x3_size @@ -20,91 +22,64 @@ module compute_grad_operator_kernel_mod_test use get_unit_test_q3x3x3_quadrature_mod, only : get_gaussian_q3x3x3_quadrature_weights_xy, & get_gaussian_q3x3x3_quadrature_weights_z use get_unit_test_3x3x3_chi_mod, only : get_w0_3x3x3_field - use funit implicit none private - public :: test_all - - @TestCase - type, extends(TestCase), public :: compute_grad_operator_test_type - private - contains - procedure setUp - procedure tearDown - procedure test_all - end type compute_grad_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, coord_space_wchi - use feign_config_mod, only : feign_finite_element_config - use sci_chi_transform_mod, only : init_chi_transforms + use sci_chi_transform_mod, only: init_chi_transforms implicit none - class(compute_grad_operator_test_type), intent(inout) :: this - - call feign_finite_element_config( & - cellshape=cellshape_quadrilateral, & - coord_order=0_i_def, & - coord_order_nonprime=1_i_def, & - coord_space=coord_space_wchi, & - 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 config_loader_mod, only: final_configuration use sci_chi_transform_mod, only: final_chi_transforms implicit none - class(compute_grad_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_compute_grad_operator_kernel_mod, & only : compute_grad_operator_code implicit none - class(compute_grad_operator_test_type), intent(inout) :: this - real(kind=r_def), dimension(12,8) :: answer real(kind=r_def), parameter :: tol = 1.0e-6_r_def real(kind=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 - integer(i_def) :: nlayers, ncells, ncell_3d, cell, icell + 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 + + integer(i_def) :: nlayers, ncells, ncell_3d, cell integer(i_def) :: ndf_w0, undf_w0, ndf_w1, undf_w1 integer(i_def) :: ndf_pid, undf_pid, ndf_wchi, undf_wchi integer(i_def) :: dim_space, dim_space_diff integer(i_def) :: nqp_h, nqp_v - integer(i_def) :: i, j, k - integer(i_def), allocatable :: map_wchi(:,:), map_pid(:,:) real(r_def), allocatable :: diff_basis_w0(:,:,:,:), & basis_w1(:,:,:,:), & @@ -167,6 +142,10 @@ contains chi2, & chi3, & panel_id, & + geometry, & + topology, & + coord_system, & + scaled_radius, & ndf_w1, & basis_w1, & ndf_w0, & diff --git a/components/science/unit-test/kernel/fem/compute_mass_matrix_kernel_w0_mod_test.pf b/components/science/unit-test/kernel/fem/compute_mass_matrix_kernel_w0_mod_test.pf index ddaaa794d..89795c1d5 100644 --- a/components/science/unit-test/kernel/fem/compute_mass_matrix_kernel_w0_mod_test.pf +++ b/components/science/unit-test/kernel/fem/compute_mass_matrix_kernel_w0_mod_test.pf @@ -8,7 +8,7 @@ !> module compute_mass_matrix_kernel_w0_mod_test - use constants_mod, only : i_def, r_def, r_single, r_double, l_def + use constants_mod, only: i_def, r_def, imdi, rmdi, r_single, r_double, l_def use get_unit_test_m3x3_q3x3x3_sizes_mod, only : get_w0_m3x3_q3x3x3_size, & get_w3_m3x3_q3x3x3_size @@ -22,72 +22,25 @@ module compute_mass_matrix_kernel_w0_mod_test use get_unit_test_3x3x3_chi_mod, only : get_w0_3x3x3_field use funit + use finite_element_config_mod, only: coord_system_xyz + implicit none private public :: test_r_single, test_mixed_precision, test_r_double - @TestCase - type, extends(TestCase), public :: compute_mass_matrix_w0_test_type - private - contains - procedure setUp - procedure tearDown - procedure test_r_single - procedure test_mixed_precision - procedure test_r_double - end type compute_mass_matrix_w0_test_type - contains - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - subroutine setUp( this ) - - use finite_element_config_mod, only : cellshape_quadrilateral, & - coord_system_xyz, coord_space_wchi - use feign_config_mod, only : feign_finite_element_config - - implicit none - - class(compute_mass_matrix_w0_test_type), intent(inout) :: this - - call feign_finite_element_config( & - cellshape = cellshape_quadrilateral, & - coord_order = 0_i_def, & - coord_order_nonprime = 1_i_def, & - coord_space = coord_space_wchi, & - 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(compute_mass_matrix_w0_test_type), intent(inout) :: this - - call final_configuration() - - end subroutine tearDown - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! TEST R_SINGLE !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! @Test - subroutine test_r_single( this ) + subroutine test_r_single() use sci_compute_mass_matrix_kernel_w_scalar_mod, only : & compute_mass_matrix_w_scalar_code implicit none - class(compute_mass_matrix_w0_test_type), intent(inout) :: this - real(r_def), parameter :: dx = 6000.0_r_def, & dy = 1000.0_r_def, & dz = 2000.0_r_def @@ -120,6 +73,10 @@ contains real(r_single), allocatable :: diff_basis_w0_r32(:,:,:,:) real(r_single), allocatable :: wh_r32(:), wv_r32(:) + 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 ! Get pre-set test data (in i_def and r_def) call get_w0_m3x3_q3x3x3_size( ndf_w0, undf_w0, ncells, & @@ -187,6 +144,9 @@ contains chi3_r32, & panel_id_r32, & extended_mesh, & + geometry, topology, & + coord_system, & + scaled_radius, & ndf_w0, & basis_w0_r32, & ndf_w0, & @@ -227,15 +187,13 @@ contains !!!!!!!!!!!!!!!!!!!!!!!!!!! TEST MIXED PRECISION !!!!!!!!!!!!!!!!!!!!!!!!!!!! @Test - subroutine test_mixed_precision( this ) + subroutine test_mixed_precision() use sci_compute_mass_matrix_kernel_w_scalar_mod, only : & compute_mass_matrix_w_scalar_code implicit none - class(compute_mass_matrix_w0_test_type), intent(inout) :: this - real(r_def), parameter :: dx = 6000.0_r_def, & dy = 1000.0_r_def, & dz = 2000.0_r_def @@ -268,6 +226,10 @@ contains real(r_double), allocatable :: diff_basis_w0_r64(:,:,:,:) real(r_double), allocatable :: wh_r64(:), wv_r64(:) + 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 ! Get pre-set test data (in i_def and r_def) call get_w0_m3x3_q3x3x3_size( ndf_w0, undf_w0, ncells, & @@ -335,6 +297,9 @@ contains chi3_r64, & panel_id_r64, & extended_mesh, & + geometry, topology, & + coord_system, & + scaled_radius, & ndf_w0, & basis_w0_r64, & ndf_w0, & @@ -375,15 +340,13 @@ contains !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! TEST R_DOUBLE !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! @Test - subroutine test_r_double( this ) + subroutine test_r_double() use sci_compute_mass_matrix_kernel_w_scalar_mod, only : & compute_mass_matrix_w_scalar_code implicit none - class(compute_mass_matrix_w0_test_type), intent(inout) :: this - real(r_def), parameter :: dx = 6000.0_r_def, & dy = 1000.0_r_def, & dz = 2000.0_r_def @@ -416,6 +379,10 @@ contains real(r_double), allocatable :: diff_basis_w0_r64(:,:,:,:) real(r_double), allocatable :: wh_r64(:), wv_r64(:) + 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 ! Get pre-set test data (in i_def and r_def) call get_w0_m3x3_q3x3x3_size( ndf_w0, undf_w0, ncells, & @@ -483,6 +450,9 @@ contains chi3_r64, & panel_id_r64, & extended_mesh, & + geometry, topology, & + coord_system, & + scaled_radius, & ndf_w0, & basis_w0_r64, & ndf_w0, & diff --git a/components/science/unit-test/kernel/fem/compute_mass_matrix_kernel_w1_mod_test.pf b/components/science/unit-test/kernel/fem/compute_mass_matrix_kernel_w1_mod_test.pf index 204233951..0e03791c7 100644 --- a/components/science/unit-test/kernel/fem/compute_mass_matrix_kernel_w1_mod_test.pf +++ b/components/science/unit-test/kernel/fem/compute_mass_matrix_kernel_w1_mod_test.pf @@ -8,94 +8,73 @@ !> module compute_mass_matrix_kernel_w1_mod_test - use constants_mod, only : i_def, r_def, imdi - use get_unit_test_m3x3_dofmap_mod, & + use constants_mod, only : i_def, r_def, imdi, rmdi + use get_unit_test_m3x3_dofmap_mod, & only : get_w0_m3x3_dofmap, get_w3_m3x3_dofmap - use get_unit_test_m3x3_q3x3x3_sizes_mod, & + use get_unit_test_m3x3_q3x3x3_sizes_mod, & only : get_w0_m3x3_q3x3x3_size, get_w1_m3x3_q3x3x3_size, & get_w3_m3x3_q3x3x3_size - use get_unit_test_q3x3x3_basis_mod, & + use get_unit_test_q3x3x3_basis_mod, & only : get_w0_q3x3x3_diff_basis, get_w1_q3x3x3_basis, & get_w0_q3x3x3_basis - use get_unit_test_q3x3x3_quadrature_mod, & + use get_unit_test_q3x3x3_quadrature_mod, & only : get_gaussian_q3x3x3_quadrature_weights_xy, & get_gaussian_q3x3x3_quadrature_weights_z - use get_unit_test_3x3x3_chi_mod, & + use get_unit_test_3x3x3_chi_mod, & only : get_w0_3x3x3_field - use funit + use funit + + use finite_element_config_mod, only: coord_system_xyz implicit none private - public :: test_all - - @TestCase - type, extends(TestCase), public :: compute_mass_matrix_w1_test_type - private - contains - procedure setUp - procedure tearDown - procedure test_all - end type compute_mass_matrix_w1_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, coord_space_wchi - use feign_config_mod, only : feign_finite_element_config - use sci_chi_transform_mod, only : init_chi_transforms + use sci_chi_transform_mod, only: init_chi_transforms implicit none - class(compute_mass_matrix_w1_test_type), intent(inout) :: this - - call feign_finite_element_config( & - cellshape=cellshape_quadrilateral, & - coord_order=0_i_def, & - coord_order_nonprime=1_i_def, & - coord_space=coord_space_wchi, & - 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 config_loader_mod, only: final_configuration use sci_chi_transform_mod, only: final_chi_transforms implicit none - class(compute_mass_matrix_w1_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_compute_mass_matrix_kernel_w1_mod, only : compute_mass_matrix_w1_code implicit none - class(compute_mass_matrix_w1_test_type), intent(inout) :: this - real(r_def), parameter :: tol = 1.0e-6_r_def 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 + integer(i_def), parameter :: nlayers = 3 + 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 integer :: i, j, k, cell @@ -159,6 +138,10 @@ contains chi2, & chi3, & panel_id, & + geometry, & + topology, & + coord_system, & + scaled_radius, & ndf_w1, & basis_w1, & ndf_w0, & diff --git a/components/science/unit-test/kernel/fem/compute_mass_matrix_kernel_w2_mod_test.pf b/components/science/unit-test/kernel/fem/compute_mass_matrix_kernel_w2_mod_test.pf index 9fe32c78b..f0ceaa119 100644 --- a/components/science/unit-test/kernel/fem/compute_mass_matrix_kernel_w2_mod_test.pf +++ b/components/science/unit-test/kernel/fem/compute_mass_matrix_kernel_w2_mod_test.pf @@ -8,7 +8,8 @@ !> module compute_mass_matrix_kernel_w2_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_dofmap_mod, & only : get_w0_m3x3_dofmap, get_w3_m3x3_dofmap use get_unit_test_m3x3_q3x3x3_sizes_mod, & @@ -24,78 +25,57 @@ module compute_mass_matrix_kernel_w2_mod_test only : get_w0_3x3x3_field use funit + use finite_element_config_mod, only: coord_system_xyz + implicit none private - public :: test_all - - @TestCase - type, extends(TestCase), public :: compute_mass_matrix_w2_test_type - private - contains - procedure setUp - procedure tearDown - procedure test_all - end type compute_mass_matrix_w2_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, coord_space_wchi - use feign_config_mod, only : feign_finite_element_config - use sci_chi_transform_mod, only : init_chi_transforms + use sci_chi_transform_mod, only: init_chi_transforms implicit none - class(compute_mass_matrix_w2_test_type), intent(inout) :: this - - call feign_finite_element_config( & - cellshape=cellshape_quadrilateral, & - coord_order=0_i_def, & - coord_order_nonprime=1_i_def, & - coord_space=coord_space_wchi, & - 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 config_loader_mod, only: final_configuration use sci_chi_transform_mod, only: final_chi_transforms implicit none - class(compute_mass_matrix_w2_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_compute_mass_matrix_kernel_w2_mod, only : compute_mass_matrix_w2_code implicit none - class(compute_mass_matrix_w2_test_type), intent(inout) :: this - real(r_def), parameter :: tol = 1.0e-6_r_def 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 + integer(i_def), parameter :: nlayers = 3 + 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 integer :: i, j, k, cell @@ -160,6 +140,10 @@ contains chi2, & chi3, & panel_id, & + geometry, & + topology, & + coord_system, & + scaled_radius, & ndf_w2, & basis_w2, & ndf_w0, & diff --git a/components/science/unit-test/kernel/fem/compute_mass_matrix_kernel_w2b_mod_test.pf b/components/science/unit-test/kernel/fem/compute_mass_matrix_kernel_w2b_mod_test.pf index fe9c58026..043fdfdf2 100644 --- a/components/science/unit-test/kernel/fem/compute_mass_matrix_kernel_w2b_mod_test.pf +++ b/components/science/unit-test/kernel/fem/compute_mass_matrix_kernel_w2b_mod_test.pf @@ -7,7 +7,8 @@ !> module compute_mass_matrix_kernel_w2b_mod_test - use constants_mod, only : i_def, r_def + 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_w2broken_m3x3_q3x3x3_size, & get_w3_m3x3_q3x3x3_size @@ -21,75 +22,33 @@ module compute_mass_matrix_kernel_w2b_mod_test use get_unit_test_3x3x3_chi_mod, only : get_w0_3x3x3_field use funit + use finite_element_config_mod, only: coord_system_xyz + implicit none private public :: test_all - @TestCase - type, extends(TestCase), public :: compute_mass_matrix_w2b_test_type - private - contains - procedure setUp - procedure tearDown - procedure test_all - end type compute_mass_matrix_w2b_test_type - contains - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - subroutine setUp( this ) - - use finite_element_config_mod, only : cellshape_quadrilateral, & - coord_system_xyz, coord_space_wchi - use feign_config_mod, only : feign_finite_element_config - - implicit none - - class(compute_mass_matrix_w2b_test_type), intent(inout) :: this - - call feign_finite_element_config( & - cellshape=cellshape_quadrilateral, & - coord_order=0_i_def, & - coord_order_nonprime=1_i_def, & - coord_space=coord_space_wchi, & - 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(compute_mass_matrix_w2b_test_type), intent(inout) :: this - - call final_configuration() - - end subroutine tearDown - - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - @Test - subroutine test_all( this ) + subroutine test_all() use sci_compute_mass_matrix_kernel_w2_mod, only : compute_mass_matrix_w2_code implicit none - class(compute_mass_matrix_w2b_test_type), intent(inout) :: this - real(kind=r_def) :: answer(6,6) real(r_def), parameter :: tol = 1.0e-6_r_def real(r_def), parameter :: dx = 6000.0_r_def, & dy = 1000.0_r_def, & dz = 2000.0_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 + integer(i_def) :: nlayers, ncells, ncell_3d, cell, icell integer(i_def) :: ndf_w0, undf_w0, ndf_w2broken, undf_w2broken, ndf_w3, undf_w3 integer(i_def) :: dim_space, dim_space_diff @@ -151,6 +110,10 @@ contains chi2, & chi3, & panel_id, & + geometry, & + topology, & + coord_system, & + scaled_radius, & ndf_w2broken, & basis_w2broken, & ndf_w0, & diff --git a/components/science/unit-test/kernel/fem/compute_mass_matrix_kernel_w3_mod_test.pf b/components/science/unit-test/kernel/fem/compute_mass_matrix_kernel_w3_mod_test.pf index 8a9cedd13..9bf127e7b 100644 --- a/components/science/unit-test/kernel/fem/compute_mass_matrix_kernel_w3_mod_test.pf +++ b/components/science/unit-test/kernel/fem/compute_mass_matrix_kernel_w3_mod_test.pf @@ -8,7 +8,7 @@ !> module compute_mass_matrix_kernel_w3_mod_test - use constants_mod, only : i_def, r_def, r_single, r_double, imdi + use constants_mod, only: i_def, r_def, l_def, r_single, r_double, imdi, rmdi use get_unit_test_m3x3_q3x3x3_sizes_mod, only : get_w0_m3x3_q3x3x3_size, & get_w3_m3x3_q3x3x3_size @@ -22,81 +22,61 @@ module compute_mass_matrix_kernel_w3_mod_test use get_unit_test_3x3x3_chi_mod, only : get_w0_3x3x3_field use funit + use finite_element_config_mod, only: coord_system_xyz + implicit none private + public :: set_up, tear_down public :: test_r_single, test_r_double, test_mixed_precision - @TestCase - type, extends(TestCase), public :: compute_mass_matrix_w3_test_type - private - contains - procedure setUp - procedure tearDown - procedure test_r_single - procedure test_r_double - procedure test_mixed_precision - end type compute_mass_matrix_w3_test_type - contains !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - subroutine setUp( this ) + @before + subroutine set_up() - use finite_element_config_mod, only : cellshape_quadrilateral, & - coord_system_xyz, coord_space_wchi - use feign_config_mod, only : feign_finite_element_config - use sci_chi_transform_mod, only : init_chi_transforms + use sci_chi_transform_mod, only: init_chi_transforms implicit none - class(compute_mass_matrix_w3_test_type), intent(inout) :: this - - call feign_finite_element_config( & - cellshape = cellshape_quadrilateral, & - coord_order = 0_i_def, & - coord_order_nonprime = 1_i_def, & - coord_space = coord_space_wchi, & - 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 config_loader_mod, only: final_configuration use sci_chi_transform_mod, only: final_chi_transforms implicit none - class(compute_mass_matrix_w3_test_type), intent(inout) :: this - - call final_configuration() call final_chi_transforms() - end subroutine tearDown + end subroutine tear_down !!!!!!!!!!!!!!!!!!!!!!!!!!!!!! TEST R_SINGLE !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! @Test - subroutine test_r_single( this ) + subroutine test_r_single() use sci_compute_mass_matrix_kernel_w3_mod, only : compute_mass_matrix_w3_code implicit none - class(compute_mass_matrix_w3_test_type), intent(inout) :: this - real(r_def), parameter :: dx = 6000.0_r_def, & dy = 1000.0_r_def, & dz = 2000.0_r_def integer(i_def), parameter :: nlayers = 3_i_def real(r_single), parameter :: tol = 1.0e-6_r_single + + 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 + logical(l_def), parameter :: rehabilitate = .true. + real(r_single) :: answer integer(i_def) :: ncells, ncell_3d, cell @@ -198,6 +178,10 @@ contains chi2_r32, & chi3_r32, & panel_id_r32, & + geometry, topology, & + coord_system, & + scaled_radius, & + rehabilitate, & ndf_w3, & basis_w3_r32, & ndf_w0, & @@ -234,19 +218,24 @@ contains !!!!!!!!!!!!!!!!!!!!!!!!!!!!!! TEST R_DOUBLE !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! @Test - subroutine test_r_double( this ) + subroutine test_r_double() use sci_compute_mass_matrix_kernel_w3_mod, only : compute_mass_matrix_w3_code implicit none - class(compute_mass_matrix_w3_test_type), intent(inout) :: this - real(r_def), parameter :: dx = 6000.0_r_def, & dy = 1000.0_r_def, & dz = 2000.0_r_def integer(i_def), parameter :: nlayers = 3_i_def real(r_double), parameter :: tol = 1.0e-6_r_double + + 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 + logical(l_def), parameter :: rehabilitate = .true. + real(r_double) :: answer integer(i_def) :: ncells, ncell_3d, cell @@ -348,6 +337,10 @@ contains chi2_r64, & chi3_r64, & panel_id_r64, & + geometry, topology, & + coord_system, & + scaled_radius, & + rehabilitate, & ndf_w3, & basis_w3_r64, & ndf_w0, & @@ -385,19 +378,24 @@ contains !!!!!!!!!!!!!!!!!!!!!!!!!!! TEST MIXED PRECISION !!!!!!!!!!!!!!!!!!!!!!!!!!!! @Test - subroutine test_mixed_precision( this ) + subroutine test_mixed_precision() use sci_compute_mass_matrix_kernel_w3_mod, only : compute_mass_matrix_w3_code implicit none - class(compute_mass_matrix_w3_test_type), intent(inout) :: this - real(r_def), parameter :: dx = 6000.0_r_def, & dy = 1000.0_r_def, & dz = 2000.0_r_def integer(i_def), parameter :: nlayers = 3_i_def real(r_single), parameter :: tol = 1.0e-6_r_single + + 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 + logical(l_def), parameter :: rehabilitate = .true. + real(r_single) :: answer integer(i_def) :: ncells, ncell_3d, cell @@ -499,6 +497,10 @@ contains chi2_r64, & chi3_r64, & panel_id_r64, & + geometry, topology, & + coord_system, & + scaled_radius, & + rehabilitate, & ndf_w3, & basis_w3_r64, & ndf_w0, & diff --git a/components/science/unit-test/kernel/fem/compute_mass_matrix_kernel_wtheta_mod_test.pf b/components/science/unit-test/kernel/fem/compute_mass_matrix_kernel_wtheta_mod_test.pf index 43334a4fe..b09876167 100644 --- a/components/science/unit-test/kernel/fem/compute_mass_matrix_kernel_wtheta_mod_test.pf +++ b/components/science/unit-test/kernel/fem/compute_mass_matrix_kernel_wtheta_mod_test.pf @@ -8,7 +8,8 @@ !> module compute_mass_matrix_kernel_wtheta_mod_test - use constants_mod, only : i_def, r_def, r_single, r_double, l_def, imdi + use constants_mod, only: i_def, r_def, r_single, r_double, & + l_def, imdi, rmdi use get_unit_test_m3x3_q3x3x3_sizes_mod, & only : get_w0_m3x3_q3x3x3_size, & @@ -28,77 +29,50 @@ module compute_mass_matrix_kernel_wtheta_mod_test only : get_w0_3x3x3_field use funit + use finite_element_config_mod, only: coord_system_xyz + implicit none private + public set_up, tear_down public :: test_r_single, test_mixed_precision, test_r_double - @TestCase - type, extends(TestCase), public :: compute_mass_matrix_wtheta_test_type - private - contains - procedure setUp - procedure tearDown - procedure test_r_single - procedure test_mixed_precision - procedure test_r_double - end type compute_mass_matrix_wtheta_test_type - contains !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - subroutine setUp( this ) + @before + subroutine set_up() - use finite_element_config_mod, only : cellshape_quadrilateral, & - coord_system_xyz, coord_space_wchi - use feign_config_mod, only : feign_finite_element_config - use sci_chi_transform_mod, only : init_chi_transforms + use sci_chi_transform_mod, only: init_chi_transforms implicit none - class(compute_mass_matrix_wtheta_test_type), intent(inout) :: this - - call feign_finite_element_config( & - cellshape = cellshape_quadrilateral, & - coord_order = 0_i_def, & - coord_order_nonprime = 1_i_def, & - coord_space = coord_space_wchi, & - 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 config_loader_mod, only: final_configuration use sci_chi_transform_mod, only: final_chi_transforms implicit none - class(compute_mass_matrix_wtheta_test_type), intent(inout) :: this - - call final_configuration() call final_chi_transforms() - end subroutine tearDown + end subroutine tear_down !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! TEST R_SINGLE !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! @Test - subroutine test_r_single( this ) + subroutine test_r_single() use sci_compute_mass_matrix_kernel_w_scalar_mod, only : & compute_mass_matrix_w_scalar_code implicit none - class(compute_mass_matrix_wtheta_test_type), intent(inout) :: this - real(r_def), parameter :: dx = 6000.0_r_def, & dy = 1000.0_r_def, & dz = 2000.0_r_def @@ -118,6 +92,10 @@ contains real(r_single), allocatable :: local_stencil_r32(:,:,:) logical(l_def), parameter :: extended_mesh = .false. + 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 ! Canned test data is kept as r_def real(r_def), allocatable :: chi1(:), chi2(:), chi3(:) @@ -211,6 +189,9 @@ contains chi3_r32, & panel_id_r32, & extended_mesh, & + geometry, topology, & + coord_system, & + scaled_radius, & ndf_wtheta, & basis_wtheta_r32, & ndf_w0, & @@ -255,15 +236,13 @@ contains !!!!!!!!!!!!!!!!!!!!!!!!!!! TEST MIXED PRECISION !!!!!!!!!!!!!!!!!!!!!!!!!!!! @Test - subroutine test_mixed_precision( this ) + subroutine test_mixed_precision() use sci_compute_mass_matrix_kernel_w_scalar_mod, only : & compute_mass_matrix_w_scalar_code implicit none - class(compute_mass_matrix_wtheta_test_type), intent(inout) :: this - real(r_def), parameter :: dx = 6000.0_r_def, & dy = 1000.0_r_def, & dz = 2000.0_r_def @@ -283,6 +262,10 @@ contains real(r_single), allocatable :: local_stencil_r32(:,:,:) logical(l_def), parameter :: extended_mesh = .false. + 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 ! Canned test data is kept as r_def real(r_def), allocatable :: chi1(:), chi2(:), chi3(:) @@ -376,6 +359,9 @@ contains chi3_r64, & panel_id_r64, & extended_mesh, & + geometry, topology, & + coord_system, & + scaled_radius, & ndf_wtheta, & basis_wtheta_r64, & ndf_w0, & @@ -420,15 +406,13 @@ contains !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! TEST R_DOUBLE !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! @Test - subroutine test_r_double( this ) + subroutine test_r_double() use sci_compute_mass_matrix_kernel_w_scalar_mod, only : & compute_mass_matrix_w_scalar_code implicit none - class(compute_mass_matrix_wtheta_test_type), intent(inout) :: this - real(r_def), parameter :: dx = 6000.0_r_def, & dy = 1000.0_r_def, & dz = 2000.0_r_def @@ -448,6 +432,10 @@ contains real(r_double), allocatable :: local_stencil_r64(:,:,:) logical(l_def), parameter :: extended_mesh = .false. + 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 ! Canned test data is kept as r_def real(r_def), allocatable :: chi1(:), chi2(:), chi3(:) @@ -541,6 +529,9 @@ contains chi3_r64, & panel_id_r64, & extended_mesh, & + geometry, topology, & + coord_system, & + scaled_radius, & ndf_wtheta, & basis_wtheta_r64, & ndf_w0, & diff --git a/components/science/unit-test/kernel/fem/gp_rhs_kernel_mod_test.pf b/components/science/unit-test/kernel/fem/gp_rhs_kernel_mod_test.pf index e3e582d9f..8c880313b 100644 --- a/components/science/unit-test/kernel/fem/gp_rhs_kernel_mod_test.pf +++ b/components/science/unit-test/kernel/fem/gp_rhs_kernel_mod_test.pf @@ -7,11 +7,12 @@ !------------------------------------------------------------------------------- module gp_rhs_kernel_mod_test - use constants_mod, only : i_def, r_def, imdi + use constants_mod, only: i_def, r_def, imdi, rmdi + use finite_element_config_mod, only: coord_system_xyz use funit - use quadrature_xyoz_mod, only: quadrature_xyoz_type, & - quadrature_xyoz_proxy_type - use quadrature_rule_gaussian_mod, only: quadrature_rule_gaussian_type + use quadrature_xyoz_mod, only: quadrature_xyoz_type, & + quadrature_xyoz_proxy_type + use quadrature_rule_gaussian_mod, only: quadrature_rule_gaussian_type implicit none @@ -40,33 +41,21 @@ module gp_rhs_kernel_mod_test contains !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + @before subroutine setUp( this ) - use finite_element_config_mod, only : cellshape_quadrilateral, & - coord_system_xyz, coord_space_wchi - use feign_config_mod, only : feign_finite_element_config - use sci_chi_transform_mod, only : init_chi_transforms + use sci_chi_transform_mod, only: init_chi_transforms implicit none class(gp_rhs_test_type), intent(inout) :: this type( quadrature_xyoz_type ) :: qr - type(quadrature_xyoz_proxy_type) :: qr_proxy + type( quadrature_xyoz_proxy_type ) :: qr_proxy type( quadrature_rule_gaussian_type ) :: quadrature_rule integer(i_def) :: nqp_h, nqp_v - call feign_finite_element_config( & - cellshape=cellshape_quadrilateral, & - coord_order=0_i_def, & - coord_order_nonprime=1_i_def, & - coord_space=coord_space_wchi, & - 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) qr = quadrature_xyoz_type(3, quadrature_rule) @@ -75,16 +64,16 @@ contains nqp_h=qr_proxy%np_xy nqp_v=qr_proxy%np_z - allocate( this%basis_f (1,ndf_f,nqp_h,nqp_v), & - this%basis (1,ndf_w0,nqp_h,nqp_v), & + allocate( this%basis_f(1,ndf_f,nqp_h,nqp_v), & + this%basis(1,ndf_w0,nqp_h,nqp_v), & this%diff_basis(3,ndf_w0,nqp_h,nqp_v) ) end subroutine setUp !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + @after subroutine tearDown( this ) - use config_loader_mod, only: final_configuration use sci_chi_transform_mod, only: final_chi_transforms implicit none @@ -93,7 +82,6 @@ contains deallocate( this%basis, this%basis_f, this%diff_basis ) - call final_configuration() call final_chi_transforms() end subroutine tearDown @@ -124,6 +112,11 @@ contains real(kind=r_def) :: fx, fy, fz integer :: df, dim, qp1, qp2 + 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 + qr = quadrature_xyoz_type(3, quadrature_rule) qr_proxy = qr%get_quadrature_proxy() @@ -150,6 +143,7 @@ contains fx = (1.0_r_def - x_vert(df,1)) + (-1.0_r_def)**(int(x_vert(df,1))+1)*xqp_h(qp1,1) fy = (1.0_r_def - x_vert(df,2)) + (-1.0_r_def)**(int(x_vert(df,2))+1)*xqp_h(qp1,2) fz = (1.0_r_def - x_vert(df,3)) + (-1.0_r_def)**(int(x_vert(df,3))+1)*xqp_v(qp2) + this%diff_basis(1,df,qp1,qp2) = fy*fz this%diff_basis(2,df,qp1,qp2) = fx*fz this%diff_basis(3,df,qp1,qp2) = fx*fy @@ -172,15 +166,16 @@ contains num_layers=1 rhs(:) = 0.0_r_def - call gp_rhs_code(num_layers, & - rhs, f_data, & - chi(:,1), chi(:,2), chi(:,3), & - panel_id, & - ndf_w0, undf, map_w0, this%basis, & - ndf_f, undf_f, map_f, this%basis_f, & - ndf_w0, undf, map_w0, & - this%basis, this%diff_basis, & - ndf_pid, undf_pid, map_pid, & + call gp_rhs_code(num_layers, & + rhs, f_data, & + chi(:,1), chi(:,2), chi(:,3), & + panel_id, geometry, topology, & + coord_system, scaled_radius, & + ndf_w0, undf, map_w0, this%basis, & + ndf_f, undf_f, map_f, this%basis_f, & + ndf_w0, undf, map_w0, & + this%basis, this%diff_basis, & + ndf_pid, undf_pid, map_pid, & nqp_h, nqp_v, wh, wv ) ! Answer for hardwired jacobian diff --git a/components/science/unit-test/kernel/fem/gp_vector_rhs_kernel_mod_test.pf b/components/science/unit-test/kernel/fem/gp_vector_rhs_kernel_mod_test.pf index 88cd570cb..54806d520 100644 --- a/components/science/unit-test/kernel/fem/gp_vector_rhs_kernel_mod_test.pf +++ b/components/science/unit-test/kernel/fem/gp_vector_rhs_kernel_mod_test.pf @@ -7,11 +7,14 @@ !------------------------------------------------------------------------------- module gp_vector_rhs_kernel_mod_test - use constants_mod, only : i_def, r_def + use constants_mod, only: i_def, r_def, rmdi + use finite_element_config_mod, only: coord_system_xyz use funit - use quadrature_xyoz_mod, only: quadrature_xyoz_type, & - quadrature_xyoz_proxy_type - use quadrature_rule_gaussian_mod, only: quadrature_rule_gaussian_type + use quadrature_xyoz_mod, only: quadrature_xyoz_type, & + quadrature_xyoz_proxy_type + use quadrature_rule_gaussian_mod, only: quadrature_rule_gaussian_type + + use mesh_mod, only: geometry_planar, topology_periodic implicit none @@ -36,44 +39,22 @@ module gp_vector_rhs_kernel_mod_test contains !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + @before subroutine setUp( this ) - 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, coord_space_wchi - 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(gp_vector_rhs_test_type), intent(inout) :: this type( quadrature_xyoz_type ) :: qr - type(quadrature_xyoz_proxy_type) :: qr_proxy + type( quadrature_xyoz_proxy_type ) :: qr_proxy type( quadrature_rule_gaussian_type ) :: quadrature_rule integer(i_def) :: nqp_h, nqp_v - 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_order_nonprime=1_i_def, & - coord_space=coord_space_wchi, & - 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) + call init_chi_transforms(geometry_planar, topology_periodic) qr = quadrature_xyoz_type(3, quadrature_rule) @@ -88,9 +69,9 @@ contains end subroutine setUp !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + @after subroutine tearDown( this ) - use config_loader_mod, only: final_configuration use sci_chi_transform_mod, only: final_chi_transforms implicit none @@ -99,7 +80,6 @@ contains deallocate( this%basis_f, this%basis, this%diff_basis ) - call final_configuration() call final_chi_transforms() end subroutine tearDown @@ -129,6 +109,10 @@ contains real(kind=r_def) :: fx, fy, fz, dx, dy, dz, t_vec(3) integer :: df, dim, qp1, qp2 + integer(i_def), parameter :: geometry = geometry_planar + integer(i_def), parameter :: topology = topology_periodic + integer(i_def), parameter :: coord_system = coord_system_xyz + real(r_def), parameter :: scaled_radius = rmdi qr = quadrature_xyoz_type(3, quadrature_rule) @@ -193,6 +177,10 @@ contains chi(:,1), chi(:,2), chi(:,3), & panel_id, & f_data, & + geometry, & + topology, & + coord_system, & + scaled_radius, & ndf_w0, undf, map_w0, this%basis, & ndf_f, undf_f,map_f, this%basis_f, & ndf_w0, undf, map_w0, this%basis, & diff --git a/components/science/unit-test/kernel/fem/mg_derham_mat_kernel_mod_test.pf b/components/science/unit-test/kernel/fem/mg_derham_mat_kernel_mod_test.pf index 3ac46475e..49dc19a8b 100644 --- a/components/science/unit-test/kernel/fem/mg_derham_mat_kernel_mod_test.pf +++ b/components/science/unit-test/kernel/fem/mg_derham_mat_kernel_mod_test.pf @@ -8,75 +8,48 @@ !> module mg_derham_mat_kernel_mod_test - use constants_mod, only : i_def, r_def, imdi + use constants_mod, only : i_def, r_def, imdi, rmdi + use finite_element_config_mod, only: coord_system_xyz + use funit implicit none private - public :: test_all - - @TestCase - type, extends(TestCase), public :: mg_derham_mat_test_type - private - contains - procedure setUp - procedure tearDown - procedure test_all - end type mg_derham_mat_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, coord_space_wchi - use feign_config_mod, only : feign_finite_element_config - use sci_chi_transform_mod, only : init_chi_transforms + use sci_chi_transform_mod, only: init_chi_transforms implicit none - class(mg_derham_mat_test_type), intent(inout) :: this - - call feign_finite_element_config( & - cellshape=cellshape_quadrilateral, & - coord_order=0_i_def, & - coord_order_nonprime=1_i_def, & - coord_space=coord_space_wchi, & - 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 ) + subroutine tear_down() - use config_loader_mod, only: final_configuration use sci_chi_transform_mod, only: final_chi_transforms implicit none - class(mg_derham_mat_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_mg_derham_mat_kernel_mod, only : mg_derham_mat_code implicit none - class(mg_derham_mat_test_type), intent(inout) :: this - real(r_def), parameter :: tol = 1.0e-12_r_def ! Operators to test @@ -105,6 +78,11 @@ contains integer(i_def) :: i real(kind=r_def) :: 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 + ! Set cell coordinates dx = 2.0_r_def dy = 3.0_r_def @@ -149,6 +127,9 @@ contains 1_i_def, div, & x, y, z, & panel_id, & + geometry, topology, & + coord_system, & + scaled_radius, & ndf, v, dv, & ndf, s, & ndf, s, & 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 26f584216..1bff56a9b 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 @@ -10,10 +10,9 @@ module chi_transform_mod_test 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 mesh_mod, only: geometry_spherical, & + topology_periodic, & + topology_non_periodic use finite_element_config_mod, only: coord_system_native, & coord_system_xyz @@ -245,7 +244,7 @@ contains select case ( this%source_coord_system ) case ( XYZ ) this%src_coord_system = coord_system_xyz - this%topology = topology_fully_periodic + this%topology = topology_periodic case ( LLH, LLH_rot ) this%src_coord_system = coord_system_native @@ -253,7 +252,7 @@ contains case ( ABH, ABH_stretch_rot ) this%src_coord_system = coord_system_native - this%topology = topology_fully_periodic + this%topology = topology_periodic end select this%geometry = geometry_spherical 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 ae1aa557a..ae63c7a4a 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 @@ -10,8 +10,7 @@ module coordinate_jacobian_alphabetaz_mod_test use funit use constants_mod, only: r_def, i_def - use base_mesh_config_mod, only: geometry_spherical, & - topology_fully_periodic + use mesh_mod, only: geometry_spherical, topology_periodic use finite_element_config_mod, only: coord_system_native implicit none @@ -30,7 +29,7 @@ contains implicit none - call init_chi_transforms(geometry_spherical, topology_fully_periodic) + call init_chi_transforms(geometry_spherical, topology_periodic) end subroutine set_up @@ -78,7 +77,7 @@ contains integer(i_def), parameter :: coord_system = coord_system_native integer(i_def), parameter :: geometry = geometry_spherical - integer(i_def), parameter :: topology = topology_fully_periodic + integer(i_def), parameter :: topology = topology_periodic real(r_def), parameter :: scaled_radius = 1000.0_r_def ! We choose a box centred on alpha = 0, beta = 0 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 e4a2c021f..b223d2a6c 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,10 +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 + use mesh_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 diff --git a/components/science/unit-test/kernel/geometry/coordinate_jacobian_stretched_mod_test.pf b/components/science/unit-test/kernel/geometry/coordinate_jacobian_stretched_mod_test.pf index abcec08a4..bc3ff0898 100644 --- a/components/science/unit-test/kernel/geometry/coordinate_jacobian_stretched_mod_test.pf +++ b/components/science/unit-test/kernel/geometry/coordinate_jacobian_stretched_mod_test.pf @@ -21,7 +21,7 @@ contains @before subroutine set_up() - use base_mesh_config_mod, only : geometry_spherical, & + use mesh_mod, only : geometry_spherical, & topology_non_periodic use extrusion_config_mod, only : method_uniform, & stretching_method_linear diff --git a/components/science/unit-test/kernel/geometry/coordinate_jacobian_xyz_mod_test.pf b/components/science/unit-test/kernel/geometry/coordinate_jacobian_xyz_mod_test.pf index bfc90c002..7f9611d9a 100644 --- a/components/science/unit-test/kernel/geometry/coordinate_jacobian_xyz_mod_test.pf +++ b/components/science/unit-test/kernel/geometry/coordinate_jacobian_xyz_mod_test.pf @@ -64,8 +64,8 @@ contains pointwise_coordinate_jacobian_inverse use finite_element_config_mod, only: coord_system_xyz - use base_mesh_config_mod, only: geometry_planar, & - topology_non_periodic + + use mesh_mod, only: geometry_planar, topology_non_periodic implicit none diff --git a/components/science/unit-test/kernel/geometry/get_dz_w3_kernel_mod_test.pf b/components/science/unit-test/kernel/geometry/get_dz_w3_kernel_mod_test.pf index f3644f666..b1bee5170 100644 --- a/components/science/unit-test/kernel/geometry/get_dz_w3_kernel_mod_test.pf +++ b/components/science/unit-test/kernel/geometry/get_dz_w3_kernel_mod_test.pf @@ -35,8 +35,8 @@ contains !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine setUp( this ) - use base_mesh_config_mod, only : geometry_planar, & - topology_fully_periodic + use mesh_mod, only : geometry_planar, & + topology_periodic use finite_element_config_mod, only : cellshape_quadrilateral, & coord_system_xyz, coord_space_wchi use feign_config_mod, only : feign_finite_element_config, & @@ -50,7 +50,7 @@ contains prime_mesh_name='unit_test', & geometry=geometry_planar, & prepartitioned=.false., & - topology=topology_fully_periodic, & + topology=topology_periodic, & fplane=.false., f_lat_deg=45.0_r_def ) call feign_finite_element_config( & diff --git a/components/science/unit-test/kernel/geometry/height_continuous_kernel_mod_test.pf b/components/science/unit-test/kernel/geometry/height_continuous_kernel_mod_test.pf index b3e3d31e6..9e41d3667 100644 --- a/components/science/unit-test/kernel/geometry/height_continuous_kernel_mod_test.pf +++ b/components/science/unit-test/kernel/geometry/height_continuous_kernel_mod_test.pf @@ -7,7 +7,7 @@ !! This test is parametrised to test both planar and spherical geometries. module height_continuous_kernel_mod_test - use base_mesh_config_mod, only : geometry_planar, & + use mesh_mod, only : geometry_planar, & geometry_spherical use constants_mod, only : i_def, r_def, str_long use funit diff --git a/components/science/unit-test/kernel/geometry/height_discontinuous_kernel_mod_test.pf b/components/science/unit-test/kernel/geometry/height_discontinuous_kernel_mod_test.pf index 400f083eb..a037734d0 100644 --- a/components/science/unit-test/kernel/geometry/height_discontinuous_kernel_mod_test.pf +++ b/components/science/unit-test/kernel/geometry/height_discontinuous_kernel_mod_test.pf @@ -20,32 +20,24 @@ module height_discontinuous_kernel_mod_test use get_unit_test_3x3x3_chi_mod, only : get_w0_3x3x3_field use get_unit_test_w3nodal_basis_mod, only : get_w0_w3nodal_basis use get_unit_test_wthetanodal_basis_mod, only : get_w0_wthetanodal_basis + use mesh_mod, only : geometry_planar + implicit none private public :: test_all - @TestCase - type, extends(TestCase), public :: height_discontinuous_kernel_test_type - private - contains - procedure test_all - end type height_discontinuous_kernel_test_type - contains !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! @Test - subroutine test_all( this ) + subroutine test_all() - use base_mesh_config_mod, only : geometry_planar use finite_element_config_mod, only : coord_system_xyz use sci_height_discontinuous_kernel_mod, only : height_discontinuous_code implicit none - class(height_discontinuous_kernel_test_type), intent(inout) :: this - real(r_def), parameter :: dx = 4.0_r_def, & dy = 3.0_r_def, & dz = 2.0_r_def, & 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 162bc1642..c1d169ea5 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 @@ -9,8 +9,8 @@ module native_jacobian_alphabetaz_mod_test use funit use constants_mod, only: r_def, i_def - use base_mesh_config_mod, only: geometry_spherical, & - topology_fully_periodic + use mesh_mod, only: geometry_spherical, & + topology_periodic use finite_element_config_mod, only: coord_system_native implicit none @@ -29,7 +29,7 @@ contains implicit none - call init_chi_transforms(geometry_spherical, topology_fully_periodic) + call init_chi_transforms(geometry_spherical, topology_periodic) end subroutine set_up @@ -73,7 +73,7 @@ contains integer(i_def), parameter :: coord_system = coord_system_native integer(i_def), parameter :: geometry = geometry_spherical - integer(i_def), parameter :: topology = topology_fully_periodic + integer(i_def), parameter :: topology = topology_periodic real(r_def), parameter :: scaled_radius = 1000.0_r_def ! We choose a box centred on alpha = 0, beta = 0 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 53ed27607..78fdc546c 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 @@ -10,10 +10,9 @@ module native_jacobian_lonlatz_mod_test use funit use constants_mod, only : r_def, i_def, PI - use base_mesh_config_mod, only: geometry_spherical, & + use mesh_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 diff --git a/components/science/unit-test/kernel/geometry/native_jacobian_stretched_mod_test.pf b/components/science/unit-test/kernel/geometry/native_jacobian_stretched_mod_test.pf index f90c9b81b..589b23af9 100644 --- a/components/science/unit-test/kernel/geometry/native_jacobian_stretched_mod_test.pf +++ b/components/science/unit-test/kernel/geometry/native_jacobian_stretched_mod_test.pf @@ -21,7 +21,7 @@ contains @before subroutine set_up() - use base_mesh_config_mod, only : geometry_spherical, & + use mesh_mod, only : geometry_spherical, & topology_non_periodic use extrusion_config_mod, only : method_uniform, & stretching_method_linear 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 5686d794b..67bb774bd 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 @@ -8,14 +8,11 @@ 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, & + use mesh_mod, only: geometry_spherical, & topology_non_periodic use finite_element_config_mod, only: coord_system_native - use funit - - implicit none private 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 53fe6db22..e31e944ae 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 @@ -6,15 +6,12 @@ module compute_sample_u_ops_kernel_mod_test - use base_mesh_config_mod, only : geometry_spherical, & + use mesh_mod, only : geometry_spherical, & topology_non_periodic 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 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 571181d93..e33781dad 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 @@ -9,10 +9,9 @@ module convert_phys_to_hdiv_kernel_mod_test use constants_mod, only: i_def, r_def, rmdi - use base_mesh_config_mod, only: geometry_planar, & - topology_fully_periodic + use mesh_mod, only: geometry_planar, & + topology_periodic use finite_element_config_mod, only: coord_system_xyz - use funit implicit none @@ -30,7 +29,7 @@ contains implicit none - call init_chi_transforms(geometry_planar, topology_fully_periodic) + call init_chi_transforms(geometry_planar, topology_periodic) end subroutine set_up @@ -77,7 +76,7 @@ contains 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 :: topology = topology_periodic integer(i_def), parameter :: coord_system = coord_system_xyz real(r_def), parameter :: scaled_radius = rmdi 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 5f1d0f284..db35d7fb6 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 @@ -23,8 +23,8 @@ 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 mesh_mod, only: geometry_planar, & + topology_periodic use finite_element_config_mod, only: coord_system_xyz implicit none @@ -42,7 +42,7 @@ contains implicit none - call init_chi_transforms(geometry_planar,topology_fully_periodic) + call init_chi_transforms(geometry_planar, topology_periodic) end subroutine set_up @@ -69,7 +69,7 @@ contains 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 :: topology = topology_periodic integer(i_def), parameter :: coord_system = coord_system_xyz real(r_def), parameter :: scaled_radius = rmdi diff --git a/components/science/unit-test/kernel/inter_function_space/sample_w3_to_wtheta_kernel_mod_test.pf b/components/science/unit-test/kernel/inter_function_space/sample_w3_to_wtheta_kernel_mod_test.pf index 9f87b1d06..af026fc74 100644 --- a/components/science/unit-test/kernel/inter_function_space/sample_w3_to_wtheta_kernel_mod_test.pf +++ b/components/science/unit-test/kernel/inter_function_space/sample_w3_to_wtheta_kernel_mod_test.pf @@ -34,10 +34,9 @@ contains !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine setUp( this ) - use base_mesh_config_mod, only : geometry_planar, & - topology_fully_periodic + use mesh_mod, only : geometry_planar, & + topology_periodic use feign_config_mod, only : feign_base_mesh_config - implicit none class(sample_w3_to_wtheta_test_type), intent(inout) :: this @@ -46,7 +45,7 @@ contains prime_mesh_name='unit_test', & geometry=geometry_planar, & prepartitioned=.false., & - topology=topology_fully_periodic, & + topology=topology_periodic, & fplane=.false., f_lat_deg=0.0_r_def ) end subroutine setUp diff --git a/components/science/unit-test/kernel/inter_function_space/sample_wtheta_to_w3_kernel_mod_test.pf b/components/science/unit-test/kernel/inter_function_space/sample_wtheta_to_w3_kernel_mod_test.pf index 1f473d77a..7cb8dc1de 100644 --- a/components/science/unit-test/kernel/inter_function_space/sample_wtheta_to_w3_kernel_mod_test.pf +++ b/components/science/unit-test/kernel/inter_function_space/sample_wtheta_to_w3_kernel_mod_test.pf @@ -32,10 +32,9 @@ contains !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine setUp( this ) - use base_mesh_config_mod, only : geometry_planar, & - topology_fully_periodic + use mesh_mod, only : geometry_planar, & + topology_periodic use feign_config_mod, only : feign_base_mesh_config - implicit none class(sample_wtheta_to_w3_test_type), intent(inout) :: this @@ -44,7 +43,7 @@ contains prime_mesh_name='unit_test', & geometry=geometry_planar, & prepartitioned=.false., & - topology=topology_fully_periodic, & + topology=topology_periodic, & fplane=.false., f_lat_deg=0.0_r_def ) end subroutine setUp 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 6d01ea321..2910ef244 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 @@ -10,13 +10,11 @@ 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 base_mesh_config_mod, only: geometry_spherical, & - topology_fully_periodic + use mesh_mod, only: geometry_spherical, & + topology_periodic use finite_element_config_mod, only: coord_system_native - use funit - implicit none private @@ -32,7 +30,7 @@ contains implicit none - call init_chi_transforms(geometry_spherical, topology_fully_periodic) + call init_chi_transforms(geometry_spherical, topology_periodic) end subroutine set_up @@ -66,7 +64,7 @@ contains 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 :: topology = topology_periodic integer(i_def), parameter :: coord_system = coord_system_native real(r_def), parameter :: scaled_radius = 1900000.0_r_def diff --git a/components/science/unit-test/kernel/inter_mesh/inject_sh_w3_to_wt_kernel_mod_test.pf b/components/science/unit-test/kernel/inter_mesh/inject_sh_w3_to_wt_kernel_mod_test.pf index 24cad0316..1b0875c44 100644 --- a/components/science/unit-test/kernel/inter_mesh/inject_sh_w3_to_wt_kernel_mod_test.pf +++ b/components/science/unit-test/kernel/inter_mesh/inject_sh_w3_to_wt_kernel_mod_test.pf @@ -32,8 +32,8 @@ contains !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine setUp( this ) - use base_mesh_config_mod, only : geometry_planar, & - topology_fully_periodic + use mesh_mod, only : geometry_planar, & + topology_periodic use feign_config_mod, only : feign_base_mesh_config implicit none @@ -44,7 +44,7 @@ contains prime_mesh_name='unit_test', & geometry=geometry_planar, & prepartitioned=.false., & - topology=topology_fully_periodic, & + topology=topology_periodic, & fplane=.false., f_lat_deg=0.0_r_def ) end subroutine setUp diff --git a/components/science/unit-test/kernel/inter_mesh/inject_wt_to_sh_w3_kernel_mod_test.pf b/components/science/unit-test/kernel/inter_mesh/inject_wt_to_sh_w3_kernel_mod_test.pf index f5467afc2..f75901ef3 100644 --- a/components/science/unit-test/kernel/inter_mesh/inject_wt_to_sh_w3_kernel_mod_test.pf +++ b/components/science/unit-test/kernel/inter_mesh/inject_wt_to_sh_w3_kernel_mod_test.pf @@ -32,8 +32,8 @@ contains !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine setUp( this ) - use base_mesh_config_mod, only : geometry_planar, & - topology_fully_periodic + use mesh_mod, only : geometry_planar, & + topology_periodic use feign_config_mod, only : feign_base_mesh_config implicit none @@ -44,7 +44,7 @@ contains prime_mesh_name='unit_test', & geometry=geometry_planar, & prepartitioned=.false., & - topology=topology_fully_periodic, & + topology=topology_periodic, & fplane=.false., f_lat_deg=0.0_r_def ) end subroutine setUp diff --git a/infrastructure/Makefile b/infrastructure/Makefile index 34a5a9cc9..163d5b302 100644 --- a/infrastructure/Makefile +++ b/infrastructure/Makefile @@ -92,6 +92,7 @@ document-api: api-documentation # Unit tests # unit-tests/%: export BIN_DIR ?= $(PROJECT_DIR)/test +unit-tests/%: export PRE_PROCESS_MACROS += INFRASTRUCTURE_UNIT_TEST unit-tests/%: export PRE_PROCESS_INCLUDE_DIRS = $(realpath unit-test/include) unit-tests/%: export PROGRAMS = infrastructure_unit_tests unit-tests/%: export PROJECT = infrastructure diff --git a/infrastructure/documentation/uml/mesh/local_mesh_mod.iuml b/infrastructure/documentation/uml/mesh/local_mesh_mod.iuml index d96f4c3b5..09fbca850 100644 --- a/infrastructure/documentation/uml/mesh/local_mesh_mod.iuml +++ b/infrastructure/documentation/uml/mesh/local_mesh_mod.iuml @@ -5,9 +5,9 @@ class local_mesh_mod::local_mesh_type { -mesh_name: string - -geometry: integer - -topology: integer - -coord_sys: integer + -mesh_geometry: integer + -mesh_topology: integer + -mesh_coord_sys: integer -coord_units_xy: string[2] -north_pole: real[2] -null_island: real[2] @@ -74,6 +74,10 @@ class local_mesh_mod::local_mesh_type { +is_coord_sys_xyz(): logical +is_coord_sys_ll(): logical + +geometry(): integer + +topology(): integer + +coord_sys(): integer + +get_north_pole(): real[2] +get_null_island(): real[2] +get_equatorial_latitude(): real diff --git a/infrastructure/documentation/uml/mesh/mesh_mod.iuml b/infrastructure/documentation/uml/mesh/mesh_mod.iuml index f22601e9f..10710f066 100644 --- a/infrastructure/documentation/uml/mesh/mesh_mod.iuml +++ b/infrastructure/documentation/uml/mesh/mesh_mod.iuml @@ -119,6 +119,10 @@ class mesh_mod::mesh_type { +get_colours( <>ncolours: integer, <>ncells_per_colour: integer[:] <>, <>colour_map: integer[::] <> ) +get_colour_map() : integer <> +is_coloured() : logical + + +geometry(): integer + +topology(): integer + +coord_sys(): integer +is_geometry_spherical(): logical +is_geometry_planar(): logical +is_topology_non_periodic(): logical diff --git a/infrastructure/source/mesh/local_mesh_mod.f90 b/infrastructure/source/mesh/local_mesh_mod.F90 similarity index 95% rename from infrastructure/source/mesh/local_mesh_mod.f90 rename to infrastructure/source/mesh/local_mesh_mod.F90 index 481d8ef3f..c7c121b4e 100644 --- a/infrastructure/source/mesh/local_mesh_mod.f90 +++ b/infrastructure/source/mesh/local_mesh_mod.F90 @@ -30,19 +30,33 @@ module local_mesh_mod LOG_LEVEL_INFO, LOG_LEVEL_DEBUG use partition_mod, only: partition_type +#if !defined(INFRASTRUCTURE_UNIT_TEST) && !defined(MESH_TOOLS) + use base_mesh_config_mod, only: & + config_geometry_spherical => geometry_spherical, & + config_geometry_planar => geometry_planar, & + config_topology_periodic => topology_fully_periodic, & + config_topology_non_periodic => topology_non_periodic +#endif + implicit none private - integer(i_def), parameter :: spherical_domain = 601 - integer(i_def), parameter :: planar_domain = 602 - - integer(i_def), parameter :: non_periodic_domain = 701 - integer(i_def), parameter :: channel_domain = 702 - integer(i_def), parameter :: periodic_domain = 703 - - integer(i_def), parameter :: lon_lat_coords = 801 - integer(i_def), parameter :: xyz_coords = 802 +#if !defined(INFRASTRUCTURE_UNIT_TEST) && !defined(MESH_TOOLS) + integer(i_def), parameter, public :: geometry_spherical = config_geometry_spherical ! 101 + integer(i_def), parameter, public :: geometry_planar = config_geometry_planar ! 202 + integer(i_def), parameter, public :: topology_non_periodic = config_topology_non_periodic ! 301 + integer(i_def), parameter, public :: topology_periodic = config_topology_periodic ! 503 +#else + integer(i_def), parameter, public :: geometry_spherical = 101_i_def + integer(i_def), parameter, public :: geometry_planar = 202_i_def + integer(i_def), parameter, public :: topology_non_periodic = 301_i_def + integer(i_def), parameter, public :: topology_periodic = 503_i_def +#endif + integer(i_def), parameter, public :: topology_channel = 402_i_def + + integer(i_def), parameter, public :: coord_sys_ll = 601_i_def + integer(i_def), parameter, public :: coord_sys_xyz = 702_i_def type, extends(linked_list_data_type), public :: local_mesh_type @@ -53,11 +67,11 @@ module local_mesh_mod ! Tag name of mesh. character(str_def) :: mesh_name ! Domain surface geometry. - integer(i_def) :: geometry = emdi + integer(i_def) :: mesh_geometry = emdi ! Domain boundaries topology. - integer(i_def) :: topology = emdi + integer(i_def) :: mesh_topology = emdi ! Co-ordinate system used to specify node locations. - integer(i_def) :: coord_sys = emdi + integer(i_def) :: mesh_coord_sys = emdi ! Co-ordinate units along xy-axes. character(str_def) :: coord_units_xy(2) = cmdi ! Marker id for cells that do not exist for mesh. @@ -226,7 +240,9 @@ module local_mesh_mod procedure, public :: get_north_pole procedure, public :: get_null_island procedure, public :: get_equatorial_latitude - + procedure, public :: geometry + procedure, public :: topology + procedure, public :: coord_sys procedure, public :: get_global_domain_extents procedure, public :: is_geometry_spherical @@ -299,23 +315,23 @@ subroutine initialise_full ( self, & ! Inherit mesh properties from the parent global mesh. if (global_mesh%is_geometry_spherical()) then - self%geometry = spherical_domain + self%mesh_geometry = geometry_spherical else if (global_mesh%is_geometry_planar()) then - self%geometry = planar_domain + self%mesh_geometry = geometry_planar end if if (global_mesh%is_topology_non_periodic()) then - self%topology = non_periodic_domain + self%mesh_topology = topology_non_periodic else if (global_mesh%is_topology_channel()) then - self%topology = channel_domain + self%mesh_topology = topology_channel else if (global_mesh%is_topology_periodic()) then - self%topology = periodic_domain + self%mesh_topology = topology_periodic end if if (global_mesh%is_coord_sys_xyz()) then - self%coord_sys = xyz_coords + self%mesh_coord_sys = coord_sys_xyz else if (global_mesh%is_coord_sys_ll()) then - self%coord_sys = lon_lat_coords + self%mesh_coord_sys = coord_sys_ll end if self%domain_extents = global_mesh%get_domain_extents() @@ -705,9 +721,9 @@ subroutine initialise_lbc ( self, & self%nverts_per_cell = local_lam_mesh%get_nverts_per_cell() self%nverts_per_edge = local_lam_mesh%get_nverts_per_edge() - self%geometry = local_lam_mesh%geometry - self%coord_sys = local_lam_mesh%coord_sys - self%topology = non_periodic_domain + self%mesh_geometry = local_lam_mesh%geometry() + self%mesh_coord_sys = local_lam_mesh%mesh_coord_sys + self%mesh_topology = topology_non_periodic self%npanels = 1_i_def self%max_stencil_depth = 0_i_def @@ -1148,17 +1164,19 @@ subroutine initialise_from_ugrid_data(self, ugrid_mesh_data) integer(i_def) :: max_face_per_node ! Only needed for global meshes that ! are to be partitioned. - logical(i_def) :: periodic_xy(2) = .false. + logical(l_def) :: periodic_xy(2) - character(str_def) :: geometry_str - character(str_def) :: topology_str - character(str_def) :: coord_sys_str + character(str_def) :: geometry_str + character(str_def) :: topology_str + character(str_def) :: coord_sys_str if (.not. ugrid_mesh_data%is_local()) then call log_event( 'Insufficient data to initialise local mesh', & LOG_LEVEL_ERROR ) end if + periodic_xy(2) = .false. + local_mesh_id_counter = local_mesh_id_counter + 1 call self%set_id( local_mesh_id_counter ) @@ -1195,14 +1213,14 @@ subroutine initialise_from_ugrid_data(self, ugrid_mesh_data) select case (trim(geometry_str)) case ('spherical') - self%geometry = spherical_domain + self%mesh_geometry = geometry_spherical case ('planar') - self%geometry = planar_domain + self%mesh_geometry = geometry_planar end select select case (trim(coord_sys_str)) case ('ll') - self%coord_sys=lon_lat_coords + self%mesh_coord_sys = coord_sys_ll ! Ensure units are in radians if ( (trim(self%coord_units_xy(1)) == 'degrees_east') .and. & @@ -1219,16 +1237,16 @@ subroutine initialise_from_ugrid_data(self, ugrid_mesh_data) end if case ('xyz') - self%coord_sys = xyz_coords + self%mesh_coord_sys = coord_sys_xyz end select select case (trim(topology_str)) case ('channel') - self%topology = channel_domain + self%mesh_topology = topology_channel case ('non_periodic') - self%topology = non_periodic_domain + self%mesh_topology = topology_non_periodic case ('periodic') - self%topology = periodic_domain + self%mesh_topology = topology_periodic end select @@ -1295,9 +1313,9 @@ subroutine initialise_unit_test ( self ) self%mesh_name = 'unit_test' - self%geometry = planar_domain - self%topology = periodic_domain - self%coord_sys = xyz_coords + self%mesh_geometry = geometry_planar + self%mesh_topology = topology_periodic + self%mesh_coord_sys = coord_sys_xyz self%void_cell = -9999_i_def local_mesh_id_counter = local_mesh_id_counter + 1 @@ -1632,7 +1650,7 @@ function is_geometry_spherical( self ) result ( answer ) logical (l_def) :: answer - answer = ( self%geometry == spherical_domain ) + answer = ( self%mesh_geometry == geometry_spherical ) end function is_geometry_spherical @@ -1650,7 +1668,7 @@ function is_geometry_planar( self ) result ( answer ) logical (l_def) :: answer - answer = ( self%geometry == planar_domain ) + answer = ( self%mesh_geometry == geometry_planar ) end function is_geometry_planar @@ -1669,7 +1687,7 @@ function is_topology_non_periodic( self ) result ( answer ) logical (l_def) :: answer - answer = ( self%topology == non_periodic_domain ) + answer = ( self%mesh_topology == topology_non_periodic ) end function is_topology_non_periodic @@ -1688,7 +1706,7 @@ function is_topology_channel( self ) result ( answer ) logical (l_def) :: answer - answer = ( self%topology == channel_domain ) + answer = ( self%mesh_topology == topology_channel ) end function is_topology_channel @@ -1707,7 +1725,7 @@ function is_topology_periodic( self ) result ( answer ) logical (l_def) :: answer - answer = ( self%topology == periodic_domain ) + answer = ( self%mesh_topology == topology_periodic ) end function is_topology_periodic @@ -1726,7 +1744,7 @@ function is_coord_sys_xyz( self ) result ( answer ) logical (l_def) :: answer - answer = ( self%coord_sys == xyz_coords ) + answer = ( self%mesh_coord_sys == coord_sys_xyz ) end function is_coord_sys_xyz @@ -1745,7 +1763,7 @@ function is_coord_sys_ll( self ) result ( answer ) logical (l_def) :: answer - answer = ( self%coord_sys == lon_lat_coords ) + answer = ( self%mesh_coord_sys == coord_sys_ll ) end function is_coord_sys_ll @@ -2505,7 +2523,54 @@ function get_mesh_maps( self ) result( local_mesh_maps ) end function get_mesh_maps + !============================================================================== + !> @brief Returns mesh geometry enumeration + !> @return geometry_enumeration Integer enumeration identifying the mesh + !> surface geometry type + !> + function geometry( self ) result( geometry_enumeration ) + + implicit none + + class(local_mesh_type), intent(in) :: self + integer(i_def) :: geometry_enumeration + + geometry_enumeration = self%mesh_geometry + + end function geometry + + !============================================================================== + !> @brief Returns mesh topology enumeration + !> @return topology_enumeration Integer enumeration identifying the mesh + !> domain boundary connectivity type + !> + function topology( self ) result( topology_enumeration ) + + implicit none + + class(local_mesh_type), intent(in) :: self + integer(i_def) :: topology_enumeration + + topology_enumeration = self%mesh_topology + + end function topology + + + !============================================================================== + !> @brief Returns mesh topology enumeration + !> @return coord_sys_enumeration Integer enumeration identifying the mesh + !> coordinate system. + !> + function coord_sys( self ) result( coord_sys_enumeration ) + + implicit none + + class(local_mesh_type), intent(in) :: self + integer(i_def) :: coord_sys_enumeration + + coord_sys_enumeration = self%mesh_coord_sys + end function coord_sys !============================================================================== !> @brief Populates a object with this local mesh object's diff --git a/infrastructure/source/mesh/mesh_mod.F90 b/infrastructure/source/mesh/mesh_mod.F90 index 58aecbdf8..492a678d5 100644 --- a/infrastructure/source/mesh/mesh_mod.F90 +++ b/infrastructure/source/mesh/mesh_mod.F90 @@ -23,7 +23,14 @@ module mesh_mod linked_list_item_type use linked_list_data_mod, only : linked_list_data_type use local_mesh_map_mod, only : local_mesh_map_type - use local_mesh_mod, only : local_mesh_type + use local_mesh_mod, only : local_mesh_type, & + geometry_spherical, & + geometry_planar, & + topology_non_periodic, & + topology_channel, & + topology_periodic, & + coord_sys_ll, & + coord_sys_xyz use log_mod, only : log_event, log_scratch_space, & LOG_LEVEL_ERROR, LOG_LEVEL_TRACE, & LOG_LEVEL_INFO, LOG_LEVEL_DEBUG @@ -44,6 +51,10 @@ module mesh_mod private + public :: geometry_spherical, geometry_planar + public :: topology_non_periodic, topology_channel, topology_periodic + public :: coord_sys_ll, coord_sys_xyz + !============================================================================ ! Declare type definitions in this module !============================================================================ @@ -174,9 +185,9 @@ module mesh_mod integer(i_def), allocatable, private :: cells_in_colour(:,:) !> integer 2-d array, how many of the first so many cells belong to each colour integer(i_def), allocatable, private :: ncells_per_colour_subset(:,:) - integer(i_def),allocatable :: last_inner_cell_per_colour(:,:) - integer(i_def),allocatable :: last_halo_cell_per_colour(:,:) - integer(i_def),allocatable :: last_edge_cell_per_colour(:) + integer(i_def), allocatable :: last_inner_cell_per_colour(:,:) + integer(i_def), allocatable :: last_halo_cell_per_colour(:,:) + integer(i_def), allocatable :: last_edge_cell_per_colour(:) !========================================================================== ! Maps that this mesh connects to ! @@ -233,6 +244,9 @@ module mesh_mod procedure, public :: get_domain procedure, public :: get_domain_top procedure, public :: get_extrusion_id + procedure, public :: geometry + procedure, public :: topology + procedure, public :: coord_sys procedure, public :: get_dz procedure, public :: get_eta procedure, public :: get_vertex_cell_owner @@ -1215,6 +1229,57 @@ function get_extrusion_id(self) result (extrusion_id) end function get_extrusion_id + + !============================================================================== + !> @brief Returns mesh geometry enumeration + !> @return geometry_enumeration Integer enumeration identifying the mesh + !> surface geometry type + !> + function geometry( self ) result( geometry_enumeration ) + + implicit none + + class(mesh_type), intent(in) :: self + integer(i_def) :: geometry_enumeration + + geometry_enumeration = self%local_mesh%geometry() + + end function geometry + + + !============================================================================== + !> @brief Returns mesh topology enumeration + !> @return topology_enumeration Integer enumeration identifying the mesh + !> domain boundary connectivity type + !> + function topology( self ) result( topology_enumeration ) + + implicit none + + class(mesh_type), intent(in) :: self + integer(i_def) :: topology_enumeration + + topology_enumeration = self%local_mesh%topology() + + end function topology + + + !> @brief Returns mesh topology enumeration + !> @return coord_sys_enumeration Integer enumeration identifying the mesh + !> coordinate system. + !> + function coord_sys( self ) result( coord_sys_enumeration ) + + implicit none + + class(mesh_type), intent(in) :: self + integer(i_def) :: coord_sys_enumeration + + coord_sys_enumeration = self%local_mesh%coord_sys() + + end function coord_sys + + !> @details This functions returns an array of 3d-layer thicknesses in !> metres !> @param[out] dz Vertical thickness of layers in [m], array of diff --git a/infrastructure/source/mesh/panel_decomposition_mod.f90 b/infrastructure/source/mesh/panel_decomposition_mod.f90 index b6713da9c..500d0a416 100644 --- a/infrastructure/source/mesh/panel_decomposition_mod.f90 +++ b/infrastructure/source/mesh/panel_decomposition_mod.f90 @@ -12,7 +12,7 @@ module panel_decomposition_mod use global_mesh_collection_mod, only: global_mesh_collection_type use constants_mod, only: i_def, l_def, r_def use log_mod, only: log_event, log_scratch_space, & - LOG_LEVEL_ERROR, LOG_LEVEL_INFO, LOG_LEVEL_DEBUG + LOG_LEVEL_ERROR, LOG_LEVEL_DEBUG implicit none @@ -166,7 +166,7 @@ subroutine get_custom_partition( self, & integer(i_def) :: num_xprocs, num_yprocs - call log_event("Using custom decomposition", LOG_LEVEL_INFO) + call log_event("Using custom decomposition", LOG_LEVEL_DEBUG) num_xprocs = self%num_xprocs num_yprocs = self%num_yprocs @@ -279,7 +279,7 @@ subroutine get_auto_partition( self, & integer(i_def) :: start_xprocs, start_width, i logical :: found_partition - call log_event("Using auto decomposition", LOG_LEVEL_INFO) + call log_event("Using auto decomposition", LOG_LEVEL_DEBUG) ! For automatic partitioning, try to partition into the squarest ! possible partitions. @@ -421,7 +421,7 @@ subroutine get_row_partition( self, & integer(i_def) :: num_xprocs, num_yprocs - call log_event("Using row decomposition", LOG_LEVEL_INFO) + call log_event("Using row decomposition", LOG_LEVEL_DEBUG) num_xprocs = panel_ranks num_yprocs = 1_i_def @@ -510,7 +510,7 @@ subroutine get_column_partition( self, & integer(i_def) :: num_xprocs, num_yprocs - call log_event("Using column decomposiiton", LOG_LEVEL_INFO) + call log_event("Using column decomposiiton", LOG_LEVEL_DEBUG) num_xprocs = 1_i_def num_yprocs = panel_ranks @@ -601,7 +601,7 @@ subroutine get_auto_nonuniform_partition( self, & integer(i_def) :: start_xprocs, start_width, i logical ::found_factors - call log_event("Using auto_nonuniform decomposition", LOG_LEVEL_INFO) + call log_event("Using auto_nonuniform decomposition", LOG_LEVEL_DEBUG) mp_num_cells_x = num_cells_x / mapping_factor mp_num_cells_y = num_cells_y / mapping_factor @@ -716,7 +716,7 @@ subroutine get_guided_nonuniform_partition( self, & integer(i_def) :: num_xprocs - call log_event("Using guided_nonuniform decomposition", LOG_LEVEL_INFO) + call log_event("Using guided_nonuniform decomposition", LOG_LEVEL_DEBUG) num_xprocs = self%num_xprocs diff --git a/infrastructure/source/mesh/partition_mod.F90 b/infrastructure/source/mesh/partition_mod.F90 index ba5ab1c11..1b10832f9 100644 --- a/infrastructure/source/mesh/partition_mod.F90 +++ b/infrastructure/source/mesh/partition_mod.F90 @@ -25,7 +25,6 @@ module partition_mod use sort_mod, only : bubble_sort use log_mod, only : log_event, & log_scratch_space, & - LOG_LEVEL_INFO, & LOG_LEVEL_ERROR, & LOG_LEVEL_DEBUG use constants_mod, only: i_def, r_def, l_def @@ -976,7 +975,7 @@ subroutine partitioner_rectangular_panels( global_mesh, & " num_y ", num_y call log_event( log_scratch_space, LOG_LEVEL_DEBUG ) write(log_scratch_space,"(a,i0,a,i0)") "Number of cells in partition ", num_x, " X ", num_y - call log_event( log_scratch_space, lOG_LEVEL_INFO ) + call log_event( log_scratch_space, lOG_LEVEL_DEBUG ) ! Create a linked list of all cells in the partition and at the same time ! create a linked-list of all edge cells known to the partition, excluding halos. diff --git a/mesh_tools/Makefile b/mesh_tools/Makefile index 6ed8916b4..2f520b0a5 100644 --- a/mesh_tools/Makefile +++ b/mesh_tools/Makefile @@ -50,6 +50,8 @@ export INTERNAL_DEPENDENCIES = $(CORE_ROOT_DIR)/infrastructure \ export SUITE_GROUP ?= developer export SUITE_GROUP_NAME ?= $(notdir $(realpath $(shell pwd)/$(CORE_ROOT_DIR)))-$(PROJECT_NAME)-.* +export PRE_PROCESS_MACROS := MESH_TOOLS + META_VN ?= HEAD META_FILE_DIR = $(PROJECT_DIR)/rose-meta/lfric-$(PROJECT_NAME)/$(META_VN) diff --git a/rose-stem/app/check_global_variables/file/dirtylist.txt b/rose-stem/app/check_global_variables/file/dirtylist.txt index 5ead624c1..e18c6aef9 100644 --- a/rose-stem/app/check_global_variables/file/dirtylist.txt +++ b/rose-stem/app/check_global_variables/file/dirtylist.txt @@ -17,7 +17,7 @@ infrastructure/source/io/io_utility_mod.f90 infrastructure/source/mesh/global_mesh_collection_mod.F90 infrastructure/source/mesh/global_mesh_mod.F90 infrastructure/source/mesh/local_mesh_collection_mod.f90 -infrastructure/source/mesh/local_mesh_mod.f90 +infrastructure/source/mesh/local_mesh_mod.F90 infrastructure/source/mesh/mesh_collection_mod.F90 infrastructure/source/mesh/mesh_mod.F90 infrastructure/source/utilities/count_mod.f90