module test_apar_shift_m !! Tests regarding apar_shift module !! Sets up an analytical expression of apar varying in time !! Solve apar_shift according to the prescribed apar !! To check the convergence, you need decrease dtau manually in !! '/src/models/braginskii/evolution_apar_shift/tests/params_apar_shift.in' use pfunit use MPI use funit use grillix_build_info_m, only : git_repo_path use comm_handler_m, only : comm_handler_t use precision_grillix_m, only : GP use status_codes_grillix_m, only : GRILLIX_ERR_NAMELIST use error_handling_grillix_m, only: handle_error use helpers_m, only : almost_equal use helpers_mesh_sample_m, only : mesh_sample_circular use screen_io_m, only : get_stdout use model_apar_shift_m, only : apar_shift_t use params_apar_shift_m, only : read_params_evol_apar_shift, & tau_filter_width, apar_shift_iselect use params_tstep_m, only : read_params_tstep, dtau, tstep_order use equilibrium_m, only: equilibrium_t use mesh_cart_m, only : mesh_cart_t use variable_m, only : variable_t use commons_misc_m, only : error_ananum use constants_m, only : two_pi implicit none contains @test(npes = [4]) subroutine test_apar_shift(this) !! Tests apar_shift module in CIRCULAR geometry class (MpiTestMethod), intent(inout) :: this integer :: comm_world, rank, ierr type(comm_handler_t) :: comm_handler class(equilibrium_t), allocatable :: equi type(mesh_cart_t) :: mesh_stag type(mesh_cart_t) :: mesh_cano type(apar_shift_t) :: ashift type(variable_t) :: apar real(GP) :: tau, tau_fin, tau_adv, tau_init_bwd, c0, c1, c2, c3, m1, m2, m3 real(GP) :: nrm2, nrmsup, err2, errsup real(GP) :: phi_stag, phi_cano, dphi integer :: l, istorage real(GP), dimension(:), allocatable :: ana_apar_shift character(len=:), allocatable :: input_params_file real(GP), dimension(:,:), allocatable :: y_firsts real(GP), dimension(:,:), allocatable :: dy_firsts character(len=22) :: fname logical :: fexist integer :: ioerr, isnaps_apar_shift character(len=256) :: io_errmsg logical :: fine logical :: write_data_on = .false. comm_world = this%getMpiCommunicator() call MPI_comm_rank(comm_world, rank, ierr) if (rank == 0) then write(*,*)'' write(*,*)'' write(*,'(A80)')'test_apar_shift ---------------------------------------------------' endif call comm_handler%init(comm_world, 4, 1) ! Create samplemesh --------------------------------------------------------------------- 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_stag, mesh_stag) call mesh_sample_circular(equi, phi_cano, mesh_cano) ! Set up parameters c0 = 10.0_GP c1 = 20.0_GP c2 = 2.0_GP c3 = 1.0_GP m1 = 50.0_GP m2 = 50.0_GP m3 = 1.0_GP tau_fin = 500.0_GP ! Read parameters input_params_file = git_repo_path // & '/src/models/braginskii/evolution_apar_shift/tests/params_apar_shift.in' call read_params_evol_apar_shift(input_params_file) call read_params_tstep(input_params_file) ! Initialise if(.not. allocated(ana_apar_shift)) allocate(ana_apar_shift(mesh_stag%get_n_points())) if(.not. allocated(y_firsts)) allocate(y_firsts(mesh_stag%get_n_points(),tstep_order-1)) if(.not. allocated(dy_firsts)) allocate(dy_firsts(mesh_stag%get_n_points(),tstep_order-1)) call apar%init(comm_handler, mesh_stag%get_n_points(), staggered=.true.) ! Initialising --------------------------------------------------------------------- tau = 0.0_GP !$omp parallel default(none) private(l) & !$omp shared(mesh_stag, apar, tau, c0, c1, c2, c3, m1, m2, m3) !$omp do do l = 1, mesh_stag%get_n_points() apar%vals(l) = time_signal(tau, c0, c1, c2, c3, m1, m2, m3) enddo !$omp end do !$omp end parallel do istorage = 1, tstep_order - 1 tau_init_bwd = tau - real(istorage, GP) * dtau !$omp parallel default(none) private(l) & !$omp shared(mesh_stag, dy_firsts, y_firsts, tau_init_bwd, istorage, tau_filter_width, & !$omp c0, c1, c2, c3, m1, m2, m3) !$omp do do l = 1, mesh_stag%get_n_points() y_firsts(l,istorage) = time_signal_shift(tau_init_bwd, c0, c1, c2, c3, m1, m2, m3, & tau_filter_width) dy_firsts(l,istorage) = 0.0_GP enddo !$omp end do !$omp end parallel enddo call ashift%init(comm_handler, equi, mesh_cano, mesh_stag, apar, tau, isnaps_apar_shift, & apar_shift_firsts=y_firsts, dapar_shift_firsts=dy_firsts) ! Timestepping --------------------------------------------------------------------- do while (tau <= tau_fin) tau_adv = tau + dtau ! Reset analytical results of apar and apar_shift !$omp parallel default(none) private(l) & !$omp shared(mesh_stag, apar, ana_apar_shift, tau_adv, tau_filter_width, & !$omp c0, c1, c2, c3, m1, m2, m3) !$omp do do l = 1, mesh_stag%get_n_points() apar%vals(l) = time_signal(tau_adv, c0, c1, c2, c3, m1, m2, m3) ana_apar_shift(l) = time_signal_shift(tau_adv, c0, c1, c2, c3, m1, m2, m3, & tau_filter_width) enddo !$omp end do !$omp end parallel ! Compute apar_shift at next timestep tau + dtau call ashift%timestep(comm_handler, mesh_stag, tau, apar) tau = tau + dtau ! Write out results for checking if ( write_data_on .and. (rank == 0) ) then fname = 'test_apar_shift.dat' inquire(file=fname, exist = fexist) if (fexist) then open(unit=56, file=fname, status='old', position='append', & action='write', iostat=ioerr, iomsg=io_errmsg) else open(unit=56, file=fname, status='new', & action='write', iostat=ioerr, iomsg=io_errmsg) endif if (ioerr /= 0) then call handle_error(io_errmsg, GRILLIX_ERR_NAMELIST, __LINE__, __FILE__) endif write(56,603)tau write(56,603)apar%vals(1) write(56,603)ashift%apar_shift%vals(1) write(56,603)ana_apar_shift(1) 603 FORMAT(ES20.12E3) close(56) endif enddo ! Compute errors between analytical and numerical solution call error_ananum(comm_handler, mesh_stag%get_n_points(), & ana_apar_shift, ashift%apar_shift%vals, & nrm2, nrmsup, err2, errsup) if (rank == 0) then write(get_stdout(),*)'test_apar_shift with ------------------------------' write(get_stdout(),202)tau_filter_width, dtau, ana_apar_shift(1), ashift%apar_shift%vals(1), err2 202 FORMAT( 3X,'tau_filter_width = ',ES14.6E3 /, & 3X,'dtau = ',ES14.6E3 /, & 3X,'apar_shift (ana) = ',ES14.6E3 /, & 3X,'apar_shift (num) = ',ES14.6E3 /, & 3X,'error = ',ES14.6E3) endif ! Check results ------------------------------------------------------- fine = almost_equal(err2, 2.760226E-3_GP, 0.0_GP, 1.0E-9_GP) @assertTrue(fine) if (rank == 0) then write(*,'(A80)')'test_apar_shift complete --------------------------------' write(*,*)'' write(*,*)'' endif contains real(GP) function time_signal(t, c0, c1, c2, c3, m1, m2, m3) !! Analytical expression of apar varyting in time real(GP), intent(in) :: t !! time real(GP), intent(in) :: c0 !! initial value real(GP), intent(in) :: c1 real(GP), intent(in) :: c2 real(GP), intent(in) :: c3 real(GP), intent(in) :: m1 real(GP), intent(in) :: m2 real(GP), intent(in) :: m3 time_signal = c0 + c1*(1.0_GP - exp(-t/m1)) & + c2*sin(t/m2) + c3*sin(t/m3) end real(GP) function time_signal_shift(t, c0, c1, c2, c3, m1, m2, m3, lambda) !! Analytical expression of apar_shift varyting in time !! Determined by time_signal real(GP), intent(in) :: t !! Time real(GP), intent(in) :: c0 !! initial value real(GP), intent(in) :: c1 real(GP), intent(in) :: c2 real(GP), intent(in) :: c3 real(GP), intent(in) :: m1 real(GP), intent(in) :: m2 real(GP), intent(in) :: m3 real(GP), intent(in) :: lambda time_signal_shift = c0 + (c1 / (m1-lambda)) * & (m1-m1*exp(-t/m1)-lambda + lambda*exp(-t/lambda)) & + (c2*m2/(m2**2+lambda**2)) * & (lambda*exp(-t/lambda)-lambda*cos(t/m2)+m2*sin(t/m2)) & + (c3*m3/(m3**2+lambda**2)) * & (lambda*exp(-t/lambda)-lambda*cos(t/m3)+m3*sin(t/m3)) end end subroutine end module