diff --git a/NavierStokes3D-Penalization/advection_mod.f90 b/NavierStokes3D-Penalization/advection_mod.f90
index 8b22ac4e4cfea2de6eb127324179b731361232a2..a9a22b23b1c4b560ff9417480e96ab982bb2d54e 100644
--- a/NavierStokes3D-Penalization/advection_mod.f90
+++ b/NavierStokes3D-Penalization/advection_mod.f90
@@ -260,7 +260,7 @@ contains
r2=kx**2+ky**2+kz**2
coef=-fac*cmplx(0.,1.)/(2.*pi*r2)
- uinfx=4./(2.4*2.4)!surface domaine physique/surface domaine numerique
+ uinfx=1/((yh-yb)*(zf-zd))!debit/surface du domaine numérique
!if ((time>=3.).and.(time<=4.)) then
! uinfy=sin(pi*(time-3.))
!else
diff --git a/NavierStokes3D-Penalization/init_mod.f90 b/NavierStokes3D-Penalization/init_mod.f90
index 2298a303d7b6e7bea2c070d65f0b0fc044625ae6..428cbf01a3db909cf2e781cc322cc2820afdeae3 100644
--- a/NavierStokes3D-Penalization/init_mod.f90
+++ b/NavierStokes3D-Penalization/init_mod.f90
@@ -16,28 +16,29 @@ contains
!---------------
! 0
!---------------
- omgx=0.
- omgy=0.
- omgz=0.
-!!$ !-------------------------
-!!$ ! test poiseuille
-!!$ !--------------------------
-!!$
-!!$
-!!$ do k=1,nz
-!!$ z=ztab(k)
-!!$ do j=1,ny
-!!$ do i=1,nx
-!!$ omgx(i,j,k)=0.
-!!$ if ( (z.ge.-1.).and.(z.le.1.)) then
-!!$ omgy(i,j,k)=-3.*z/(1.)**2
-!!$ else
-!!$ omgy(i,j,k)=0.
-!!$ end if
-!!$ omgz(i,j,k)=0.
-!!$ end do
-!!$ end do
-!!$ end do
+!!$ omgx=0.
+!!$ omgy=0.
+!!$ omgz=0.
+ !-------------------------
+ ! test poiseuille
+ !--------------------------
+
+
+ do k=1,nz
+ z=ztab(k)
+ do j=1,ny
+ do i=1,nx
+ omgx(i,j,k)=0.
+ omgy(i,j,k)=0.
+ if ( (z.ge.-1.).and.(z.le.1.)) then
+ omgy(i,j,k)=-3.*z/(4.*0.5*(yh-yb))
+ else
+ omgy(i,j,k)=0.
+ end if
+ omgz(i,j,k)=0.
+ end do
+ end do
+ end do
!!$
!!$
@@ -103,47 +104,47 @@ contains
- if (num_suite==1) then
- open(unit=25,file=trim(name_relance),form="formatted")
-
- read(25,'(A26)')text_tmp
- read(25,'(A7)')text_tmp
- read(25,'(A5)')text_tmp
- read(25,'(A25)')text_tmp
- read(25,'(A10,3(i4,1x))')text_tmp,i_tmp,i_tmp,i_tmp
- read(25,'(a6,3(f10.5))')text_tmp,r_tmp,r_tmp,r_tmp
- read(25,'(A7,3(f10.5))')text_tmp,r_tmp,r_tmp,r_tmp
-
-
- read(25,'(A10,i10)') text_tmp,i_tmp
- read(25,'(A21)') text_tmp
-
- do k=1,nz
- do j=1,ny
- do i=1,nx
- read(25,'(3(f20.9))') romgx,romgy,romgz
- omgx(i,j,k)=romgx
- omgy(i,j,k)=romgy
- omgz(i,j,k)=romgz
- end do
- end do
- end do
- close(25)
- end if
+!!$ if (num_suite==1) then
+!!$ open(unit=25,file=trim(name_relance),form="formatted")
+!!$
+!!$ read(25,'(A26)')text_tmp
+!!$ read(25,'(A7)')text_tmp
+!!$ read(25,'(A5)')text_tmp
+!!$ read(25,'(A25)')text_tmp
+!!$ read(25,'(A10,3(i4,1x))')text_tmp,i_tmp,i_tmp,i_tmp
+!!$ read(25,'(a6,3(f10.5))')text_tmp,r_tmp,r_tmp,r_tmp
+!!$ read(25,'(A7,3(f10.5))')text_tmp,r_tmp,r_tmp,r_tmp
+!!$
+!!$
+!!$ read(25,'(A10,i10)') text_tmp,i_tmp
+!!$ read(25,'(A21)') text_tmp
+!!$
+!!$ do k=1,nz
+!!$ do j=1,ny
+!!$ do i=1,nx
+!!$ read(25,'(3(f20.9))') romgx,romgy,romgz
+!!$ omgx(i,j,k)=romgx
+!!$ omgy(i,j,k)=romgy
+!!$ omgz(i,j,k)=romgz
+!!$ end do
+!!$ end do
+!!$ end do
+!!$ close(25)
+!!$ end if
-
- do k=1,nz
- do j=1,ny
- do i=1,nx
- write(11,'(3(f20.9))') real(omgx(i,j,k)),real(omgy(i,j,k)),real(omgz(i,j,k))
- end do
- end do
- end do
+!!$
+!!$ do k=1,nz
+!!$ do j=1,ny
+!!$ do i=1,nx
+!!$ write(11,'(3(f20.9))') real(omgx(i,j,k)),real(omgy(i,j,k)),real(omgz(i,j,k))
+!!$ end do
+!!$ end do
+!!$ end do
diff --git a/NavierStokes3D-Penalization/main.f90 b/NavierStokes3D-Penalization/main.f90
index f7227078d2c2abc914c25e7ba51725d30fbbda1f..e08bd390d2d574282a2d67c72f77104dced582e1 100644
--- a/NavierStokes3D-Penalization/main.f90
+++ b/NavierStokes3D-Penalization/main.f90
@@ -58,8 +58,8 @@ program NS
!maxv=1.D0
- maxv=2.D0
- cfl=0.4D0
+ !maxv=2.D0
+ !cfl=0.4D0
!dt=cfl*dx/maxv
@@ -117,6 +117,14 @@ program NS
allocate(stx(1:npart),sty(1:npart),stz(1:npart))
open(unit=13,file=trim(nom_fich_drag),form="formatted")
+ open(unit=14,file="RES/test/fich1.deb",form="formatted")
+ open(unit=15,file="RES/test/fich2.deb",form="formatted")
+ open(unit=16,file="RES/test/fich3.deb",form="formatted")
+ open(unit=17,file="RES/test/fich4.deb",form="formatted")
+ open(unit=18,file="RES/test/fich5.deb",form="formatted")
+ open(unit=19,file="RES/test/fich6.deb",form="formatted")
+ open(unit=20,file="RES/test/fich7.deb",form="formatted")
+
@@ -181,6 +189,7 @@ program NS
!!$ print*,"dx",dx
!!$ stop
+
@@ -200,26 +209,47 @@ program NS
call cpu_time(t1)
call vitesse_fft_uinf
call cpu_time(t2)
-
-
- !call imposer_debit(vxg,5.76d0,nx,ny,nz,dx,dy,dz)
-
+
+ !call imposer_debit(vxg,1.d0,nx,ny,nz,dx,dy,dz)
+ !vxg=vxg+1./(4.*1.*1.5) !debit=debit voulu +c*surface du domaine num, c etant la constante à ajouter au champ de vitesse.
+
!test la valeur du débit
!-----------------------
- t4=0.
- t5=0.
- t6=0.
- do k=1,nz
- do j=1,ny
- t4=t4+vxg(1,j,k)
- t5=t5+vxg(10,j,k)
- t6=t6+vxg(10,j,k)*(1.-chi(10,j,k))
- end do
- end do
- print*, "debit i=1 ",t4*dy*dz
- print*, "debit i=10 ",t5*dy*dz
- print*, "debit i=10 domaine penal ",t6*dy*dz
+!!$ t4=0.
+!!$ t5=0.
+!!$ t6=0.
+!!$ do k=1,nz
+!!$ do j=1,ny
+!!$ t4=t4+vxg(1,j,k)
+!!$ t5=t5+vxg(10,j,k)
+!!$ t6=t6+vxg(10,j,k)*(1.-chi(10,j,k))
+!!$ end do
+!!$ end do
+!!$ print*, "debit i=1 ",t4*dy*dz
+!!$ print*, "debit i=10 ",t5*dy*dz
+!!$ print*, "debit i=10 domaine penal ",t6*dy*dz
+
+
+
+!-------------------------test calcul vitesse avec omg exact------------------------------------------------------
+!!$ indi1=-xg/dx+1
+!!$ indi2=(0.7-xg)/dx+1
+!!$ indj=-yb/dy+1
+!!$
+!!$ do indk=1,nz
+!!$ if ( ((zd+(indk-1)*dz).ge.-1.).and.((zd+(indk-1)*dz).le.1.)) then
+!!$ vexa=1.5*(1.-(zd+(indk-1)*dz)**2)/(4.*0.5*abs(yh-yb))
+!!$ omexa=-3.*(zd+(indk-1)*dz)/(4.*0.5*abs(yh-yb))
+!!$ else
+!!$ vexa=0.
+!!$ omexa=0.
+!!$ end if
+!!$ write(14,'(8(e25.15,2x))') time,zd+(indk-1)*dz,vxg(indi1,indj,indk),omgy(indi1,indj,indk),vexa,omexa,vxg(indi2,indj,indk),omgy(indi2,indj,indk)
+!!$ end do
+!!$!stop
+!-----------------------------------------------------------------------------------------
+
!penalisation: nouveau calcul de v
!---------------------------------
@@ -237,6 +267,23 @@ program NS
end do
end do
+!-------------------------test ----------------------------------------------------------------
+!!$ indi1=-xg/dx+1
+!!$ indi2=(0.7-xg)/dx+1
+!!$ indj=-yb/dy+1
+!!$
+!!$ do indk=1,nz
+!!$ if ( ((zd+(indk-1)*dz).ge.-1.).and.((zd+(indk-1)*dz).le.1.)) then
+!!$ vexa=1.5*(1.-(zd+(indk-1)*dz)**2)/4.
+!!$ omexa=-3.*(zd+(indk-1)*dz)/4.
+!!$ else
+!!$ vexa=0.
+!!$ omexa=0.
+!!$ end if
+!!$ write(15,'(8(e25.15,2x))') time,zd+(indk-1)*dz,vxg(indi1,indj,indk),omgy(indi1,indj,indk),vexa,omexa,vxg(indi2,indj,indk),omgy(indi2,indj,indk)
+!!$ end do
+!!$!stop
+!-----------------------------------------------------------------------------------------
call cpu_time(t3)
@@ -249,8 +296,8 @@ program NS
if (mod(cpt_ite,5)==0) then
do indk=1,nz
if ( ((zd+(indk-1)*dz).ge.-1.).and.((zd+(indk-1)*dz).le.1.)) then
- vexa=1.5*(1.-(zd+(indk-1)*dz)**2)
- omexa=-3.*(zd+(indk-1)*dz)/(1.)**2
+ vexa=1.5*(1.-(zd+(indk-1)*dz)**2)/(4.*0.5*abs(yh-yb))
+ omexa=-3.*(zd+(indk-1)*dz)/(4.*0.5*abs(yh-yb))
else
vexa=0.
omexa=0.
@@ -291,7 +338,25 @@ program NS
!call penal_fft !calcul omg=rot(v)
call rotv_df4
+ !call rotv_df2
+!-------------------------test ----------------------------------------------------------------
+!!$ indi1=-xg/dx+1
+!!$ indi2=(0.7-xg)/dx+1
+!!$ indj=-yb/dy+1
+!!$
+!!$ do indk=1,nz
+!!$ if ( ((zd+(indk-1)*dz).ge.-1.).and.((zd+(indk-1)*dz).le.1.)) then
+!!$ vexa=1.5*(1.-(zd+(indk-1)*dz)**2)/4.
+!!$ omexa=-3.*(zd+(indk-1)*dz)/4.
+!!$ else
+!!$ vexa=0.
+!!$ omexa=0.
+!!$ end if
+!!$ write(16,'(8(e25.15,2x))') time,zd+(indk-1)*dz,vxg(indi1,indj,indk),omgy(indi1,indj,indk),vexa,omexa,vxg(indi2,indj,indk),omgy(indi2,indj,indk)
+!!$ end do
+!!$!stop
+!-----------------------------------------------------------------------------------------
call strech_diff_penal() ! penal commentee
@@ -307,6 +372,23 @@ program NS
call cpu_time(t5)
+!-------------------------test ----------------------------------------------------------------
+!!$ indi1=-xg/dx+1
+!!$ indi2=(0.7-xg)/dx+1
+!!$ indj=-yb/dy+1
+!!$
+!!$ do indk=1,nz
+!!$ if ( ((zd+(indk-1)*dz).ge.-1.).and.((zd+(indk-1)*dz).le.1.)) then
+!!$ vexa=1.5*(1.-(zd+(indk-1)*dz)**2)/4.
+!!$ omexa=-3.*(zd+(indk-1)*dz)/4.
+!!$ else
+!!$ vexa=0.
+!!$ omexa=0.
+!!$ end if
+!!$ write(17,'(8(e25.15,2x))') time,zd+(indk-1)*dz,vxg(indi1,indj,indk),omgy(indi1,indj,indk),vexa,omexa,vxg(indi2,indj,indk),omgy(indi2,indj,indk)
+!!$ end do
+!!$!stop
+!-----------------------------------------------------------------------------------------
!!$ !======================================
@@ -497,6 +579,13 @@ program NS
deallocate(stxg,styg,stzg,stx,sty,stz)
close(13)
+ close(14)
+ close(15)
+ close(16)
+ close(17)
+ close(18)
+ close(19)
+ close(20)
end program NS
diff --git a/NavierStokes3D-Penalization/old_mod.f90 b/NavierStokes3D-Penalization/old_mod.f90
index ee33c696b78069440f365521fdfb451bd8601515..aacee781f0c22bf5db63d0d2b4027c47c3a74eae 100644
--- a/NavierStokes3D-Penalization/old_mod.f90
+++ b/NavierStokes3D-Penalization/old_mod.f90
@@ -754,6 +754,46 @@ contains
end subroutine rotv_df4
+ subroutine rotv_df2
+ implicit none
+ integer :: l,m,n,i,j,k
+
+ integer :: ip1,ip2,im1,im2,jp1,jp2,jm1,jm2,kp1,kp2,km1,km2
+ real(kind=8) :: facx,facy,facz,dxsca,dysca,dzsca
+
+ !omgx = dy vz - dz vy
+ !omgy = dz vx - dx vz
+ !omgz = dx vy - dy vx
+
+ facx=1.D0/(2.D0*dx)
+ facy=1.D0/(2.D0*dy)
+ facz=1.D0/(2.D0*dz)
+
+ do k=1,nz
+ do j=1,ny
+ do i=1,nx
+
+ ip1=mod(i+1+nx-1,nx)+1
+ im1=mod(i-1+nx-1,nx)+1
+
+ jp1=mod(j+1+ny-1,ny)+1
+ jm1=mod(j-1+ny-1,ny)+1
+
+ kp1=mod(k+1+nz-1,nz)+1
+ km1=mod(k-1+nz-1,nz)+1
+
+ omgx(i,j,k)=facy*(-vzg(i,jm1,k)+vzg(i,jp1,k)) &
+ -facz*(-vyg(i,j,km1)+vyg(i,j,kp1))
+ omgy(i,j,k)=facz*(-vxg(i,j,km1)+vxg(i,j,kp1)) &
+ -facx*(-vzg(im1,j,k)+vzg(ip1,j,k))
+ omgz(i,j,k)=facx*(-vyg(im1,j,k)+vyg(ip1,j,k)) &
+ -facy*(-vxg(i,jm1,k)+vxg(i,jp1,k))
+ end do
+ end do
+ end do
+
+ end subroutine rotv_df2
+
diff --git a/NavierStokes3D-Penalization/parameter b/NavierStokes3D-Penalization/parameter
index 587d7fefd497f7f51b045f3a9e66d08c8eee7e08..04b6519931d089a35c34ac54ba9661caac762b77 100644
--- a/NavierStokes3D-Penalization/parameter
+++ b/NavierStokes3D-Penalization/parameter
@@ -1,12 +1,12 @@
--1.2 1.2 -1.2 1.2 -1.2 1.2 ! xg xd yb yh zd zf
+-1 1 -1 1 -1.2 1.2 ! xg xd yb yh zd zf
33 33 33 33 65 129 257!nx ny nz
-50. ! TFin
-1.D-6 !cutoff
+10. ! TFin
+1.D-8 !cutoff
4 ! type bloc
2 !longeur bloc
"test.vtk" ! sortie fin vtk
5.5D5 !lambda
-0.1 7.496251874D-03 (Re=133.4) 4.74833808D-03 (Re=210.6) !nu=1/Re
+0.01 7.496251874D-03 (Re=133.4) 4.74833808D-03 (Re=210.6) !nu=1/Re
"RES/test/omg_" "RES/test/vg_" 20000 !vtk
-"RES/test/coupe1.res" !drag/lift ou coupe vitesse
+"RES/test/pois_33_df4.res" !drag/lift ou coupe vitesse
0 200000 "RES/test/sauv.vtk" !num_suite (0non 1 oui) num_ite name_relance
\ No newline at end of file
diff --git a/NavierStokes3D-Penalization/utile_mod.f90 b/NavierStokes3D-Penalization/utile_mod.f90
index bad09481c9a36ad7b12ce37718c76b55dfe31fd3..89b057611d505564c4f018f925e676d15ba63b6e 100644
--- a/NavierStokes3D-Penalization/utile_mod.f90
+++ b/NavierStokes3D-Penalization/utile_mod.f90
@@ -359,6 +359,7 @@ contains
+
end module utile_mod