submodule (multistep_m) karniadakis_s implicit none contains pure integer module function get_ncount(self) class(karniadakis_t), intent(in) :: self get_ncount = self%ncount end function pure integer module function get_order(self) class(karniadakis_t), intent(in) :: self get_order = self%order end function pure integer module function get_current_order(self) class(karniadakis_t), intent(in) :: self get_current_order = self%current_order end function real(GP) module function get_dtau(self) class(karniadakis_t), intent(in) :: self get_dtau = self%dtau end function real(GP) module function get_dtau_bdf(self) class(karniadakis_t), intent(in) :: self select case (self%current_order) case(1) get_dtau_bdf = self%dtau case(2) get_dtau_bdf = 2.0_GP/3.0_GP*self%dtau case(3) get_dtau_bdf = 6.0_GP/11.0_GP*self%dtau case(4) get_dtau_bdf = 12.0_GP/25.0_GP*self%dtau case(5) get_dtau_bdf = 60.0_GP/137.0_GP*self%dtau case(6) get_dtau_bdf = 60.0_GP/147.0_GP*self%dtau case default call handle_error('Karniadakis order not valid', & GRILLIX_ERR_OTHER, __LINE__, __FILE__, & error_info_t('Order: ',[self%current_order])) end select end function module subroutine init_karniadakis(self, ndim, order, dtau, y_firsts, dy_firsts) class(karniadakis_t), intent(inout) :: self integer, intent(in) :: ndim integer, intent(in) :: order real(GP), intent(in) :: dtau real(GP), dimension(ndim, order-1), intent(in), optional :: y_firsts real(GP), dimension(ndim, order-1), intent(in), optional :: dy_firsts integer :: i, iorder real(GP), dimension(ndim, order-1, 2) :: wrk !G put yfirsts and dyfirsts in self%ncount = 0 self%order = order self%dtau = dtau if ( (present(y_firsts)) .and. (present(dy_firsts)) ) then self%current_order = order !$OMP PARALLEL PRIVATE(iorder, i) do iorder = 1, order-1 !$OMP DO do i = 1, ndim wrk(i, iorder, 1) = y_firsts(i, iorder) wrk(i, iorder, 2) = dy_firsts(i, iorder) enddo !$OMP END DO enddo !$OMP END PARALLEL call self%init_storage(ndim, order-1, 2, wrk) else self%current_order = 0 call self%init_storage(ndim, order-1, 2) endif end subroutine module subroutine advance(self, dy, y, yn, dbgout) class(karniadakis_t), intent(inout) :: self real(GP), dimension(self%ndim), intent(in) :: dy real(GP), dimension(self%ndim), intent(in) :: y real(GP), dimension(self%ndim), intent(out) :: yn integer, intent(in), optional :: dbgout integer :: outi, i, iorder, shft_bwd, shft_fwd ! Debug output level if (present(dbgout)) then outi = dbgout else outi = 0 endif ! Increase self%current_order = min(self%current_order+1, self%order) if ( (outi == 1) .and. (is_rank_info_writer) ) then write(get_stdout(),*)'karniadakis%advance' write(get_stdout(),*)'ncount = ', self%ncount write(get_stdout(),*)'actual order employed = ', self%current_order endif ! Advance in time to valnew --------------------------------------------------------------- !$OMP PARALLEL PRIVATE(i, iorder, shft_bwd, shft_fwd) select case (self%current_order) case(1) !$OMP DO SIMD do i = 1, self%ndim yn(i) = y(i) + self%dtau * dy(i) enddo !$OMP END DO SIMD case(2) !$OMP DO SIMD do i = 1, self%ndim yn(i) = ( 4*y(i) & - self%vstore(i,1,1) ) / 3.0_GP & + 2.0_GP/3.0_GP*self%dtau * ( & + 2*dy(i) & - self%vstore(i,1,2) ) enddo !$OMP END DO SIMD case(3) !$OMP DO SIMD do i = 1, self%ndim yn(i) = ( 18*y(i) & - 9 * self%vstore(i,1,1) & + 2 * self%vstore(i,2,1) ) / 11.0_GP & + 6.0_GP/11.0_GP*self%dtau * ( & + 3 * dy(i) & - 3 * self%vstore(i,1,2) & + self%vstore(i,2,2) ) enddo !$OMP END DO SIMD case(4) !$OMP DO SIMD do i = 1, self%ndim yn(i) = ( 48*y(i) & - 36 * self%vstore(i,1,1) & + 16 * self%vstore(i,2,1) & - 3 * self%vstore(i,3,1) ) / 25.0_GP & + 12.0_GP/25.0_GP*self%dtau * ( & + 4*dy(i) & - 6*self%vstore(i,1,2) & + 4*self%vstore(i,2,2) & - self%vstore(i,3,2) ) enddo !$OMP END DO SIMD case(5) !$OMP DO SIMD do i = 1, self%ndim yn(i) = ( 300 * y(i) & - 300 * self%vstore(i,1,1) & + 200 * self%vstore(i,2,1) & - 75 * self%vstore(i,3,1) & + 12 * self%vstore(i,4,1)) / 137.0_GP & + 60.0_GP/137.0_GP*self%dtau * ( & + 5 * dy(i) & - 10 * self%vstore(i,1,2) & + 10 * self%vstore(i,2,2) & - 5 * self%vstore(i,3,2) & + self%vstore(i,4,2) ) enddo !$OMP END DO SIMD case(6) !$OMP DO SIMD do i = 1, self%ndim yn(i) = ( 360 * y(i) & - 450 * self%vstore(i,1,1) & + 400 * self%vstore(i,2,1) & - 225 * self%vstore(i,3,1) & + 72 * self%vstore(i,4,1) & - 10 * self%vstore(i,5,1)) / 147.0_GP & + 60.0_GP/147.0_GP*self%dtau * ( & + 6 * dy(i) & - 15 * self%vstore(i,1,2) & + 20 * self%vstore(i,2,2) & - 15 * self%vstore(i,3,2) & + 6 * self%vstore(i,4,2) & - self%vstore(i,5,2) ) enddo !$OMP END DO SIMD case default !$omp critical call handle_error('Karniadakis order not valid', & GRILLIX_ERR_OTHER, __LINE__, __FILE__, & error_info_t('Order: ',[self%current_order])) !$omp end critical end select !$OMP END PARALLEL self%ncount = self%ncount + 1 end subroutine module subroutine display_karniadakis(self) class(karniadakis_t), intent(in) :: self if (.not. is_rank_info_writer) then return endif write(get_stdout(),*)'Information on karniadakis ----------------------------' write(get_stdout(),101)self%ncount, self%ndim, self%order, self%current_order, self%dtau 101 FORMAT( 3X,'ncount = ',I10 /, & 3X,'ndim = ',I10 /, & 3X,'order = ',I3 /, & 3X,'current_order = ',I3 /, & 3X,'dtau = ',Es14.6E3 ) write(get_stdout(),*)'-------------------------------------------------------' end subroutine module subroutine destructor_karniadakis(self) type(karniadakis_t), intent(inout) :: self self%ncount = 0 self%order = 0 self%current_order = 0 self%dtau = 0.0_GP end subroutine end submodule