module test_polar_operators_m use pfunit use MPI implicit none contains @test(npes = [4]) subroutine polar_operators_circular(this) !! Test of polar operators in circular geometry use precision_grillix_m, only : GP use comm_handler_m, only : comm_handler_t use helpers_mesh_sample_m, only : mesh_sample_circular use helpers_m, only : almost_equal use polars_m, only : polars_t use polar_operators_m, only : map_cart_to_polar, map_polar_to_cart, & surface_average, flux_surface_average, & drift_surface_integral use equilibrium_storage_m, only : equilibrium_storage_t use parallel_map_m, only : parallel_map_t use constants_m, only : two_pi use equilibrium_m, only : equilibrium_t use circular_equilibrium_m, only : circular_t use mesh_cart_m, only : mesh_cart_t use coords_polar_m, only : polar_to_cart class (MpiTestMethod), intent(inout) :: this integer :: comm_world type(comm_handler_t) :: comm_handler class(equilibrium_t), allocatable :: equi type(equilibrium_storage_t) :: equi_storage type(mesh_cart_t) :: mesh type(mesh_cart_t) :: mesh_stag type(parallel_map_t) :: map type(polars_t) :: polars integer :: rank, nprocs, ierr logical :: fine integer :: l, nrho, ntheta, ip, jp, ic, jc, lc real(GP) :: x, y, phi, dphi, phi_stag, rho, theta real(GP), allocatable, dimension(:) :: u, v real(GP), allocatable, dimension(:,:) :: upol real(GP), allocatable, dimension(:) :: ucartmap real(GP), allocatable, dimension(:) :: usrf, uflsrf, udrf ! reference values real(GP), dimension(0:3) :: up_ref, ucm_ref real(GP), dimension(20) :: uzon_ref real(GP) :: zonint_ref real(GP), dimension(20) :: drift_flux_ref comm_world = this%getMpiCommunicator() call MPI_comm_rank(comm_world, rank, ierr) call MPI_comm_size(comm_world, nprocs, ierr) if (rank == 0) then write(*,*)'' write(*,*)'' write(*,'(A80)')'test_polar_operators_circular ---------------------------------------------' endif call comm_handler%init(comm_world, 4, 1) dphi = TWO_PI / comm_handler%get_nplanes() phi = dphi * comm_handler%get_plane() phi_stag = phi + dphi / 2.0_GP ! Create samplemesh --------------------------------------------------------------------------- call mesh_sample_circular(equi, phi, mesh) call mesh_sample_circular(equi, phi_stag, mesh_stag) ! Create equi_storage ------------------------------------------------------------------------- call equi_storage%fill_storage(equi, mesh) ! 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, mesh_stag) ! Create polars ------------------------------------------------------------------------------- call polars%set_parameters(nrho_in = 20, & ntheta_in = 100, & intorder_in = 3 ) call polars%build(equi, mesh) nrho = polars%grid%get_nrho() ntheta = polars%grid%get_ntheta() ! Check xpol and ypol ------------------------------------------------------------------------- do ip = 1, nrho rho = polars%grid%get_rho(ip) do jp = 1, ntheta theta = polars%grid%get_theta(jp) call polar_to_cart(equi, rho, theta, phi, x, y) fine = almost_equal(polars%xpol(jp, ip), x , rtol=0.0_GP, atol=1.0E-10_GP) @assertTrue(fine) fine = almost_equal(polars%ypol(jp, ip), y , rtol=0.0_GP, atol=1.0E-10_GP) @assertTrue(fine) enddo enddo ! Create some state --------------------------------------------------------------------------- allocate(u(mesh%get_n_points())) allocate(v(mesh%get_n_points())) !$omp parallel default(none) private (l, x, y) shared(mesh, u, v, phi) !$omp do do l = 1, mesh%get_n_points() x = mesh%get_x(l) y = mesh%get_y(l) u(l) = ((x-0.05_GP)*(x+0.1_GP)-(y-0.1)*(y-0.15_GP))*sin(phi) & + 0.8_GP*cos(30*x)*cos(60*y) * (x*(x-0.15_GP)-(y-0.2_GP)) v(l) = ((x-0.1_GP) * (x+0.05_GP) - (y-0.05_GP) * (y-0.1_GP)) * cos(phi) & + 0.6_GP * sin(20.0_GP*x) * sin(50.0_GP*y) * (x*(x-0.1_GP)-(y-0.15_GP)) enddo !$omp end do !$omp end parallel ! Reference values ----------------------------------------------------------------------------- up_ref = [ 0.205873640022745_GP, 0.149052075118978_GP, & 0.205873640022745_GP, 0.262695204926513_GP ] ucm_ref = [ -0.23123346478239221_GP, -0.27211358271335950_GP, & -0.23123346478239221_GP, -0.19035334685142502_GP ] uzon_ref = [ -3.170009530211E-02_GP, -3.297795781436385E-02_GP, -2.771959180182548E-02_GP, & -1.710546798000904E-02_GP, -3.357540300754342E-03_GP, 1.081099583085063E-02_GP, & 2.260493548883027E-02_GP, 2.975316048358748E-02_GP, 3.095151236802976E-02_GP, & 2.603102297108491E-02_GP, 1.604157664448725E-02_GP, 3.063021726150079E-03_GP, & -1.036745063434765E-02_GP, -2.157334237052562E-02_GP, -2.839188574491934E-02_GP, & -2.952442586529498E-02_GP, -2.481493549655764E-02_GP, -1.528624085726353E-02_GP, & -2.830958671176336E-03_GP, 9.976806998472E-003_GP ] drift_flux_ref = [ 0.806452290989191_GP, 0.783730443947067_GP, 0.468444941897647_GP, & -0.109060511021405_GP, -0.798873131210421_GP, -1.38373915767915_GP, & -1.68342177119950_GP, -1.64592021130817_GP, -1.37537250434989_GP, & -1.07679180520316_GP, -0.942791234820728_GP, -1.04718438275051_GP, & -1.29279101817654_GP, -1.46181148266171_GP, -1.34215823012756_GP, & -0.867090625675885_GP, -0.176801025818685_GP, 0.423204668701239_GP, & 0.607194583493583_GP, 0.191117511721084_GP ] zonint_ref = -2.380074589475478E-004_GP ! test map_cart_to_polar ----------------------------------------------------------------------- allocate(upol(ntheta, nrho)) call map_cart_to_polar(mesh, polars, u, upol) fine = almost_equal(upol(89,13), up_ref(rank) , rtol = 0.0_GP, atol = 1.0E-10_GP) @assertTrue(fine) ! test map_polar_to_cart ---------------------------------------------------------------------- allocate(ucartmap(mesh%get_n_points())) call map_polar_to_cart(equi, mesh, polars, upol, ucartmap) ic = 24 jc = 33 lc = mesh%cart_to_mesh_index(ic, jc) fine = almost_equal(ucartmap(lc), ucm_ref(rank), rtol = 0.0_GP, atol = 1.0E-8_GP) @assertTrue(fine) ! test surface average ------------------------------------------------------------------------ allocate(usrf(nrho)) call surface_average(comm_world, mesh, polars, u, usrf) do ip = 1, nrho fine = almost_equal(usrf(ip), uzon_ref(ip) , rtol = 0.0_GP, atol = 1.0E-10_GP) @assertTrue(fine) enddo ! test flux surface average ------------------------------------------------------------------- allocate(uflsrf(nrho)) call flux_surface_average(comm_world, mesh, polars, u, uflsrf) do ip = 1, nrho fine = almost_equal(uflsrf(ip), uzon_ref(ip) , rtol = 0.0_GP, atol = 1.0E-10_GP) @assertTrue(fine) enddo ! Test drift_surface_integral ------------------------------------------------------------- allocate(udrf(nrho)) call drift_surface_integral(comm_handler%get_comm_planes(), equi, mesh, polars, & map%get_dphi(), u, udrf) ! For homogeneous magenetic field (circular), the result should be zero do ip = 1, nrho fine = almost_equal(udrf(ip), 0.0_GP , rtol = 0.0_GP, atol = 1.0E-10_GP) @assertTrue(fine) enddo call drift_surface_integral(comm_handler%get_comm_planes(), equi, mesh, polars, & map%get_dphi(), u, udrf, v) do ip = 1, nrho fine = almost_equal(udrf(ip), drift_flux_ref(ip) , rtol = 0.0_GP, atol = 1.0E-8_GP) @assertTrue(fine) enddo if (rank == 0) then write(*,'(A80)')'test_polar_operators_circular complete ------------------------------------' write(*,*)'' write(*,*)'' endif end subroutine @test(npes = [4]) subroutine polar_operators_cerfons(this) !! Test of polar operators in circular geometry use precision_grillix_m, only : GP use comm_handler_m, only : comm_handler_t use helpers_mesh_sample_m, only : mesh_sample_cerfons use helpers_m, only : almost_equal use polars_m, only : polars_t use polar_operators_m, only : map_cart_to_polar, map_polar_to_cart, & surface_average, flux_surface_average, & drift_surface_integral use equilibrium_storage_m, only : equilibrium_storage_t use parallel_map_m, only : parallel_map_t use constants_m, only : two_pi use equilibrium_m, only : equilibrium_t use cerfons_equilibrium_m, only : cerfons_t use mesh_cart_m, only : mesh_cart_t use coords_polar_m, only : polar_to_cart class (MpiTestMethod), intent(inout) :: this integer :: comm_world type(comm_handler_t) :: comm_handler class(equilibrium_t), allocatable :: equi type(equilibrium_storage_t) :: equi_storage type(mesh_cart_t) :: mesh type(mesh_cart_t) :: mesh_stag type(parallel_map_t) :: map type(polars_t) :: polars integer :: rank, nprocs, ierr logical :: fine integer :: l, nrho, ntheta, ip, jp, ic, jc, lc real(GP) :: x, y, phi, dphi, phi_stag, rho, theta real(GP), allocatable, dimension(:) :: u, v real(GP), allocatable, dimension(:) :: ucartmap real(GP), allocatable, dimension(:,:) :: upol real(GP), allocatable, dimension(:) :: usrf, uflsrf, udrf ! reference values real(GP), dimension(0:3) :: up_ref, ucm_ref real(GP), dimension(20) :: usrf_ref, uflsrf_ref, drift_flux_ref real(GP) :: zonint_ref comm_world = this%getMpiCommunicator() call MPI_comm_rank(comm_world, rank, ierr) call MPI_comm_size(comm_world, nprocs, ierr) if (rank == 0) then write(*,*)'' write(*,*)'' write(*,'(A80)')'test_polar_operators_cerfons ----------------------------------------------' endif call comm_handler%init(comm_world, 4, 1) dphi = TWO_PI / comm_handler%get_nplanes() phi = dphi * comm_handler%get_plane() phi_stag = phi + dphi / 2.0_GP ! Create samplemesh --------------------------------------------------------------------------- call mesh_sample_cerfons(equi, phi, mesh) call mesh_sample_cerfons(equi, phi_stag, mesh_stag) ! Create equi_storage ------------------------------------------------------------------------- call equi_storage%fill_storage(equi, mesh) ! 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, mesh_stag) ! Create polars ------------------------------------------------------------------------------- call polars%set_parameters(nrho_in = 20, & ntheta_in = 100, & intorder_in = 3 ) call polars%build(equi, mesh) nrho = polars%grid%get_nrho() ntheta = polars%grid%get_ntheta() ! Check xpol and ypol ------------------------------------------------------------------------- do ip = 1, nrho rho = polars%grid%get_rho(ip) do jp = 1, ntheta theta = polars%grid%get_theta(jp) call polar_to_cart(equi, rho, theta, phi, x, y) fine = almost_equal(polars%xpol(jp, ip), x , rtol=0.0_GP, atol=1.0E-10_GP) @assertTrue(fine) fine = almost_equal(polars%ypol(jp, ip), y , rtol=0.0_GP, atol=1.0E-10_GP) @assertTrue(fine) enddo enddo ! Create some state --------------------------------------------------------------------------- allocate(u(mesh%get_n_points())) allocate(v(mesh%get_n_points())) !$omp parallel default(none) private (l, x, y) shared(mesh, u, v, phi) !$omp do do l = 1, mesh%get_n_points() x = mesh%get_x(l) y = mesh%get_y(l) u(l) = ((x-1.05_GP)*(x+1.1_GP)-(y-1.1)*(y-1.15_GP))*sin(phi) & + 0.8_GP*cos(30*x)*cos(60*y) * (x*(x-1.15_GP)-(y-1.2_GP)) v(l) = ((x-1.10_GP)*(x+1.05_GP)-(y-1.15)*(y-1.0_GP))*cos(phi) & + 0.6_GP*sin(30.0_GP*x)*sin(50.0_GP*y) * (x*(x-1.05_GP)-(y-1.1_GP)) enddo !$omp end do !$omp end parallel ! Reference values ----------------------------------------------------------------------------- up_ref = [ 0.272682526866696_GP, -1.10535710939789_GP, & 0.272682526866695_GP, 1.65072216313128_GP ] ucm_ref = [ 0.26217561082406066_GP, -2.1648215526130934_GP, & 0.26217561082406043_GP, 2.6891727742612148_GP ] usrf_ref = [ 9.665525387754E-002_GP, 7.757342574733E-003_GP, -9.039100247919E-002_GP, & -1.215628662821E-001_GP, -6.250181425138E-002_GP, 3.396619817826E-002_GP, & 8.836215385409E-002_GP, 6.999031115968E-002_GP, 1.862688988754E-002_GP, & -1.360672195000E-002_GP, -2.636295827749E-002_GP, -4.844905619327E-002_GP, & -6.224993416605E-002_GP, -1.362950268601E-002_GP, 6.960785452847E-002_GP, & 7.393012742691E-002_GP, -4.681504709458E-003_GP, -3.853136750694E-002_GP, & -2.885548825614E-002_GP, 4.879389699515E-003_GP ] uflsrf_ref = [ 1.180944635568E-001_GP, 1.271762080217E-002_GP, -1.073043424058E-001_GP, & -1.432826101975E-001_GP, -6.574432770571E-002_GP, 5.285481259821E-002_GP, & 1.068145973500E-001_GP, 6.456891103536E-002_GP, -3.600027576022E-003_GP, & -2.101890761877E-002_GP, -7.808934329393E-003_GP, -3.536420389551E-002_GP, & -8.024819979330E-002_GP, -2.950108742164E-002_GP, 9.447643363413E-002_GP, & 9.515932480939E-002_GP, -3.233136294540E-002_GP, -4.983260974696E-002_GP, & -2.864963074305E-002_GP, 2.352492046971E-002_GP ] drift_flux_ref = [ -24.2483111593422_GP, -29.5292992027886_GP, -24.6126046458981_GP, & -19.7218113935363_GP, -18.6762529659512_GP, -12.4035426442304_GP, & -2.475512504820415E-002_GP, 6.28455667933672_GP, 7.85321010347493_GP, & 17.3257328247358_GP, 24.4998029789076_GP, 18.4661640921833_GP, & 14.3615681453487_GP, 6.16974204946794_GP, -10.7643014258170_GP, & -16.9901361295376_GP, -26.1160621477588_GP, -23.2510862556703_GP, & -20.4161558112943_GP, -11.1086660103638_GP] zonint_ref = -7.795814956601165E-003 ! test map_cart_to_polar ----------------------------------------------------------------------- allocate(upol(ntheta, nrho)) call map_cart_to_polar(mesh, polars, u, upol) fine = almost_equal(upol(89,13), up_ref(rank) , rtol = 0.0_GP, atol = 1.0E-10_GP) @assertTrue(fine) ! test map_polar_to_cart ---------------------------------------------------------------------- allocate(ucartmap(mesh%get_n_points())) call map_polar_to_cart(equi, mesh, polars, upol, ucartmap) ic = 67 jc = 57 lc = mesh%cart_to_mesh_index(ic, jc) fine = almost_equal(ucartmap(lc), ucm_ref(rank), rtol = 0.0_GP, atol = 1.0E-8_GP) @assertTrue(fine) ! test surface average ------------------------------------------------------------------------ allocate(usrf(nrho)) call surface_average(comm_world, mesh, polars, u, usrf) do ip = 1, nrho fine = almost_equal(usrf(ip), usrf_ref(ip) , rtol = 0.0_GP, atol = 1.0E-8_GP) @assertTrue(fine) enddo ! test flux surface average ------------------------------------------------------------------- allocate(uflsrf(nrho)) call flux_surface_average(comm_world, mesh, polars, u, uflsrf) do ip = 1, nrho fine = almost_equal(uflsrf(ip), uflsrf_ref(ip) , rtol = 0.0_GP, atol = 1.0E-8_GP) @assertTrue(fine) enddo ! Test drift_surface_integral ------------------------------------------------------------- allocate(udrf(nrho)) call drift_surface_integral(comm_handler%get_comm_planes(), equi, mesh, polars, & map%get_dphi(), u, udrf, v) do ip = 1, nrho fine = almost_equal(udrf(ip), drift_flux_ref(ip) , rtol = 0.0_GP, atol = 1.0E-8_GP) @assertTrue(fine) enddo if (rank == 0) then write(*,'(A80)')'test_polar_operators_cerfons complete -------------------------------------' write(*,*)'' write(*,*)'' endif end subroutine end module