TABLE OF CONTENTS
- 1. NUM/rk2
NUM/rk2 [ Modules ]
NOM
rk2(calc)
DESCRIPTION
Integration en temps, par la methode de Runge-Kutta d ordre 1 ou 2 des bilans des flux.
Integration in time, by the Runge-Kutta method of order 1 or 2 of flux balances.
ENTREES / INPUT
calc : objet calcul / calculation object
SORTIES / OUTPUT
calc : objet calcul mis a jour / updated calculation object
SOURCE
! Ce logiciel est regi par la licence [CeCILL-B] ! This software is governed by the [CeCILL-B] license !=========================== DEBUT DES DECLARATIONS ==================== !.1----- Implicit, Use use num_typ use num ,only: dtloc, synchro_mpi, u_penalisation,display,calc_flux, calc_source,synchro_w,sharp IMPLICIT NONE !.2----- Declaration type(calcul),intent(inout) :: calc !! Objet calcul / Calcul object integer :: i,nitmax real(kind=kind(0.d0)) :: pm,e,vol integer :: lev,ic1,ic2 !=========================== DEBUT DU CODE EXECUTABLE ================== ! level: niveau de raffinement a considerer. Ici lev=grand. lev=1000 ! Temps courant tc=calc%tmin ! Nbre de pas de temps max nitmax=100000 ! Boucle en temps DO WHILE (tc.lt.calc%tmax) ! Determination du pas de temps du calcul CALL dtloc(calc) ! Ajustement du pas de temps si necessaire au dernier pas de temps IF (tc+calc%dt.gt.calc%tmax) calc%dt=calc%tmax-tc ! Synchronisation du pas de temps et des cellules aux frontieres du domaine CALL synchro_mpi(calc) ! Calcul des gradients si necessaire IF (calc%iordre_e.eq.2) then CALL gradprim(calc%msh) SELECT CASE (calc%ilimiteur) CASE (0) ! barth CALL gradlim(calc%msh) CASE (1) ! cartesien CALL gradlimxy(calc%msh) END SELECT CALL synchro_w(calc) END IF ! Calcul de l entropie a l instant t DO i=1,calc%msh%nb_cell CALL entropie(calc%msh%list_cell(i)%w,calc%msh%list_cell(i)%w%enti) END DO !------- Initialisation du stockage des flux et sources DO i=1,calc%msh%nb_cell calc%msh%list_cell(i)%balance%v=0.d0 END DO !------- Calcul des flux aux faces (et eventuellement stockage) CALL calc_flux(calc,lev) !------- Calcul des termes source (et eventuellement stockage) CALL calc_source(calc,lev) tc=tc+calc%dt calc%ipas=calc%ipas+1 IF (calc%msh%numdom.eq.1) write( *,1000) calc%ipas,tc,calc%dt,calc%vmax 1000 format(' Ipas: ',i7,' T=',e12.5,' dt=',e12.5,' vmax=',e12.5) ! Mise a jour des variables conservatives pm=0.d0 DO i=1,calc%msh%nb_cell vol=calc%msh%list_cell(i)%vol calc%msh%list_cell(i)%w%vbal%v=calc%msh%list_cell(i)%w%vbal%v & - calc%dt/calc%iordre_t/vol * calc%msh%list_cell(i)%balance%v ! point milieu CALL Bal2Prim(calc%msh%list_cell(i)%w) IF (dabs(calc%dt/vol * calc%msh%list_cell(i)%balance%v(5)).gt.pm& !!!!!! PEUT-ON LE METTRE AILLEURS ? .and. (imodel_phy.eq.2)) pm=dabs(calc%dt/vol * calc%msh%list_cell(i)%balance%v(5)) END DO if(imodel_phy.eq.1)then calc%dt=calc%dt/calc%iordre_t call friction(calc) calc%dt=calc%dt*calc%iordre_t endif ! Correction par le point milieu ou heun si ordre 2 en temps IF (calc%iordre_t.eq.2) THEN IF (calc%iordre_e.eq.2) THEN CALL synchro_w(calc) CALL gradprim(calc%msh) SELECT CASE(calc%ilimiteur) CASE(0) ! barth CALL gradlim(calc%msh) CASE(1) ! cartesien CALL gradlimxy(calc%msh) END SELECT CALL synchro_w(calc) END IF ! On revient en arriere d un demi pas de temps (en variables conservatives) DO i=1,calc%msh%nb_cell calc%msh%list_cell(i)%balance%v=-0.5d0*calc%msh%list_cell(i)%balance%v IF (ive_phy.ne.0) calc%msh%list_cell(i)%balance%v(ive_phy)=0.d0 END DO !------- Calcul des flux aux faces (et eventuellement stockage) CALL calc_flux(calc,lev) !------- Calcul des termes source (et eventuellement stockage) CALL calc_source(calc,lev) pm=0.d0 DO i=1,calc%msh%nb_cell vol=calc%msh%list_cell(i)%vol calc%msh%list_cell(i)%w%vbal%v=calc%msh%list_cell(i)%w%vbal%v & - calc%dt/vol * calc%msh%list_cell(i)%balance%v CALL Bal2Prim(calc%msh%list_cell(i)%w) IF (dabs(calc%dt/vol * calc%msh%list_cell(i)%balance%v(5)).gt.pm& .and. imodel_phy.eq.2) pm=dabs(calc%dt/vol *& calc%msh%list_cell(i)%balance%v(5)) END DO if(imodel_phy.eq.1)then calc%dt=calc%dt/calc%iordre_t call friction(calc) calc%dt=calc%dt*calc%iordre_t endif END IF !------- Calcul de l entropie a l instant t+dt DO i=1,calc%msh%nb_cell CALL entropie(calc%msh%list_cell(i)%w,e) calc%msh%list_cell(i)%w%entf=e IF (ive_phy.ne.0) THEN calc%msh%list_cell(i)%w%pe=(e-calc%msh%list_cell(i)%w%enti)/calc%dt +& calc%msh%list_cell(i)%balance%v(ive_phy)/calc%msh%list_cell(i)%vol calc%msh%list_cell(i)%w%vprim%v(ive_phy)=calc%msh%list_cell(i)%w%pe END IF END DO ! dans le cas de saint venant determination de la shoreline et calcul de la friction if(imodel_phy.eq.1)then DO i=1,calc%msh%nb_cell calc%msh%list_cell(i)%w%chi=0.d0 END DO DO i=1,calc%msh%nb_face ic1=calc%msh%list_face(i)%ic1 ic2=calc%msh%list_face(i)%ic2 if(ic2.gt.0)then if(calc%msh%list_cell(ic1)%w%vprim%v(1).gt.hzero_phy.and.& calc%msh%list_cell(ic2)%w%vprim%v(1).lt.hzero_phy)calc%msh%list_cell(ic1)%w%chi=1.d0 if(calc%msh%list_cell(ic1)%w%vprim%v(1).lt.hzero_phy.and.& calc%msh%list_cell(ic2)%w%vprim%v(1).gt.hzero_phy)calc%msh%list_cell(ic2)%w%chi=1.d0 if(calc%msh%list_cell(ic1)%w%vprim%v(1).gt.hzero_phy.and.& calc%msh%list_cell(ic2)%w%vprim%v(1).lt.hzero_phy)calc%msh%list_cell(ic2)%w%chi=1.d0 if(calc%msh%list_cell(ic1)%w%vprim%v(1).lt.hzero_phy.and.& calc%msh%list_cell(ic2)%w%vprim%v(1).gt.hzero_phy)calc%msh%list_cell(ic1)%w%chi=1.d0 endif END DO !call friction(calc) endif ! Compression de l interface IF (imodel_phy.eq.2.and.sharp_phy.ne.0.d0) CALL sharp(calc,pm) ! Penalisation si obstacle IF (calc%nb_sol .NE. 0) CALL u_penalisation(calc) ! Verification conservation !call info_cons(calc) END DO !if(calc%msh%numdom.eq.1.and.calc%nb_sol .NE. 0)write(13,'(7e13.5)')tc,calc%list_sol(1)%xg%v CALL synchro_w(calc) !call display(calc%msh) !=========================== FIN DE LA ROUTINE ==================== END SUBROUTINE rk2