TABLE OF CONTENTS
GEO/gradlim.f90 [ Modules ]
NOM
gradlim(msh)
DESCRIPTION
Limite le gradient des variables primitives par la methode de Barth 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 a MPI pour mettre a jours TOUS les limiteurs Limits the gradient of primitive variables using Barth s method Therefore assumes that the gradient has already been calculated (call gradprim) Warning: the gradient of 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 :: dlim ! tableau des limiteurs de pente dans chaque cellule integer :: i,iv,ic1,ic2,idom1,idom2,idom,ii,j real(kind=kind(0.d0)),dimension(nvarmax) :: odg real(kind=kind(0.d0)) :: g,r,dlim2 type(vect_3) :: vcf type(state) :: w1 type(state) :: w2 !=========================== DEBUT DU CODE EXECUTABLE ================== allocate(dlim(3,nvar,msh%nb_cell)) ! valeur min, valeur max, limiteur do i=1, msh%nb_cell do j=1, nvar dlim(1,j,i)=1.d20 !valeur min sur la cellule i de la variable j dlim(2,j,i)=-1.d20 !valeur max sur la cellule i de la variable j dlim(3,j,i)=1.d0 !limiteur de la variable j de la cellule i enddo enddo ! 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 !------------------------------------------------------------------------- ! Determination du min et du max entourant une cellule !------------------------------------------------------------------------- ! boucle sur les faces du maillage DO i=1,msh%nb_face 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 DO iv=1,nvar if(w2%vprim%v(iv).lt.dlim(1,iv,ic1))dlim(1,iv,ic1)=w2%vprim%v(iv) if(w2%vprim%v(iv).gt.dlim(2,iv,ic1))dlim(2,iv,ic1)=w2%vprim%v(iv) IF (ic2.gt.0)then if(w1%vprim%v(iv).lt.dlim(1,iv,ic2))dlim(1,iv,ic2)=w1%vprim%v(iv) if(w1%vprim%v(iv).gt.dlim(2,iv,ic2))dlim(2,iv,ic2)=w1%vprim%v(iv) endif END DO END DO DO i=1,msh%nb_cell w1=msh%list_cell(i)%w DO iv=1,nvar if(w1%vprim%v(iv).lt.dlim(1,iv,i))dlim(1,iv,i)=w1%vprim%v(iv) if(w1%vprim%v(iv).gt.dlim(2,iv,i))dlim(2,iv,i)=w1%vprim%v(iv) END DO END DO !------------------------------------------------------------------------- ! calcul du limiteur de pente de Barth-Jespersen !------------------------------------------------------------------------- ! boucle sur les faces du maillage DO i=1,msh%nb_face ! reperage 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 ! determination 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 r=0.d0 IF (g.ne.0.d0) THEN ! moins compliqué et plus stable if(g.gt.0.d0)r=(dlim(2,iv,ic1)-w1%vprim%v(iv))/g if(g.lt.0.d0)r=(dlim(1,iv,ic1)-w1%vprim%v(iv))/g IF (r.lt.0.d0) THEN dlim2=0.d0 !write(*,*)'r negatif',ic1,r ELSE IF (r.gt.1.d0) THEN dlim2=1.d0 !write(*,*)'r sup a 1',ic1,r ELSE dlim2=r END IF ELSE dlim2=1.d0 END IF IF (dlim(3,iv,ic1).gt.dlim2) dlim(3,iv,ic1)=dlim2 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) r=0.d0 IF (g.ne.0.d0) THEN ! moins compliqué et plus stable ! if (dabs(g/odg(iv)).gt.1.d-6) then ! si trop limité if(g.ge.0.d0)r=(dlim(2,iv,ic2)-w2%vprim%v(iv))/g if(g.lt.0.d0)r=(dlim(1,iv,ic2)-w2%vprim%v(iv))/g IF (r.lt.0.d0) THEN dlim2=0.d0 !write(*,*)'r negatif',ic2,r ELSE IF (r.gt.1.d0) THEN dlim2=1.d0 !write(*,*)'r sup a 1',ic2,r ELSE dlim2=r END IF ELSE dlim2=1.d0 END IF IF (dlim(3,iv,ic2).gt.dlim2) dlim(3,iv,ic2)=dlim2 END DO END IF END IF END DO ! Mise a 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 DO ii=1,msh%ndim msh%list_cell(i)%w%gprim(ii)%v(iv)=0.9d0*dlim(3,iv,i)*msh%list_cell(i)%w%gprim(ii)%v(iv) ! précaution en plus... if(msh%list_cell(i)%w%chi.ne.0.d0)msh%list_cell(i)%w%gprim(ii)%v(iv)=0.d0 END DO END DO END IF END DO deallocate(dlim) !=========================== FIN DE LA ROUTINE ==================== END SUBROUTINE GRADLIM