module test_parallel_map_operators_m !! Test of parallel operators within map module use pfunit use grillix_build_info_m, only : git_repo_path use MPI use precision_grillix_m, only : GP use constants_m, only : TWO_PI use screen_io_m, only : get_stdout use elementary_functions_m, only : gaussian use comm_handler_m, only : comm_handler_t use csrmat_m, only : csrmat_t use equilibrium_m, only : equilibrium_t use equilibrium_factory_m, only : create_equilibrium, DOMMASCHK use mesh_cart_m, only : mesh_cart_t use parallel_map_m, only : parallel_map_t use helpers_m, only : almost_equal use helpers_mesh_sample_m, only : mesh_sample_circular, mesh_sample_cerfons use mpi_mapping_auxiliaries_m, only : getdata_fwdbwdplane implicit none contains @test(npes = [4]) subroutine test_parallel_map_operators_circular(this) !! Test of parallel operators in circular geometry class (MpiTestMethod), intent(inout) :: this ! Reference data real(GP), dimension(4), parameter :: ref_par_grad = & [-6.223015336415E-002_GP, -5.164509702544E-002_GP, -1.345652334683E-003_GP, 4.319200809822E-001_GP] real(GP), dimension(4), parameter :: ref_par_divb = & [-7.362700505646E-001_GP, -2.088306581708E-001_GP, -2.706162886153E-003_GP, 4.701891193702E-002_GP] real(GP), dimension(4), parameter :: ref_lift_cano_to_stag = & [ 2.471933301369E-001_GP, 4.555677185862E-002_GP, 1.069731575332E-003_GP, 3.405262724933E-001_GP] real(GP), dimension(4), parameter :: ref_flat_stag_to_cano = & [ 4.755317446589E-001_GP, 1.418992425092E-001_GP, 5.011460665839E-003_GP, 2.909123914153E-002_GP] real(GP), dimension(4), parameter :: ref_upstream_cano_from_cano_fwd = & [3.383005384639E-002_GP, 8.354773206648E-004_GP, 1.749798174382E-006_GP, 1.161694759674E-001_GP] real(GP), dimension(4), parameter :: ref_upstream_cano_from_cano_bwd = & [3.418780954873E-007_GP, 2.269736006068E-002_GP, 6.609764799487E-003_GP, 1.632367660416E-004_GP] real(GP), dimension(4), parameter :: ref_upstream_stag_from_stag_fwd = & [9.918286308906E-002_GP, 2.449045564990E-003_GP, 2.494517314496E-003_GP, 3.412488281323E-001_GP] real(GP), dimension(4), parameter :: ref_upstream_stag_from_stag_bwd = & [1.480229771901E-001_GP, 4.333953060744E-002_GP, 1.132511585313E-003_GP, 3.907645481862E-006_GP] real(GP), dimension(4), parameter :: ref_upstream_cano_from_cano_fwd2 = & [4.377349081126231E-007_GP, 9.167786176042117E-010_GP, 6.086512898711020E-005_GP, 1.772471274278609E-005_GP] real(GP), dimension(4), parameter :: ref_upstream_cano_from_cano_bwd2 = & [2.172802448042987E-008_GP, 4.550651062380860E-011_GP, 3.021186997259530E-006_GP, 8.798087272600082E-007_GP] integer :: comm_world type(comm_handler_t) :: comm_handler class(equilibrium_t), allocatable :: equi type(mesh_cart_t) :: mesh_cano type(mesh_cart_t) :: mesh_stag type(parallel_map_t) :: map integer :: rank, ierr real(GP) :: phi_cano, phi_stag, dphi, x, y integer :: ic, jc, lc, l_cano, nrecv logical :: fine real(GP), allocatable, dimension(:) :: u_cano , wrk_halo, & par_grad_u_stag, ulift, & uret comm_world = this%getMpiCommunicator() call MPI_comm_rank(comm_world, rank, ierr) if (rank == 0) then write(*,*)'' write(*,*)'' write(get_stdout(),'(A80)')'test_parallel_map_operators_circular'& //repeat('-', 80) endif call comm_handler%init(comm_world, 4, 1) dphi = TWO_PI / comm_handler%get_nplanes() phi_cano = dphi * comm_handler%get_plane() phi_stag = phi_cano + dphi / 2.0_GP ! Create samplemeshes call mesh_sample_circular(equi, phi_cano, mesh_cano) call mesh_sample_circular(equi, phi_stag, mesh_stag) ! Create map call map%set_parameters(comm_handler, & intorder_in=3, & xorder_in=1, & use_fixed_stencil_in=.false., & use_gauss_quadrature_in=.false.) call map%build(comm_handler, equi, mesh_cano, mesh_stag) ! Select test point ic = 88 jc = 55 lc = mesh_stag%cart_to_mesh_index(ic, jc) ! Create some variable allocate(u_cano(mesh_cano%get_n_points())) do l_cano = 1, mesh_cano%get_n_points() x = mesh_cano%get_x(l_cano) y = mesh_cano%get_y(l_cano) u_cano(l_cano) = gaussian(0.3_GP, sqrt(2.0_GP)*0.05_GP, x) * & gaussian(0.05_GP, sqrt(2.0_GP)*0.05_GP, y) * & gaussian(0.0_GP, sqrt(2.0_GP), phi_cano) enddo ! Parallel gradient allocate(par_grad_u_stag(mesh_stag%get_n_points())) call getdata_fwdbwdplane(comm_world, 1, mesh_cano%get_n_points(), & u_cano, nrecv, wrk_halo) !$omp parallel call map%par_grad(u_cano, wrk_halo, par_grad_u_stag) !$omp end parallel deallocate(wrk_halo) fine = almost_equal(par_grad_u_stag(lc), ref_par_grad(rank+1), & rtol=0.0_GP, atol=1.0E-8_GP) @assertTrue(fine) ! Parallel divergence of parallel gradient (=parallel diffusion) allocate(uret(mesh_cano%get_n_points())) call getdata_fwdbwdplane(comm_world, -1, mesh_stag%get_n_points(), & par_grad_u_stag, nrecv, wrk_halo) !$omp parallel call map%par_divb(par_grad_u_stag, wrk_halo, uret) !$omp end parallel deallocate(wrk_halo) fine = almost_equal(uret(lc), ref_par_divb(rank+1), & rtol=0.0_GP, atol=1.0E-8_GP) @assertTrue(fine) deallocate(uret) ! Lift from cano to stag allocate(ulift(mesh_stag%get_n_points())) call getdata_fwdbwdplane(comm_world, 1, mesh_cano%get_n_points(), & u_cano, nrecv, wrk_halo) !$omp parallel call map%lift_cano_to_stag(u_cano, wrk_halo, ulift) !$omp end parallel deallocate(wrk_halo) fine = almost_equal(ulift(lc), ref_lift_cano_to_stag(rank+1), & rtol=0.0_GP, atol=1.0E-8_GP) @assertTrue(fine) ! Flat from stag to cano allocate(uret(mesh_cano%get_n_points())) call getdata_fwdbwdplane(comm_world, -1, mesh_stag%get_n_points(), & ulift, nrecv, wrk_halo) !$omp parallel call map%flat_stag_to_cano(ulift, wrk_halo, uret) !$omp end parallel deallocate(wrk_halo) fine = almost_equal(uret(lc), ref_flat_stag_to_cano(rank+1), & rtol=0.0_GP, atol=1.0E-8_GP) @assertTrue(fine) deallocate(uret) ! Upstream cano from cano_fwd allocate(uret(mesh_cano%get_n_points())) call getdata_fwdbwdplane(comm_world, 1, mesh_cano%get_n_points(), & u_cano, nrecv, wrk_halo) !$omp parallel call map%upstream_cano_from_cano_fwd(wrk_halo, uret) !$omp end parallel deallocate(wrk_halo) fine = almost_equal(uret(lc), ref_upstream_cano_from_cano_fwd(rank+1), & rtol=0.0_GP, atol=1.0E-8_GP) @assertTrue(fine) deallocate(uret) ! Upstream cano from cano_bwd allocate(uret(mesh_cano%get_n_points())) call getdata_fwdbwdplane(comm_world, -1, mesh_cano%get_n_points(), & u_cano, nrecv, wrk_halo) !$omp parallel call map%upstream_cano_from_cano_bwd(wrk_halo, uret) !$omp end parallel deallocate(wrk_halo) fine = almost_equal(uret(lc), ref_upstream_cano_from_cano_bwd(rank+1), & rtol=0.0_GP, atol=1.0E-8_GP) @assertTrue(fine) deallocate(uret) ! Upstream stag from stag_fwd allocate(uret(mesh_stag%get_n_points())) call getdata_fwdbwdplane(comm_world, 1, mesh_stag%get_n_points(), & ulift, nrecv, wrk_halo) !$omp parallel call map%upstream_stag_from_stag_fwd(wrk_halo, uret) !$omp end parallel deallocate(wrk_halo) fine = almost_equal(uret(lc), ref_upstream_stag_from_stag_fwd(rank+1), & rtol=0.0_GP, atol=1.0E-8_GP) @assertTrue(fine) deallocate(uret) ! Upstream stag from stag_bwd allocate(uret(mesh_stag%get_n_points())) call getdata_fwdbwdplane(comm_world, -1, mesh_stag%get_n_points(), & ulift, nrecv, wrk_halo) !$omp parallel call map%upstream_stag_from_stag_bwd(wrk_halo, uret) !$omp end parallel deallocate(wrk_halo) fine = almost_equal(uret(lc), ref_upstream_stag_from_stag_bwd(rank+1), & rtol=0.0_GP, atol=1.0E-8_GP) @assertTrue(fine) deallocate(uret) ! Upstream cano from cano_fwd2 allocate(uret(mesh_cano%get_n_points())) call getdata_fwdbwdplane(comm_world, 2, mesh_cano%get_n_points(), & u_cano, nrecv, wrk_halo) !$omp parallel call map%upstream_cano_from_cano_fwd2(wrk_halo, uret) !$omp end parallel deallocate(wrk_halo) fine = almost_equal(uret(lc), ref_upstream_cano_from_cano_fwd2(rank+1), & rtol=0.0_GP, atol=1.0E-8_GP) @assertTrue(fine) deallocate(uret) ! Upstream cano from cano_bwd2 allocate(uret(mesh_cano%get_n_points())) call getdata_fwdbwdplane(comm_world, -2, mesh_cano%get_n_points(), & u_cano, nrecv, wrk_halo) !$omp parallel call map%upstream_cano_from_cano_bwd2(wrk_halo, uret) !$omp end parallel deallocate(wrk_halo) fine = almost_equal(uret(lc), ref_upstream_cano_from_cano_bwd2(rank+1), & rtol=0.0_GP, atol=1.0E-8_GP) @assertTrue(fine) deallocate(uret) if (rank == 0) then write(get_stdout(),'(A80)')'test_parallel_map_operators_circular complete' & //repeat('-', 80) write(*,*)'' write(*,*)'' endif end subroutine @test(npes = [4]) subroutine test_parallel_map_operators_cerfons(this) !! Test of parallel operators in CERFONS geometry class (MpiTestMethod), intent(inout) :: this ! Reference data real(GP), dimension(4), parameter :: ref_par_grad = & [-7.864002961751E-001_GP,-2.742990794921E-001_GP,-6.869025865711E-003_GP, 5.835101409647E-001_GP] real(GP), dimension(4), parameter :: ref_par_divb = & [-1.942608737908E+000_GP, 5.410583256838E-002_GP, 1.704791771843E-001_GP, 1.970177633976E-001_GP] real(GP), dimension(4), parameter :: ref_lift_cano_to_stag = & [5.438147275781E-001_GP, 1.362451646519E-001_GP, 3.318423157202E-003_GP, 2.850196565647E-001_GP] real(GP), dimension(4), parameter :: ref_flat_stag_to_cano = & [4.872730022456E-001_GP, 2.858838761626E-001_GP, 4.637469026462E-002_GP, 4.688100590737E-002_GP] real(GP), dimension(4), parameter :: ref_upstream_cano_from_cano_fwd = & [5.275286994380E-002_GP, 1.302800954386E-003_GP, 2.728546514889E-006_GP, 1.811487881448E-001_GP] real(GP), dimension(4), parameter :: ref_upstream_cano_from_cano_bwd = & [ 9.043355549649E-006_GP, 6.003903139060E-001_GP, 1.748414243858E-001_GP, 4.317937105580E-003_GP] real(GP), dimension(4), parameter :: ref_upstream_stag_from_stag_fwd = & [8.229912253626E-002_GP, 2.030093081135E-003_GP, 1.471434305373E-002_GP, 2.865285455473E-001_GP] real(GP), dimension(4), parameter :: ref_upstream_stag_from_stag_bwd = & [4.706223163375E-001_GP, 2.782732478967E-001_GP, 4.451048207284E-002_GP, 1.022744892520E-003_GP] real(GP), dimension(4), parameter :: ref_upstream_cano_from_cano_fwd2 = & [1.733348253823518E-005_GP, 3.630271624483158E-008_GP, 2.410145115075145E-003_GP, 7.018654284324773E-004_GP] real(GP), dimension(4), parameter :: ref_upstream_cano_from_cano_bwd2 = & [8.025716044014930E-004_GP, 1.680881447595878E-006_GP, 0.111594079699764_GP, 3.249763927868744E-002_GP] integer :: comm_world type(comm_handler_t) :: comm_handler class(equilibrium_t), allocatable :: equi type(mesh_cart_t) :: mesh_cano type(mesh_cart_t) :: mesh_stag type(parallel_map_t) :: map integer :: rank, ierr real(GP) :: phi_cano, phi_stag, dphi, x, y integer :: ic, jc, lc, l_cano, nrecv logical :: fine real(GP), allocatable, dimension(:) :: u_cano , wrk_halo, & par_grad_u_stag, ulift, & uret comm_world = this%getMpiCommunicator() call MPI_comm_rank(comm_world, rank, ierr) if (rank == 0) then write(*,*)'' write(*,*)'' write(get_stdout(),'(A80)')'test_parallel_map_operators_cerfons'& //repeat('-', 80) endif call comm_handler%init(comm_world, 4, 1) dphi = TWO_PI / comm_handler%get_nplanes() phi_cano = dphi * comm_handler%get_plane() phi_stag = phi_cano + dphi / 2.0_GP ! Create samplemeshes call mesh_sample_cerfons(equi, phi_cano, mesh_cano) call mesh_sample_cerfons(equi, phi_stag, mesh_stag) ! Create map call map%set_parameters(comm_handler, & intorder_in=3, & xorder_in=1, & use_fixed_stencil_in=.false., & use_gauss_quadrature_in=.false.) call map%build(comm_handler, equi, mesh_cano, mesh_stag) ! Select test point ic = 22 jc = 130 lc = mesh_stag%cart_to_mesh_index(ic, jc) ! Create some variable allocate(u_cano(mesh_cano%get_n_points())) do l_cano = 1, mesh_cano%get_n_points() x = mesh_cano%get_x(l_cano) y = mesh_cano%get_y(l_cano) u_cano(l_cano) = gaussian(0.625_GP, sqrt(2.0_GP)*0.1_GP, x) * & gaussian(0.32_GP, sqrt(2.0_GP)*0.1_GP, y) * & gaussian(0.0_GP, sqrt(2.0_GP), phi_cano) enddo ! Parallel gradient allocate(par_grad_u_stag(mesh_stag%get_n_points())) call getdata_fwdbwdplane(comm_world, 1, mesh_cano%get_n_points(), & u_cano, nrecv, wrk_halo) !$omp parallel call map%par_grad(u_cano, wrk_halo, par_grad_u_stag) !$omp end parallel deallocate(wrk_halo) fine = almost_equal(par_grad_u_stag(lc), ref_par_grad(rank+1), & rtol=0.0_GP, atol=1.0E-8_GP) @assertTrue(fine) ! Parallel divergence of parallel gradient (=parallel diffusion) allocate(uret(mesh_cano%get_n_points())) call getdata_fwdbwdplane(comm_world, -1, mesh_stag%get_n_points(), & par_grad_u_stag, nrecv, wrk_halo) !$omp parallel call map%par_divb(par_grad_u_stag, wrk_halo, uret) !$omp end parallel deallocate(wrk_halo) fine = almost_equal(uret(lc), ref_par_divb(rank+1), & rtol=0.0_GP, atol=1.0E-8_GP) @assertTrue(fine) deallocate(uret) ! Lift from cano to stag allocate(ulift(mesh_stag%get_n_points())) call getdata_fwdbwdplane(comm_world, 1, mesh_cano%get_n_points(), & u_cano, nrecv, wrk_halo) !$omp parallel call map%lift_cano_to_stag(u_cano, wrk_halo, ulift) !$omp end parallel deallocate(wrk_halo) fine = almost_equal(ulift(lc), ref_lift_cano_to_stag(rank+1), & rtol=0.0_GP, atol=1.0E-8_GP) @assertTrue(fine) ! Flat from stag to cano allocate(uret(mesh_cano%get_n_points())) call getdata_fwdbwdplane(comm_world, -1, mesh_stag%get_n_points(), & ulift, nrecv, wrk_halo) !$omp parallel call map%flat_stag_to_cano(ulift, wrk_halo, uret) !$omp end parallel deallocate(wrk_halo) fine = almost_equal(uret(lc), ref_flat_stag_to_cano(rank+1), & rtol=0.0_GP, atol=1.0E-8_GP) @assertTrue(fine) deallocate(uret) ! Upstream cano from cano_fwd allocate(uret(mesh_cano%get_n_points())) call getdata_fwdbwdplane(comm_world, 1, mesh_cano%get_n_points(), & u_cano, nrecv, wrk_halo) !$omp parallel call map%upstream_cano_from_cano_fwd(wrk_halo, uret) !$omp end parallel deallocate(wrk_halo) fine = almost_equal(uret(lc), ref_upstream_cano_from_cano_fwd(rank+1), & rtol=0.0_GP, atol=1.0E-8_GP) @assertTrue(fine) deallocate(uret) ! Upstream cano from cano_bwd allocate(uret(mesh_cano%get_n_points())) call getdata_fwdbwdplane(comm_world, -1, mesh_cano%get_n_points(), & u_cano, nrecv, wrk_halo) !$omp parallel call map%upstream_cano_from_cano_bwd(wrk_halo, uret) !$omp end parallel deallocate(wrk_halo) fine = almost_equal(uret(lc), ref_upstream_cano_from_cano_bwd(rank+1), & rtol=0.0_GP, atol=1.0E-8_GP) @assertTrue(fine) deallocate(uret) ! Upstream stag from stag_fwd allocate(uret(mesh_stag%get_n_points())) call getdata_fwdbwdplane(comm_world, 1, mesh_stag%get_n_points(), & ulift, nrecv, wrk_halo) !$omp parallel call map%upstream_stag_from_stag_fwd(wrk_halo, uret) !$omp end parallel deallocate(wrk_halo) fine = almost_equal(uret(lc), ref_upstream_stag_from_stag_fwd(rank+1), & rtol=0.0_GP, atol=1.0E-8_GP) @assertTrue(fine) deallocate(uret) ! Upstream stag from stag_bwd allocate(uret(mesh_stag%get_n_points())) call getdata_fwdbwdplane(comm_world, -1, mesh_stag%get_n_points(), & ulift, nrecv, wrk_halo) !$omp parallel call map%upstream_stag_from_stag_bwd(wrk_halo, uret) !$omp end parallel deallocate(wrk_halo) fine = almost_equal(uret(lc), ref_upstream_stag_from_stag_bwd(rank+1), & rtol=0.0_GP, atol=1.0E-8_GP) @assertTrue(fine) deallocate(uret) ! Upstream cano from cano_fwd2 allocate(uret(mesh_cano%get_n_points())) call getdata_fwdbwdplane(comm_world, 2, mesh_cano%get_n_points(), & u_cano, nrecv, wrk_halo) !$omp parallel call map%upstream_cano_from_cano_fwd2(wrk_halo, uret) !$omp end parallel deallocate(wrk_halo) fine = almost_equal(uret(lc), ref_upstream_cano_from_cano_fwd2(rank+1), & rtol=0.0_GP, atol=1.0E-8_GP) @assertTrue(fine) deallocate(uret) ! Upstream cano from cano_bwd2 allocate(uret(mesh_cano%get_n_points())) call getdata_fwdbwdplane(comm_world, -2, mesh_cano%get_n_points(), & u_cano, nrecv, wrk_halo) !$omp parallel call map%upstream_cano_from_cano_bwd2(wrk_halo, uret) !$omp end parallel deallocate(wrk_halo) fine = almost_equal(uret(lc), ref_upstream_cano_from_cano_bwd2(rank+1), & rtol=0.0_GP, atol=1.0E-8_GP) @assertTrue(fine) deallocate(uret) if (rank == 0) then write(get_stdout(),'(A80)')'test_parallel_map_operators_cerfons complete' & //repeat('-', 80) write(*,*)'' write(*,*)'' endif end subroutine @test(npes = [4]) subroutine test_parallel_map_operators_dommaschk(this) !! Test of parallel operators in DOMMASCHK geometry class (MpiTestMethod), intent(inout) :: this ! Reference data real(GP), dimension(4), parameter :: ref_par_grad = & [2.779905629320E-001_GP,-5.255974317293E-001_GP,-6.118377346607E-001_GP, 2.925139349134E-001_GP] real(GP), dimension(4), parameter :: ref_par_divb = & [2.866733465135E-002_GP,-6.905372725838E-001_GP,-1.312856536729E-002_GP, 6.926204204244E-001_GP] real(GP), dimension(4), parameter :: ref_lift_cano_to_stag = & [ 2.316111848482E-001_GP, 4.320751533931E-001_GP,-5.004551140213E-001_GP,-2.408941951855E-001_GP] real(GP), dimension(4), parameter :: ref_flat_stag_to_cano = & [ 1.936385391946E-002_GP, 4.675727074930E-001_GP,-1.467780101449E-002_GP,-4.662794678176E-001_GP] real(GP), dimension(4), parameter :: ref_upstream_cano_from_cano_fwd = & [ 8.974177312118E-001_GP, 1.079195331548E-016_GP,-7.981296894031E-001_GP, 0.000000000000E+000_GP] real(GP), dimension(4), parameter :: ref_upstream_cano_from_cano_bwd = & [-8.294652731052E-001_GP, 0.000000000000E+000_GP, 7.443812163106E-001_GP, 9.093783009664E-017_GP] real(GP), dimension(4), parameter :: ref_upstream_stag_from_stag_fwd = & [ 2.313171262900E-001_GP,-2.211530441758E-001_GP,-4.928008603293E-001_GP, 3.374181768854E-001_GP] real(GP), dimension(4), parameter :: ref_upstream_stag_from_stag_bwd = & [-1.664854227309E-001_GP, 4.250984632892E-001_GP, 3.839619462922E-001_GP,-2.433039336376E-001_GP] real(GP), dimension(4), parameter :: ref_upstream_cano_from_cano_fwd2 = & [ 7.598169429512891E-017_GP, -0.607795252631153_GP,0.000000000000000E+000_GP, 0.578996700617533_GP] real(GP), dimension(4), parameter :: ref_upstream_cano_from_cano_bwd2 = & [ 7.749685229419076E-017_GP, -0.611528140543295_GP, 0.000000000000000E+000_GP, 0.470069721699250_GP] integer :: comm_world type(comm_handler_t) :: comm_handler class(equilibrium_t), allocatable :: equi character(len=:), allocatable :: equi_params_file type(mesh_cart_t) :: mesh_cano type(mesh_cart_t) :: mesh_stag type(parallel_map_t) :: map integer :: rank, ierr real(GP) :: spacing_f, phi_cano, phi_stag, dphi, x, y integer :: ic, jc, lc_cano, lc_stag, l_cano, nrecv logical :: fine real(GP), allocatable, dimension(:) :: u_cano , wrk_halo, & par_grad_u_stag, ulift, & uret comm_world = this%getMpiCommunicator() call MPI_comm_rank(comm_world, rank, ierr) if (rank == 0) then write(*,*)'' write(*,*)'' write(get_stdout(),'(A80)')'test_parallel_map_operators_dommaschk'& //repeat('-', 80) endif call comm_handler%init(comm_world, 4, 1) dphi = TWO_PI / comm_handler%get_nplanes() phi_cano = dphi * comm_handler%get_plane() phi_stag = phi_cano + dphi / 2.0_GP ! Create equilibrium equi_params_file = git_repo_path // & '/cicd/res/params_equi_dommaschk.in' call create_equilibrium(equi, DOMMASCHK, equi_params_file) ! Create mesh spacing_f = 1E-2 call mesh_cano%init(equi, phi_cano, lvl=1, spacing_f=spacing_f, & size_neighbor=2, size_ghost_layer=2, & extend_beyond_wall=.false., dbgout=0) call mesh_cano%reorder(grid_reorder_block=8, dbgout=0) call mesh_stag%init(equi, phi_stag, lvl=1, spacing_f=spacing_f, & size_neighbor=2, size_ghost_layer=2, & extend_beyond_wall=.false., dbgout=0) call mesh_stag%reorder(grid_reorder_block=8, dbgout=0) ! Create map call map%set_parameters(comm_handler, & intorder_in=3, & xorder_in=1, & use_fixed_stencil_in=.false., & use_gauss_quadrature_in=.false.) call map%build(comm_handler, equi, mesh_cano, mesh_stag) ! Select test point ic = 25 jc = 22 lc_cano = mesh_cano%cart_to_mesh_index(ic, jc) lc_stag = mesh_stag%cart_to_mesh_index(ic, jc) ! Create some variable allocate(u_cano(mesh_cano%get_n_points())) do l_cano = 1, mesh_cano%get_n_points() x = mesh_cano%get_x(l_cano) y = mesh_cano%get_y(l_cano) u_cano(l_cano) = gaussian(1.025_GP, 0.06_GP, x) * & gaussian(0.01_GP, 0.08_GP, y) * & sin(phi_cano) enddo ! Parallel gradient allocate(par_grad_u_stag(mesh_stag%get_n_points())) call getdata_fwdbwdplane(comm_world, 1, mesh_cano%get_n_points(), & u_cano, nrecv, wrk_halo) !$omp parallel call map%par_grad(u_cano, wrk_halo, par_grad_u_stag) !$omp end parallel deallocate(wrk_halo) fine = almost_equal(par_grad_u_stag(lc_stag), ref_par_grad(rank+1), & rtol=0.0_GP, atol=1.0E-8_GP) @assertTrue(fine) ! Parallel divergence of parallel gradient (=parallel diffusion) allocate(uret(mesh_cano%get_n_points())) call getdata_fwdbwdplane(comm_world, -1, mesh_stag%get_n_points(), & par_grad_u_stag, nrecv, wrk_halo) !$omp parallel call map%par_divb(par_grad_u_stag, wrk_halo, uret) !$omp end parallel deallocate(wrk_halo) fine = almost_equal(uret(lc_cano), ref_par_divb(rank+1), & rtol=0.0_GP, atol=1.0E-8_GP) @assertTrue(fine) deallocate(uret) ! Lift from cano to stag allocate(ulift(mesh_stag%get_n_points())) call getdata_fwdbwdplane(comm_world, 1, mesh_cano%get_n_points(), & u_cano, nrecv, wrk_halo) !$omp parallel call map%lift_cano_to_stag(u_cano, wrk_halo, ulift) !$omp end parallel deallocate(wrk_halo) fine = almost_equal(ulift(lc_stag), ref_lift_cano_to_stag(rank+1), & rtol=0.0_GP, atol=1.0E-8_GP) @assertTrue(fine) ! Flat from stag to cano allocate(uret(mesh_cano%get_n_points())) call getdata_fwdbwdplane(comm_world, -1, mesh_stag%get_n_points(), & ulift, nrecv, wrk_halo) !$omp parallel call map%flat_stag_to_cano(ulift, wrk_halo, uret) !$omp end parallel deallocate(wrk_halo) fine = almost_equal(uret(lc_cano), ref_flat_stag_to_cano(rank+1), & rtol=0.0_GP, atol=1.0E-8_GP) @assertTrue(fine) deallocate(uret) ! Upstream cano from cano_fwd allocate(uret(mesh_cano%get_n_points())) call getdata_fwdbwdplane(comm_world, 1, mesh_cano%get_n_points(), & u_cano, nrecv, wrk_halo) !$omp parallel call map%upstream_cano_from_cano_fwd(wrk_halo, uret) !$omp end parallel deallocate(wrk_halo) fine = almost_equal(uret(lc_cano), ref_upstream_cano_from_cano_fwd(rank+1), & rtol=0.0_GP, atol=1.0E-8_GP) @assertTrue(fine) deallocate(uret) ! Upstream cano from cano_bwd allocate(uret(mesh_cano%get_n_points())) call getdata_fwdbwdplane(comm_world, -1, mesh_cano%get_n_points(), & u_cano, nrecv, wrk_halo) !$omp parallel call map%upstream_cano_from_cano_bwd(wrk_halo, uret) !$omp end parallel deallocate(wrk_halo) fine = almost_equal(uret(lc_cano), ref_upstream_cano_from_cano_bwd(rank+1), & rtol=0.0_GP, atol=1.0E-8_GP) @assertTrue(fine) deallocate(uret) ! Upstream stag from stag_fwd allocate(uret(mesh_stag%get_n_points())) call getdata_fwdbwdplane(comm_world, 1, mesh_stag%get_n_points(), & ulift, nrecv, wrk_halo) !$omp parallel call map%upstream_stag_from_stag_fwd(wrk_halo, uret) !$omp end parallel deallocate(wrk_halo) fine = almost_equal(uret(lc_stag), ref_upstream_stag_from_stag_fwd(rank+1), & rtol=0.0_GP, atol=1.0E-8_GP) @assertTrue(fine) deallocate(uret) ! Upstream stag from stag_bwd allocate(uret(mesh_stag%get_n_points())) call getdata_fwdbwdplane(comm_world, -1, mesh_stag%get_n_points(), & ulift, nrecv, wrk_halo) !$omp parallel call map%upstream_stag_from_stag_bwd(wrk_halo, uret) !$omp end parallel deallocate(wrk_halo) fine = almost_equal(uret(lc_stag), ref_upstream_stag_from_stag_bwd(rank+1), & rtol=0.0_GP, atol=1.0E-8_GP) @assertTrue(fine) deallocate(uret) ! Upstream cano from cano_fwd2 allocate(uret(mesh_cano%get_n_points())) call getdata_fwdbwdplane(comm_world, 2, mesh_cano%get_n_points(), & u_cano, nrecv, wrk_halo) !$omp parallel call map%upstream_cano_from_cano_fwd2(wrk_halo, uret) !$omp end parallel deallocate(wrk_halo) fine = almost_equal(uret(lc_cano), ref_upstream_cano_from_cano_fwd2(rank+1), & rtol=0.0_GP, atol=1.0E-8_GP) @assertTrue(fine) deallocate(uret) ! Upstream cano from cano_bwd2 allocate(uret(mesh_cano%get_n_points())) call getdata_fwdbwdplane(comm_world, -2, mesh_cano%get_n_points(), & u_cano, nrecv, wrk_halo) !$omp parallel call map%upstream_cano_from_cano_bwd2(wrk_halo, uret) !$omp end parallel deallocate(wrk_halo) fine = almost_equal(uret(lc_cano), ref_upstream_cano_from_cano_bwd2(rank+1), & rtol=0.0_GP, atol=1.0E-8_GP) @assertTrue(fine) deallocate(uret) if (rank == 0) then write(get_stdout(),'(A80)')'test_parallel_map_operators_dommaschk complete' & //repeat('-', 80) write(*,*)'' write(*,*)'' endif end subroutine end module