module timestep_template_m !! Advances template model by one timestep use perf_m, only : perf_start, perf_stop use precision_grillix_m, only : GP use error_handling_grillix_m, only : handle_error use status_codes_grillix_m, only : GRILLIX_ERR_TSTEP use comm_handler_m, only : comm_handler_t use static_data_m, only : mesh_cano, masks_cano, mesh_stag, masks_stag, map use state_template_m, only : np_max, dens, parmom, upar, pion, tion, update_derived_fields use params_template_m, only : tstep_method, tstep_order, dtau, mms_on, perptransp_on, partransp_on, pardiff_on, & parbnd_method use runge_kutta_m, only : runge_kutta_t use multistep_m, only : karniadakis_t use immersed_template_m, only : apply_immersed_sources use taylor_template_m, only : apply_taylor_parbnds use, intrinsic :: iso_c_binding implicit none type, bind(C) :: static_data_struct !! This struct shall hold algorithm-relevant data structures !! that are fixed for a fixed numerical experiment (i.e. numerical solution). !! (e.g. mesh) type(c_ptr) :: mesh_cano_data end type static_data_struct type, bind(C) :: runtime_data_struct !! This struct shall hold fixed runtime data. !! mpi communicator. integer in fortran, MPI_Fint in C !! use MPI_Comm_f2c in C/C++ to get an MPI_Comm integer :: mpi_communicator end type runtime_data_struct type, bind(C) :: state_data_struct !! This struct shall hold data needed in a checkpoint !! (e.g. dynamical system state variables) integer(c_int64_t) :: npoints type(c_ptr) :: density real(GP) :: tau end type state_data_struct real(GP), private, allocatable, dimension(:) :: ddens !! Changerate for density real(GP), private, allocatable, dimension(:) :: dparmom !! Changerate for density real(GP), private, allocatable, dimension(:) :: dpion !! Changerate for ion pressure integer, parameter, private :: neq = 3 !! Number of equations type(runge_kutta_t), private :: tstepper_rk !! Runge-Kutta timestepper type(comm_handler_t), pointer, private :: ptr_comm_handler => null() !! Pointer to communication handler, needed for interface handling in Runge-Kutta time-step integrator type(karniadakis_t), private :: tstepper_karks !! Karniadakis timestepper public :: timestep_template #ifdef ENABLE_GRILLACCX interface cxx_timestep integer(c_int32_t) function GrillAccX_timestep( & static_object, & runtime_object, & state_object) & bind (c, name='timestep') use, intrinsic :: iso_c_binding import static_data_struct import runtime_data_struct import state_data_struct type(static_data_struct), intent(inout) :: static_object type(runtime_data_struct), intent(inout) :: runtime_object type(state_data_struct), intent(inout) :: state_object end function GrillAccX_timestep end interface #endif interface ! Definition of routines that compute changerate from individual model comonents module subroutine apply_partransp(comm_handler) !! Applies changerates due to parallel transport model type(comm_handler_t), intent(in) :: comm_handler !! Communication handler end subroutine module subroutine apply_perptransp(comm_handler) !! Applies changerates due to perpendicular diffusive transport model type(comm_handler_t), intent(in) :: comm_handler !! Communication handler end subroutine module subroutine apply_pardiff(comm_handler) !! Applies changerates due to parallel diffusion type(comm_handler_t), intent(in) :: comm_handler !! Communication handler end subroutine module subroutine apply_mms_sources(comm_handler, tau) !! Applies changerates due to MMS sources type(comm_handler_t), intent(in) :: comm_handler !! Communication handler real(GP), intent(in) :: tau !! Current time end subroutine module subroutine apply_mms_boundaries(comm_handler, tau) !! Applies boundaries according MMS solutions type(comm_handler_t), intent(in) :: comm_handler !! Communication handler real(GP), intent(in) :: tau !! Current time end subroutine end interface contains integer(c_int32_t) function fortran_timestep( & comm_handler, & tau) type(comm_handler_t), target :: comm_handler !! Communication handler real(GP), intent(inout) :: tau !! Current time logical, save :: first_call = .true. real(GP) :: tau_adv real(GP), dimension(neq*np_max) :: wrk, dwrk, wrk_adv integer :: l tau_adv = tau + dtau ! Initialize time-steppers if (first_call) then first_call = .false. ptr_comm_handler => comm_handler select case(tstep_method) case('Runge-Kutta') call tstepper_rk%init(neq*np_max, tstep_order) case('Karniadakis') call tstepper_karks%init(neq*np_max, tstep_order, dtau) case default call handle_error('tstep_method not valid: ' // tstep_method, & GRILLIX_ERR_TSTEP, __LINE__, __FILE__) end select endif ! Perform explicit time stepping part select case(tstep_method) case('Runge-Kutta') !$omp parallel default(none) private(l) shared(np_max, dens, parmom, pion, wrk) !$omp do do l = 1, np_max wrk(l) = dens(l) wrk(np_max + l) = parmom(l) wrk(2*np_max + l) = pion(l) enddo !$omp end do !$omp end parallel call tstepper_rk%advance(dtau, tau, wrk, rhs_fun_for_rk) !$omp parallel default(none) shared(np_max, dens, parmom, pion, wrk) !$omp do do l = 1, np_max dens(l) = wrk(l) parmom(l) = wrk(np_max + l) pion(l) = wrk(2*np_max + l) enddo !$omp end do !$omp end parallel case('Karniadakis') call compute_rhs_explicit(comm_handler, tau) !$omp parallel default(none) private(l) & !$omp shared(np_max, dens, ddens, parmom, dparmom, pion, dpion, wrk, dwrk) !$omp do do l = 1, np_max wrk(l) = dens(l) dwrk(l) = ddens(l) wrk(np_max + l) = parmom(l) dwrk(np_max + l) = dparmom(l) wrk(2*np_max + l) = pion(l) dwrk(2*np_max + l) = dpion(l) enddo !$omp end do !$omp end parallel call tstepper_karks%advance(dwrk, wrk, wrk_adv) call tstepper_karks%shift_storage((/wrk, dwrk/)) !$omp parallel default(none) shared(np_max, dens, parmom, pion, wrk_adv) !$omp do do l = 1, np_max dens(l) = wrk_adv(l) parmom(l) = wrk_adv(np_max+l) pion(l) = wrk_adv(2*np_max+l) enddo !$omp end do !$omp end parallel case default call handle_error('tstep_method not valid: ' // tstep_method, & GRILLIX_ERR_TSTEP, __LINE__, __FILE__) end select ! Apply boundary conditions call apply_boundaries(comm_handler, tau_adv) ! Compute derived variables call update_derived_fields(comm_handler) fortran_timestep = 0 end function subroutine timestep_template(comm_handler, tau) !! Model template subroutine !! Yet a simple in-plane diffusion equation with first order explicit Euler scheme is implemented type(comm_handler_t), target :: comm_handler !! Communication handler real(GP), intent(inout) :: tau !! Current time integer(c_int32_t) :: ccall_result type(static_data_struct) :: static_object type(runtime_data_struct) :: runtime_object type(state_data_struct) :: state_object call perf_start('timestep_template') #ifdef ENABLE_GRILLACCX runtime_object%mpi_communicator = comm_handler%get_comm_cart() state_object%npoints = np_max state_object%density = c_loc(dens) ccall_result = cxx_timestep(static_object, runtime_object, state_object) #else ccall_result = fortran_timestep(comm_handler, tau) #endif call perf_stop('timestep_template') end subroutine subroutine compute_rhs_explicit(comm_handler, tau) !! Computes explicit changerate at time tau type(comm_handler_t), intent(in) :: comm_handler !! Communication handler real(GP), intent(inout) :: tau !! Time-point integer :: l ! Initialise changerates to zero if (.not.allocated(ddens)) allocate(ddens(np_max)) if (.not.allocated(dparmom)) allocate(dparmom(np_max)) if (.not.allocated(dpion)) allocate(dpion(np_max)) !$omp parallel default(none) private(l) shared(np_max, ddens, dparmom, dpion) !$omp do do l = 1, np_max ddens(l) = 0.0_GP dparmom(l) = 0.0_GP dpion(l) = 0.0_GP enddo !$omp end do !$omp end parallel if (partransp_on) then call apply_partransp(comm_handler) endif if (perptransp_on) then call apply_perptransp(comm_handler) endif if (pardiff_on) then call apply_pardiff(comm_handler) endif if (mms_on) then call apply_mms_sources(comm_handler, tau) endif if (parbnd_method == 'Immersed') then call apply_immersed_sources(comm_handler, ddens, dparmom, dpion) endif end subroutine subroutine apply_boundaries(comm_handler, tau) !! Applies boundaries type(comm_handler_t), intent(in) :: comm_handler !! Communication handler real(GP), intent(in) :: tau !! Time-point ! Set boundary conditions if (mms_on) then call apply_mms_boundaries(comm_handler, tau) return endif if (parbnd_method == 'Taylor') then call apply_taylor_parbnds(comm_handler) endif end subroutine subroutine rhs_fun_for_rk(ndim, tau, y, dy) !! Evaluates right hand side (explicit terms) for use in Runge-Kutta time-step integrator !! Be aware, that this function changes the state of the template model. !! For interface, see the corresponding Runge-Kutta module integer, intent(in) :: ndim real(GP), intent(in) :: tau real(GP), dimension(ndim), intent(in) :: y real(GP), dimension(ndim), intent(out) :: dy integer :: l real(GP) :: tau_wrk tau_wrk = tau !$omp parallel default(none) private(l) shared(np_max, dens, parmom, pion, y) !$omp do do l = 1, np_max dens(l) = y(l) parmom(l) = y(np_max+l) pion(l) = y(2*np_max+l) enddo !$omp end do !$omp end parallel call update_derived_fields(ptr_comm_handler) call apply_boundaries(ptr_comm_handler, tau_wrk) call compute_rhs_explicit(ptr_comm_handler, tau_wrk) !$omp parallel default(none) private(l) shared(np_max, ddens, dparmom, dpion, dy) !$omp do do l = 1, np_max dy(l) = ddens(l) dy(np_max+l) = dparmom(l) dy(2*np_max+l) = dpion(l) enddo !$omp end do !$omp end parallel end subroutine end module