TABLE OF CONTENTS
GEO/gradlimxy.f90 [ Modules ]
NOM
gradlimxy(msh)
DESCRIPTION
Limite le gradient des variables primitives par la méthode de Barth La limitation est faite dans chaque direction x y indépendamment Suppose donc que le maillage est régulier en x y et 2D Suppose donc que le gradient a deja ete calcule (appel a gradprim) Attention: le gradient des cellules frontiere n est pas limite Il faut donc un appel à MPI pour mettre à jours TOUS les limiteurs Limit the gradient of primitive variables using Barth s method Limitation is done in each x y direction independently Assumes that the mesh is regular in x y and 2D. Assumes that the gradient has already been calculated (call gradprim). Warning: the gradient of the border cells is not limited A call to MPI is therefore required to update ALL the limiters. ENTREES / INPUT msh : maillage / mesh SORTIES / OUTPUT msh : maillage / mesh
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 geo_typ IMPLICIT NONE !.2----- Declaration type(mesh), intent(inout) :: msh !! Maillage / Mesh real(kind=kind(0.d0)),dimension(:,:),allocatable :: dlimx , dlimy ! tableau des limiteurs de pente dans chaque cellule real(kind=kind(0.d0)),dimension(nvarmax) :: odg integer :: i,iv,ic1,ic2,idom1,idom2,idom real(kind=kind(0.d0)) :: g,r,dlim2 type(vect_3) :: vcf type(state) :: w1,w2 !=========================== DEBUT DU CODE EXECUTABLE ================== allocate(dlimx(nvar,msh%nb_cell)) allocate(dlimy(nvar,msh%nb_cell)) dlimx=1 dlimy=1 ! vecteur des ordres de grandeur (permet d avoir un limiteur adimensionne) DO iv=1,nvar odg(iv)=1.d0 END DO odg(5)=1.d5 ! odg(1)=1.d0 ! boucle sur les faces du maillage DO i=1,msh%nb_face ! repérage des voisins ic1=msh%list_face(i)%ic1 ic2=msh%list_face(i)%ic2 w1=msh%list_cell(ic1)%w IF (ic2.gt.0) THEN w2=msh%list_cell(ic2)%w ELSE CALL Wext(w1,w2,msh%list_face(i)%vnorm,ic2, & msh%list_face(i)%center%x , & msh%list_face(i)%center%y , & msh%list_face(i)%center%z , & msh%list_cell(ic1)%center%x, & msh%list_cell(ic1)%center%y, & msh%list_cell(ic1)%center%z ) END IF ! numero de domaine à gauche idom1=msh%list_cell(ic1)%idom IF (idom1.ne.msh%numdom) THEN write(*,*) 'erreur maillage domaine' STOP END IF IF (ic2.ne.-999) THEN ! ne pas boucler sur les faces fictives ! numero de domaine à droite idom2=0 IF (ic2.gt.0) idom2=msh%list_cell(ic2)%idom IF (idom2.ne.msh%numdom) idom2=0 ! limitation cellule de gauche ! détermination du vecteur centre de la cellule - centre de la face vcf%v(1)=msh%list_face(i)%center%x-msh%list_cell(ic1)%center%x vcf%v(2)=msh%list_face(i)%center%y-msh%list_cell(ic1)%center%y vcf%v(3)=msh%list_face(i)%center%z-msh%list_cell(ic1)%center%z DO iv=1,nvar dlim2=1.d0 g=w1%gprim(1)%v(iv)*vcf%v(1)+& & w1%gprim(2)%v(iv)*vcf%v(2)+& & w1%gprim(3)%v(iv)*vcf%v(3) !IF (dabs(g/odg(iv)).gt.1.d-6) THEN IF (g.ne.0.d0) THEN ! moins compliqué et plus stable r=(w2%vprim%v(iv)-w1%vprim%v(iv))/g IF (r.lt.0.d0) THEN dlim2=0.d0 ELSE IF (r.gt.1.d0) THEN dlim2=1.d0 ELSE dlim2=r END IF ELSE dlim2=1.d0 END IF IF (dabs(msh%list_face(i)%vnorm%v(2)).lt.1.d-10) THEN ! direction x IF (dlimx(iv,ic1).gt.dlim2) dlimx(iv,ic1)=dlim2 END IF IF (dabs(msh%list_face(i)%vnorm%v(1)).lt.1.d-10) THEN ! direction y IF (dlimy(iv,ic1).gt.dlim2) dlimy(iv,ic1)=dlim2 END IF END DO ! limitation cellule de droite IF (idom2.ne.0) THEN ! détermination du vecteur centre de la cellule - centre de la face vcf%v(1)=msh%list_face(i)%center%x-msh%list_cell(ic2)%center%x vcf%v(2)=msh%list_face(i)%center%y-msh%list_cell(ic2)%center%y vcf%v(3)=msh%list_face(i)%center%z-msh%list_cell(ic2)%center%z DO iv=1,nvar dlim2=1.d0 g=w2%gprim(1)%v(iv)*vcf%v(1)+& & w2%gprim(2)%v(iv)*vcf%v(2)+& & w2%gprim(3)%v(iv)*vcf%v(3) IF (g.ne.0.d0) THEN ! moins compliqué et plus stable ! if (dabs(g/odg(iv)).gt.1.d-6) then ! si trop limité r=(w1%vprim%v(iv)-w2%vprim%v(iv))/g IF (r.lt.0.d0) THEN dlim2=0.d0 ELSE IF (r.gt.1.d0) THEN dlim2=1.d0 ELSE dlim2=r END IF ELSE dlim2=1.d0 END IF IF (dabs(msh%list_face(i)%vnorm%v(2)).lt.1.d-10) THEN ! direction x IF (dlimx(iv,ic2).gt.dlim2) dlimx(iv,ic2)=dlim2 END IF IF (dabs(msh%list_face(i)%vnorm%v(1)).lt.1.d-10) THEN ! direction y IF (dlimy(iv,ic2).gt.dlim2) dlimy(iv,ic2)=dlim2 END IF END DO END IF END IF END DO ! mise à jour du gradient dans les cellules DO i=1,msh%nb_cell idom=msh%list_cell(i)%idom IF (idom.eq.msh%numdom) THEN DO iv=1,nvar msh%list_cell(i)%w%gprim(1)%v(iv)=0.95*dlimx(iv,i)*msh%list_cell(i)%w%gprim(1)%v(iv) msh%list_cell(i)%w%gprim(2)%v(iv)=0.95*dlimy(iv,i)*msh%list_cell(i)%w%gprim(2)%v(iv) END DO END IF END DO deallocate(dlimx) deallocate(dlimy) !=========================== FIN DE LA ROUTINE ==================== END SUBROUTINE gradlimxy