! Copyright (c) 2012-2014 Jens Zudrop <j.zudrop@grs-sim.de> ! Copyright (c) 2013 Verena Krupp <verena.krupp@uni-siegen.de> ! Copyright (c) 2013-2015 Harald Klimach <harald.klimach@uni-siegen.de> ! Copyright (c) 2013, 2016-2018 Peter Vitt <peter.vitt2@uni-siegen.de> ! Copyright (c) 2016 Tobias Girresser <tobias.girresser@student.uni-siegen.de> ! Copyright (c) 2019 Daniel Fleischer <daniel.fleischer@student.uni-siegen.de> ! ! Permission to use, copy, modify, and distribute this software for any ! purpose with or without fee is hereby granted, provided that the above ! copyright notice and this permission notice appear in all copies. ! ! THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHORS DISCLAIM ALL WARRANTIES ! WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF ! MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHORS BE LIABLE FOR ! ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES ! WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN ! ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF ! OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE. ! **************************************************************************** ! !> author: Jens Zudrop !! Routines and datatypes realted to specification and definition !! of the 2D MODG scheme. module atl_modg_2d_scheme_module use env_module, only: rk, labelLen use aotus_module, only: flu_State, aot_get_val use aot_top_module, only: aoterr_Fatal, aoterr_NonExistent use aot_table_module, only: aot_table_open, aot_table_close, aot_table_length use tem_aux_module, only: tem_abort use tem_tools_module, only: upper_to_lower use tem_logging_module, only: logUnit use ply_dof_module, only: Q_space,P_space implicit none private !> Datatype for the modg scheme type atl_modg_2d_scheme_type !> The maximal polynomial degree per spatial direction. integer :: maxPolyDegree !> Polynomial basis type. !! !! 3D Monomials have the form x^i * y^j * z^k !! - Q_space: quadratic polynomial space (i,j,k) <= maxPolyDegree !! - P_space: polynomial space i+j+k <= maxPolyDegree integer :: basisType !> Facial Chebyshev nodes (reference element) for 2 spatial direction !! and left and right face. These points are necessary to transfer !! boundary conditions given in physical space to modal space by !! means of an FPT. end type atl_modg_2d_scheme_type public :: atl_modg_2d_scheme_type public :: atl_modg_2d_scheme_init contains subroutine atl_modg_2d_scheme_init(me, nDofs, nFaceDofs, & & conf, thandle, currentLevel, maxLevel ) type(atl_modg_2d_scheme_type), intent(out) :: me integer, intent(out) :: nDofs integer, intent(out) :: nFaceDofs type(flu_State) :: conf integer, intent(in) :: thandle !> The current level of the mesh integer, intent(in) :: currentLevel !> The maximal level of the mesh integer, intent(in) :: maxLevel character(len=labelLen) :: modg_2d_space character(len=labelLen) :: polyspace integer :: degree integer :: iError ! Find the maximal polynomial degree for the current level. call atl_modg_2d_scheme_load_polyDegree(degree, conf, thandle, & & currentlevel, maxlevel) ! Read the chosen polynomial space call aot_get_val(L = conf, thandle = thandle, & & key = 'modg_2d_space', & & val = modg_2d_space, & & ErrCode = iError, & & default = 'Q' ) polyspace = upper_to_lower(modg_2d_space) polyspace = adjustl(polyspace) select case(trim(polyspace)) case('q') me%basisType = Q_space nDofs = (degree+1)**2 nFaceDofs = (degree+1)**1 me%maxPolyDegree = degree case('p') me%basisType = P_space nDofs = ((degree+1)*(degree+2))/2 nFacedofs = degree+1 me%maxPolyDegree = degree case default write(logUnit(1),*) 'ERROR in init_kernel_state: unknown modg_2d_space: ' & & // trim(modg_2d_space) // ' !' write(logUnit(1),*) 'Supported are:' write(logUnit(1),*) '* Q (quadratic with i,j,k <= maxDegree' write(logUnit(1),*) 'Stopping....' call tem_abort() end select write(logUnit(1),*) 'MODG scheme with '//trim(modg_2d_space)//' basis on level ', & & currentLevel write(logUnit(1),*) '... maximal polynomial degree: ', degree end subroutine atl_modg_2d_scheme_init subroutine atl_modg_2d_scheme_load_polyDegree(degree, conf, thandle, & & currentlevel, maxlevel) integer, intent(out) :: degree type(flu_state) :: conf integer, intent(in) :: thandle integer, intent(in) :: currentLevel integer, intent(in) :: maxLevel character(len=labelLen) :: predefDegrees character(len=labelLen) :: predeg integer :: mtable integer :: degTable integer :: iError integer :: nDegTables integer :: iTable integer :: tableLevel logical :: foundLevel real(kind=rk) :: levelFact, base_order ! Get the degree of the scheme call aot_get_val(L = conf, thandle = thandle, & & key = 'm', & & val = degree, & & ErrCode = iError) ! if no error occurs, we have a global polynomial degree if (iError.ne.0) then ! open the spatial subtable call aot_table_open(L=conf, parent=thandle, & & thandle = mtable, key = 'm') call aot_get_val(L = conf, thandle = mtable, key = 'predefined', & & val = predefDegrees, ErrCode = iError) predef: if (btest(iError, aoterr_nonExistent)) then ! No predefined set of polynomial degrees, try to get them ! individually for each level. ! The number of tables in the m table nDegTables = aot_table_length(L = conf, thandle = mTable) ! Iterate over all the sub tables and check if one has ! the current mesh level foundLevel = .false. do iTable = 1, nDegTables call aot_table_open(L = conf, parent = mTable, thandle = degTable, & & pos = iTable) ! get the level of the current table call aot_get_val( L = conf, thandle = degTable, key = 'level', & & val = tableLevel, ErrCode = iError ) if (iError.ne.0) then write(logUnit(1),*) 'ERROR in atl_modg_2d_scheme_init: not able to read' & & //' level for ' write(logUnit(1),*) ' the level specific polynomial degree!' write(logUnit(1),*) 'STOPPING ...' call tem_abort() end if if (tableLevel.eq.currentLevel) then ! read the degree for the level call aot_get_val( L = conf, thandle = degTable, key = 'm', & & val = degree, ErrCode = iError ) if (iError.ne.0) then write(logUnit(1),*) 'ERROR in atl_modg_2d_scheme_init: not able to ' & & //' read degree for ' write(logUnit(1),*) ' the level specific polynomial degree!' write(logUnit(1),*) ' Failing level: ', currentLevel write(logUnit(1),*) 'STOPPING ...' call tem_abort() end if foundLevel = .true. end if call aot_table_close(L = conf, thandle = degtable) end do ! Check if we found a degree for the current level of the mesh if (.not.foundlevel) then write(logUnit(1),*) 'ERROR in atl_modg_2d_scheme_init:' write(logUnit(1),*) ' level not found for the level specific' & & // ' polynomial degree!' write(logUnit(1),*) ' Failing level: ', currentLevel write(logUnit(1),*) 'STOPPING ...' call tem_abort() end if else predef if (btest(iError, aoterr_Fatal)) then write(logUnit(1),*) 'ERROR in atl_modg_2d_scheme_init: ' & & // 'could not read the predefined polynomial' write(logUnit(1),*) ' degrees for different levels!' write(logUnit(1),*) 'ABORTING...' call tem_abort() end if predeg = upper_to_lower(predefDegrees) predeg = adjustl(predeg) select case(trim(predeg)) case ('fixedfact') ! This set of degrees increases the polynomial degrees by a fixed ! factor from level to level according to the following formula: ! m(iLevel) = nint(baseOrder*(factor**(maxLevel-iLevel)) - 1) ! Get the factor to increase the order level by level. ! Defaults to sqrt(2), as this maintains approximately the same ! time step restriction throughout the domain. call aot_get_val(L = conf, thandle = mtable, key = 'factor', & & val = levelFact, default = sqrt(2.0_rk), & & ErrCode = iError) if (btest(iError, aoterr_Fatal)) then write(logUnit(1),*) 'ERROR in atl_modg_2d_scheme_init: ' & & // 'could not read the factor for the' write(logUnit(1),*) ' predefined fixedfact set of degrees!' write(logUnit(1),*) 'ABORTING...' call tem_abort() end if ! Get the base order to use on the maximal level. ! Defaults to 1. call aot_get_val(L = conf, thandle = mtable, key = 'base_order', & & val = base_order, default = 1.0_rk, & & ErrCode = iError) if (btest(iError, aoterr_Fatal)) then write(logUnit(1),*) 'ERROR in atl_modg_2d_scheme_init: ' & & // 'could not read the base_order for the' write(logUnit(1),*) ' predefined fixedfact set of degrees!' write(logUnit(1),*) 'ABORTING...' call tem_abort() end if degree = nint( base_order*(levelFact**(maxLevel-currentLevel)) & & - 1.0_rk ) case default write(logUnit(1),*) 'ERROR in atl_modg_2d_scheme_init: ' & & // 'unknoen predefined polynomial degrees:' write(logUnit(1),*) trim(predefDegrees) write(logUnit(1),*) 'Currently supported predefined sets are:' write(logUnit(1),*) ' * fixedFact (use a fixed factor between levels)' write(logUnit(1),*) 'ABORTING...' call tem_abort() end select end if predef ! Close the m table call aot_table_close(L = conf, thandle = mtable) end if end subroutine atl_modg_2d_scheme_load_polyDegree end module atl_modg_2d_scheme_module