! Copyright (c) 2011-2013 Manuel Hasert <m.hasert@grs-sim.de> ! Copyright (c) 2011-2017 Jiaxing Qi <jiaxing.qi@uni-siegen.de> ! Copyright (c) 2011-2021 Kannan Masilamani <kannan.masilamani@uni-siegen.de> ! Copyright (c) 2011-2012, 2021 Harald Klimach <harald.klimach@uni-siegen.de> ! Copyright (c) 2011-2014 Simon Zimny <s.zimny@grs-sim.de> ! Copyright (c) 2011 Jan Hueckelheim <j.hueckelheim@grs-sim.de> ! Copyright (c) 2012-2014 Kartik Jain <kartik.jain@uni-siegen.de> ! Copyright (c) 2016 Tobias Schneider <tobias1.schneider@student.uni-siegen.de> ! Copyright (c) 2016-2017 Raphael Haupt <raphael.haupt@uni-siegen.de> ! Copyright (c) 2021 Gregorio Gerardo Spinelli <gregoriogerardo.spinelli@dlr.de> ! Copyright (c) 2022 Kannan Masilamani <kannan.masilamani@dlr.de> ! ! Redistribution and use in source and binary forms, with or without ! modification, are permitted provided that the following conditions are met: ! ! 1. Redistributions of source code must retain the above copyright notice, ! this list of conditions and the following disclaimer. ! ! 2. Redistributions in binary form must reproduce the above copyright notice, ! this list of conditions and the following disclaimer in the documentation ! and/or other materials provided with the distribution. ! ! THIS SOFTWARE IS PROVIDED BY THE UNIVERSITY OF SIEGEN “AS IS” AND ANY EXPRESS ! OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES ! OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. ! IN NO EVENT SHALL UNIVERSITY OF SIEGEN OR CONTRIBUTORS BE LIABLE FOR ANY ! DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES ! (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; ! LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ! ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT ! (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS ! SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. ! **************************************************************************** ! !> In this module, the control structure for computing each time step is !! established. !! !! Various different control structures can be defined and set by function !! pointers. !! The current main routine is recursive_multilevel which recursively computes !! the new time step for all the schemes defined on all levels, starting from !! the coarsest. !! module mus_control_module ! include treelm modules use env_module, only: rk, labelLen use tem_timer_module, only: tem_startTimer, tem_stopTimer use tem_general_module, only: tem_general_type use tem_debug_module, only: main_debug, dbgUnit use tem_time_module, only: tem_time_advance, tem_time_dump use tem_logging_module, only: logUnit use tem_varSys_module, only: tem_varSys_type use tem_stencil_module, only: tem_stencilHeader_type ! include musubi modules use mus_pdf_module, only: pdf_data_type, mus_swap_now_next use mus_param_module, only: mus_param_type use mus_geom_module, only: mus_geom_type use mus_bc_general_module, only: set_boundary use mus_aux_module, only: mus_update_relaxParams use mus_scheme_type_module, only: mus_scheme_type use mus_scheme_header_module, only: mus_scheme_header_type use mus_field_module, only: mus_field_type use mus_source_module, only: mus_apply_sourceTerms, & & mus_update_sourceVars use mus_debug_module, only: dump_debug_Info use mus_IBM_module, only: mus_inamuro_IBM, mus_buildBuffIBM use mus_timer_module, only: mus_timerHandles, nStages use mus_physics_module, only: mus_convertFac_type use mus_auxField_module, only: mus_calcAuxFieldAndExchange, & & mus_intpAuxFieldCoarserAndExchange use mus_derVarPos_module, only: mus_derVarPos_type implicit none private public :: mus_control_type public :: mus_init_control !> Datatype containing mapping of control routines to function pointers type mus_control_type procedure( computation ), pointer :: do_computation => null() end type mus_control_type abstract interface !> Interface describes the main control routine which does computation !! set boundary and check flow status subroutine computation( me, scheme, geometry, params, iLevel) import :: mus_scheme_type, mus_geom_type, mus_param_type, & & mus_control_type !> self control type class( mus_control_type ) :: me !> container for the scheme type( mus_scheme_type ), intent(inout) :: scheme !> geometry infomation type( mus_geom_type ), intent(inout) :: geometry !> global parameters type( mus_param_type ), intent(inout) :: params !> Level counter variable integer, intent(in) :: iLevel end subroutine computation end interface integer, save :: iStage = 0 logical, save :: running = .false. contains ! ------------------------------------------------------------------------ ! !> This routines sets the function pointer to main control routine !! !! control routine is chosen based on the type of simulation. !! like single level, multi-level because multilevel requires !! recursive routine !! !! in the lua file you can now select !! `control_routine = '...'` !! where possible values are currently !! !! - `benchmark`: strongly reduced control routine with only single level, !! no sources, etc. !! mainly for benchmarking !! !! - `multiLevel`: full multilevel, multiLevel routine !! - if nothing is given, the full multilevel, multiLevel routine is chosen subroutine mus_init_control( controlRoutine, me, minLevel, maxLevel ) ! -------------------------------------------------------------------- ! character(len=labelLen), intent(in) :: controlRoutine !> contains function pointer to point control routine type( mus_control_type ), intent(out) :: me integer, intent(in) :: minLevel, maxLevel ! -------------------------------------------------------------------- ! ! Select according to special need if ( trim( controlRoutine ) == 'benchmark' ) then me%do_computation => do_benchmark write(logUnit(5),"(A)") "Select benchmark control routine." else if ( trim( controlRoutine ) == 'debug' ) then ! me%do_computation => do_recursive_debug write(logUnit(5),"(A)") "Select debug recursive control routine." else ! No special requirement from user, then select by single or multi levels if ( minLevel == maxLevel ) then me%do_computation => do_fast_singleLevel write(logUnit(5),"(A)") "Select fast single level control routine." else me%do_computation => do_recursive_multiLevel write(logUnit(5),"(A)") "Select recursive multi level control routine." end if end if end subroutine mus_init_control ! ------------------------------------------------------------------------ ! ! ------------------------------------------------------------------------ ! !> Main control routine: Update the time step for all levels. !! Main steps: !! * if iLevel < maxLevel do recursive at iLevel+1 !! * do BC at iLevel !! * do auxField calculation at iLevel !! * do compute kernel at iLevel !! * do apply source at iLevel !! * do do_IntpFinerAndExchange at iLevel if iLevel < maxLevel !! * intp My Coarser ghost (iLevel) from Finer (iLevel+1) !! * do exchange bufferFromFiner at iLevel !! * exchange buffer at iLevel !! * exchange bufferFromCoarser at iLevel if iLevel > minLevel !! * do do_intpCoarserAndExchange at iLevel if iLevel < maxLevel !! * intp Finer Ghost (iLevel+1) from my coarser (iLevel) !! * exchange bufferFromCoarser at iLevel+1 !! recursive subroutine do_recursive_multiLevel( me, scheme, geometry, params, & & iLevel) ! -------------------------------------------------------------------- ! !> self control type class( mus_control_type ) :: me !> container for the scheme type( mus_scheme_type ), intent(inout) :: scheme !> geometry infomation type( mus_geom_type ), intent(inout) :: geometry !> global parameters type( mus_param_type ), intent(inout) :: params !> the current level integer, intent(in) :: iLevel ! -------------------------------------------------------------------- ! integer :: now, next, iNestingLoop ! -------------------------------------------------------------------- ! ! ------------------------------------------------------------------------- write(logUnit(10), "(A,I0)") 'Entering do_computation with iLevel: ', iLevel ! Update auxField dependent source fields before adding source term to state ! and auxField such that both auxField and apply_source uses same source in ! one multilevel cycle. write(logUnit(10), "(A)") 'Update source variables which depend on auxField' call mus_update_sourceVars( nFields = scheme%nFields, & & field = scheme%field, & & globSrc = scheme%globSrc, & & varSys = scheme%varSys, & & iLevel = iLevel, & & auxField = scheme%auxField(iLevel)%val, & & phyConvFac = params%physics%fac(iLevel), & & derVarPos = scheme%derVarPos ) ! when not on finest level, go to next level if( iLevel < geometry%tree%global%maxLevel ) then ! Perform the number of nested time steps on the finer level L+1 ! according to scaling type : ! diffusive: nNesting = 4 ! acoustic: nNesting = 2 do iNestingLoop = 1, params%nNesting write(logUnit(10), "(A,I0)") 'Nesting loop ', iNestingloop call me%do_computation( scheme, geometry, params, iLevel+1 ) end do end if write(logUnit(10), "(A,I0)") 'Compute on level: ', iLevel call start_stageTimer() ! update the time counters. MH: checked. Please dont move ! Increasing with the smallest time step (maxLevel) ! KM: time is advanced here since new time is required to update sources. if( iLevel == geometry%tree%global%maxLevel ) then write(logUnit(10), "(A)") 'Advance time t+dt_maxLevel' call tem_time_advance( me = params%general%simControl%now, & & sim_dt = params%physics%dtLvl( geometry & & %tree%global%maxLevel ) ) endif write(logUnit(10), "(A)") 'Set boundary condition' ! set boundary for each field in current scheme call set_boundary( field = scheme%field, & & pdf = scheme%pdf(iLevel), & & state = scheme%state(iLevel)%val, & & levelDesc = scheme%levelDesc(iLevel), & & tree = geometry%tree, & & iLevel = iLevel, & & nBCs = geometry%boundary%nBCtypes, & & params = params, & & layout = scheme%layout, & & physics = params%physics, & & varSys = scheme%varSys, & & mixture = scheme%mixture, & & derVarPos = scheme%derVarPos, & & globBC = scheme%globBC ) ! ------------------------------------------------------------------------- write(logUnit(10), "(A)") 'Swap now and next' ! swap double buffer index for current level call mus_swap_now_next( scheme%pdf( iLevel ) ) now = scheme%pdf(iLevel)%nNow next = scheme%pdf(iLevel)%nNext ! -------------------------------------------------------------------------- ! Compute auxField from pre-collision state for fluid and ghostFromCoarser ! and exchange them if turbulence is active call tem_startTimer( timerHandle = mus_timerHandles%aux(iLevel) ) write(logUnit(10), "(A)") 'Calculate auxField' call mus_calcAuxFieldAndExchange( & & auxField = scheme%auxField(iLevel), & & calcAuxField = scheme%calcAuxField, & & state = scheme%state(iLevel)%val(:, now), & & pdfData = scheme%pdf(iLevel), & & nFields = scheme%nFields, & & field = scheme%field(:), & & globSrc = scheme%globSrc, & & stencil = scheme%layout%fStencil, & & varSys = scheme%varSys, & & derVarPos = scheme%derVarPos, & & general = params%general, & & phyConvFac = params%physics%fac(iLevel), & & iLevel = iLevel, & & minLevel = geometry%tree%global%minLevel, & & schemeHeader = scheme%header ) if (iLevel < geometry%tree%global%maxLevel) then write(logUnit(10), "(A)") 'Interpolate and exchange auxField in ' & & //'ghostFromFiner' call mus_intpAuxFieldCoarserAndExchange( & & intp = scheme%intp, & & tAuxField = scheme%auxField(iLevel), & & sAuxField = scheme%auxField(iLevel+1), & & tLevelDesc = scheme%levelDesc(iLevel), & & stencil = scheme%layout%fStencil, & & iLevel = iLevel, & & nAuxScalars = scheme%varSys%nAuxScalars, & & general = params%general ) end if call tem_stopTimer( timerHandle = mus_timerHandles%aux(iLevel) ) ! -------------------------------------------------------------------------- write(logUnit(10), "(A)") 'Update relaxparams' ! -------------------------------------------------------------------------- ! Update parameters, relaxation time .etc call tem_startTimer( timerHandle = mus_timerHandles%relax(iLevel) ) call mus_update_relaxParams( scheme = scheme, & & iLevel = iLevel, & & tNow = params%general%simControl%now, & & physics = params%physics, & & lattice = params%lattice, & & nBCs = geometry%boundary%nBCtypes ) call tem_stopTimer( timerHandle = mus_timerHandles%relax(iLevel) ) ! -------------------------------------------------------------------------- ! ------------------------------------------------------------------------- write(logUnit(10), "(A)") 'Stream and collide' ! Compute current scheme of current level call tem_startTimer( timerHandle = mus_timerHandles%compute(iLevel) ) !$omp parallel call scheme%compute( & & fieldProp = scheme%field(:)%fieldProp, & & inState = scheme%state(iLevel)%val(:,Now), & & outState = scheme%state(iLevel)%val(:,Next), & & auxField = scheme%auxField(ilevel)%val(:), & & neigh = scheme%pdf(iLevel)%neigh(:), & & nElems = scheme%pdf(iLevel)%nSize, & & nSolve = scheme%pdf(iLevel)%nElems_solve, & & level = iLevel, & & layout = scheme%layout, & & params = params, & & derVarPos = scheme%derVarPos, & & varSys = scheme%varSys ) !$omp end parallel call tem_stopTimer( timerHandle = mus_timerHandles%compute(iLevel) ) ! ------------------------------------------------------------------------- ! -------------------------------------------------------------------------- write(logUnit(10), "(A)") 'Apply source' call mus_apply_sourceTerms( field = scheme%field(:), & & nFields = scheme%nFields, & & globSrc = scheme%globSrc, & & pdf = scheme%pdf(iLevel), & & varSys = scheme%varSys, & & iLevel = iLevel, & & time = params%general%simControl%now, & & state = scheme%state(iLevel)%val, & & auxField = scheme%auxField(iLevel)%val(:), & & derVarPos = scheme%derVarPos(:), & & phyConvFac = params%physics%fac(iLevel) ) ! -------------------------------------------------------------------------- ! ------------------------------------------------------------------------- write(logUnit(10), "(A)") 'Communicate fluids' ! Communicate the halo elements of each scheme on current level call tem_startTimer( timerHandle = mus_timerHandles%comm(iLevel) ) ! communicate halo elements for Next call params%general%commPattern%exchange_real( & & send = scheme%levelDesc(iLevel)%sendbuffer, & & recv = scheme%levelDesc(iLevel)%recvbuffer, & & state = scheme%state(iLevel)%val(:,Next), & & message_flag = iLevel, & & comm = params%general%proc%comm ) ! communicate turbulent viscosity, required for interpolation if (trim(scheme%header%kind) == 'fluid' .or. & & trim(scheme%header%kind) == 'fluid_incompressible') then if (scheme%field(1)%fieldProp%fluid%turbulence%active) then call params%general%commPattern%exchange_real( & & recv = scheme%field(1)%fieldProp%fluid%turbulence & & %dataOnLvl(iLevel)%recvbuffer, & & send = scheme%field(1)%fieldProp%fluid%turbulence & & %dataOnLvl(iLevel)%sendbuffer, & & state = scheme%field(1)%fieldProp%fluid%turbulence & & %dataOnLvl(iLevel)%visc(:), & & message_flag = iLevel+100, & & comm = params%general%proc%comm ) end if end if call tem_stopTimer( timerHandle = mus_timerHandles%comm(iLevel) ) ! ------------------------------------------------------------------------- ! ------------------------------------------------------------------------- ! communicate my finer ghost element from coarse level if ( iLevel > geometry%tree%global%minLevel ) then write(logUnit(10), "(A)") 'Communicate ghostFromCoarser' call tem_startTimer( timerHandle & & = mus_timerHandles%commFromCoarser(iLevel) ) call params%general%commPattern%exchange_real( & & send = scheme%levelDesc(iLevel)%sendbufferFromCoarser, & & recv = scheme%levelDesc(iLevel)%recvbufferFromCoarser, & & state = scheme%state(iLevel)%val(:, Next), & & message_flag = iLevel, & & comm = params%general%proc%comm ) ! communicate turbulent viscosity, required for interpolation if (trim(scheme%header%kind) == 'fluid' .or. & & trim(scheme%header%kind) == 'fluid_incompressible') then if (scheme%field(1)%fieldProp%fluid%turbulence%active) then call params%general%commPattern%exchange_real( & & recv = scheme%field(1)%fieldProp%fluid%turbulence & & %dataOnLvl(iLevel)%recvBufferFromCoarser, & & send = scheme%field(1)%fieldProp%fluid%turbulence & & %dataOnLvl(iLevel)%sendBufferFromCoarser, & & state = scheme%field(1)%fieldProp%fluid%turbulence & & %dataOnLvl(iLevel)%visc(:), & & message_flag = iLevel+200, & & comm = params%general%proc%comm ) end if end if call tem_stopTimer( timerHandle & & = mus_timerHandles%commFromCoarser(iLevel) ) end if ! -------------------------------------------------------------------------- call stop_stageTimer() ! Interpolate ghost elements if( iLevel < geometry%tree%global%maxLevel ) then ! Fill my coarser element (L) from finer (L+1) call do_intpFinerAndExchange( scheme, params, iLevel ) ! Interpolate the ghost elements on the finer level(L+1) with data provided ! from current level(L). call do_intpCoarserAndExchange( scheme, params, iLevel ) end if ! if not on finest level ! -------------------------------------------------------------------------- ! -------------------------------------------------------------------------- if( iLevel == geometry%tree%global%minLevel ) then iStage = 0 running = .false. end if ! -------------------------------------------------------------------------- end subroutine do_recursive_multiLevel ! ------------------------------------------------------------------------ ! ! ------------------------------------------------------------------------ ! !> This routine does: !! 1. interpolate my coarse ghost element (iLevel) from finer level (iLevel+1) !! 2. exchange the data of my coarse ghost elements between process subroutine do_intpFinerAndExchange( scheme, params, iLevel ) ! -------------------------------------------------------------------- ! !> containers for the different schemes type( mus_scheme_type ), intent(inout), target :: scheme !> global parameters type( mus_param_type ), intent(inout) :: params integer, intent(in) :: iLevel !< Level counter variable ! -------------------------------------------------------------------- ! integer :: thisLevelNext, nextLevelNext ! -------------------------------------------------------------------- ! ! Interpolate my coarse ghost elements from the finer ones thisLevelNext = scheme%pdf(iLevel)%nNext nextLevelNext = scheme%pdf(iLevel+1)%nNext call start_stageTimer() ! write(logUnit(5), "(A,I0)") 'intpFromFiner on level: ', iLevel ! ------------------------------------------------------------------------- call tem_startTimer( timerHandle = mus_timerHandles%intpFromFiner(iLevel) ) ! ------------------------------------------------------------------------- write(logUnit(10), "(A)") 'Interpolate ghostFromFiner' !$omp parallel call scheme%intp%fillMineFromFiner%do_intp( & & fieldProp = scheme%field(:)%fieldProp, & & sState = scheme%state(iLevel+1)%val(:,nextLevelNext), & & sNeigh = scheme%pdf(iLevel+1)%Neigh, & & snSize = scheme%pdf(iLevel+1)%nSize, & & sAuxField = scheme%auxField(iLevel+1)%val(:), & & tnSize = scheme%pdf(iLevel)%nSize, & & tState = scheme%state(iLevel)%val(:,thisLevelNext), & & tNeigh = scheme%pdf(iLevel)%Neigh, & & tLevelDesc = scheme%levelDesc(iLevel), & & level = iLevel, & & nTargets = scheme%levelDesc(iLevel)%intpFromFiner%nVals, & & targetList = scheme%levelDesc(iLevel)%intpFromFiner%val, & & layout = scheme%layout, & & physics = params%physics, & & varSys = scheme%varSys, & & derVarPos = scheme%derVarPos(:), & & time = params%general%simControl%now ) !$omp end parallel call tem_stopTimer( timerHandle = mus_timerHandles%intpFromFiner(iLevel) ) ! ------------------------------------------------------------------------- ! if ( main_debug%checkEachAlgorithmicStep ) then ! buffer=' after fill Mine FromFiner ' ! call dump_debug_info( buffer, scheme, params, iLevel, 1, & ! & pdf=scheme%pdf(iLevel) ) ! end if call stop_stageTimer() write(logUnit(10), "(A)") 'Communicate ghostFromFiner' ! ------------------------------------------------------------------------- ! Exchange the coarse ghost elements call tem_startTimer( timerHandle = mus_timerHandles%commFromFiner(iLevel) ) call params%general%commPattern%exchange_real( & & send = scheme%levelDesc(iLevel)%sendbufferFromFiner, & & recv = scheme%levelDesc(iLevel)%recvbufferFromFiner, & & state = scheme%state(iLevel)%val(:,thisLevelNext), & & message_flag = iLevel, & & comm = params%general%proc%comm ) if (trim(scheme%header%kind) == 'fluid' .or. & & trim(scheme%header%kind) == 'fluid_incompressible') then if (scheme%field(1)%fieldProp%fluid%turbulence%active) then ! communicate turbulent viscosity, required for interpolation call params%general%commPattern%exchange_real( & & recv = scheme%field(1)%fieldProp%fluid%turbulence & & %dataOnLvl(iLevel)%recvBufferFromFiner, & & send = scheme%field(1)%fieldProp%fluid%turbulence & & %dataOnLvl(iLevel)%sendBufferFromFiner, & & state = scheme%field(1)%fieldProp%fluid%turbulence & & %dataOnLvl(iLevel)%visc(:), & & message_flag = iLevel+200, & & comm = params%general%proc%comm ) end if end if call tem_stopTimer( timerHandle = mus_timerHandles%commFromFiner(iLevel) ) ! ------------------------------------------------------------------------- ! if ( main_debug%checkEachAlgorithmicStep ) then ! buffer = ' after exch Mine FromFiner' ! call dump_debug_info( buffer, scheme, params, iLevel, 1, & ! & pdf=scheme%pdf(iLevel) ) ! end if end subroutine do_intpFinerandExchange ! ------------------------------------------------------------------------ ! ! ------------------------------------------------------------------------ ! !> This routine utilizes fluid elements on my level (L) to fill finer !! ghost elements on next level (L+1). !! Then it exchanges the datas of finer ghost elements (L+1) between process. subroutine do_intpCoarserAndExchange( scheme, params, iLevel ) ! -------------------------------------------------------------------- ! !> containers for the different schemes type( mus_scheme_type ), target, intent(inout) :: scheme !> global parameters type( mus_param_type ), intent(in) :: params !> Level counter variable integer, intent(in) :: iLevel ! -------------------------------------------------------------------- ! integer :: thisLevelNext, nextLevelNext integer :: iOrder ! -------------------------------------------------------------------- ! thisLevelNext = scheme%pdf(iLevel)%nNext nextLevelNext = scheme%pdf(iLevel+1)%nNext call start_stageTimer() ! write(logUnit(5), "(A,I0)") 'intpFromCoarser on level: ', iLevel+1 ! Fill finer ghost elements (L+1) by low order interpolation call tem_startTimer( & & timerHandle = mus_timerHandles%intpFromCoarser(iLevel+1) ) ! Fill finer ghost elements (L+1) by starting from lower order ! to higher order interpolation write(logUnit(10), "(A)") 'Interpolate ghostFromCoarser iL+1' !$omp parallel do iOrder = 0, scheme%intp%config%order call scheme%intp%fillFinerFromMe(iOrder)%do_intp( & & fieldProp = scheme%field(:)%fieldProp, & & sState = scheme%state(iLevel)%val(:,thisLevelNext), & & sNeigh = scheme%pdf(iLevel)%Neigh, & & snSize = scheme%pdf(iLevel)%nSize, & & sAuxField = scheme%auxField(iLevel)%val(:), & & tnSize = scheme%pdf(iLevel+1)%nSize, & & tState = scheme%state(iLevel+1)%val(:,nextLevelNext), & & tNeigh = scheme%pdf(iLevel+1)%Neigh, & & tLevelDesc = scheme%levelDesc(iLevel+1), & & level = iLevel, & & nTargets = scheme%levelDesc(iLevel+1) & & %intpFromCoarser(iOrder)%nVals, & & targetList = scheme%levelDesc(iLevel+1) & & %intpFromCoarser(iOrder)%Val, & & layout = scheme%layout, & & physics = params%physics, & & varSys = scheme%varSys, & & derVarPos = scheme%derVarPos(:), & & time = params%general%simControl%now ) end do !$omp end parallel call tem_stopTimer( timerHandle & & = mus_timerHandles%intpFromCoarser(iLevel+1) ) ! Debug output ! if ( main_debug%checkEachAlgorithmicStep ) then ! buffer = ' after fillFinerFromMe' ! call dump_debug_info( buffer, scheme, params, iLevel, 1, & ! & pdf = scheme%pdf(iLevel)) ! end if call stop_stageTimer() ! Exchange the fine ghost elements write(logUnit(10), "(A)") 'Communicate ghostFromCoarser iL+1' call tem_startTimer( & & timerHandle = mus_timerHandles%commFromCoarser(iLevel+1) ) call params%general%commPattern%exchange_real( & & send = scheme%levelDesc(iLevel+1)%SendBufferFromCoarser, & & recv = scheme%levelDesc(iLevel+1)%RecvBufferFromCoarser, & & state = scheme%state(iLevel+1)%val(:,nextLevelNext), & & message_flag = iLevel, & & comm = params%general%proc%comm ) if (trim(scheme%header%kind) == 'fluid' .or. & & trim(scheme%header%kind) == 'fluid_incompressible') then if (scheme%field(1)%fieldProp%fluid%turbulence%active) then ! communicate turbulent viscosity, required for interpolation call params%general%commPattern%exchange_real( & & recv = scheme%field(1)%fieldProp%fluid%turbulence & & %dataOnLvl(iLevel+1)%recvBufferFromCoarser, & & send = scheme%field(1)%fieldProp%fluid%turbulence & & %dataOnLvl(iLevel+1)%sendBufferFromCoarser, & & state = scheme%field(1)%fieldProp%fluid%turbulence & & %dataOnLvl(iLevel+1)%visc(:), & & message_flag = iLevel+200, & & comm = params%general%proc%comm ) end if end if call tem_stopTimer( & & timerHandle = mus_timerHandles%commFromCoarser(iLevel+1) ) ! Debug output ! if ( main_debug%checkEachAlgorithmicStep ) then ! buffer = ' after exchFinerFromMe' ! call dump_debug_info( buffer, scheme, params, iLevel, 1, & ! & pdf = scheme%pdf(iLevel) ) ! end if end subroutine do_intpCoarserandExchange ! ------------------------------------------------------------------------ ! ! ------------------------------------------------------------------------ ! !> Control routine for an optimized workflow with reduced functionality. !! !! No sources, no multilevel, no multiLevel. !! Use for benchmarking !! subroutine do_fast_singleLevel( me, scheme, geometry, params, iLevel ) ! -------------------------------------------------------------------- ! !> self control type !! dummy variable in this routine, required by interface class( mus_control_type ) :: me !> container for the scheme type( mus_scheme_type ), intent(inout) :: scheme !> geometry infomation type( mus_geom_type ), intent(inout) :: geometry !> global parameters type( mus_param_type ), intent(inout) :: params !> Level counter variable integer, intent(in) :: iLevel ! -------------------------------------------------------------------- ! integer :: now, next ! -------------------------------------------------------------------- ! ! Update auxField dependent source fields before adding source term to state ! and auxField such that both auxField and apply_source uses same source. call mus_update_sourceVars( nFields = scheme%nFields, & & field = scheme%field, & & globSrc = scheme%globSrc, & & varSys = scheme%varSys, & & iLevel = iLevel, & & auxField = scheme%auxField(iLevel)%val, & & phyConvFac = params%physics%fac(iLevel), & & derVarPos = scheme%derVarPos ) ! ------------------------------------------------------------------------- ! Increasing with the smallest time step (maxLevel) ! KM: time is advanced here since new time is required to update sources and BCs call tem_time_advance( me = params%general%simControl%now, & & sim_dt = params%physics%dtLvl(iLevel )) ! -------------------------------------------------------------------------- !set boundary for each field in current scheme call set_boundary( field = scheme%field, & & pdf = scheme%pdf(iLevel), & & state = scheme%state(iLevel)%val, & & levelDesc = scheme%levelDesc(iLevel), & & tree = geometry%tree, & & iLevel = iLevel, & & nBCs = geometry%boundary%nBCtypes, & & params = params, & & layout = scheme%layout, & & physics = params%physics, & & varSys = scheme%varSys, & & mixture = scheme%mixture, & & derVarPos = scheme%derVarPos, & & globBC = scheme%globBC ) ! -------------------------------------------------------------------------- ! swap double buffer index for current level call mus_swap_now_next( scheme%pdf( iLevel ) ) now = scheme%pdf(iLevel)%nNow next = scheme%pdf(iLevel)%nNext ! -------------------------------------------------------------------------- ! Compute auxField from pre-collision state for fluid and ghostFromCoarser ! and exchange them if turbulence is active call tem_startTimer( timerHandle = mus_timerHandles%aux(iLevel) ) call mus_calcAuxFieldAndExchange( & & auxField = scheme%auxField(iLevel), & & calcAuxField = scheme%calcAuxField, & & state = scheme%state(iLevel)%val(:, now), & & pdfData = scheme%pdf(iLevel), & & nFields = scheme%nFields, & & field = scheme%field(:), & & globSrc = scheme%globSrc, & & stencil = scheme%layout%fStencil, & & varSys = scheme%varSys, & & derVarPos = scheme%derVarPos, & & general = params%general, & & phyConvFac = params%physics%fac(iLevel), & & iLevel = iLevel, & & minLevel = geometry%tree%global%minLevel, & & schemeHeader = scheme%header ) call tem_stopTimer( timerHandle = mus_timerHandles%aux(iLevel) ) ! -------------------------------------------------------------------------- ! -------------------------------------------------------------------------- ! Update parameters, relaxation time .etc call tem_startTimer( timerHandle = mus_timerHandles%relax(iLevel) ) call mus_update_relaxParams( scheme = scheme, & & iLevel = iLevel, & & tNow = params%general%simControl%now, & & physics = params%physics, & & lattice = params%lattice, & & nBCs = geometry%boundary%nBCtypes ) call tem_stopTimer( timerHandle = mus_timerHandles%relax(iLevel) ) ! -------------------------------------------------------------------------- ! ------------------------------------------------------------------------- ! Compute current scheme of current level call tem_startTimer( timerHandle = mus_timerHandles%compute(iLevel) ) !$omp parallel call scheme%compute( & & fieldProp = scheme%field(:)%fieldProp, & & inState = scheme%state(iLevel)%val(:,Now), & & outState = scheme%state(iLevel)%val(:,Next), & & auxField = scheme%auxField(ilevel)%val(:), & & neigh = scheme%pdf(iLevel)%neigh(:), & & nElems = scheme%pdf(iLevel)%nSize, & & nSolve = scheme%pdf(iLevel)%nElems_solve, & & level = iLevel, & & layout = scheme%layout, & & params = params, & & derVarPos = scheme%derVarPos, & & varSys = scheme%varSys ) !$omp end parallel call tem_stopTimer( timerHandle = mus_timerHandles%compute(iLevel) ) ! -------------------------------------------------------------------------- ! -------------------------------------------------------------------------- call mus_apply_sourceTerms( field = scheme%field(:), & & nFields = scheme%nFields, & & globSrc = scheme%globSrc, & & pdf = scheme%pdf(iLevel), & & varSys = scheme%varSys, & & iLevel = iLevel, & & time = params%general%simControl%now, & & state = scheme%state(iLevel)%val, & & auxField = scheme%auxField(iLevel)%val(:), & & derVarPos = scheme%derVarPos(:), & & phyConvFac = params%physics%fac(iLevel) ) ! ------------------------------------------------------------------------- ! Communicate the halo elements of each scheme on current level ! KM: Communicate post-collision before set_boundary because nonEq_expol ! BC depends on post-collision from neighbor at next time step call tem_startTimer( timerHandle = mus_timerHandles%comm(iLevel) ) call params%general%commPattern%exchange_real( & & send = scheme%levelDesc(iLevel)%sendbuffer, & & recv = scheme%levelDesc(iLevel)%recvbuffer, & & state = scheme%state(iLevel)%val(:,Next), & & message_flag = iLevel, & & comm = params%general%proc%comm ) ! communicate turbulent viscosity, required for interpolation if (trim(scheme%header%kind) == 'fluid' .or. & & trim(scheme%header%kind) == 'fluid_incompressible') then if (scheme%field(1)%fieldProp%fluid%turbulence%active) then call params%general%commPattern%exchange_real( & & recv = scheme%field(1)%fieldProp%fluid%turbulence & & %dataOnLvl(iLevel)%recvbuffer, & & send = scheme%field(1)%fieldProp%fluid%turbulence & & %dataOnLvl(iLevel)%sendbuffer, & & state = scheme%field(1)%fieldProp%fluid%turbulence & & %dataOnLvl(iLevel)%visc(:), & & message_flag = iLevel+100, & & comm = params%general%proc%comm ) end if end if call tem_stopTimer( timerHandle = mus_timerHandles%comm(iLevel) ) ! ... check if at least one of the IBMs is active if ( geometry%globIBM%nIBMs > 0 ) then call mus_buildBuffIBM( & & me = geometry%globIBM%IBM, & & commPattern = params%general%commPattern, & & globTree = geometry%tree, & & params = params, & & layout = scheme%layout, & & levelDesc = scheme%levelDesc(iLevel), & & iLevel = iLevel ) end if ! update the immersed boundaries if available ! ... and over the schemes ! ... check if at least one of the IBMs is active if( geometry%globIBM%nIBMs > 0 )then call mus_inamuro_IBM( & & me = geometry%globIBM%IBM, & & commPattern = params%general%commPattern, & & globTree = geometry%tree, & & general = params%general, & & pdf = scheme%pdf(iLevel), & & layout = scheme%layout, & & levelDesc = scheme%levelDesc(iLevel), & & globSys = scheme%varSys, & & stateVarMap = scheme%stateVarMap%varPos%val(:), & & convFac = params%physics%fac(iLevel), & & iField = 1, & & state = scheme%state(iLevel)%val, & & iLevel = iLevel ) end if ! ------------------------------------------------------------------------- end subroutine do_fast_singleLevel ! ------------------------------------------------------------------------ ! ! ------------------------------------------------------------------------ ! subroutine do_benchmark( me, scheme, geometry, params, iLevel ) ! -------------------------------------------------------------------- ! !> self control type !! dummy variable in this routine, required by interface class( mus_control_type ) :: me !> containers for the different schemes type( mus_scheme_type ), intent(inout) :: scheme !> geometry infomation type( mus_geom_type ), intent(inout) :: geometry !> global parameters type( mus_param_type ), intent(inout) :: params !> Level counter variable integer, intent(in) :: iLevel ! -------------------------------------------------------------------- ! integer :: now, next ! -------------------------------------------------------------------- ! ! Update auxField dependent source fields before adding source term to state ! and auxField such that both auxField and apply_source uses same source. call mus_update_sourceVars( nFields = scheme%nFields, & & field = scheme%field, & & globSrc = scheme%globSrc, & & varSys = scheme%varSys, & & iLevel = iLevel, & & auxField = scheme%auxField(iLevel)%val, & & phyConvFac = params%physics%fac(iLevel), & & derVarPos = scheme%derVarPos ) ! Increasing with the smallest time step (maxLevel) call tem_time_advance( me = params%general%simControl%now, & & sim_dt = params%physics%dtLvl(iLevel )) ! -------------------------------------------------------------------------- !set boundary for each field in current scheme call set_boundary( field = scheme%field, & & pdf = scheme%pdf(iLevel), & & state = scheme%state(iLevel)%val, & & levelDesc = scheme%levelDesc(iLevel), & & tree = geometry%tree, & & iLevel = iLevel, & & nBCs = geometry%boundary%nBCtypes, & & params = params, & & layout = scheme%layout, & & physics = params%physics, & & varSys = scheme%varSys, & & mixture = scheme%mixture, & & derVarPos = scheme%derVarPos, & & globBC = scheme%globBC ) ! -------------------------------------------------------------------------- ! swap double buffer index for current level call mus_swap_now_next( scheme%pdf( iLevel ) ) now = scheme%pdf(iLevel)%nNow next = scheme%pdf(iLevel)%nNext ! -------------------------------------------------------------------------- ! Compute auxField from pre-collision state for fluid and ghostFromCoarser ! and exchange them if turbulence is active call tem_startTimer( timerHandle = mus_timerHandles%aux(iLevel) ) call mus_calcAuxFieldAndExchange( & & auxField = scheme%auxField(iLevel), & & calcAuxField = scheme%calcAuxField, & & state = scheme%state(iLevel)%val(:, now), & & pdfData = scheme%pdf(iLevel), & & nFields = scheme%nFields, & & field = scheme%field(:), & & globSrc = scheme%globSrc, & & stencil = scheme%layout%fStencil, & & varSys = scheme%varSys, & & derVarPos = scheme%derVarPos, & & general = params%general, & & phyConvFac = params%physics%fac(iLevel), & & iLevel = iLevel, & & minLevel = geometry%tree%global%minLevel, & & schemeHeader = scheme%header ) call tem_stopTimer( timerHandle = mus_timerHandles%aux(iLevel) ) ! -------------------------------------------------------------------------- ! -------------------------------------------------------------------------- ! Update parameters, relaxation time .etc call mus_update_relaxParams( scheme = scheme, & & iLevel = iLevel, & & tNow = params%general%simControl%now, & & physics = params%physics, & & lattice = params%lattice, & & nBCs = geometry%boundary%nBCtypes ) ! -------------------------------------------------------------------------- ! -------------------------------------------------------------------------- ! Compute current scheme of current level call tem_startTimer( timerHandle = mus_timerHandles%compute(iLevel) ) !$omp parallel call scheme%compute( & & fieldProp = scheme%field(:)%fieldProp, & & inState = scheme%state(iLevel)%val(:,Now), & & outState = scheme%state(iLevel)%val(:,Next), & & auxField = scheme%auxField(ilevel)%val(:), & & neigh = scheme%pdf(iLevel)%neigh(:), & & nElems = scheme%pdf(iLevel)%nSize, & & nSolve = scheme%pdf(iLevel)%nElems_solve, & & level = iLevel, & & layout = scheme%layout, & & params = params, & & derVarPos = scheme%derVarPos, & & varSys = scheme%varSys ) !$omp end parallel call tem_stopTimer( timerHandle = mus_timerHandles%compute(iLevel) ) ! -------------------------------------------------------------------------- ! -------------------------------------------------------------------------- call mus_apply_sourceTerms( field = scheme%field(:), & & nFields = scheme%nFields, & & globSrc = scheme%globSrc, & & pdf = scheme%pdf(iLevel), & & varSys = scheme%varSys, & & iLevel = iLevel, & & time = params%general%simControl%now, & & state = scheme%state(iLevel)%val, & & auxField = scheme%auxField(iLevel)%val(:), & & derVarPos = scheme%derVarPos(:), & & phyConvFac = params%physics%fac(iLevel) ) ! -------------------------------------------------------------------------- ! Communicate the halo elements of each scheme on current level call tem_startTimer( timerHandle = mus_timerHandles%comm(iLevel) ) call params%general%commPattern%exchange_real( & & send = scheme%levelDesc(iLevel)%sendbuffer, & & recv = scheme%levelDesc(iLevel)%recvbuffer, & & state = scheme%state(iLevel)%val(:,Next), & & message_flag = iLevel, & & comm = params%general%proc%comm ) ! communicate turbulent viscosity, required for interpolation if (trim(scheme%header%kind) == 'fluid' .or. & & trim(scheme%header%kind) == 'fluid_incompressible') then if (scheme%field(1)%fieldProp%fluid%turbulence%active) then call params%general%commPattern%exchange_real( & & recv = scheme%field(1)%fieldProp%fluid%turbulence & & %dataOnLvl(iLevel)%recvbuffer, & & send = scheme%field(1)%fieldProp%fluid%turbulence & & %dataOnLvl(iLevel)%sendbuffer, & & state = scheme%field(1)%fieldProp%fluid%turbulence & & %dataOnLvl(iLevel)%visc(:), & & message_flag = iLevel+100, & & comm = params%general%proc%comm ) end if end if call tem_stopTimer( timerHandle = mus_timerHandles%comm(iLevel) ) ! ------------------------------------------------------------------------- end subroutine do_benchmark ! ------------------------------------------------------------------------ ! ! ------------------------------------------------------------------------ ! subroutine start_stageTimer( ) if ( .not. running ) then iStage = mod( iStage, nStages ) + 1 call tem_startTimer( timerHandle = mus_timerHandles%stage(iStage) ) running = .true. end if end subroutine start_stageTimer ! ------------------------------------------------------------------------ ! ! ------------------------------------------------------------------------ ! subroutine stop_stageTimer( ) if ( running ) then call tem_stopTimer( timerHandle = mus_timerHandles%stage(iStage) ) running = .false. end if end subroutine stop_stageTimer ! ------------------------------------------------------------------------ ! end module mus_control_module ! ***************************************************************************** !