diff --git a/Remesh/split_2d/Makefile b/Remesh/split_2d/Makefile
index ed4a869ecc136905b1844188f7e95c99b827b410..b2ad1583c132376fed0cfea4fd236b723449aa5b 100644
--- a/Remesh/split_2d/Makefile
+++ b/Remesh/split_2d/Makefile
@@ -1,8 +1,8 @@
#
# Makefile
#
-#Compilateur = ifort
-Compilateur = ifort -O0 -pg
+#Compilateur = ifort
+Compilateur = ifort -pg
#Compilateur = ifort -O0 -debug extended
Programme = split_2d
FILES = donnees_mod.f90 donnees_limit_mod.f90 tab_mod.f90 init_mod.f90 remaillage_mod.f90 interpolation_mod.f90 advection_mod.f90 utile_mod.f90 resultats_mod.f90 main.f90
diff --git a/Remesh/split_2d/advection_mod.f90 b/Remesh/split_2d/advection_mod.f90
index e8e5d425c5debcd4c02a949aa46cb297861f6f72..de520926217234862dba09e64c6fdacf70a913d3 100644
--- a/Remesh/split_2d/advection_mod.f90
+++ b/Remesh/split_2d/advection_mod.f90
@@ -20,8 +20,8 @@ contains
y=yb+(j-1)*dy_gro
- !rotation avec chgmt signe
- !-------------------------
+!!$ !rotation avec chgmt signe
+!!$ !-------------------------
!!$ r2=x**2+y**2
!!$ !tmp=1
!!$ tmp=cos(3.*0.5*pi*sqrt(r2))
@@ -51,7 +51,7 @@ contains
!test instationaire
!--------------------
-!!$
+
!!$ t1=1.
!!$ t4=0.3
!!$ t5=1.
@@ -66,9 +66,9 @@ contains
vitx(i,j)=x/r
vity(i,j)=y/r
- if (r<0.45) then
- vitx(i,j)=x/0.45
- vity(i,j)=y/0.45
+ if (r<0.4) then
+ vitx(i,j)=x/0.4
+ vity(i,j)=y/0.4
end if
@@ -190,8 +190,11 @@ end subroutine evalv
yp1(ib)=yp(ib)+0.5*dt*vy(ib)
end do
- call interpo_l3_2d(vxg1,xp1,yp1,vx)
- call interpo_l3_2d(vyg1,xp1,yp1,vy)
+ call interpo_l2_2d(vxg,xp1,yp1,vx)
+ call interpo_l2_2d(vyg,xp1,yp1,vy)
+
+ !call interpo_l3_2d(vxg1,xp1,yp1,vx)
+ !call interpo_l3_2d(vyg1,xp1,yp1,vy)
do ib=1,npart
xp(ib)=xp(ib)+dt*vx(ib)
@@ -220,16 +223,16 @@ end subroutine evalv
yp1(ib)=yp(ib)+0.5*dt*vy(ib)
end do
- call interpo_l4_2d(vxg1,xp1,yp1,vx1)
- call interpo_l4_2d(vyg1,xp1,yp1,vy1)
+ call interpo_l3_2d(vxg1,xp1,yp1,vx1)
+ call interpo_l3_2d(vyg1,xp1,yp1,vy1)
do ib=1,npart
xp2(ib)=xp(ib)-dt*vx(ib)+2.*dt*vx1(ib)
yp2(ib)=yp(ib)-dt*vy(ib)+2.*dt*vy1(ib)
end do
- call interpo_l4_2d(vxg2,xp2,yp2,vx2)
- call interpo_l4_2d(vyg2,xp2,yp2,vy2)
+ call interpo_l3_2d(vxg2,xp2,yp2,vx2)
+ call interpo_l3_2d(vyg2,xp2,yp2,vy2)
do ib=1,npart
xp(ib)=xp(ib)+dt*(vx(ib)/6.+2.*vx1(ib)/3.+vx2(ib)/6.)
@@ -256,24 +259,24 @@ end subroutine evalv
yp1(ib)=yp(ib)+0.5*dt*vy(ib)
end do
- call interpo_l5_2d(vxg1,xp1,yp1,vx1)
- call interpo_l5_2d(vyg1,xp1,yp1,vy1)
+ call interpo_l4_2d(vxg1,xp1,yp1,vx1)
+ call interpo_l4_2d(vyg1,xp1,yp1,vy1)
do ib=1,npart
xp2(ib)=xp(ib)+0.5*dt*vx1(ib)
yp2(ib)=yp(ib)+0.5*dt*vy1(ib)
end do
- call interpo_l5_2d(vxg1,xp2,yp2,vx2)
- call interpo_l5_2d(vyg1,xp2,yp2,vy2)
+ call interpo_l4_2d(vxg1,xp2,yp2,vx2)
+ call interpo_l4_2d(vyg1,xp2,yp2,vy2)
do ib=1,npart
xp3(ib)=xp(ib)+dt*vx2(ib)
yp3(ib)=yp(ib)+dt*vy2(ib)
end do
- call interpo_l5_2d(vxg2,xp3,yp3,vx3)
- call interpo_l5_2d(vyg2,xp3,yp3,vy3)
+ call interpo_l4_2d(vxg2,xp3,yp3,vx3)
+ call interpo_l4_2d(vyg2,xp3,yp3,vy3)
do ib=1,npart
xp(ib)=xp(ib)+dt*(vx(ib)+2.*vx1(ib)+2.*vx2(ib)+vx3(ib))/6.
@@ -330,7 +333,8 @@ end subroutine evalv
xp1(ib)=xp(ib)+0.5*dt*vx(ib)
yp1(ib)=yp(ib)+0.5*dt*vy(ib)
end do
- call interpo_l3_2d(vxg1,xp1,yp1,vx)
+ !call interpo_l2_2d(vxg1,xp1,yp1,vx)
+ call interpo_l2_2d(vxg,xp1,yp1,vx)
deallocate(xp1,yp1)
end subroutine update_vx_2
@@ -344,7 +348,8 @@ end subroutine evalv
xp1(ib)=xp(ib)-0.5*dt*vx(ib)
yp1(ib)=yp(ib)+0.5*dt*vy(ib)
end do
- call interpo_l3_2d(vyg1,xp1,yp1,vy)
+ !call interpo_l2_2d(vyg1,xp1,yp1,vy)
+ call interpo_l2_2d(vyg,xp1,yp1,vy)
deallocate(xp1,yp1)
end subroutine update_vy_2
@@ -364,15 +369,19 @@ end subroutine evalv
yp1(ib)=yp(ib)+2.*dt*vy(ib)/3.
end do
- call interpo_l4_2d(vxg2,xp1,yp1,vx1)
- call interpo_l3_2d(vyg2,xp1,yp1,vy1)
+ !call interpo_l3_2d(vxg2,xp1,yp1,vx1)
+ !call interpo_l3_2d(vyg2,xp1,yp1,vy1)
+ call interpo_l3_2d(vxg,xp1,yp1,vx1)
+ call interpo_l3_2d(vyg,xp1,yp1,vy1)
do ib=1,npart
xp2(ib)=xp(ib)+dt*(-vx(ib)/4.+vx1(ib)/4.)
yp2(ib)=yp(ib)+dt*(-vy(ib)+vy1(ib))
end do
- call interpo_l4_2d(vxg,xp2,yp,vx2)
+ !call interpo_l3_2d(vxg,xp2,yp,vx2)
+ !call interpo_l3_2d(vxg,xp,yp2,vy2)
+ call interpo_l3_2d(vxg,xp2,yp,vx2)
call interpo_l3_2d(vxg,xp,yp2,vy2)
do ib=1,npart
@@ -399,16 +408,20 @@ end subroutine evalv
xp0(ib)=xp(ib)-dt*vx(ib)/6.
end do
- call interpo_l4_2d(vxg1,xp0,yp1,vx1)
- call interpo_l3_2d(vyg1,xp1,yp1,vy1)
+ !call interpo_l3_2d(vxg1,xp0,yp1,vx1)
+ !call interpo_l3_2d(vyg1,xp1,yp1,vy1)
+ call interpo_l3_2d(vxg,xp0,yp1,vx1)
+ call interpo_l3_2d(vyg,xp1,yp1,vy1)
do ib=1,npart
xp2(ib)=xp(ib)+dt*(vx(ib)-vx1(ib))
yp2(ib)=yp(ib)+dt*(-vy(ib)+2.*vy1(ib))
end do
- call interpo_l4_2d(vyg,xp2,yp,vy2)
- call interpo_l3_2d(vyg3,xp,yp2,vy3)
+ !call interpo_l3_2d(vyg,xp2,yp,vy2)
+ !call interpo_l3_2d(vyg3,xp,yp2,vy3)
+ call interpo_l3_2d(vyg,xp2,yp,vy2)
+ call interpo_l3_2d(vyg,xp,yp2,vy3)
do ib=1,npart
vy(ib)=-vy(ib)+3.*vy1(ib)/4.+vy2(ib)+vy3(ib)/4.
diff --git a/Remesh/split_2d/init_mod.f90 b/Remesh/split_2d/init_mod.f90
index d12c5fe1bd5ef95b76ffddc9958a052856b775f1..a3e1ce4e8d8b82944d1a261c6d16ded7cfeee87d 100644
--- a/Remesh/split_2d/init_mod.f90
+++ b/Remesh/split_2d/init_mod.f90
@@ -107,19 +107,36 @@ contains
!---anneau-----------
!-----------------------
+!!$ do i=1,nx
+!!$ do j=1,ny
+!!$ x=xtab(i,j)
+!!$ y=ytab(i,j)
+!!$ r=sqrt(x**2+y**2)
+!!$ qg(i,j)=0.
+!!$ !if ( (r>0.5).and.(r<1.4) ) qg(i,j)=1.
+!!$ !if ( (r>0.).and.(r<0.8) ) qg(i,j)=1.
+!!$ !if ( (r>=1.).and.(r<1.5) ) qg(i,j)=1.
+!!$ if ( (r>0.5).and.(r<1.) ) qg(i,j)=1.
+!!$ end do
+!!$ end do
+
+ !---anneau- regulier----------
+ !-----------------------
do i=1,nx
do j=1,ny
x=xtab(i,j)
y=ytab(i,j)
r=sqrt(x**2+y**2)
qg(i,j)=0.
- !if ( (r>0.5).and.(r<1.4) ) qg(i,j)=1.
- !if ( (r>0.).and.(r<0.8) ) qg(i,j)=1.
- !if ( (r>=1.).and.(r<1.5) ) qg(i,j)=1.
- if ( (r>0.5).and.(r<1.) ) qg(i,j)=1.
+ if ( (r>0.4).and.(r<0.6) ) qg(i,j)=-512.+5400.*r-22500.*r**2+46250*r**3-46875*r**4+18750*r**5
+ if ( (r>=0.6).and.(r<=0.8) ) qg(i,j)=1.
+ if ( (r>0.8).and.(r<1.) ) qg(i,j)=10625.-60000*r+135000*r**2-151250*r**3+84375*r**4-18750*r**5
+
+ !if ( (r>0.4).and.(r<1.) ) qg(i,j)=1.
end do
end do
+
end subroutine init_grille
diff --git a/Remesh/split_2d/main.f90 b/Remesh/split_2d/main.f90
index fee6ec6cb657921e1341f2f303fd0e71348df10e..84d6afbd7bd6ccbc9e6f22acaf7e31c7b7c592ef 100644
--- a/Remesh/split_2d/main.f90
+++ b/Remesh/split_2d/main.f90
@@ -13,10 +13,10 @@ program split
use utile_mod ! fonctions en tout genre
real(kind=8) :: x,y,dt_deb,dt_boucle,dtsauv,maxv
- character(len=50) :: name_err
+ character(len=50) :: name_err,name_coupe
integer :: i,j,cpt_ite,ib,jb,nx_boucle
- real(kind=8) :: deb,masse
+ real(kind=8) :: deb,masse,t1,t2,t3
!lecture données
@@ -29,14 +29,13 @@ program split
read(10,*) long_bloc
read(10,*) int_limit
read(10,*) nom_fich_vtk
- read(10,*) name_err
+ read(10,*) name_err,name_coupe
close(10)
-
- !erreur:do nx_boucle=10,2000,30
- erreur:do nx_boucle=300,300
+ erreur:do nx_boucle=201,201
!erreur:do nx_boucle=10,5000,30
+ !erreur:do nx_boucle=10,5000,100
!nx=(long_bloc+1)*k+1 => bloc bord de bonne taille
open(unit=10,file="parameter",form="formatted")
@@ -58,12 +57,7 @@ program split
maxv=1.
!maxv=2.
- cfl=10
- !cfl=12
- !cfl=24
- !cfl=3.
- !cfl=20.
- !cfl=10.
+ cfl=3.
dt=cfl*dx/maxv
!dt=0.0001
@@ -141,8 +135,10 @@ program split
allocate (xp(1:npart),qp(1:npart),vx(1:npart),yp(1:npart),vy(1:npart))
allocate (blocg(0:npart),blocd(0:npart),Nbloc(0:npart))
- open(unit=29,file=name_err,form="formatted")
+ open(unit=29,file=name_err,form="formatted")
+ open(unit=30,file=name_coupe,form="formatted")
+ call cpu_time(t1)
!
!initialisation
@@ -151,7 +147,7 @@ program split
call init_grille()
masse=dx*dy*sum(qg)
-
+ call def_v_advec(vxg,vyg,time)
!
!boucle temps
@@ -164,7 +160,8 @@ program split
if (time<=2.*dt) dt_deb=dt
- call def_v_advec(vxg,vyg,time)
+
+ !call def_v_advec(vxg,vyg,time) !inutile de le garder dans la boucle quand depend pas du temps.
call crea_part_x()
@@ -177,23 +174,23 @@ program split
!!$ !champ de vitesse
!!$ !-----------------
-!!$ allocate(vxg1(1:nx_gro,1:ny_gro),vyg1(1:nx_gro,1:ny_gro))
+!!$ !allocate(vxg1(1:nx_gro,1:ny_gro),vyg1(1:nx_gro,1:ny_gro))
!!$ !rk2
!!$ !call def_v_advec(vxg1,vyg1,time+0.5*dt)
-!!$ vxg1=vxg
-!!$ vyg1=vyg
+!!$ !vxg1=vxg
+!!$ !vyg1=vyg
!!$
!!$ !rk3
-!!$ allocate(vxg2(1:nx_gro,1:ny_gro),vyg2(1:nx_gro,1:ny_gro))
+!!$ !allocate(vxg2(1:nx_gro,1:ny_gro),vyg2(1:nx_gro,1:ny_gro))
!!$ !call def_v_advec(vxg1,vyg1,time+0.5*dt)
!!$ !call def_v_advec(vxg2,vyg2,time+dt)
-!!$ vxg2=vxg
-!!$ vyg2=vyg
+!!$ !vxg2=vxg
+!!$ !vyg2=vyg
!!$
!!$ !rk4
-!!$ allocate(vxg3(1:nx_gro,1:ny_gro),vyg3(1:nx_gro,1:ny_gro))
-!!$ vxg3=vxg
-!!$ vyg3=vyg
+!!$ !allocate(vxg3(1:nx_gro,1:ny_gro),vyg3(1:nx_gro,1:ny_gro))
+!!$ !vxg3=vxg
+!!$ !vyg3=vyg
!!$
!!$
!!$ !t source:rk4
@@ -203,24 +200,24 @@ program split
!!$ !advection
!!$ !---------
!!$ !call ad_tenso_euler
-!!$ !call ad_tenso_rk2
+!!$ call ad_tenso_rk2
!!$ !call ad_tenso_rk3
-!!$ call ad_tenso_rk4
+!!$ !call ad_tenso_rk4
!!$
!!$
!!$
!!$ !remaillage
!!$ !----------
!!$ !call remaill_l2(qp,xp,yp,qg)
-!!$ !call remaill_l3(qp,xp,yp,qg)
+!!$ call remaill_l3(qp,xp,yp,qg)
!!$ !call remaill_l4(qp,xp,yp,qg)
!!$ !call remaill_l6(qp,xp,yp,qg)
!!$ !call remaill_mp4(qp,xp,yp,qg)
-!!$ call remaill_mppp6(qp,xp,yp,qg)
+!!$ !call remaill_mppp6(qp,xp,yp,qg)
!!$
-!!$ deallocate(vxg1,vyg1)
-!!$ deallocate(vxg2,vyg2)
-!!$ deallocate(vxg3,vyg3)
+!!$ !deallocate(vxg1,vyg1)
+!!$ !deallocate(vxg2,vyg2)
+!!$ !deallocate(vxg3,vyg3)
!=====================
!splitting en espace
@@ -228,11 +225,11 @@ program split
!champ de vitesse
!-----------------
- allocate(vxg1(1:nx_gro,1:ny_gro),vyg1(1:nx_gro,1:ny_gro))
+ !allocate(vxg1(1:nx_gro,1:ny_gro),vyg1(1:nx_gro,1:ny_gro))
!rk2
!call def_v_advec(vxg1,vyg1,time+0.5*dt)
- vxg1=vxg
- vyg1=vyg ! attention sinon modifier dans les update vitesses d'advection
+ !vxg1=vxg
+ !vyg1=vyg ! attention sinon modifier dans les update vitesses d'advection
!rk3
!allocate(vxg2(1:nx_gro,1:ny_gro),vxg3(1:nx_gro,1:ny_gro),vyg2(1:nx_gro,1:ny_gro),vyg3(1:nx_gro,1:ny_gro))
@@ -250,7 +247,7 @@ program split
!update v
!--------
!call update_vx_2
- !call update_vx_3
+ call update_vx_3
!test cfl
@@ -260,23 +257,23 @@ program split
!print*,"max vx",maxval(vx)
- ainf=0.
- do j=1,ny
- do i=1,nx-1
- if ( (numpg(i+1,j)/=0).and.(numpg(i,j)/=0) ) ainf=max(ainf,abs(vx(numpg(i+1,j))-vx(numpg(i,j))))
- end do
- end do
- if (dt>(0.5*dx/(ainf))) then
- print*,""
- print*,"cfl l2 viole" ,dt,0.5*dx/(ainf),0.5*dx/(2.*ainf)
- print*,""
- end if
- if (dt>(0.5*dx/(2.*ainf))) then
- print*,""
- print*,"cfl l4 viole" ,dt,0.5*dx/(ainf),0.5*dx/(2.*ainf)
- print*,"dx,dt",dx,dt
- print*,""
- end if
+! ainf=0.
+! do j=1,ny
+! do i=1,nx-1
+! if ( (numpg(i+1,j)/=0).and.(numpg(i,j)/=0) ) ainf=max(ainf,abs(vx(numpg(i+1,j))-vx(numpg(i,j))))
+! end do
+! end do
+! if (dt>(0.5*dx/(ainf))) then
+! print*,""
+! print*,"cfl l2 viole" ,dt,0.5*dx/(ainf),0.5*dx/(2.*ainf)
+! print*,""
+! end if
+! if (dt>(0.5*dx/(2.*ainf))) then
+! print*,""
+! print*,"cfl l4 viole" ,dt,0.5*dx/(ainf),0.5*dx/(2.*ainf)
+! print*,"dx,dt",dx,dt
+! print*,""
+! end if
! if (dt>(dx/(2.*ainf/(sqrt(3.)-1.)))) then
! print*,""
! print*,"TSC gauche pas TVD (l2)" ,dt,dx/(2.*ainf/(sqrt(3.)-1.))
@@ -315,6 +312,7 @@ program split
if (type_b==20) call remaill_l2_x(qp,xp,yp,qg)
if (type_b==50)call remaill_l5_x(qp,xp,yp,qg)
if (type_b==60)call remaill_l6_x(qp,xp,yp,qg)
+ if (type_b==70)call remaill_mp4_x(qp,xp,yp,qg)
!if (type_b==2) call remaill_l2_bloc_x(qp,xp,yp,qg)
!if (type_b==4) call remaill_l4_bloc_x(qp,xp,yp,qg)
if (type_b==2) then
@@ -326,8 +324,8 @@ program split
if (int_limit==1)call remaill_l4_bloc_v2_limit_x(qp,xp,yp,qg)
end if
- if (type_b==200)call remaill_l2_bloc_v2_limit_x(qp,xp,yp,qg)
- if (type_b==400)call remaill_l4_bloc_v2_limit_x(qp,xp,yp,qg)
+ if (type_b==200)call remaill_l2_limit_x(qp,xp,yp,qg)
+ if (type_b==400)call remaill_l4_limit_x(qp,xp,yp,qg)
!y
!--
@@ -336,7 +334,7 @@ program split
!update v
!--------
!call update_vy_2
- !call update_vy_3
+ call update_vy_3
!test cfl
@@ -379,34 +377,34 @@ program split
!-------test----------------------------------------
-!!$ if (cpt_ite <=1) then
-!!$
-!!$ !trace la position des particules
-!!$ open(unit=28,file="RES/pos_part.deb",form="formatted")
-!!$ open(unit=29,file="RES/pos_part_2.deb",form="formatted")
-!!$ do ib=1,nx
-!!$ x=xg+(ib-1)*dx
-!!$ if ( (x>0.8).and.(x<0.81) ) then
-!!$ do jb=1,(ny+10)
-!!$ y=yb+(jb-1)*dy
-!!$ write(28,'(2(f35.20,2x))') y,0.5
-!!$ end do
-!!$ end if
-!!$ end do
-!!$
-!!$ do ib=1,npart
-!!$ if ((xp(ib)>0.8).and.(xp(ib)<0.81)) then
-!!$ write(29,'(4(f35.20,2x))') yp(ib),0.5,blocg(ib)*0.1 ,blocd(ib)*0.1
-!!$ end if
-!!$ end do
-!!$
-!!$ close (28)
-!!$ close (29)
-!!$
-!!$ end if
+! if (cpt_ite <=1) then
+!
+! !trace la position des particules
+! open(unit=28,file="RES/pos_part.deb",form="formatted")
+! open(unit=29,file="RES/pos_part_2.deb",form="formatted")
+! do ib=1,nx
+! x=xg+(ib-1)*dx
+! if ( (x>0.8).and.(x<0.81) ) then
+! do jb=1,(ny+10)
+! y=yb+(jb-1)*dy
+! write(28,'(2(f35.20,2x))') y,0.5
+! end do
+! end if
+! end do
+!
+! do ib=1,npart
+! if ((xp(ib)>0.8).and.(xp(ib)<0.81)) then
+! write(29,'(4(f35.20,2x))') yp(ib),0.5,blocg(ib)*0.1 ,blocd(ib)*0.1
+! end if
+! end do
+!
+! close (28)
+! close (29)
+!
+! end if
!-------test----------------------------------------
-
+!!$
!remaillage
!----------
if (type_b==30)call remaill_l3_y(qp,xp,yp,qg)
@@ -414,6 +412,7 @@ program split
if (type_b==20)call remaill_l2_y(qp,xp,yp,qg)
if (type_b==50)call remaill_l5_y(qp,xp,yp,qg)
if (type_b==60)call remaill_l6_y(qp,xp,yp,qg)
+ if (type_b==70)call remaill_mp4_y(qp,xp,yp,qg)
!if (type_b==2)call remaill_l2_bloc_y(qp,xp,yp,qg)
!if (type_b==4)call remaill_l4_bloc_y(qp,xp,yp,qg)
if (type_b==2) then
@@ -425,15 +424,15 @@ program split
if (int_limit==1)call remaill_l4_bloc_v2_limit_y(qp,xp,yp,qg)
end if
- if (type_b==200)call remaill_l2_bloc_v2_limit_y(qp,xp,yp,qg)
- if (type_b==400)call remaill_l4_bloc_v2_limit_y(qp,xp,yp,qg)
+ if (type_b==200)call remaill_l2_limit_y(qp,xp,yp,qg)
+ if (type_b==400)call remaill_l4_limit_y(qp,xp,yp,qg)
!t source
!---------
!call source_split_2
- deallocate(vxg1,vyg1)
+ !deallocate(vxg1,vyg1)
!deallocate(vxg2,vxg3,vyg2,vyg3)
@@ -1147,24 +1146,24 @@ program split
!erreur pour sol init gaussienne
!--------------------------------
- do i=1,nx
- do j=1,ny
- x=xtab(i,j)
- y=ytab(i,j)
- exa(i,j)=0.
-
- !exa(i,j)=(1.-((x-0.5)**2+(y-0.5)**2))**6
- !if (((x-0.5)**2+(y-0.5)**2)>1.) exa(i,j)=0.
-
- exa(i,j)=(1.-((x)**2+(y)**2))**6
- if (((x)**2+(y)**2)>1.) exa(i,j)=0.
-
- !exa(i,j)=1.
-
-
- end do
- end do
-
+!!$ do i=1,nx
+!!$ do j=1,ny
+!!$ x=xtab(i,j)
+!!$ y=ytab(i,j)
+!!$ exa(i,j)=0.
+!!$
+!!$ !exa(i,j)=(1.-((x-0.5)**2+(y-0.5)**2))**6
+!!$ !if (((x-0.5)**2+(y-0.5)**2)>1.) exa(i,j)=0.
+!!$
+!!$ exa(i,j)=(1.-((x)**2+(y)**2))**6
+!!$ if (((x)**2+(y)**2)>1.) exa(i,j)=0.
+!!$
+!!$ !exa(i,j)=1.
+!!$
+!!$
+!!$ end do
+!!$ end do
+!!$
!!$ !-----------Analytique en cos--------------------------
!!$ !------------------------------------------------------
!!$ t1=1.
@@ -1192,17 +1191,22 @@ program split
!!$ end do
!!$ end do
-
-
- write(29,'(4(e18.10,2x))') dt_deb,dx,dx*dy*sum(abs(exa-qg)),sqrt(dx*dy*sum((exa-qg)**2))
-
+ !write(29,'(4(e18.10,2x))') dt_deb,dx,dx*dy*sum(abs(exa-qg)),sqrt(dx*dy*sum((exa-qg)**2))
+
+ !COUPE en y=1.5
+ !--------------
+ j=(1.5-yb)/dy+1
+ do i=1,nx
+ write(30,'(2(e18.10,2x))') xg+(i-1)*dx,qg(i,j)
+ end do
+
end do erreur
-
+ call cpu_time(t2)
!resultat final
- !--------------
+ !--------------
print*,"tps final",time
print*,"nbre d'ite",cpt_ite
print*,"nart final",npart
@@ -1210,6 +1214,7 @@ program split
print*,"max vx",maxval(abs(vxg)),maxloc(abs(vxg))
print*,"max vy",maxval(abs(vyg)),maxloc(abs(vxg))
print*,"perte masse",masse-dx*dy*sum(qg)
+ print*,"tps dexecution",t2-t1
@@ -1226,6 +1231,8 @@ program split
close(53)
close(66)
close(67)
+ close(29)
+ close(30)
deallocate (xp,qp,vx,vy,yp)
deallocate (qg,vxg,vyg,exa)
diff --git a/Remesh/split_2d/parameter b/Remesh/split_2d/parameter
index d1f8979835ee73278ceb83ad3489f6d38628e027..1998801abbe062f49397728b71bd4d0046f2ba0e 100644
--- a/Remesh/split_2d/parameter
+++ b/Remesh/split_2d/parameter
@@ -1,8 +1,8 @@
--1.5 1.5 -1.5 1.5 -1 1 0 1 0 1 -2 2 -2 2 ! xg xd yb yh
-0.4 2. 5. 0.05 ! TFin
-0. 0.00001 !cutoff
-2 !type bloc
+-2 2 -2 2 ! xg xd yb yh
+0.8 2. 5. 0.05 ! TFin
+0.00001 !cutoff
+70 !type bloc
1 !longeur des blocs
-1 ! limitation 0=non 1=oui
-"test_radial.vtk" ! sortie fin vtk
-"RES/test.er" ! fichier erreur
\ No newline at end of file
+0 ! limitation 0=non 1=oui
+"anregu_o3_cfl3_200_t08_cutm5_mp4.vtk" ! sortie fin vtk
+"RES/test.er" "RES/anregu_o3_cfl3_200_t08_cutm5_mp4.co" ! fichier erreur /coupe
\ No newline at end of file
diff --git a/Remesh/split_2d/remaillage_mod.f90 b/Remesh/split_2d/remaillage_mod.f90
index 6a214f0ea782ec880aec6a2525b31f4160cbfacf..bbbd60860539dbf4c3ab034a2f5bcaf876c71c02 100644
--- a/Remesh/split_2d/remaillage_mod.f90
+++ b/Remesh/split_2d/remaillage_mod.f90
@@ -1010,6 +1010,144 @@ contains
end subroutine remaill_l3_y
+
+
+
+
+
+ subroutine remaill_mp4_x(donne,posx,posy,remaille)
+ implicit none
+ real(kind=8),dimension(:),intent(in) :: donne
+ real(kind=8),dimension(:),intent(in) :: posx,posy
+ real(kind=8),dimension(1:nx,1:ny),intent(out) :: remaille
+ integer :: i,c,d,ib,j,ivar,k
+ integer,dimension(0:3) :: ip
+ real(kind=8),dimension(0:3) :: poidx
+ real(kind=8) :: xx1,tmp_pos
+ real(kind=8) :: x2,x3,x4,x5
+
+
+ remaille=0.
+
+
+ !interpole en 2d la vorticite sur la grille en fonction de la position des particules et de leur vorticité
+ !---------------------------------------------------------------------------------------------------------
+
+ do i=1,npart
+
+ j=nint((posy(i)-yb)/dy)+1
+
+ !numero des points sur le maillage
+ !--------------------------------
+ ip(1) = floor((posx(i)-xg)/dx)
+ ip(0) = ip(1) - 1
+ ip(2) = ip(1) + 1
+ ip(3) = ip(1) + 2
+
+
+ !distance de la particule à remailler au second point (des quatres utilisés pour le remaillage)
+ !----------------------------------------------------
+ xx1 = (posx(i) - real(ip(1),kind=8)*dx-xg)/dx !relatif
+
+ !conditions au bord
+ !------------------
+ !periodique:
+ do c=0,3
+ ip(c)=mod(ip(c)+nx,nx)
+ end do
+
+
+ !calcul des poids
+ !----------------
+ x2=xx1**2
+ x3=xx1**3
+ x4=xx1**4
+ x5=xx1**5
+ poidx(0)=-0.5*xx1+x2-0.5*x3
+ poidx(1)=1-5./2.*x2+3.*0.5*x3
+ poidx(2)=0.5*xx1+2.*x2-3.*0.5*x3
+ poidx(3)=-0.5*x2+0.5*x3
+
+
+ !remaillage à l' interrieur domaine
+ !---------------------------------
+ do c=0,3
+ remaille(ip(c)+1,j)=remaille(ip(c)+1,j)+donne(i)*poidx(c)
+ end do
+
+
+ end do
+
+ end subroutine remaill_mp4_x
+
+
+ subroutine remaill_mp4_y(donne,posx,posy,remaille)
+ implicit none
+ real(kind=8),dimension(:),intent(in) :: donne
+ real(kind=8),dimension(:),intent(in) :: posx,posy
+ real(kind=8),dimension(1:nx,1:ny),intent(out) :: remaille
+ integer :: i,c,d,ib,j,ivar,k
+ integer,dimension(0:3) :: ip
+ real(kind=8),dimension(0:3) :: poidx
+ real(kind=8) :: xx1,tmp_pos
+ real(kind=8) :: x2,x3,x4,x5
+
+ remaille=0.
+
+
+ !interpole en 2d la vorticite sur la grille en fonction de la position des particules et de leur vorticité
+ !---------------------------------------------------------------------------------------------------------
+
+ do i=1,npart
+
+ ivar=nint((posx(i)-xg)/dx)+1
+
+ !numero des points sur le maillage
+ !--------------------------------
+ ip(1) = floor((posy(i)-yb)/dy)
+ ip(0) = ip(1) - 1
+ ip(2) = ip(1) + 1
+ ip(3) = ip(1) + 2
+
+
+ !distance de la particule à remailler au second point (des quatres utilisés pour le remaillage)
+ !----------------------------------------------------
+ xx1 = (posy(i) - real(ip(1),kind=8)*dy-yb)/dy !relatif
+
+ !conditions au bord
+ !------------------
+ !periodique:
+ do c=0,3
+ ip(c)=mod(ip(c)+ny,ny)
+ end do
+
+ !calcul des poids
+ !----------------
+ x2=xx1**2
+ x3=xx1**3
+ x4=xx1**4
+ x5=xx1**5
+ poidx(0)=-0.5*xx1+x2-0.5*x3
+ poidx(1)=1-5./2.*x2+3.*0.5*x3
+ poidx(2)=0.5*xx1+2.*x2-3.*0.5*x3
+ poidx(3)=-0.5*x2+0.5*x3
+
+ !remaillage à l' interrieur domaine
+ !---------------------------------
+ do c=0,3
+ remaille(ivar,ip(c)+1)=remaille(ivar,ip(c)+1)+donne(i)*poidx(c)
+ end do
+
+
+ end do
+
+ end subroutine remaill_mp4_y
+
+
+
+
+
+
subroutine remaill_l4_x(donne,posx,posy,remaille)
implicit none
real(kind=8),dimension(:),intent(in) :: donne
@@ -1782,224 +1920,84 @@ contains
!!$ end subroutine remaill_m6kernel_y
- subroutine remaill_mp4_x(donne,posx,posy,remaille)
- implicit none
- real(kind=8),dimension(:),intent(in) :: donne
- real(kind=8),dimension(:),intent(in) :: posx,posy
- real(kind=8),dimension(1:nx,1:ny),intent(out) :: remaille
- integer :: i,c,d,ib,j,ivar,k
- integer,dimension(0:5) :: ip
- real(kind=8),dimension(0:5) :: poidx
- real(kind=8) :: xx1,tmp_pos,x2,x3,x4,x5
- remaille=0.
- !interpole en 2d la vorticite sur la grille en fonction de la position des particules et de leur vorticité
- !---------------------------------------------------------------------------------------------------------
- do i=1,npart
- j=nint((posy(i)-yb)/dy)+1
- !numero des points sur le maillage
- !--------------------------------
- ip(2) = int((posx(i)-xg)/dx)
- ip(0) = ip(2) - 2
- ip(1) = ip(2) - 1
- ip(3) = ip(2) + 1
- ip(4) = ip(2) + 2
- ip(5) = ip(2) + 3
- if ((posx(i)-xg)<0) then
- do c=0,5
- ip(c)=ip(c)-1
- end do
- end if
- !distance de la particule à remailler au second point (des quatres utilisés pour le remaillage)
- !----------------------------------------------------
- xx1 = (posx(i) - real(ip(2),kind=8)*dx-xg)/dx !relatif
- do c=0,5
- ip(c)=mod(ip(c)+nx,nx)
- end do
+ subroutine remaill_l2_bloc_x (donne,posx,posy,remaille)
+ implicit none
+ real(kind=8),dimension(:),intent(in) :: donne
+ real(kind=8),dimension(:),intent(in) :: posx,posy
+ real(kind=8),dimension(1:nx,1:ny),intent(out) :: remaille
+ integer :: i,c,d,ib,j,per
+ integer,dimension(-2:2) :: ip
+ real(kind=8),dimension(-2:2) :: poids
+ real(kind=8) :: xx,tmp_pos
- !calcul des poids
- !----------------
- x2=xx1**2
- x3=xx1**3
- x4=xx1**4
- x5=xx1**5
- poidx(0)=0.
- poidx(1)=-0.5*xx1+x2-0.5*x3
- poidx(2)=1-5./2.*x2+3.*0.5*x3
- poidx(3)=0.5*xx1+2.*x2-3.*0.5*x3
- poidx(4)=-0.5*x2+0.5*x3
- poidx(5)=0.
+ remaille=0.
- !remaillage à l' interrieur domaine
- !---------------------------------
- do c=0,5
- remaille(ip(c)+1,j)=remaille(ip(c)+1,j)+donne(i)*poidx(c)
- end do
+ do i=1,npart
+ j=nint((posy(i)-yb)/dy)+1
+ poids=0.
- end do
+ !numero des points sur le maillage
+ !--------------------------------
+ ip(0) = floor((posx(i)-xg)/dx)
+ ip(-2) = ip(0) - 2
+ ip(-1) = ip(0) - 1
+ ip(1) = ip(0) + 1
+ ip(2) = ip(0) + 2
- end subroutine remaill_mp4_x
- subroutine remaill_mp4_y(donne,posx,posy,remaille)
- implicit none
- real(kind=8),dimension(:),intent(in) :: donne
- real(kind=8),dimension(:),intent(in) :: posx,posy
- real(kind=8),dimension(1:nx,1:ny),intent(out) :: remaille
- integer :: i,c,d,ib,j,ivar,k
- integer,dimension(0:5) :: ip
- real(kind=8),dimension(0:5) :: poidx
- real(kind=8) :: xx1,tmp_pos,x2,x3,x4,x5
+ !distance de la particule à remailler au point de grille de gauche
+ !-----------------------------------------------------------------
+ xx = (posx(i) - real(ip(0) ,kind=8)*dx-xg)/dx
- remaille=0.
+ !conditions au bord
+ !------------------
+ !periodique:
+ do c=-2,2
+ ip(c)=mod(ip(c)+nx,nx)
+ end do
- !interpole en 2d la vorticite sur la grille en fonction de la position des particules et de leur vorticité
- !---------------------------------------------------------------------------------------------------------
+ !calcul des poids
+ !----------------
- do i=1,npart
+ select case (blocg(i))
- ivar=nint((posx(i)-xg)/dx)+1
+ case(0)
+ if (xx<=0.5) then
+ poids(-1)=0.5*xx*(xx-1.)
+ poids(0)=1.-xx**2
+ else
+ poids(0)=0.5*(1.-xx)*(2.-xx)
+ end if
+ case(1)
+ poids(-1)=0.5*xx*(xx-1.)
+ poids(0)=1.-xx**2
+ case(2)
+ poids(0)=1.-0.5*xx*(1.+xx)
+ case(3)
+ poids(-1)=0.5*xx*(xx-1.)
+ poids(0)=1.-xx
+ case(4)
+ poids(-2)=0.5*xx*(1.+xx)
+ poids(-1)=-xx
+ poids(0)=1.-xx**2
+ case(5)
+ poids(-1)=-0.5*xx+0.5*xx**2
+ poids(0)=1.-poids(-1)
- !numero des points sur le maillage
- !--------------------------------
- ip(2) = int((posy(i)-yb)/dy)
- ip(0) = ip(2) - 2
- ip(1) = ip(2) - 1
- ip(3) = ip(2) + 1
- ip(4) = ip(2) + 2
- ip(5) = ip(2) + 3
-
- if ((posy(i)-yb)<0) then
- do c=0,5
- ip(c)=ip(c)-1
- end do
- end if
-
- !distance de la particule à remailler au second point (des quatres utilisés pour le remaillage)
- !----------------------------------------------------
- xx1 = (posy(i) - real(ip(2),kind=8)*dy-yb)/dy !relatif
-
- do c=0,5
- ip(c)=mod(ip(c)+ny,ny)
- end do
-
-
- !calcul des poids
- !----------------
- x2=xx1**2
- x3=xx1**3
- x4=xx1**4
- x5=xx1**5
-
- poidx(0)=0.
- poidx(1)=-0.5*xx1+x2-0.5*x3
- poidx(2)=1-5./2.*x2+3.*0.5*x3
- poidx(3)=0.5*xx1+2.*x2-3.*0.5*x3
- poidx(4)=-0.5*x2+0.5*x3
- poidx(5)=0.
-
-
- !remaillage à l' interrieur domaine
- !---------------------------------
- do c=0,5
- remaille(ivar,ip(c)+1)=remaille(ivar,ip(c)+1)+donne(i)*poidx(c)
- end do
-
-
- end do
-
- end subroutine remaill_mp4_y
-
-
-
-
-
-
-
-
-
-
-
-
- subroutine remaill_l2_bloc_x (donne,posx,posy,remaille)
- implicit none
- real(kind=8),dimension(:),intent(in) :: donne
- real(kind=8),dimension(:),intent(in) :: posx,posy
- real(kind=8),dimension(1:nx,1:ny),intent(out) :: remaille
- integer :: i,c,d,ib,j,per
- integer,dimension(-2:2) :: ip
- real(kind=8),dimension(-2:2) :: poids
- real(kind=8) :: xx,tmp_pos
-
- remaille=0.
-
- do i=1,npart
-
- j=nint((posy(i)-yb)/dy)+1
- poids=0.
-
- !numero des points sur le maillage
- !--------------------------------
- ip(0) = floor((posx(i)-xg)/dx)
- ip(-2) = ip(0) - 2
- ip(-1) = ip(0) - 1
- ip(1) = ip(0) + 1
- ip(2) = ip(0) + 2
-
-
-
- !distance de la particule à remailler au point de grille de gauche
- !-----------------------------------------------------------------
- xx = (posx(i) - real(ip(0) ,kind=8)*dx-xg)/dx
-
-
- !conditions au bord
- !------------------
- !periodique:
- do c=-2,2
- ip(c)=mod(ip(c)+nx,nx)
- end do
-
- !calcul des poids
- !----------------
-
- select case (blocg(i))
-
- case(0)
- if (xx<=0.5) then
- poids(-1)=0.5*xx*(xx-1.)
- poids(0)=1.-xx**2
- else
- poids(0)=0.5*(1.-xx)*(2.-xx)
- end if
- case(1)
- poids(-1)=0.5*xx*(xx-1.)
- poids(0)=1.-xx**2
- case(2)
- poids(0)=1.-0.5*xx*(1.+xx)
- case(3)
- poids(-1)=0.5*xx*(xx-1.)
- poids(0)=1.-xx
- case(4)
- poids(-2)=0.5*xx*(1.+xx)
- poids(-1)=-xx
- poids(0)=1.-xx**2
- case(5)
- poids(-1)=-0.5*xx+0.5*xx**2
- poids(0)=1.-poids(-1)
-
- end select
+ end select
select case (blocd(i))
@@ -3280,10 +3278,10 @@ contains
r=(u(ib+1)-u(ib))/(u(ib)-u(ib-1))
!anti diffusif
- phi_out(ib)=max(0.,min(tp3(ib),tp1(ib)*r))
+ !phi_out(ib)=max(0.,min(tp3(ib),tp1(ib)*r))
!du type Minmod (attention non justifié pour y>=0.8)
- !phi_out(ib)=max(0.,min(1.,4.*r))
+ phi_out(ib)=max(0.,min(1.,4.*r))
!bornes dépendent de la condition cfl
!phi_out(ib)=max(0.,min(1.,tp3(ib),tp1(ib)*r))
@@ -3307,8 +3305,8 @@ contains
do ib=0,1
r=(u(ib-1)-u(ib-2))/(u(ib)-u(ib-1))
- phi_out(ib)=max(0.,min(tp2(ib-1),tp1(ib-1)*r))
- !phi_out(ib)=max(0.,min(1.,4.*r))
+ !phi_out(ib)=max(0.,min(tp2(ib-1),tp1(ib-1)*r))
+ phi_out(ib)=max(0.,min(1.,4.*r))
!phi_out(ib)=max(0.,min(1.,tp2(ib-1),tp1(ib-1)*r))
if(abs(u(ib)-u(ib-1))<0.0000001) phi_out(ib)=0.
end do
@@ -3579,10 +3577,10 @@ contains
r=(u(ib+1)-u(ib))/(u(ib)-u(ib-1))
!anti diffusif
- phi_out(ib)=max(0.,min(tp3(ib),tp1(ib)*r))
+ !phi_out(ib)=max(0.,min(tp3(ib),tp1(ib)*r))
!du type Minmod (attention non justifié pour y>=0.8)
- !phi_out(ib)=max(0.,min(1.,4.*r))
+ phi_out(ib)=max(0.,min(1.,4.*r))
!bornes dépendent de la condition cfl
!phi_out(ib)=max(0.,min(1.,tp3(ib),tp1(ib)*r))
@@ -3606,8 +3604,8 @@ contains
do ib=0,1
r=(u(ib-1)-u(ib-2))/(u(ib)-u(ib-1))
- phi_out(ib)=max(0.,min(tp2(ib-1),tp1(ib-1)*r))
- !phi_out(ib)=max(0.,min(1.,4.*r))
+ !phi_out(ib)=max(0.,min(tp2(ib-1),tp1(ib-1)*r))
+ phi_out(ib)=max(0.,min(1.,4.*r))
!phi_out(ib)=max(0.,min(1.,tp2(ib-1),tp1(ib-1)*r))
if(abs(u(ib)-u(ib-1))<0.0000001) phi_out(ib)=0.
end do
@@ -3620,185 +3618,1307 @@ contains
end subroutine remaill_l2_bloc_v2_limit_y
- subroutine remaill_l4_bloc_v2_limit_x (donne,posx,posy,remaille)
+ subroutine remaill_l2_limit_x (donne,posx,posy,remaille)
implicit none
real(kind=8),dimension(:),intent(in) :: donne
real(kind=8),dimension(:),intent(in) :: posx,posy
real(kind=8),dimension(1:nx,1:ny),intent(out) :: remaille
integer :: i,c,d,ib,j,per
- integer,dimension(-3:4) :: ip
- real(kind=8),dimension(-3:4) :: poids
- real(kind=8) :: xx,tmp_pos,t1,t2,t3
- real(kind=8) :: phi_m3demi,phi_mdemi,phi_pdemi,phi_p3demi
- real(kind=8),dimension(-1:2) :: phi
-
-
- remaille=0.
+ integer,dimension(-2:2) :: ip
+ real(kind=8),dimension(-2:2) :: poids
+ real(kind=8) :: xx,tmp_pos
+ real(kind=8),dimension(0:2) :: phi_mdemi,phi_pdemi
+ real(kind=8),dimension(0:1) :: phi
+ remaille=0.
do i=1,npart
- j=nint((posy(i)-yb)/dy)+1
+ j=nint((posy(i)-yb)/dy)+1
poids=0.
!numero des points sur le maillage
!--------------------------------
ip(0) = floor((posx(i)-xg)/dx)
- ip(-3) = ip(0) - 3
ip(-2) = ip(0) - 2
ip(-1) = ip(0) - 1
ip(1) = ip(0) + 1
ip(2) = ip(0) + 2
- ip(3) = ip(0) + 3
- ip(4) = ip(0) + 4
!distance de la particule à remailler au point de grille de gauche
!-----------------------------------------------------------------
- xx = (posx(i) - real(ip(0) ,kind=8)*dx-xg)/dx
+ xx = (posx(i) - real(ip(0) ,kind=8)*dx-xg)/dx
!conditions au bord
- !------------------
- do c=-3,4
+ !------------------
+ !periodique:
+ do c=-2,2
ip(c)=mod(ip(c)+nx,nx)
end do
!calcul des poids
- !----------------
+ !----------------
+ !calcul du limiteur pour cas de bloc centre avec xx>0.5
+ !--------------------------------------------------------
+ call calcul_phi_tvd(phi,1)
+ phi_mdemi(1)=phi(0)
+ phi_pdemi(1)=phi(1)
- call calcul_phi_minmax(phi)
- phi_m3demi=phi(-1)
- phi_mdemi=phi(0)
- phi_pdemi=phi(1)
- phi_p3demi=phi(2)
-
- select case (blocg(i))
+ !les autres cas ...
+ !---------------------
+ call calcul_phi_tvd(phi,0)
+ phi_mdemi(0)=phi(0)
+ phi_pdemi(0)=phi(1)
- case(0)
- if (xx<=0.5) then
- poids(-3)=0.
- poids(-2)=ca(xx)
- poids(-1)=cb(xx)
- poids(0)=cc(xx)
- else
- poids(-3)=0.
- poids(-2)=0.
- poids(-1)=ca(xx-1.)
- poids(0)=cb(xx-1.)
- poids(1)=cc(xx-1)
+ phi_mdemi(2)=0.
+ phi_pdemi(2)=0.
+
+ if (xx<=0.5) then
+
+ poids(-1)=ca(xx,0)
+ if (test (i,1) ) then
+ poids(0)=cb(xx,0,1)
+ poids(1)=cc(xx,1)
+ else
+ poids(0)=cb(xx,0,0)
+ poids(1)=cc(xx,0)
+ end if
+
+ else
+
+ poids(0)=ca(xx-1.,1)
+ if (test (i,1) ) then
+ poids(1)=cb(xx-1.,1,1)
+ poids(2)=cc(xx-1.,1)
+ else
+ poids(1)=cb(xx-1.,1,0)
+ poids(2)=cc(xx-1.,0)
+ end if
+
end if
- case(1)
- poids(-3)=0.
- poids(-2)=ca(xx)
- poids(-1)=cb(xx)
- poids(0)=cc(xx)
- case(2)
- poids(-3)=0.
- poids(-2)=0.
- poids(-1)=ca(xx-1.)
- poids(0)=cb(xx-1.)
- poids(1)=cc(xx-1.)+cd(xx-1.)+ce(xx-1.)-ce(xx)
- case(3)
- poids(-3)=0.
- poids(-2)=ca(xx)
- poids(-1)=cb(xx)
- poids(0)=cc(xx)+cd(xx)+ce(xx)-ce(xx+1.)
- case(4)
- poids(-3)=0.
- poids(-2)=0.
- poids(-1)=ca(xx-1.)
- poids(0)=ca(xx)-ca(xx-1.)+cb(xx)+cc(xx)
- case(5)
- poids(-3)=0.
- poids(-2)=0.
- poids(-1)=0.
- poids(0)=ca(xx-1.)+cb(xx-1.)
- poids(1)=cc(xx-1.)
- case(6)
- poids(-3)=0.
- poids(-2)=0.
- poids(-1)=ca(xx)+cb(xx)
- poids(0)=cc(xx)
- case(7)
- poids(-3)=0.
- poids(-2)=ca(xx)
- poids(-1)=cb(xx)
- poids(0)=ca(xx-1.)-ca(xx)+cb(xx-1.)-cb(xx)
- poids(1)=cc(xx-1)
- case(8)
- poids(-3)=ca(xx+1.)
- poids(-2)=cb(xx+1.)
- poids(-1)=ca(xx)-ca(xx+1.)+cb(xx)-cb(xx+1.)
- poids(0)=cc(xx)
- case(9)
- poids(-3)=0.
- poids(-2)=ca(xx)
- poids(-1)=ca(xx-1.)-ca(xx)
- poids(0)=cb(xx-1.)
- poids(1)=cc(xx-1)
- case(10)
- poids(-3)=ca(xx+1.)
- poids(-2)=ca(xx)-ca(xx+1.)
- poids(-1)=cb(xx)
- poids(0)=cc(xx)
-
- end select
- select case (blocd(i))
- case(0)
+
+ !remaillage à l' interrieur domaine
+ !---------------------------------
+ do c=-2,2
+ remaille(ip(c)+1,j)=remaille(ip(c)+1,j)+donne(i)*poids(c)
+ end do
+
+
+ end do
+
+ contains
+
+ function coeff_l2g (num,x) result(c)
+ implicit none
+ integer,intent(in) :: num
+ real(kind=8),intent(in) :: x
+ real(kind=8) :: c
+
+ select case(num)
+ case(-1)
+ c=0.5*x*(x-1.)
+ case(0)
+ c=1.-x**2
+ case(1)
+ c=0.5*x*(1.+x)
+ end select
+
+ end function coeff_l2g
+
+ function coeff_tscg (num,x) result(c)
+ implicit none
+ integer,intent(in) :: num
+ real(kind=8),intent(in) :: x
+ real(kind=8) :: c
+
+ select case(num)
+ case(-1)
+ c=0.5*x*(x-1.) +1./8. ! +1./6.
+ case(0)
+ c=-x**2+1. -2./8. !-2./6.
+ case(1)
+ c=0.5*x*(1.+x) +1./8. ! +1./6.
+ end select
+
+ end function coeff_tscg
+ function ca(x,num)
+ integer :: num
+ real(kind=8) :: x
+ real(kind=8) :: ca
+
+ ca=coeff_tscg(-1,x)+phi_mdemi(num)*(coeff_l2g(-1,x)-coeff_tscg(-1,x))
+
+ end function ca
+
+ function cb(x,num,num2)
+ integer :: num
+ real(kind=8) :: x
+ integer,optional :: num2
+ real(kind=8) :: cb
+
+ if (present (num2)) then
+ cb=1.-(coeff_tscg(-1,x)+coeff_tscg(1,x))-phi_pdemi(num2)*(coeff_l2g(1,x)-coeff_tscg(1,x))-phi_mdemi(num)*(coeff_l2g(-1,x)-coeff_tscg(-1,x))
+ else
+ cb=1.-(coeff_tscg(-1,x)+coeff_tscg(1,x))-phi_pdemi(num)*(coeff_l2g(1,x)-coeff_tscg(1,x))-phi_mdemi(num)*(coeff_l2g(-1,x)-coeff_tscg(-1,x))
+ end if
+
+ end function cb
+
+ function cc(x,num)
+ integer :: num
+ real(kind=8) :: x
+ real(kind=8) :: cc
+
+ cc=coeff_tscg(1,x)+phi_pdemi(num)*(coeff_l2g(1,x)-coeff_tscg(1,x))
+
+
+ end function cc
+
+ function test (ipart,vois)
+ implicit none
+ integer :: ipart,vois
+ logical :: test
+ real(kind=8) :: x
+
+ test=.false.
+ x=(posx(mod(ipart-1+vois+npart,npart)+1) - real( floor((posx(mod(ipart-1+vois+npart,npart)+1)-xg)/dx) ,kind=8)*dx-xg)/dx
+ if ( ((blocg(mod(ipart-1+vois+npart,npart)+1)==0).or.(blocg(mod(ipart-1+vois+npart,npart)+1)==2)).and.(x>0.5)) test=.true.
+
+ end function test
+
+
+
+ subroutine calcul_phi_tvd(phi_out,num)
+ implicit none
+ integer :: num
+ real(kind=8),dimension(0:1),intent(out):: phi_out
+ real(kind=8),dimension(-5:5) :: u
+ real(kind=8) :: rp,rm,diff1,diff2,diff3,diff0,rmm,y,r
+ integer :: ib,jb
+ real(kind=8),dimension(-5:5) :: xx_vois,x,tp1,tp2,tp3
+
+ !particules voisines
+ !-------------------
+ do ib=-2,2
+ u(ib)=donne(mod(i-1+ib+npart,npart)+1)
+ end do
+
+ do ib=-1,2
+ x(ib)=(posx(mod(i-1+ib+npart,npart)+1) - real( floor((posx(mod(i-1+ib+npart,npart)+1)-xg)/dx) ,kind=8)*dx-xg)/dx
+ end do
+
+
+ !cas bloc centre xx>0.5 (psi)
+ !-----------------------
+ if (num==1) then
+
+ do ib=0,1
+ x(ib)=x(ib)-1.
+ tp1(ib)=6.-8.*x(ib)*x(ib)
+ tp3(ib)=4.*(x(ib)-0.5)**2
+ end do
+
+ do ib=0,1
+ r=(u(ib+1)-u(ib))/(u(ib)-u(ib-1))
+
+ !anti diffusif
+ !phi_out(ib)=max(0.,min(tp3(ib),tp1(ib)*r))
+
+ !du type Minmod (attention non justifié pour y>=0.8)
+ phi_out(ib)=max(0.,min(1.,4.*r))
+
+ !bornes dépendent de la condition cfl
+ !phi_out(ib)=max(0.,min(1.,tp3(ib),tp1(ib)*r))
+
+ if(abs(u(ib)-u(ib-1))<0.0000001) phi_out(ib)=0.
+ end do
+
+
+
+
+ else
+ !sinon (phi)
+ !-----------------------
+
+
+ do ib=-1,0
+ tp1(ib)=6.-8.*x(ib)*x(ib)
+ tp2(ib)=4.*(x(ib)+0.5)**2
+ end do
+
+
+ do ib=0,1
+ r=(u(ib-1)-u(ib-2))/(u(ib)-u(ib-1))
+ !phi_out(ib)=max(0.,min(tp2(ib-1),tp1(ib-1)*r))
+ phi_out(ib)=max(0.,min(1.,4.*r))
+ !phi_out(ib)=max(0.,min(1.,tp2(ib-1),tp1(ib-1)*r))
+ if(abs(u(ib)-u(ib-1))<0.0000001) phi_out(ib)=0.
+ end do
+
+ end if
+
+
+ end subroutine calcul_phi_tvd
+
+
+
+ end subroutine remaill_l2_limit_x
+
+ subroutine remaill_l2_limit_y (donne,posx,posy,remaille)
+ implicit none
+ real(kind=8),dimension(:),intent(in) :: donne
+ real(kind=8),dimension(:),intent(in) :: posx,posy
+ real(kind=8),dimension(1:nx,1:ny),intent(out) :: remaille
+ integer :: i,c,d,ib,ivar,per
+ integer,dimension(-2:2) :: ip
+ real(kind=8),dimension(-2:2) :: poids
+ real(kind=8) :: xx,tmp_pos
+ real(kind=8),dimension(0:2):: phi_mdemi,phi_pdemi
+ real(kind=8),dimension(0:1):: phi
+
+
+ remaille=0.
+
+ do i=1,npart
+
+ ivar=nint((posx(i)-xg)/dx)+1
+ poids=0.
+
+ !numero des points sur le maillage
+ !--------------------------------
+ ip(0) = floor((posy(i)-yb)/dy)
+ ip(-2) = ip(0) - 2
+ ip(-1) = ip(0) - 1
+ ip(1) = ip(0) + 1
+ ip(2) = ip(0) + 2
+
+
+
+ !distance de la particule à remailler au point de grille de gauche
+ !-----------------------------------------------------------------
+ xx = (posy(i) - real(ip(0) ,kind=8)*dy-yb)/dy
+
+ !conditions au bord
+ !------------------
+ !periodique:
+ do c=-2,2
+ ip(c)=mod(ip(c)+ny,ny)
+ end do
+
+ !calcul des poids
+ !----------------
+ !calcul du limiteur pour cas de bloc centre avec xx>0.5
+ !--------------------------------------------------------
+ call calcul_phi_tvd(phi,1)
+ phi_mdemi(1)=phi(0)
+ phi_pdemi(1)=phi(1)
+
+ !les autres cas ...
+ !---------------------
+ call calcul_phi_tvd(phi,0)
+ phi_mdemi(0)=phi(0)
+ phi_pdemi(0)=phi(1)
+
+ phi_mdemi(2)=0.
+ phi_pdemi(2)=0.
+
+
+ if (xx<=0.5) then
+
+ poids(-1)=ca(xx,0)
+ if (test (i,1) ) then
+ poids(0)=cb(xx,0,1)
+ poids(1)=cc(xx,1)
+ else
+ poids(0)=cb(xx,0,0)
+ poids(1)=cc(xx,0)
+ end if
+
+ else
+
+ poids(0)=ca(xx-1.,1)
+ if (test (i,1) ) then
+ poids(1)=cb(xx-1.,1,1)
+ poids(2)=cc(xx-1.,1)
+ else
+ poids(1)=cb(xx-1.,1,0)
+ poids(2)=cc(xx-1.,0)
+ end if
+
+ end if
+
+
+
+ !remaillage à l' interrieur domaine
+ !---------------------------------
+ do c=-2,2
+ remaille(ivar,ip(c)+1)=remaille(ivar,ip(c)+1)+donne(i)*poids(c)
+ end do
+
+
+ end do
+ contains
+
+ function coeff_l2g (num,x) result(c)
+ implicit none
+ integer,intent(in) :: num
+ real(kind=8),intent(in) :: x
+ real(kind=8) :: c
+
+ select case(num)
+ case(-1)
+ c=0.5*x*(x-1.)
+ case(0)
+ c=1.-x**2
+ case(1)
+ c=0.5*x*(1.+x)
+ end select
+
+ end function coeff_l2g
+
+ function coeff_tscg (num,x) result(c)
+ implicit none
+ integer,intent(in) :: num
+ real(kind=8),intent(in) :: x
+ real(kind=8) :: c
+
+ select case(num)
+ case(-1)
+ c=0.5*x*(x-1.) +1./8. ! +1./6.
+ case(0)
+ c=-x**2+1. -2./8. !-2./6.
+ case(1)
+ c=0.5*x*(1.+x) +1./8. ! +1./6.
+ end select
+
+ end function coeff_tscg
+
+ function ca(x,num)
+ integer :: num
+ real(kind=8) :: x
+ real(kind=8) :: ca
+
+ ca=coeff_tscg(-1,x)+phi_mdemi(num)*(coeff_l2g(-1,x)-coeff_tscg(-1,x))
+
+ end function ca
+
+ function cb(x,num,num2)
+ integer :: num
+ real(kind=8) :: x
+ integer,optional :: num2
+ real(kind=8) :: cb
+
+ if (present (num2)) then
+ cb=1.-(coeff_tscg(-1,x)+coeff_tscg(1,x))-phi_pdemi(num2)*(coeff_l2g(1,x)-coeff_tscg(1,x))-phi_mdemi(num)*(coeff_l2g(-1,x)-coeff_tscg(-1,x))
+ else
+ cb=1.-(coeff_tscg(-1,x)+coeff_tscg(1,x))-phi_pdemi(num)*(coeff_l2g(1,x)-coeff_tscg(1,x))-phi_mdemi(num)*(coeff_l2g(-1,x)-coeff_tscg(-1,x))
+ end if
+
+ end function cb
+
+ function cc(x,num)
+ integer :: num
+ real(kind=8) :: x
+ real(kind=8) :: cc
+
+ cc=coeff_tscg(1,x)+phi_pdemi(num)*(coeff_l2g(1,x)-coeff_tscg(1,x))
+
+
+ end function cc
+
+ function test (ipart,vois)
+ implicit none
+ integer :: ipart,vois
+ logical :: test
+ real(kind=8) :: x
+
+ test=.false.
+ x=(posy(mod(ipart-1+vois+npart,npart)+1) - real( floor((posy(mod(ipart-1+vois+npart,npart)+1)-yb)/dy) ,kind=8)*dy-yb)/dy
+ if ( ((blocg(mod(ipart-1+vois+npart,npart)+1)==0).or.(blocg(mod(ipart-1+vois+npart,npart)+1)==2)).and.(x>0.5)) test=.true.
+
+ end function test
+
+
+
+ subroutine calcul_phi_tvd(phi_out,num)
+ implicit none
+ integer :: num
+ real(kind=8),dimension(0:1),intent(out):: phi_out
+ real(kind=8),dimension(-5:5) :: u
+ real(kind=8) :: rp,rm,diff1,diff2,diff3,diff0,rmm,y,r
+ integer :: ib,jb
+ real(kind=8),dimension(-5:5) :: xx_vois,x,tp1,tp2,tp3
+
+ !particules voisines
+ !-------------------
+ do ib=-2,2
+ u(ib)=donne(mod(i-1+ib+npart,npart)+1)
+ end do
+
+ do ib=-1,2
+ x(ib)=(posx(mod(i-1+ib+npart,npart)+1) - real( floor((posx(mod(i-1+ib+npart,npart)+1)-xg)/dx) ,kind=8)*dx-xg)/dx
+ end do
+
+
+ !cas bloc centre xx>0.5 (psi)
+ !-----------------------
+ if (num==1) then
+
+ do ib=0,1
+ x(ib)=x(ib)-1.
+ tp1(ib)=6.-8.*x(ib)*x(ib)
+ tp3(ib)=4.*(x(ib)-0.5)**2
+ end do
+
+ do ib=0,1
+ r=(u(ib+1)-u(ib))/(u(ib)-u(ib-1))
+
+ !anti diffusif
+ !phi_out(ib)=max(0.,min(tp3(ib),tp1(ib)*r))
+
+ !du type Minmod (attention non justifié pour y>=0.8)
+ phi_out(ib)=max(0.,min(1.,4.*r))
+
+ !bornes dépendent de la condition cfl
+ !phi_out(ib)=max(0.,min(1.,tp3(ib),tp1(ib)*r))
+
+ if(abs(u(ib)-u(ib-1))<0.0000001) phi_out(ib)=0.
+ end do
+
+
+
+
+ else
+ !sinon (phi)
+ !-----------------------
+
+
+ do ib=-1,0
+ tp1(ib)=6.-8.*x(ib)*x(ib)
+ tp2(ib)=4.*(x(ib)+0.5)**2
+ end do
+
+
+ do ib=0,1
+ r=(u(ib-1)-u(ib-2))/(u(ib)-u(ib-1))
+ !phi_out(ib)=max(0.,min(tp2(ib-1),tp1(ib-1)*r))
+ phi_out(ib)=max(0.,min(1.,4.*r))
+ !phi_out(ib)=max(0.,min(1.,tp2(ib-1),tp1(ib-1)*r))
+ if(abs(u(ib)-u(ib-1))<0.0000001) phi_out(ib)=0.
+ end do
+
+ end if
+
+
+ end subroutine calcul_phi_tvd
+
+ end subroutine remaill_l2_limit_y
+
+
+
+
+
+ subroutine remaill_l4_bloc_v2_limit_x (donne,posx,posy,remaille)
+ implicit none
+ real(kind=8),dimension(:),intent(in) :: donne
+ real(kind=8),dimension(:),intent(in) :: posx,posy
+ real(kind=8),dimension(1:nx,1:ny),intent(out) :: remaille
+ integer :: i,c,d,ib,j,per
+ integer,dimension(-3:4) :: ip
+ real(kind=8),dimension(-3:4) :: poids
+ real(kind=8) :: xx,tmp_pos,t1,t2,t3
+ real(kind=8) :: phi_m3demi,phi_mdemi,phi_pdemi,phi_p3demi
+ real(kind=8),dimension(-1:2) :: phi
+
+
+ remaille=0.
+
+
+ do i=1,npart
+
+ j=nint((posy(i)-yb)/dy)+1
+ poids=0.
+
+ !numero des points sur le maillage
+ !--------------------------------
+ ip(0) = floor((posx(i)-xg)/dx)
+ ip(-3) = ip(0) - 3
+ ip(-2) = ip(0) - 2
+ ip(-1) = ip(0) - 1
+ ip(1) = ip(0) + 1
+ ip(2) = ip(0) + 2
+ ip(3) = ip(0) + 3
+ ip(4) = ip(0) + 4
+
+
+
+ !distance de la particule à remailler au point de grille de gauche
+ !-----------------------------------------------------------------
+ xx = (posx(i) - real(ip(0) ,kind=8)*dx-xg)/dx
+
+
+ !conditions au bord
+ !------------------
+ do c=-3,4
+ ip(c)=mod(ip(c)+nx,nx)
+ end do
+
+ !calcul des poids
+ !----------------
+
+ call calcul_phi_minmax(phi)
+ phi_m3demi=phi(-1)
+ phi_mdemi=phi(0)
+ phi_pdemi=phi(1)
+ phi_p3demi=phi(2)
+
+ select case (blocg(i))
+
+ case(0)
+ if (xx<=0.5) then
+ poids(-3)=0.
+ poids(-2)=ca(xx)
+ poids(-1)=cb(xx)
+ poids(0)=cc(xx)
+ else
+ poids(-3)=0.
+ poids(-2)=0.
+ poids(-1)=ca(xx-1.)
+ poids(0)=cb(xx-1.)
+ poids(1)=cc(xx-1)
+ end if
+ case(1)
+ poids(-3)=0.
+ poids(-2)=ca(xx)
+ poids(-1)=cb(xx)
+ poids(0)=cc(xx)
+ case(2)
+ poids(-3)=0.
+ poids(-2)=0.
+ poids(-1)=ca(xx-1.)
+ poids(0)=cb(xx-1.)
+ poids(1)=cc(xx-1.)+cd(xx-1.)+ce(xx-1.)-ce(xx)
+ case(3)
+ poids(-3)=0.
+ poids(-2)=ca(xx)
+ poids(-1)=cb(xx)
+ poids(0)=cc(xx)+cd(xx)+ce(xx)-ce(xx+1.)
+ case(4)
+ poids(-3)=0.
+ poids(-2)=0.
+ poids(-1)=ca(xx-1.)
+ poids(0)=ca(xx)-ca(xx-1.)+cb(xx)+cc(xx)
+ case(5)
+ poids(-3)=0.
+ poids(-2)=0.
+ poids(-1)=0.
+ poids(0)=ca(xx-1.)+cb(xx-1.)
+ poids(1)=cc(xx-1.)
+ case(6)
+ poids(-3)=0.
+ poids(-2)=0.
+ poids(-1)=ca(xx)+cb(xx)
+ poids(0)=cc(xx)
+ case(7)
+ poids(-3)=0.
+ poids(-2)=ca(xx)
+ poids(-1)=cb(xx)
+ poids(0)=ca(xx-1.)-ca(xx)+cb(xx-1.)-cb(xx)
+ poids(1)=cc(xx-1)
+ case(8)
+ poids(-3)=ca(xx+1.)
+ poids(-2)=cb(xx+1.)
+ poids(-1)=ca(xx)-ca(xx+1.)+cb(xx)-cb(xx+1.)
+ poids(0)=cc(xx)
+ case(9)
+ poids(-3)=0.
+ poids(-2)=ca(xx)
+ poids(-1)=ca(xx-1.)-ca(xx)
+ poids(0)=cb(xx-1.)
+ poids(1)=cc(xx-1)
+ case(10)
+ poids(-3)=ca(xx+1.)
+ poids(-2)=ca(xx)-ca(xx+1.)
+ poids(-1)=cb(xx)
+ poids(0)=cc(xx)
+
+ end select
+
+ select case (blocd(i))
+ case(0)
+ if (xx<=0.5) then
+ poids(1)=cd(xx)
+ poids(2)=ce(xx)
+ poids(3)=0.
+ poids(4)=0.
+ else
+ poids(2)=cd(xx-1.)
+ poids(3)=ce(xx-1.)
+ poids(4)=0.
+ end if
+ case(1)
+ poids(1)=cd(xx)
+ poids(2)=ce(xx)
+ poids(3)=0.
+ poids(4)=0.
+ case(2)
+ poids(2)=cd(xx-1.)+ce(xx-1.)
+ poids(3)=0.
+ poids(4)=0.
+ case(3)
+ poids(1)=cd(xx)+ce(xx)
+ poids(2)=0.
+ poids(3)=0.
+ poids(4)=0.
+ case(4)
+ poids(2)=ce(xx)
+ poids(3)=0.
+ poids(4)=0.
+ case(5)
+ poids(1)=ce(xx+1.)
+ poids(2)=0.
+ poids(3)=0.
+ poids(4)=0.
+ case(6)
+ poids(2)=cd(xx-1.)
+ poids(3)=ce(xx-1.)-ce(xx-2.)
+ poids(4)=ce(xx-2.)
+ case(7)
+ poids(1)=cd(xx)
+ poids(2)=ce(xx)-ce(xx-1.)
+ poids(3)=ce(xx-1.)
+ poids(4)=0.
+ case(8)
+ poids(1)=cd(xx)-cd(xx-1.)+ce(xx)-ce(xx-1.)
+ poids(2)=cd(xx-1.)
+ poids(3)=ce(xx-1.)
+ poids(4)=0.
+
+
+ end select
+
+
+ !remaillage à l' interrieur domaine
+ !---------------------------------
+ do c=-3,4
+ remaille(ip(c)+1,j)=remaille(ip(c)+1,j)+donne(i)*poids(c)
+ end do
+
+
+ end do
+
+ contains
+
+ function ca(x)
+ implicit none
+ real(kind=8),intent(in) :: x
+ real(kind=8) :: ca
+ !ca=x*(x-1.)*(x-2.)*(x+1.)/24.
+ ca=coeff_M5g(-2,x)+phi_m3demi*coeff_l4mM5g(-2,x)
+ end function ca
+ function cb(x)
+ implicit none
+ real(kind=8),intent(in) :: x
+ real(kind=8) :: cb
+ !cb=-x*(x-1.)*(x**2-4.)/6.
+ cb=coeff_M5g(-1,x)-phi_m3demi*coeff_l4mM5g(-2,x)+phi_mdemi*(coeff_l4mM5g(-2,x)+coeff_l4mM5g(-1,x))
+ end function cb
+ function cc(x)
+ implicit none
+ real(kind=8),intent(in) :: x
+ real(kind=8) :: cc
+ !cc=(x**2-4.)*(x**2-1.)/4.
+ cc=coeff_M5g(0,x)-phi_mdemi*(coeff_l4mM5g(-2,x)+coeff_l4mM5g(-1,x))-phi_pdemi*(coeff_l4mM5g(1,x)+coeff_l4mM5g(2,x))
+ end function cc
+ function cd(x)
+ implicit none
+ real(kind=8),intent(in) :: x
+ real(kind=8) :: cd
+ !cd=-x*(x+1.)*(x+2.)*(x-2.)/6.
+ cd=coeff_M5g(1,x)+phi_pdemi*(coeff_l4mM5g(1,x)+coeff_l4mM5g(2,x))-phi_p3demi*coeff_l4mM5g(2,x)
+ end function cd
+ function ce(x)
+ implicit none
+ real(kind=8),intent(in) :: x
+ real(kind=8) :: ce
+ !ce=x*(x+1.)*(x+2.)*(x-1.)/24.
+ ce=coeff_M5g(2,x)+phi_p3demi*coeff_l4mM5g(2,x)
+ end function ce
+
+
+ function coeff_l4g (num,x) result(c)
+ implicit none
+ integer,intent(in) :: num
+ real(kind=8),intent(in) :: x
+ real(kind=8) :: c
+
+ select case(num)
+ case(-2)
+ c=x*(x-1)*(x-2)*(x+1)/24.
+ case(-1)
+ c=-x*(x-1)*(x-2)*(x+2)/6.
+ case(0)
+ c=(x-2)*(x-1)*(x+1)*(x+2)/4.
+ case(1)
+ c=-x*(x-2)*(x+1)*(x+2)/6.
+ case(2)
+ c=x*(x-1)*(x+1)*(x+2)/24.
+ end select
+
+ end function coeff_l4g
+
+ function coeff_M5g (num,x) result(c)
+ implicit none
+ integer,intent(in) :: num
+ real(kind=8),intent(in) :: x
+ real(kind=8) :: c
+
+ select case(num)
+ case(-2)
+ c=(2.*x-1.)**4/384.
+ case(-1)
+ c=19./96.-11/24.*x+x**2/4.+x**3/6.-x**4/6.
+ case(0)
+ c=115./192.-5./8.*x**2+x**4/4.
+ case(1)
+ c=19./96.+11/24.*x-x**3/6.+x**2/4.-x**4/6.
+ case(2)
+ c=(2.*x+1.)**4/384.
+ end select
+
+ end function coeff_M5g
+
+ function coeff_l4mM5g (num,x) result(c)
+ !lambda4 - M5 gauche
+ implicit none
+ integer,intent(in) :: num
+ real(kind=8),intent(in) :: x
+ real(kind=8) :: c
+
+ select case(num)
+ case(-2)
+ c=5.*x/48.*(1.-x)-1./384.
+ case(-1)
+ c=-5.*x/24.*(1.-2.*x)-19/96.
+ case(0)
+ c=-5.*x**2/8.+77./192.
+ case(1)
+ c=5.*x/24.*(1.+2.*x)-19./96.
+ case(2)
+ c=-5.*x/48.*(1.+x)-1./384.
+ end select
+
+ end function coeff_l4mM5g
+
+
+
+
+ subroutine calcul_phi_minmax(phi_out)
+ implicit none
+ real(kind=8),dimension(-1:2),intent(out):: phi_out
+ real(kind=8),dimension(-5:5) :: u,bmin,bmax,x,xx_vois
+ real(kind=8),dimension(-5:5,-5:5) :: ej,fj
+ real(kind=8) :: tp1,tp2,phi1,phi2
+ integer :: ib,jb
+
+
+ !particules voisines
+ !-------------------
+ do ib=-4,5
+ u(ib)=donne(mod(i-1+ib+npart,npart)+1)
+ end do
+
+ !coeff decentre droite = coeff decentre gauche -1
+ !-------------------------------------------------
+ do ib=-4,3
+ xx_vois(ib)=(posx(mod(i-1+ib+npart,npart)+1) - real( floor((posx(mod(i-1+ib+npart,npart)+1)-xg)/dx) ,kind=8)*dx-xg)/dx
+
+ if ((xx_vois(ib)>0.5).and.((blocg(mod(i-1+ib+npart,npart)+1)==0).or.(blocd(mod(i-1+ib+npart,npart)+1)==0))) then
+ x(ib)=xx_vois(ib)-1.
+ else
+ x(ib)=xx_vois(ib)
+ end if
+ end do
+
+
+ do ib=-4,3
+
+ bmin(ib)=min(u(ib-2),u(ib-1),u(ib),u(ib+1),u(ib+2))
+ bmax(ib)=max(u(ib-2),u(ib-1),u(ib),u(ib+1),u(ib+2))
+
+
+ do jb=-4,3
+ ej(ib,jb)=coeff_M5g(-2,x(jb))*u(ib+2)+coeff_M5g(-1,x(jb))*u(ib+1)+coeff_M5g(0,x(jb))*u(ib)+coeff_M5g(1,x(jb))*u(ib-1)+coeff_M5g(2,x(jb))*u(ib-2)
+ fj(ib,jb)=coeff_l4mM5g(-2,x(jb))*u(ib+1)+(coeff_l4mM5g(-2,x(jb))+coeff_l4mM5g(-1,x(jb)))*u(ib)-(coeff_l4mM5g(1,x(jb))+coeff_l4mM5g(2,x(jb)))*u(ib-1)-coeff_l4mM5g(2,x(jb))*u(ib-2)
+ end do
+
+ end do
+
+
+ do ib=-1,2
+
+ tp1=(bmin(ib-1)-ej(ib-1,ib)+max(0.,fj(ib-1,ib)))/fj(ib,ib)
+ tp2=(bmax(ib-1)-ej(ib-1,ib)+min(0.,fj(ib-1,ib)))/fj(ib,ib)
+ !phi1=max(0.,min(tp1,tp2),min(1.,max(tp1,tp2)))
+ phi1=max(0.,min(1.,max(tp1,tp2)))
+
+ tp1=(bmin(ib)-ej(ib,ib-1)+max(0.,-fj(ib+1,ib-1)))/(-fj(ib,ib-1))
+ tp2=(bmax(ib)-ej(ib,ib-1)+min(0.,-fj(ib+1,ib-1)))/(-fj(ib,ib-1))
+ !phi2=max(0.,min(tp1,tp2),min(1.,max(tp1,tp2)))
+ phi2=max(0.,min(1.,max(tp1,tp2)))
+
+ phi_out(ib)=min(phi1,phi2)
+
+ if(abs(fj(ib,ib))<0.0000001) phi_out(ib)=0.
+ if(abs(fj(ib,ib-1))<0.0000001) phi_out(ib)=0.
+
+ end do
+
+
+ end subroutine calcul_phi_minmax
+
+
+ end subroutine remaill_l4_bloc_v2_limit_x
+
+
+
+
+ subroutine remaill_l4_bloc_v2_limit_y (donne,posx,posy,remaille)
+ implicit none
+ real(kind=8),dimension(:),intent(in) :: donne
+ real(kind=8),dimension(:),intent(in) :: posx,posy
+ real(kind=8),dimension(1:nx,1:ny),intent(out) :: remaille
+ integer :: i,c,d,ib,ivar,per
+ integer,dimension(-3:4) :: ip
+ real(kind=8),dimension(-3:4) :: poids
+ real(kind=8) :: xx,tmp_pos,t1,t2,t3
+ real(kind=8) :: phi_m3demi,phi_mdemi,phi_pdemi,phi_p3demi
+ real(kind=8),dimension(-1:2) :: phi
+
+
+ remaille=0.
+
+ do i=1,npart
+
+ poids=0.
+ ivar=nint((posx(i)-xg)/dx)+1
+
+ !numero des points sur le maillage
+ !--------------------------------
+ ip(0) = floor((posy(i)-yb)/dy)
+ ip(-3) = ip(0) - 3
+ ip(-2) = ip(0) - 2
+ ip(-1) = ip(0) - 1
+ ip(1) = ip(0) + 1
+ ip(2) = ip(0) + 2
+ ip(3) = ip(0) + 3
+ ip(4) = ip(0) + 4
+
+
+
+ !distance de la particule à remailler au point de grille de gauche
+ !-----------------------------------------------------------------
+ xx = (posy(i) - real(ip(0) ,kind=8)*dy-yb)/dy
+
+ !conditions au bord
+ !------------------
+ do c=-3,4
+ ip(c)=mod(ip(c)+ny,ny)
+ end do
+
+ !calcul des poidss
+ !----------------
+ call calcul_phi_minmax(phi)
+ phi_m3demi=phi(-1)
+ phi_mdemi=phi(0)
+ phi_pdemi=phi(1)
+ phi_p3demi=phi(2)
+ select case (blocg(i))
+
+ case(0)
+ if (xx<=0.5) then
+ poids(-3)=0.
+ poids(-2)=ca(xx)
+ poids(-1)=cb(xx)
+ poids(0)=cc(xx)
+ else
+ poids(-3)=0.
+ poids(-2)=0.
+ poids(-1)=ca(xx-1.)
+ poids(0)=cb(xx-1.)
+ poids(1)=cc(xx-1)
+ end if
+ case(1)
+ poids(-3)=0.
+ poids(-2)=ca(xx)
+ poids(-1)=cb(xx)
+ poids(0)=cc(xx)
+ case(2)
+ poids(-3)=0.
+ poids(-2)=0.
+ poids(-1)=ca(xx-1.)
+ poids(0)=cb(xx-1.)
+ poids(1)=cc(xx-1.)+cd(xx-1.)+ce(xx-1.)-ce(xx)
+ case(3)
+ poids(-3)=0.
+ poids(-2)=ca(xx)
+ poids(-1)=cb(xx)
+ poids(0)=cc(xx)+cd(xx)+ce(xx)-ce(xx+1.)
+ case(4)
+ poids(-3)=0.
+ poids(-2)=0.
+ poids(-1)=ca(xx-1.)
+ poids(0)=ca(xx)-ca(xx-1.)+cb(xx)+cc(xx)
+ case(5)
+ poids(-3)=0.
+ poids(-2)=0.
+ poids(-1)=0.
+ poids(0)=ca(xx-1.)+cb(xx-1.)
+ poids(1)=cc(xx-1.)
+ case(6)
+ poids(-3)=0.
+ poids(-2)=0.
+ poids(-1)=ca(xx)+cb(xx)
+ poids(0)=cc(xx)
+ case(7)
+ poids(-3)=0.
+ poids(-2)=ca(xx)
+ poids(-1)=cb(xx)
+ poids(0)=ca(xx-1.)-ca(xx)+cb(xx-1.)-cb(xx)
+ poids(1)=cc(xx-1)
+ case(8)
+ poids(-3)=ca(xx+1.)
+ poids(-2)=cb(xx+1.)
+ poids(-1)=ca(xx)-ca(xx+1.)+cb(xx)-cb(xx+1.)
+ poids(0)=cc(xx)
+ case(9)
+ poids(-3)=0.
+ poids(-2)=ca(xx)
+ poids(-1)=ca(xx-1.)-ca(xx)
+ poids(0)=cb(xx-1.)
+ poids(1)=cc(xx-1)
+ case(10)
+ poids(-3)=ca(xx+1.)
+ poids(-2)=ca(xx)-ca(xx+1.)
+ poids(-1)=cb(xx)
+ poids(0)=cc(xx)
+
+ end select
+
+ select case (blocd(i))
+ case(0)
if (xx<=0.5) then
poids(1)=cd(xx)
poids(2)=ce(xx)
poids(3)=0.
- poids(4)=0.
- else
+ poids(4)=0.
+ else
+ poids(2)=cd(xx-1.)
+ poids(3)=ce(xx-1.)
+ poids(4)=0.
+ end if
+ case(1)
+ poids(1)=cd(xx)
+ poids(2)=ce(xx)
+ poids(3)=0.
+ poids(4)=0.
+ case(2)
+ poids(2)=cd(xx-1.)+ce(xx-1.)
+ poids(3)=0.
+ poids(4)=0.
+ case(3)
+ poids(1)=cd(xx)+ce(xx)
+ poids(2)=0.
+ poids(3)=0.
+ poids(4)=0.
+ case(4)
+ poids(2)=ce(xx)
+ poids(3)=0.
+ poids(4)=0.
+ case(5)
+ poids(1)=ce(xx+1.)
+ poids(2)=0.
+ poids(3)=0.
+ poids(4)=0.
+ case(6)
+ poids(2)=cd(xx-1.)
+ poids(3)=ce(xx-1.)-ce(xx-2.)
+ poids(4)=ce(xx-2.)
+ case(7)
+ poids(1)=cd(xx)
+ poids(2)=ce(xx)-ce(xx-1.)
+ poids(3)=ce(xx-1.)
+ poids(4)=0.
+ case(8)
+ poids(1)=cd(xx)-cd(xx-1.)+ce(xx)-ce(xx-1.)
+ poids(2)=cd(xx-1.)
+ poids(3)=ce(xx-1.)
+ poids(4)=0.
+
+
+ end select
+ !remaillage à l' interrieur domaine
+ !---------------------------------
+ do c=-3,4
+ remaille(ivar,ip(c)+1)=remaille(ivar,ip(c)+1)+donne(i)*poids(c)
+ end do
+
+
+ end do
+ contains
+
+ function ca(x)
+ implicit none
+ real(kind=8),intent(in) :: x
+ real(kind=8) :: ca
+ !ca=x*(x-1.)*(x-2.)*(x+1.)/24.
+ ca=coeff_M5g(-2,x)+phi_m3demi*coeff_l4mM5g(-2,x)
+ end function ca
+ function cb(x)
+ implicit none
+ real(kind=8),intent(in) :: x
+ real(kind=8) :: cb
+ !cb=-x*(x-1.)*(x**2-4.)/6.
+ cb=coeff_M5g(-1,x)-phi_m3demi*coeff_l4mM5g(-2,x)+phi_mdemi*(coeff_l4mM5g(-2,x)+coeff_l4mM5g(-1,x))
+ end function cb
+ function cc(x)
+ implicit none
+ real(kind=8),intent(in) :: x
+ real(kind=8) :: cc
+ !cc=(x**2-4.)*(x**2-1.)/4.
+ cc=coeff_M5g(0,x)-phi_mdemi*(coeff_l4mM5g(-2,x)+coeff_l4mM5g(-1,x))-phi_pdemi*(coeff_l4mM5g(1,x)+coeff_l4mM5g(2,x))
+ end function cc
+ function cd(x)
+ implicit none
+ real(kind=8),intent(in) :: x
+ real(kind=8) :: cd
+ !cd=-x*(x+1.)*(x+2.)*(x-2.)/6.
+ cd=coeff_M5g(1,x)+phi_pdemi*(coeff_l4mM5g(1,x)+coeff_l4mM5g(2,x))-phi_p3demi*coeff_l4mM5g(2,x)
+ end function cd
+ function ce(x)
+ implicit none
+ real(kind=8),intent(in) :: x
+ real(kind=8) :: ce
+ !ce=x*(x+1.)*(x+2.)*(x-1.)/24.
+ ce=coeff_M5g(2,x)+phi_p3demi*coeff_l4mM5g(2,x)
+ end function ce
+
+
+ function coeff_l4g (num,x) result(c)
+ implicit none
+ integer,intent(in) :: num
+ real(kind=8),intent(in) :: x
+ real(kind=8) :: c
+
+ select case(num)
+ case(-2)
+ c=x*(x-1)*(x-2)*(x+1)/24.
+ case(-1)
+ c=-x*(x-1)*(x-2)*(x+2)/6.
+ case(0)
+ c=(x-2)*(x-1)*(x+1)*(x+2)/4.
+ case(1)
+ c=-x*(x-2)*(x+1)*(x+2)/6.
+ case(2)
+ c=x*(x-1)*(x+1)*(x+2)/24.
+ end select
+
+ end function coeff_l4g
+
+ function coeff_M5g (num,x) result(c)
+ implicit none
+ integer,intent(in) :: num
+ real(kind=8),intent(in) :: x
+ real(kind=8) :: c
+
+ select case(num)
+ case(-2)
+ c=(2.*x-1.)**4/384.
+ case(-1)
+ c=19./96.-11/24.*x+x**2/4.+x**3/6.-x**4/6.
+ case(0)
+ c=115./192.-5./8.*x**2+x**4/4.
+ case(1)
+ c=19./96.+11/24.*x-x**3/6.+x**2/4.-x**4/6.
+ case(2)
+ c=(2.*x+1.)**4/384.
+ end select
+
+ end function coeff_M5g
+
+ function coeff_l4mM5g (num,x) result(c)
+ !lambda4 - M5 gauche
+ implicit none
+ integer,intent(in) :: num
+ real(kind=8),intent(in) :: x
+ real(kind=8) :: c
+
+ select case(num)
+ case(-2)
+ c=5.*x/48.*(1.-x)-1./384.
+ case(-1)
+ c=-5.*x/24.*(1.-2.*x)-19/96.
+ case(0)
+ c=-5.*x**2/8.+77./192.
+ case(1)
+ c=5.*x/24.*(1.+2.*x)-19./96.
+ case(2)
+ c=-5.*x/48.*(1.+x)-1./384.
+ end select
+
+ end function coeff_l4mM5g
+
+
+
+
+ subroutine calcul_phi_minmax(phi_out)
+ implicit none
+ real(kind=8),dimension(-1:2),intent(out):: phi_out
+ real(kind=8),dimension(-5:5) :: u,bmin,bmax,x,xx_vois
+ real(kind=8),dimension(-5:5,-5:5) :: ej,fj
+ real(kind=8) :: tp1,tp2,phi1,phi2
+ integer :: ib,jb
+
+
+ !particules voisines
+ !-------------------
+ do ib=-4,5
+ u(ib)=donne(mod(i-1+ib+npart,npart)+1)
+ end do
+
+ !coeff decentre droite = coeff decentre gauche -1
+ !-------------------------------------------------
+ do ib=-4,3
+ xx_vois(ib)=(posx(mod(i-1+ib+npart,npart)+1) - real( floor((posx(mod(i-1+ib+npart,npart)+1)-xg)/dx) ,kind=8)*dx-xg)/dx
+
+ if ((xx_vois(ib)>0.5).and.((blocg(mod(i-1+ib+npart,npart)+1)==0).or.(blocd(mod(i-1+ib+npart,npart)+1)==0))) then
+ x(ib)=xx_vois(ib)-1.
+ else
+ x(ib)=xx_vois(ib)
+ end if
+ end do
+
+
+ do ib=-4,3
+
+ bmin(ib)=min(u(ib-2),u(ib-1),u(ib),u(ib+1),u(ib+2))
+ bmax(ib)=max(u(ib-2),u(ib-1),u(ib),u(ib+1),u(ib+2))
+
+
+ do jb=-4,3
+ ej(ib,jb)=coeff_M5g(-2,x(jb))*u(ib+2)+coeff_M5g(-1,x(jb))*u(ib+1)+coeff_M5g(0,x(jb))*u(ib)+coeff_M5g(1,x(jb))*u(ib-1)+coeff_M5g(2,x(jb))*u(ib-2)
+ fj(ib,jb)=coeff_l4mM5g(-2,x(jb))*u(ib+1)+(coeff_l4mM5g(-2,x(jb))+coeff_l4mM5g(-1,x(jb)))*u(ib)-(coeff_l4mM5g(1,x(jb))+coeff_l4mM5g(2,x(jb)))*u(ib-1)-coeff_l4mM5g(2,x(jb))*u(ib-2)
+ end do
+
+ end do
+
+
+ do ib=-1,2
+
+ tp1=(bmin(ib-1)-ej(ib-1,ib)+max(0.,fj(ib-1,ib)))/fj(ib,ib)
+ tp2=(bmax(ib-1)-ej(ib-1,ib)+min(0.,fj(ib-1,ib)))/fj(ib,ib)
+ !phi1=max(0.,min(tp1,tp2),min(1.,max(tp1,tp2)))
+ phi1=max(0.,min(1.,max(tp1,tp2)))
+
+ tp1=(bmin(ib)-ej(ib,ib-1)+max(0.,-fj(ib+1,ib-1)))/(-fj(ib,ib-1))
+ tp2=(bmax(ib)-ej(ib,ib-1)+min(0.,-fj(ib+1,ib-1)))/(-fj(ib,ib-1))
+ !phi2=max(0.,min(tp1,tp2),min(1.,max(tp1,tp2)))
+ phi2=max(0.,min(1.,max(tp1,tp2)))
+
+ phi_out(ib)=min(phi1,phi2)
+
+ if(abs(fj(ib,ib))<0.0000001) phi_out(ib)=0.
+ if(abs(fj(ib,ib-1))<0.0000001) phi_out(ib)=0.
+
+ end do
+
+
+ end subroutine calcul_phi_minmax
+
+ end subroutine remaill_l4_bloc_v2_limit_y
+
+
+
+
+ subroutine remaill_l4_limit_x (donne,posx,posy,remaille)
+ implicit none
+ real(kind=8),dimension(:),intent(in) :: donne
+ real(kind=8),dimension(:),intent(in) :: posx,posy
+ real(kind=8),dimension(1:nx,1:ny),intent(out) :: remaille
+ integer :: i,c,d,ib,j,per
+ integer,dimension(-3:4) :: ip
+ real(kind=8),dimension(-3:4) :: poids
+ real(kind=8) :: xx,tmp_pos,t1,t2,t3
+ real(kind=8) :: phi_m3demi,phi_mdemi,phi_pdemi,phi_p3demi
+ real(kind=8),dimension(-1:2) :: phi
+
+
+ remaille=0.
+
+
+ do i=1,npart
+
+ j=nint((posy(i)-yb)/dy)+1
+ poids=0.
+
+ !numero des points sur le maillage
+ !--------------------------------
+ ip(0) = floor((posx(i)-xg)/dx)
+ ip(-3) = ip(0) - 3
+ ip(-2) = ip(0) - 2
+ ip(-1) = ip(0) - 1
+ ip(1) = ip(0) + 1
+ ip(2) = ip(0) + 2
+ ip(3) = ip(0) + 3
+ ip(4) = ip(0) + 4
+
+
+
+ !distance de la particule à remailler au point de grille de gauche
+ !-----------------------------------------------------------------
+ xx = (posx(i) - real(ip(0) ,kind=8)*dx-xg)/dx
+
+
+ !conditions au bord
+ !------------------
+ do c=-3,4
+ ip(c)=mod(ip(c)+nx,nx)
+ end do
+
+ !calcul des poids
+ !----------------
+
+ call calcul_phi_minmax(phi)
+ phi_m3demi=phi(-1)
+ phi_mdemi=phi(0)
+ phi_pdemi=phi(1)
+ phi_p3demi=phi(2)
+
+ if (xx<=0.5) then
+ poids(-3)=0.
+ poids(-2)=ca(xx)
+ poids(-1)=cb(xx)
+ poids(0)=cc(xx)
+ poids(1)=cd(xx)
+ poids(2)=ce(xx)
+ poids(3)=0.
+ poids(4)=0.
+ else
+ poids(-3)=0.
+ poids(-2)=0.
+ poids(-1)=ca(xx-1.)
+ poids(0)=cb(xx-1.)
+ poids(1)=cc(xx-1)
poids(2)=cd(xx-1.)
poids(3)=ce(xx-1.)
poids(4)=0.
end if
- case(1)
- poids(1)=cd(xx)
- poids(2)=ce(xx)
- poids(3)=0.
- poids(4)=0.
- case(2)
- poids(2)=cd(xx-1.)+ce(xx-1.)
- poids(3)=0.
- poids(4)=0.
- case(3)
- poids(1)=cd(xx)+ce(xx)
- poids(2)=0.
- poids(3)=0.
- poids(4)=0.
- case(4)
- poids(2)=ce(xx)
- poids(3)=0.
- poids(4)=0.
- case(5)
- poids(1)=ce(xx+1.)
- poids(2)=0.
- poids(3)=0.
- poids(4)=0.
- case(6)
- poids(2)=cd(xx-1.)
- poids(3)=ce(xx-1.)-ce(xx-2.)
- poids(4)=ce(xx-2.)
- case(7)
- poids(1)=cd(xx)
- poids(2)=ce(xx)-ce(xx-1.)
- poids(3)=ce(xx-1.)
- poids(4)=0.
- case(8)
- poids(1)=cd(xx)-cd(xx-1.)+ce(xx)-ce(xx-1.)
- poids(2)=cd(xx-1.)
- poids(3)=ce(xx-1.)
- poids(4)=0.
-
- end select
-
+
!remaillage à l' interrieur domaine
!---------------------------------
@@ -3980,12 +5100,12 @@ contains
end subroutine calcul_phi_minmax
- end subroutine remaill_l4_bloc_v2_limit_x
+ end subroutine remaill_l4_limit_x
- subroutine remaill_l4_bloc_v2_limit_y (donne,posx,posy,remaille)
+ subroutine remaill_l4_limit_y (donne,posx,posy,remaille)
implicit none
real(kind=8),dimension(:),intent(in) :: donne
real(kind=8),dimension(:),intent(in) :: posx,posy
@@ -4035,130 +5155,28 @@ contains
phi_mdemi=phi(0)
phi_pdemi=phi(1)
phi_p3demi=phi(2)
- select case (blocg(i))
-
- case(0)
+
if (xx<=0.5) then
poids(-3)=0.
poids(-2)=ca(xx)
poids(-1)=cb(xx)
- poids(0)=cc(xx)
+ poids(0)=cc(xx)
+ poids(1)=cd(xx)
+ poids(2)=ce(xx)
+ poids(3)=0.
+ poids(4)=0.
else
poids(-3)=0.
poids(-2)=0.
poids(-1)=ca(xx-1.)
poids(0)=cb(xx-1.)
poids(1)=cc(xx-1)
- end if
- case(1)
- poids(-3)=0.
- poids(-2)=ca(xx)
- poids(-1)=cb(xx)
- poids(0)=cc(xx)
- case(2)
- poids(-3)=0.
- poids(-2)=0.
- poids(-1)=ca(xx-1.)
- poids(0)=cb(xx-1.)
- poids(1)=cc(xx-1.)+cd(xx-1.)+ce(xx-1.)-ce(xx)
- case(3)
- poids(-3)=0.
- poids(-2)=ca(xx)
- poids(-1)=cb(xx)
- poids(0)=cc(xx)+cd(xx)+ce(xx)-ce(xx+1.)
- case(4)
- poids(-3)=0.
- poids(-2)=0.
- poids(-1)=ca(xx-1.)
- poids(0)=ca(xx)-ca(xx-1.)+cb(xx)+cc(xx)
- case(5)
- poids(-3)=0.
- poids(-2)=0.
- poids(-1)=0.
- poids(0)=ca(xx-1.)+cb(xx-1.)
- poids(1)=cc(xx-1.)
- case(6)
- poids(-3)=0.
- poids(-2)=0.
- poids(-1)=ca(xx)+cb(xx)
- poids(0)=cc(xx)
- case(7)
- poids(-3)=0.
- poids(-2)=ca(xx)
- poids(-1)=cb(xx)
- poids(0)=ca(xx-1.)-ca(xx)+cb(xx-1.)-cb(xx)
- poids(1)=cc(xx-1)
- case(8)
- poids(-3)=ca(xx+1.)
- poids(-2)=cb(xx+1.)
- poids(-1)=ca(xx)-ca(xx+1.)+cb(xx)-cb(xx+1.)
- poids(0)=cc(xx)
- case(9)
- poids(-3)=0.
- poids(-2)=ca(xx)
- poids(-1)=ca(xx-1.)-ca(xx)
- poids(0)=cb(xx-1.)
- poids(1)=cc(xx-1)
- case(10)
- poids(-3)=ca(xx+1.)
- poids(-2)=ca(xx)-ca(xx+1.)
- poids(-1)=cb(xx)
- poids(0)=cc(xx)
-
- end select
-
- select case (blocd(i))
- case(0)
- if (xx<=0.5) then
- poids(1)=cd(xx)
- poids(2)=ce(xx)
- poids(3)=0.
- poids(4)=0.
- else
poids(2)=cd(xx-1.)
poids(3)=ce(xx-1.)
poids(4)=0.
end if
- case(1)
- poids(1)=cd(xx)
- poids(2)=ce(xx)
- poids(3)=0.
- poids(4)=0.
- case(2)
- poids(2)=cd(xx-1.)+ce(xx-1.)
- poids(3)=0.
- poids(4)=0.
- case(3)
- poids(1)=cd(xx)+ce(xx)
- poids(2)=0.
- poids(3)=0.
- poids(4)=0.
- case(4)
- poids(2)=ce(xx)
- poids(3)=0.
- poids(4)=0.
- case(5)
- poids(1)=ce(xx+1.)
- poids(2)=0.
- poids(3)=0.
- poids(4)=0.
- case(6)
- poids(2)=cd(xx-1.)
- poids(3)=ce(xx-1.)-ce(xx-2.)
- poids(4)=ce(xx-2.)
- case(7)
- poids(1)=cd(xx)
- poids(2)=ce(xx)-ce(xx-1.)
- poids(3)=ce(xx-1.)
- poids(4)=0.
- case(8)
- poids(1)=cd(xx)-cd(xx-1.)+ce(xx)-ce(xx-1.)
- poids(2)=cd(xx-1.)
- poids(3)=ce(xx-1.)
- poids(4)=0.
-
-
- end select
+
+
!remaillage à l' interrieur domaine
!---------------------------------
do c=-3,4
@@ -4337,7 +5355,7 @@ contains
end subroutine calcul_phi_minmax
- end subroutine remaill_l4_bloc_v2_limit_y
+ end subroutine remaill_l4_limit_y
diff --git a/Remesh/split_3d_rapide/Makefile b/Remesh/split_3d_rapide/Makefile
index fd7983ca53ebb5137028df1dd4605202fe892e45..17d0e77ef14580711262137f79ab141d79e4f9e3 100644
--- a/Remesh/split_3d_rapide/Makefile
+++ b/Remesh/split_3d_rapide/Makefile
@@ -2,7 +2,7 @@
# Makefile
#
#Compilateur = ifort
-Compilateur = ifort -O0 -pg -debug
+Compilateur = ifort -pg
#Compilateur = ifort -O0 -debug
Programme = split_3d
FILES = donnees_mod.f90 tab_mod.f90 init_mod.f90 remaillage_mod.f90 interpolation_mod.f90 advection_mod.f90 resultats_mod.f90 utile_mod.f90 main.f90
diff --git a/Remesh/split_3d_rapide/advection_mod.f90 b/Remesh/split_3d_rapide/advection_mod.f90
index 714093445b3b3e478ccf1a8d5dd24bbe9d6ce6f9..94667e3ffab99f119114099e266a13982d723233 100644
--- a/Remesh/split_3d_rapide/advection_mod.f90
+++ b/Remesh/split_3d_rapide/advection_mod.f90
@@ -18,30 +18,30 @@ contains
!-------------------
!deformation sphere:
!--------------------
- do k=1,nz_gro
- z=zd+(k-1)*dz_gro
- do j=1,ny_gro
- y=yb+(j-1)*dy_gro
- do i=1,nx_gro
- x=xg+(i-1)*dx_gro
-
- !t3=0.3
- !t3=3.
- !t3=1.
-
- if (tps<=1.) then
- vitx(i,j,k)=2.*sin(pi*x)**2*sin(2.*pi*y)*sin(2.*pi*z) !*cos(pi*tps/t3)
- vity(i,j,k)=-sin(2.*pi*x)*sin(pi*y)**2*sin(2.*pi*z) !*cos(pi*tps/t3)
- vitz(i,j,k)=-sin(2.*pi*x)*sin(2.*pi*y)*sin(pi*z)**2 !*cos(pi*tps/t3)
- else
- vitx(i,j,k)=-2.*sin(pi*x)**2*sin(2.*pi*y)*sin(2.*pi*z)
- vity(i,j,k)=sin(2.*pi*x)*sin(pi*y)**2*sin(2.*pi*z)
- vitz(i,j,k)=sin(2.*pi*x)*sin(2.*pi*y)*sin(pi*z)**2
- end if
-
- end do
- end do
- end do
+!!$ do k=1,nz_gro
+!!$ z=zd+(k-1)*dz_gro
+!!$ do j=1,ny_gro
+!!$ y=yb+(j-1)*dy_gro
+!!$ do i=1,nx_gro
+!!$ x=xg+(i-1)*dx_gro
+!!$
+!!$ !t3=0.3
+!!$ !t3=3.
+!!$ !t3=1.
+!!$
+!!$ if (tps<=1.) then
+!!$ vitx(i,j,k)=2.*sin(pi*x)**2*sin(2.*pi*y)*sin(2.*pi*z) !*cos(pi*tps/t3)
+!!$ vity(i,j,k)=-sin(2.*pi*x)*sin(pi*y)**2*sin(2.*pi*z) !*cos(pi*tps/t3)
+!!$ vitz(i,j,k)=-sin(2.*pi*x)*sin(2.*pi*y)*sin(pi*z)**2 !*cos(pi*tps/t3)
+!!$ else
+!!$ vitx(i,j,k)=-2.*sin(pi*x)**2*sin(2.*pi*y)*sin(2.*pi*z)
+!!$ vity(i,j,k)=sin(2.*pi*x)*sin(pi*y)**2*sin(2.*pi*z)
+!!$ vitz(i,j,k)=sin(2.*pi*x)*sin(2.*pi*y)*sin(pi*z)**2
+!!$ end if
+!!$
+!!$ end do
+!!$ end do
+!!$ end do
!----------------------------------
@@ -99,22 +99,21 @@ contains
!----------------------------------
-!!$ open(20,file='datavelx_adr',form='unformatted',convert='big_endian',status='unknown')
-!!$ read(20) (((tvitx(i,j,k),i=1,nx_gro),j=1,nx_gro),k=1,nx_gro)
-!!$ close(20)
-!!$
-!!$ open(20,file='datavely_adr',form='unformatted',convert='big_endian',status='unknown')
-!!$ read(20) (((tvity(i,j,k),i=1,nx_gro),j=1,nx_gro),k=1,nx_gro)
-!!$ close(20)
-!!$
-!!$ open(20,file='datavelz_adr',form='unformatted',convert='big_endian',status='unknown')
-!!$ read(20) (((tvitz(i,j,k),i=1,nx_gro),j=1,nx_gro),k=1,nx_gro)
-!!$ close(20)
-!!$
-!!$ vitx=tvitx
-!!$ vity=tvity
-!!$ vitz=tvitz
+ open(20,file='datavelx',form='unformatted',convert='big_endian',status='unknown')
+ read(20) (((tvitx(i,j,k),i=1,nx_gro),j=1,nx_gro),k=1,nx_gro)
+ close(20)
+
+ open(20,file='datavely',form='unformatted',convert='big_endian',status='unknown')
+ read(20) (((tvity(i,j,k),i=1,nx_gro),j=1,nx_gro),k=1,nx_gro)
+ close(20)
+
+ open(20,file='datavelz',form='unformatted',convert='big_endian',status='unknown')
+ read(20) (((tvitz(i,j,k),i=1,nx_gro),j=1,nx_gro),k=1,nx_gro)
+ close(20)
+ vitx=tvitx
+ vity=tvity
+ vitz=tvitz
@@ -152,7 +151,7 @@ contains
integer :: i,j,k
real(kind=8) :: x,y,z,m
- m=cutoff*maxval(abs(qg))
+ m=cutoff!*maxval(abs(qg))
npart=0
numpg=0
do k=1,nz
@@ -187,7 +186,7 @@ contains
integer :: i,j,k
real(kind=8) :: x,y,z,m
- m=cutoff*maxval(abs(qg))
+ m=cutoff!*maxval(abs(qg))
npart=0
numpg=0
do k=1,nz
@@ -222,7 +221,7 @@ contains
integer :: i,j,k
real(kind=8) :: x,y,z,m
- m=cutoff*maxval(abs(qg))
+ m=cutoff!*maxval(abs(qg))
npart=0
numpg=0
diff --git a/Remesh/split_3d_rapide/init_mod.f90 b/Remesh/split_3d_rapide/init_mod.f90
index 9d652be312e4a366b0f588736c634ab3ddb98ab1..7eb85c85d9463ec656704362a5bace45bbdb12f4 100644
--- a/Remesh/split_3d_rapide/init_mod.f90
+++ b/Remesh/split_3d_rapide/init_mod.f90
@@ -13,20 +13,20 @@ contains
!------------
!sphere
!------------
- do k=1,nz
- z=ztab(k)
- do j=1,ny
- y=ytab(j)
- do i=1,nx
- x=xtab(i)
-
- r=sqrt((x-0.35)**2+(y-0.35)**2+(z-0.35)**2)
- qg(i,j,k)=0.
- if (r<=0.15) qg(i,j,k)=1.
-
- end do
- end do
- end do
+!!$ do k=1,nz
+!!$ z=ztab(k)
+!!$ do j=1,ny
+!!$ y=ytab(j)
+!!$ do i=1,nx
+!!$ x=xtab(i)
+!!$
+!!$ r=sqrt((x-0.35)**2+(y-0.35)**2+(z-0.35)**2)
+!!$ qg(i,j,k)=0.
+!!$ if (r<=0.15) qg(i,j,k)=1.
+!!$
+!!$ end do
+!!$ end do
+!!$ end do
!!$ !---gaussienne sur [0,1]^3-------------------------
@@ -68,19 +68,19 @@ contains
!------------
!test
-!!$ !------------
-!!$ do k=1,nz
-!!$ z=ztab(k)
-!!$ do j=1,ny
-!!$ y=ytab(j)
-!!$ do i=1,nx
-!!$ x=xtab(i)
-!!$
-!!$ qg(i,j,k)=1.
-!!$
-!!$ end do
-!!$ end do
-!!$ end do
+ !------------
+ do k=1,nz
+ z=ztab(k)
+ do j=1,ny
+ y=ytab(j)
+ do i=1,nx
+ x=xtab(i)
+
+ qg(i,j,k)=1.
+
+ end do
+ end do
+ end do
end subroutine init_grille
diff --git a/Remesh/split_3d_rapide/main.f90 b/Remesh/split_3d_rapide/main.f90
index 4ff00918f63e8f388bdd3019c33e6694f823d613..87be5adbead22acd19b7dd1a4993f975fca8b804 100644
--- a/Remesh/split_3d_rapide/main.f90
+++ b/Remesh/split_3d_rapide/main.f90
@@ -36,9 +36,12 @@ program split
read(10,*) nom_fich_tps_cpu
close(10)
+
+ name_err=trim(name_err)
+
call cpu_time(t1)
- erreur:do nx_boucle=100,100
+ erreur:do nx_boucle=129,129
open(unit=10,file="parameter",form="formatted")
read(10,*) xg,xd,yb,yh,zd,zf
@@ -64,14 +67,15 @@ program split
maxv=2.
!maxv=3.3
!maxv=1.
- !cfl=0.4
- !cfl=7.
+ cfl=0.4
+ !cfl=5.
!cfl=14.
- cfl=4.
+ !cfl=4.
! cfl=8.
!cfl=2.*0.025252525252525/dx! cfl5
!cfl=2.*0.02424242520573163/dx !cfl 4.8
+
!print*,"cfl=",cfl
@@ -79,6 +83,8 @@ program split
dt=cfl*dx/maxv
! dt=2./150.
+ dt=0.049
+
print*,"dt",dt
!
@@ -120,6 +126,8 @@ program split
call make_grille
call init_grille()
+ cutoff=cutoff*maxval(abs(qg))
+
!=============================
!champ de vitesse
@@ -134,7 +142,7 @@ program split
!print*,"cflmax:",dt*max(maxval(vxg),maxval(vyg),maxval(vzg))/dx
- if ((time+dt)>tfin) dt=tfin-time
+ !if ((time+dt)>tfin) dt=tfin-time
if (time<=2.*dt) dt_deb=dt
@@ -148,7 +156,7 @@ program split
!================================
!si remaillage tenso en espace
!================================
-!!$
+
!!$ !RK2
!!$ !allocate(vxg1(1:nx_gro,1:ny_gro,1:nz_gro),vyg1(1:nx_gro,1:ny_gro,1:nz_gro),vzg1(1:nx_gro,1:ny_gro,1:nz_gro))
!!$ !call def_v_advec (vxg1,vyg1,vzg1,time+0.5*dt)
@@ -183,11 +191,11 @@ program split
!!$ !deallocate(vxg1,vyg1,vzg1)
!!$ !deallocate(vxg2,vyg2,vzg2)
!!$
-!!$
-!!$ !=====================
-!!$ !splitting en espace
-!!$ !=====================
-
+
+ !=====================
+ !splitting en espace
+ !=====================
+!!$
!RK2
!allocate(vxg1(1:nx_gro,1:ny_gro,1:nz_gro),vyg1(1:nx_gro,1:ny_gro,1:nz_gro),vzg1(1:nx_gro,1:ny_gro,1:nz_gro))
!call def_v_advec (vxg1,vyg1,vzg1,time +0.5*dt)
@@ -214,7 +222,7 @@ program split
!update v
!--------
- call update_vx_2
+ !call update_vx_2
!call update_vx_3
!print*,"vmax interpole rk2",max(maxval(abs(vx)),maxval(abs(vy)),maxval(abs(vz)))
@@ -265,13 +273,14 @@ program split
!remaillage
!----------
- !call remaill_l2_x(qp,xp,yp,zp,qg)
- !call remaill_l4_x(qp,xp,yp,zp,qg)
+ if (type_b==200) call remaill_l2_x(qp,xp,yp,zp,qg)
+ if (type_b==400) call remaill_l4_x(qp,xp,yp,zp,qg)
!call remaill_l5_x(qp,xp,yp,zp,qg)
- !call remaill_l6_x(qp,xp,yp,zp,qg)
+ if (type_b==600) call remaill_l6_x(qp,xp,yp,zp,qg)
!call remaill_l8_x(qp,xp,yp,zp,qg)
if (type_b==2) call remaill_l2_bloc_x_v2(qp,xp,yp,zp,qg)
if (type_b==4) call remaill_l4_bloc_x_v2(qp,xp,yp,zp,qg)
+ if (type_b==20) call limit_l2_bloc_x(qp,xp,yp,zp,qg)
!contrainte cfl
!--------------
@@ -290,11 +299,12 @@ program split
!y
!--
call crea_part_y()
+ !call crea_part_x()
!update v
!--------
- call update_vy_2
+ !call update_vy_2
!call update_vy_3
!bloc y
@@ -309,13 +319,14 @@ program split
call ad_euler_y
!remaillage
!----------
- !call remaill_l2_y(qp,xp,yp,zp,qg)
- !call remaill_l4_y(qp,xp,yp,zp,qg)
+ if (type_b==200) call remaill_l2_y(qp,xp,yp,zp,qg)
+ if (type_b==400) call remaill_l4_y(qp,xp,yp,zp,qg)
!call remaill_l5_y(qp,xp,yp,zp,qg)
- !call remaill_l6_y(qp,xp,yp,zp,qg)
+ if (type_b==600) call remaill_l6_y(qp,xp,yp,zp,qg)
!call remaill_l8_y(qp,xp,yp,zp,qg)
if (type_b==2)call remaill_l2_bloc_y_v2(qp,xp,yp,zp,qg)
- if (type_b==4)call remaill_l4_bloc_y_v2(qp,xp,yp,zp,qg)
+ if (type_b==4)call remaill_l4_bloc_y_v2(qp,xp,yp,zp,qg)
+ if (type_b==20) call limit_l2_bloc_y(qp,xp,yp,zp,qg)
!contrainte cfl
!--------------
@@ -332,10 +343,11 @@ program split
!z
!--
call crea_part_z()
+ !call crea_part_x()
!update v
!--------
- call update_vz_2
+ !call update_vz_2
!call update_vz_3
!call def_v_advec (vxg1,vyg1,vzg1,time+dt/3.)
@@ -356,13 +368,14 @@ program split
!remaillage
!----------
- !call remaill_l2_z(qp,xp,yp,zp,qg)
- !call remaill_l4_z(qp,xp,yp,zp,qg)
+ if (type_b==200) call remaill_l2_z(qp,xp,yp,zp,qg)
+ if (type_b==400) call remaill_l4_z(qp,xp,yp,zp,qg)
!call remaill_l5_z(qp,xp,yp,zp,qg)
- !call remaill_l6_z(qp,xp,yp,zp,qg)
+ if (type_b==600) call remaill_l6_z(qp,xp,yp,zp,qg)
!call remaill_l8_z(qp,xp,yp,zp,qg)
if (type_b==2)call remaill_l2_bloc_z_v2(qp,xp,yp,zp,qg)
- if (type_b==4)call remaill_l4_bloc_z_v2(qp,xp,yp,zp,qg)
+ if (type_b==4)call remaill_l4_bloc_z_v2(qp,xp,yp,zp,qg)
+ if (type_b==20) call limit_l2_bloc_z(qp,xp,yp,zp,qg)
!contrainte cfl
!--------------
@@ -427,10 +440,10 @@ program split
!!$ !if (type_b==2)call remaill_l2_bloc_x_v2(qp,xp,yp,zp,qg)
!!$ !if (type_b==4)call remaill_l4_bloc_x_v2(qp,xp,yp,zp,qg)
!!$ !call remaill_l2_x(qp,xp,yp,zp,qg)
-!!$ !call remaill_l4_x(qp,xp,yp,zp,qg)
+!!$ call remaill_l4_x(qp,xp,yp,zp,qg)
!!$ !call remaill_l5_x(qp,xp,yp,zp,qg)
!!$ !call remaill_l6_x(qp,xp,yp,zp,qg)
-!!$ call remaill_l8_x(qp,xp,yp,zp,qg)
+!!$ !call remaill_l8_x(qp,xp,yp,zp,qg)
!!$
!!$ !y
!!$ !--
@@ -455,10 +468,10 @@ program split
!!$ !if (type_b==2)call remaill_l2_bloc_y_v2(qp,xp,yp,zp,qg)
!!$ !if (type_b==4)call remaill_l4_bloc_y_v2(qp,xp,yp,zp,qg)
!!$ !call remaill_l2_y(qp,xp,yp,zp,qg)
-!!$ !call remaill_l4_y(qp,xp,yp,zp,qg)
+!!$ call remaill_l4_y(qp,xp,yp,zp,qg)
!!$ !call remaill_l5_y(qp,xp,yp,zp,qg)
!!$ !call remaill_l6_y(qp,xp,yp,zp,qg)
-!!$ call remaill_l8_y(qp,xp,yp,zp,qg)
+!!$ !call remaill_l8_y(qp,xp,yp,zp,qg)
!!$
!!$
!!$ !z
@@ -487,10 +500,10 @@ program split
!!$ !if (type_b==2)call remaill_l2_bloc_z_v2(qp,xp,yp,zp,qg)
!!$ !if (type_b==4)call remaill_l4_bloc_z_v2(qp,xp,yp,zp,qg)
!!$ !call remaill_l2_z(qp,xp,yp,zp,qg)
-!!$ !call remaill_l4_z(qp,xp,yp,zp,qg)
+!!$ call remaill_l4_z(qp,xp,yp,zp,qg)
!!$ !call remaill_l5_z(qp,xp,yp,zp,qg)
!!$ !call remaill_l6_z(qp,xp,yp,zp,qg)
-!!$ call remaill_l8_z(qp,xp,yp,zp,qg)
+!!$ !call remaill_l8_z(qp,xp,yp,zp,qg)
!!$
!!$
!!$ !y
@@ -517,10 +530,10 @@ program split
!!$ !if (type_b==2)call remaill_l2_bloc_y_v2(qp,xp,yp,zp,qg)
!!$ !if (type_b==4)call remaill_l4_bloc_y_v2(qp,xp,yp,zp,qg)
!!$ !call remaill_l2_y(qp,xp,yp,zp,qg)
-!!$ !call remaill_l4_y(qp,xp,yp,zp,qg)
+!!$ call remaill_l4_y(qp,xp,yp,zp,qg)
!!$ !call remaill_l5_y(qp,xp,yp,zp,qg)
!!$ !call remaill_l6_y(qp,xp,yp,zp,qg)
-!!$ call remaill_l8_y(qp,xp,yp,zp,qg)
+!!$ !call remaill_l8_y(qp,xp,yp,zp,qg)
!!$
!!$
!!$ !x
@@ -547,10 +560,10 @@ program split
!!$ !if (type_b==2)call remaill_l2_bloc_x_v2(qp,xp,yp,zp,qg)
!!$ !if (type_b==4)call remaill_l4_bloc_x_v2(qp,xp,yp,zp,qg)
!!$ !call remaill_l2_x(qp,xp,yp,zp,qg)
-!!$ !call remaill_l4_x(qp,xp,yp,zp,qg)
+!!$ call remaill_l4_x(qp,xp,yp,zp,qg)
!!$ !call remaill_l5_x(qp,xp,yp,zp,qg)
!!$ !call remaill_l6_x(qp,xp,yp,zp,qg)
-!!$ call remaill_l8_x(qp,xp,yp,zp,qg)
+!!$ !call remaill_l8_x(qp,xp,yp,zp,qg)
!!$
!!$
!!$ dt=dt_sauv
@@ -626,7 +639,7 @@ program split
-
+ ! write(29,'(3(e18.10,2x))') time,minval(qg),maxval(qg)
@@ -675,6 +688,7 @@ program split
close(66)
close(67)
close(10)
+ close(29)
deallocate (xp,qp,vx,vy,yp,zp,vz)
deallocate (qg,vxg,vyg,vzg)
diff --git a/Remesh/split_3d_rapide/parameter b/Remesh/split_3d_rapide/parameter
index 5cd38c67cf01ddf1a6549bb649a57360bf4692c4..07a0d9f7f905c1d9d7de64e1859611921dc1ae1d 100644
--- a/Remesh/split_3d_rapide/parameter
+++ b/Remesh/split_3d_rapide/parameter
@@ -1,9 +1,9 @@
0. 1. 0. 1. 0. 1. ! xg xd yb yh zd zf
-1. ! TFin
-0.0001 0.1 !cutoff
-2 ! type bloc
-1 !longeur bloc
-"test.vtk" ! sortie fin vtk
+0.0001 ! TFin
+0. 0.0001 0.0001 0.1 !cutoff
+4 ! type bloc
+3 !longeur bloc
+"vtk/test_o1_long_b_3.vtk" ! sortie fin vtk
"RES/test.er" ! fichier erreur
-"/users/these/magni/forge_limit/code-these/split_3d_rapide/RES/" "/home/magni/split_3d_rapide/" "RES/vtk/" !dossier
+"/users/these/magni/parmes/trunk/Remesh/split_3d_rapide/RES/" "/home/magni/split_3d_rapide/" "RES/vtk/" !dossier
"test.tps" ! fichier tps cpu
\ No newline at end of file
diff --git a/Remesh/split_3d_rapide/remaillage_mod.f90 b/Remesh/split_3d_rapide/remaillage_mod.f90
index 79ceec71c323f02803a2a45fba2d6ce26c7670b3..41f2f0647170e147fa20034feb93726b3eb3af21 100644
--- a/Remesh/split_3d_rapide/remaillage_mod.f90
+++ b/Remesh/split_3d_rapide/remaillage_mod.f90
@@ -1244,29 +1244,43 @@ contains
!calcul des poids
!----------------
- x2=xx**2
- x3=xx**3
- x4=xx**4
-
- if (xx<=0.5) then
- poidx(-2)=(2.*xx-x2-2*x3+x4)/24.
- poidx(-1)=(-4.*xx+4.*x2+x3-x4)/6.
- poidx(0)=1.+(-5.*x2+x4)/4.
- poidx(1)=(4.*xx+4.*x2-x3-x4)/6.
- poidx(2)=(-2.*xx-x2+2*x3+x4)/24.
+!!$ x2=xx**2
+!!$ x3=xx**3
+!!$ x4=xx**4
+!!$
+!!$ if (xx<=0.5) then
+!!$ poidx(-2)=(2.*xx-x2-2*x3+x4)/24.
+!!$ poidx(-1)=(-4.*xx+4.*x2+x3-x4)/6.
+!!$ poidx(0)=1.+(-5.*x2+x4)/4.
+!!$ poidx(1)=(4.*xx+4.*x2-x3-x4)/6.
+!!$ poidx(2)=(-2.*xx-x2+2*x3+x4)/24.
+!!$ poidx(3)=0.
+!!$ else
+!!$ poidx(-2)=0.
+!!$ poidx(-1)=(-6.*xx+11.*x2-6.*x3+x4)/24.
+!!$ poidx(0)=1.+(-5.*xx-5.*x2+5.*x3-x4)/6.
+!!$ poidx(1)=(6.*xx+x2-4.*x3+x4)/4.
+!!$ poidx(2)=(-3.*xx+x2+3.*x3-x4)/6.
+!!$ poidx(3)=(2.*xx-x2-2.*x3+x4)/24.
+!!$ end if
+
+ if (xx<=0.5) then
+ poidx(-2)=xx*(xx-1.)*(xx-2.)*(xx+1.)/24. !a
+ poidx(-1)=-xx*(xx-1.)*(xx**2-4.)/6. !b
+ poidx(0)=(xx**2-4.)*(xx**2-1.)/4. !c
+ poidx(1)=-xx*(xx+1.)*(xx+2.)*(xx-2.)/6. !d
+ poidx(2)=xx*(xx+1.)*(xx+2.)*(xx-1.)/24. !e
poidx(3)=0.
- else
- poidx(-2)=0.
- poidx(-1)=(-6.*xx+11.*x2-6.*x3+x4)/24.
- poidx(0)=1.+(-5.*xx-5.*x2+5.*x3-x4)/6.
- poidx(1)=(6.*xx+x2-4.*x3+x4)/4.
- poidx(2)=(-3.*xx+x2+3.*x3-x4)/6.
- poidx(3)=(2.*xx-x2-2.*x3+x4)/24.
+ else
+ poidx(-2)=0.
+ poidx(-1)=(xx-1.)*(xx-2.)*(xx-3.)*xx/24. !a'
+ poidx(0)=-(xx-1.)*(xx-2.)*(xx-3.)*(xx+1.)/6. !b'
+ poidx(1)=xx*(xx+1.)*(xx-3.)*(xx-2.)/4. !c'
+ poidx(2)=-xx*(xx+1.)*(xx-3.)*(xx-1.)/6. !d'
+ poidx(3)=xx*(xx-2.)*(xx**2-1.)/24. !e'
end if
-
-
!remaillage à l' interrieur domaine
!---------------------------------
do c=-2,3
@@ -1333,30 +1347,45 @@ contains
!calcul des poids
!----------------
- x2=xx**2
- x3=xx**3
- x4=xx**4
-
- if (xx<=0.5) then
- poidx(-2)=(2.*xx-x2-2*x3+x4)/24.
- poidx(-1)=(-4.*xx+4.*x2+x3-x4)/6.
- poidx(0)=1.+(-5.*x2+x4)/4.
- poidx(1)=(4.*xx+4.*x2-x3-x4)/6.
- poidx(2)=(-2.*xx-x2+2*x3+x4)/24.
+!!$ x2=xx**2
+!!$ x3=xx**3
+!!$ x4=xx**4
+!!$
+!!$ if (xx<=0.5) then
+!!$ poidx(-2)=(2.*xx-x2-2*x3+x4)/24.
+!!$ poidx(-1)=(-4.*xx+4.*x2+x3-x4)/6.
+!!$ poidx(0)=1.+(-5.*x2+x4)/4.
+!!$ poidx(1)=(4.*xx+4.*x2-x3-x4)/6.
+!!$ poidx(2)=(-2.*xx-x2+2*x3+x4)/24.
+!!$ poidx(3)=0.
+!!$ else
+!!$ poidx(-2)=0.
+!!$ poidx(-1)=(-6.*xx+11.*x2-6.*x3+x4)/24.
+!!$ poidx(0)=1.+(-5.*xx-5.*x2+5.*x3-x4)/6.
+!!$ poidx(1)=(6.*xx+x2-4.*x3+x4)/4.
+!!$ poidx(2)=(-3.*xx+x2+3.*x3-x4)/6.
+!!$ poidx(3)=(2.*xx-x2-2.*x3+x4)/24.
+!!$ end if
+
+ if (xx<=0.5) then
+ poidx(-2)=xx*(xx-1.)*(xx-2.)*(xx+1.)/24. !a
+ poidx(-1)=-xx*(xx-1.)*(xx**2-4.)/6. !b
+ poidx(0)=(xx**2-4.)*(xx**2-1.)/4. !c
+ poidx(1)=-xx*(xx+1.)*(xx+2.)*(xx-2.)/6. !d
+ poidx(2)=xx*(xx+1.)*(xx+2.)*(xx-1.)/24. !e
poidx(3)=0.
- else
- poidx(-2)=0.
- poidx(-1)=(-6.*xx+11.*x2-6.*x3+x4)/24.
- poidx(0)=1.+(-5.*xx-5.*x2+5.*x3-x4)/6.
- poidx(1)=(6.*xx+x2-4.*x3+x4)/4.
- poidx(2)=(-3.*xx+x2+3.*x3-x4)/6.
- poidx(3)=(2.*xx-x2-2.*x3+x4)/24.
+ else
+ poidx(-2)=0.
+ poidx(-1)=(xx-1.)*(xx-2.)*(xx-3.)*xx/24. !a'
+ poidx(0)=-(xx-1.)*(xx-2.)*(xx-3.)*(xx+1.)/6. !b'
+ poidx(1)=xx*(xx+1.)*(xx-3.)*(xx-2.)/4. !c'
+ poidx(2)=-xx*(xx+1.)*(xx-3.)*(xx-1.)/6. !d'
+ poidx(3)=xx*(xx-2.)*(xx**2-1.)/24. !e'
end if
-
!remaillage à l' interrieur domaine
!---------------------------------
@@ -1418,30 +1447,45 @@ contains
!calcul des poids
!----------------
- x2=xx**2
- x3=xx**3
- x4=xx**4
-
- if (xx<=0.5) then
- poidx(-2)=(2.*xx-x2-2*x3+x4)/24.
- poidx(-1)=(-4.*xx+4.*x2+x3-x4)/6.
- poidx(0)=1.+(-5.*x2+x4)/4.
- poidx(1)=(4.*xx+4.*x2-x3-x4)/6.
- poidx(2)=(-2.*xx-x2+2*x3+x4)/24.
+!!$ x2=xx**2
+!!$ x3=xx**3
+!!$ x4=xx**4
+!!$
+!!$ if (xx<=0.5) then
+!!$ poidx(-2)=(2.*xx-x2-2*x3+x4)/24.
+!!$ poidx(-1)=(-4.*xx+4.*x2+x3-x4)/6.
+!!$ poidx(0)=1.+(-5.*x2+x4)/4.
+!!$ poidx(1)=(4.*xx+4.*x2-x3-x4)/6.
+!!$ poidx(2)=(-2.*xx-x2+2*x3+x4)/24.
+!!$ poidx(3)=0.
+!!$ else
+!!$ poidx(-2)=0.
+!!$ poidx(-1)=(-6.*xx+11.*x2-6.*x3+x4)/24.
+!!$ poidx(0)=1.+(-5.*xx-5.*x2+5.*x3-x4)/6.
+!!$ poidx(1)=(6.*xx+x2-4.*x3+x4)/4.
+!!$ poidx(2)=(-3.*xx+x2+3.*x3-x4)/6.
+!!$ poidx(3)=(2.*xx-x2-2.*x3+x4)/24.
+!!$ end if
+
+ if (xx<=0.5) then
+ poidx(-2)=xx*(xx-1.)*(xx-2.)*(xx+1.)/24. !a
+ poidx(-1)=-xx*(xx-1.)*(xx**2-4.)/6. !b
+ poidx(0)=(xx**2-4.)*(xx**2-1.)/4. !c
+ poidx(1)=-xx*(xx+1.)*(xx+2.)*(xx-2.)/6. !d
+ poidx(2)=xx*(xx+1.)*(xx+2.)*(xx-1.)/24. !e
poidx(3)=0.
- else
- poidx(-2)=0.
- poidx(-1)=(-6.*xx+11.*x2-6.*x3+x4)/24.
- poidx(0)=1.+(-5.*xx-5.*x2+5.*x3-x4)/6.
- poidx(1)=(6.*xx+x2-4.*x3+x4)/4.
- poidx(2)=(-3.*xx+x2+3.*x3-x4)/6.
- poidx(3)=(2.*xx-x2-2.*x3+x4)/24.
+ else
+ poidx(-2)=0.
+ poidx(-1)=(xx-1.)*(xx-2.)*(xx-3.)*xx/24. !a'
+ poidx(0)=-(xx-1.)*(xx-2.)*(xx-3.)*(xx+1.)/6. !b'
+ poidx(1)=xx*(xx+1.)*(xx-3.)*(xx-2.)/4. !c'
+ poidx(2)=-xx*(xx+1.)*(xx-3.)*(xx-1.)/6. !d'
+ poidx(3)=xx*(xx-2.)*(xx**2-1.)/24. !e'
end if
-
!remaillage à l' interrieur domaine
!---------------------------------
@@ -3482,9 +3526,7 @@ subroutine remaill_l8_z (donne,posx,posy,posz,remaille)
!calcul des poids
!----------------
-
- select case (blocg(i))
-
+ select case (blocg(i))
case(0)
if (xx<=0.5) then
poids(-3)=0.
@@ -3506,57 +3548,52 @@ subroutine remaill_l8_z (donne,posx,posy,posz,remaille)
case(2)
poids(-3)=0.
poids(-2)=0.
- poids(-1)=(xx-1.)*(xx-2.)*(xx-3.)*xx/24. !a'
- poids(0)=-(xx-1.)*(xx-2.)*(xx-3.)*(xx+1.)/6. !b'
- poids(1)=xx*(xx+1.)*(xx-3.)*(xx-2.)/4. !c'
- case(3)
- poids(-3)=(xx+2.)*xx*(xx**2-1.)/24.
- poids(-2)=-xx*(xx+3.)*(xx**2-1.)/6.
- poids(-1)=0.5*xx**3+0.5*xx**2-xx
- case(4)
- poids(-3)=0.
- poids(-2)=xx*(xx-2.)*(xx**2-1.)/24.
- poids(-1)=-xx*(xx-1.)*(xx**2-4.)/6.
- poids(0)=0.5*xx**3-xx**2-0.5*xx+1.
- case(5)
- poids(-3)=(xx+2.)*xx*(xx**2-1.)/24.
- poids(-2)=xx*(1.-xx**2)/6.
- poids(-1)=-xx*(xx-1.)*(xx**2-4.)/6.
- poids(0)=(xx**2-4.)*(xx**2-1.)/4.
- case(6)
- poids(-3)=0.
- poids(-2)=xx*(xx-2.)*(xx**2-1.)/24.
- poids(-1)=-xx**3/6.+0.5*xx**2-xx/3.
+ poids(-1)=(xx-1.)*(xx-2.)*(xx-3.)*xx/24.
poids(0)=-(xx-1.)*(xx-2.)*(xx-3.)*(xx+1.)/6.
- poids(1)=xx*(xx+1.)*(xx-3.)*(xx-2.)/4.
- case(7)
+ poids(1)=xx**4/12.-2.*xx**3/3.+5.*xx**2/12.+7.*xx/6.
+ case(3)
poids(-3)=0.
poids(-2)=xx*(xx-1.)*(xx-2.)*(xx+1.)/24.
poids(-1)=-xx*(xx-1.)*(xx**2-4.)/6.
- poids(0)=1.+xx**4/12.-13.*xx**2/12.-xx**3/3.+xx/3.
- case(8)
- poids(-3)=0.
- poids(-2)=0.
- poids(-1)=(xx-1.)*(xx-2.)*(xx-3.)*xx/24.
- poids(0)=-(xx-1.)*(xx-2.)*(xx-3.)*(xx+1.)/6.
- poids(1)=xx**4/12.-2.*xx**3/3.+5.*xx**2/12.+7.*xx/6.
- case(9)
+ poids(0)=1.+xx**4/12.-13.*xx**2/12.-xx**3/3.+xx/3.
+ case(4)
poids(-3)=0.
poids(-2)=0.
poids(-1)=xx*(xx-1.)*(xx-2.)*(xx-3.)/24.
poids(0)=xx**4/12.+xx**3/3.-13.*xx**2/12.-xx/3.+1.
- poids(1)=-xx*(xx+1.)*(xx+2.)*(xx-2.)/6.
- case(10)
- poids(-3)=0.
- poids(-2)=0.
- poids(-1)=-xx**4/8.+xx**3/12.+5.*xx**2/8.-7.*xx/12.
- poids(0)=(xx**2-4.)*(xx**2-1.)/4.
- case(11)
+ case(5)
poids(-3)=0.
poids(-2)=0.
poids(-1)=0.
poids(0)=1.-xx**4/8.+7.*xx**3/12.-3.*xx**2/8.-13.*xx/12.
poids(1)=xx*(xx+1.)*(xx-3.)*(xx-2.)/4.
+ case(6)
+ poids(-3)=0.
+ poids(-2)=0.
+ poids(-1)=-xx**4/8.+xx**3/12.+5.*xx**2/8.-7.*xx/12.
+ poids(0)=(xx**2-4.)*(xx**2-1.)/4.
+ case(7)
+ poids(-3)=0.
+ poids(-2)=xx*(xx-2.)*(xx**2-1.)/24.
+ poids(-1)=-xx*(xx-1.)*(xx**2-4.)/6.
+ poids(0)=0.5*xx**3-xx**2-0.5*xx+1.
+ poids(1)=xx*(xx+1.)*(xx-3.)*(xx-2.)/4.
+ case(8)
+ poids(-3)=(xx+2.)*xx*(xx**2-1.)/24.
+ poids(-2)=-xx*(xx+3.)*(xx**2-1.)/6.
+ poids(-1)=0.5*xx**3+0.5*xx**2-xx
+ poids(0)=(xx**2-4.)*(xx**2-1.)/4.
+ case(9)
+ poids(-3)=0.
+ poids(-2)=xx*(xx-2.)*(xx**2-1.)/24.
+ poids(-1)=-xx**3/6.+0.5*xx**2-xx/3.
+ poids(0)=-(xx-1.)*(xx-2.)*(xx-3.)*(xx+1.)/6.
+ poids(1)=xx*(xx+1.)*(xx-3.)*(xx-2.)/4.
+ case(10)
+ poids(-3)=(xx+2.)*xx*(xx**2-1.)/24.
+ poids(-2)=xx*(1.-xx**2)/6.
+ poids(-1)=-xx*(xx-1.)*(xx**2-4.)/6.
+ poids(0)=(xx**2-4.)*(xx**2-1.)/4.
end select
select case (blocd(i))
@@ -3572,53 +3609,44 @@ subroutine remaill_l8_z (donne,posx,posy,posz,remaille)
poids(4)=0.
end if
case(1)
- poids(0)=(xx**2-4.)*(xx**2-1.)/4. !c
poids(1)=-xx*(xx+1.)*(xx+2.)*(xx-2.)/6. !d
poids(2)=xx*(xx+1.)*(xx+2.)*(xx-1.)/24. !e
poids(3)=0.
poids(4)=0.
case(2)
- poids(1)=xx*(xx+1.)*(xx-3.)*(xx-2.)/4. !c'
- poids(2)=-xx*(xx+1.)*(xx-3.)*(xx-1.)/6.
- poids(3)=xx*(xx-2.)*(xx**2-1.)/24.
+ poids(2)=-xx**4/8.+5.*xx**3/12.+xx**2/8.-5.*xx/12.
+ poids(3)=0.
poids(4)=0.
case(3)
- poids(1)=-0.5*xx**3+0.5*xx**2+xx
- poids(2)=-xx*(xx-3)*(xx**2-1.)/6.
- poids(3)=xx*(xx-2.)*(xx**2-1.)/24.
+ poids(1)=-xx**4/8.-xx**3/12.+5.*xx**2/8.+7.*xx/12.
+ poids(2)=0.
+ poids(3)=0.
poids(4)=0.
case(4)
- poids(1)=-xx*(xx+1.)*(xx+2.)*(xx-2.)/6.
- poids(2)=xx*(xx**2-1.)/6.
- poids(3)=xx*(xx-2.)*(xx**2-1.)/24.
+ poids(2)=xx*(xx+2.)*(xx**2-1.)/24.
+ poids(3)=0.
poids(4)=0.
case(5)
- poids(2)=-xx*(xx+1.)*(xx-3.)*(xx-1.)/6.
- poids(3)=xx**3/6.-0.5*xx*2+xx/3.
- poids(4)=xx*(xx-1)*(xx-2.)*(xx-3.)/24.
- case(6)
- poids(1)=-xx**4/8.-xx**3/12.+5.*xx**2/8.+7.*xx/12.
+ poids(1)=xx*(xx+1.)*(xx+2.)*(xx+3.)/24.
poids(2)=0.
poids(3)=0.
poids(4)=0.
+ case(6)
+ poids(2)=-xx*(xx+1.)*(xx-3.)*(xx-1.)/6.
+ poids(3)=xx**3/6.-0.5*xx*2+xx/3.
+ poids(4)=xx*(xx-1)*(xx-2.)*(xx-3.)/24.
case(7)
- poids(2)=-xx**4/8.+5.*xx**3/12.+xx**2/8.-5.*xx/12.
- poids(3)=0.
+ poids(1)=-xx*(xx+1.)*(xx+2.)*(xx-2.)/6.
+ poids(2)=xx*(xx**2-1.)/6.
+ poids(3)=xx*(xx-2.)*(xx**2-1.)/24.
poids(4)=0.
case(8)
- poids(1)=xx*(xx+1.)*(xx+2.)*(xx+3.)/24.
- poids(2)=0.
- poids(3)=0.
- poids(4)=0.
- case(9)
- poids(2)=xx*(xx+2.)*(xx**2-1.)/24.
- poids(3)=0.
- poids(4)=0.
- case(10)
- poids(2)=xx*(xx+1.)*(xx+2.)*(xx-1.)/24.
- poids(3)=0.
+ poids(1)=-0.5*xx**3+0.5*xx**2+xx
+ poids(2)=-xx*(xx-3)*(xx**2-1.)/6.
+ poids(3)=xx*(xx-2.)*(xx**2-1.)/24.
poids(4)=0.
end select
+
@@ -3677,9 +3705,7 @@ subroutine remaill_l8_z (donne,posx,posy,posz,remaille)
!calcul des poidss
!----------------
-
- select case (blocg(i))
-
+ select case (blocg(i))
case(0)
if (xx<=0.5) then
poids(-3)=0.
@@ -3701,57 +3727,52 @@ subroutine remaill_l8_z (donne,posx,posy,posz,remaille)
case(2)
poids(-3)=0.
poids(-2)=0.
- poids(-1)=(xx-1.)*(xx-2.)*(xx-3.)*xx/24. !a'
- poids(0)=-(xx-1.)*(xx-2.)*(xx-3.)*(xx+1.)/6. !b'
- poids(1)=xx*(xx+1.)*(xx-3.)*(xx-2.)/4. !c'
- case(3)
- poids(-3)=(xx+2.)*xx*(xx**2-1.)/24.
- poids(-2)=-xx*(xx+3.)*(xx**2-1.)/6.
- poids(-1)=0.5*xx**3+0.5*xx**2-xx
- case(4)
- poids(-3)=0.
- poids(-2)=xx*(xx-2.)*(xx**2-1.)/24.
- poids(-1)=-xx*(xx-1.)*(xx**2-4.)/6.
- poids(0)=0.5*xx**3-xx**2-0.5*xx+1.
- case(5)
- poids(-3)=(xx+2.)*xx*(xx**2-1.)/24.
- poids(-2)=xx*(1.-xx**2)/6.
- poids(-1)=-xx*(xx-1.)*(xx**2-4.)/6.
- poids(0)=(xx**2-4.)*(xx**2-1.)/4.
- case(6)
- poids(-3)=0.
- poids(-2)=xx*(xx-2.)*(xx**2-1.)/24.
- poids(-1)=-xx**3/6.+0.5*xx**2-xx/3.
+ poids(-1)=(xx-1.)*(xx-2.)*(xx-3.)*xx/24.
poids(0)=-(xx-1.)*(xx-2.)*(xx-3.)*(xx+1.)/6.
- poids(1)=xx*(xx+1.)*(xx-3.)*(xx-2.)/4.
- case(7)
+ poids(1)=xx**4/12.-2.*xx**3/3.+5.*xx**2/12.+7.*xx/6.
+ case(3)
poids(-3)=0.
poids(-2)=xx*(xx-1.)*(xx-2.)*(xx+1.)/24.
poids(-1)=-xx*(xx-1.)*(xx**2-4.)/6.
- poids(0)=1.+xx**4/12.-13.*xx**2/12.-xx**3/3.+xx/3.
- case(8)
- poids(-3)=0.
- poids(-2)=0.
- poids(-1)=(xx-1.)*(xx-2.)*(xx-3.)*xx/24.
- poids(0)=-(xx-1.)*(xx-2.)*(xx-3.)*(xx+1.)/6.
- poids(1)=xx**4/12.-2.*xx**3/3.+5.*xx**2/12.+7.*xx/6.
- case(9)
+ poids(0)=1.+xx**4/12.-13.*xx**2/12.-xx**3/3.+xx/3.
+ case(4)
poids(-3)=0.
poids(-2)=0.
poids(-1)=xx*(xx-1.)*(xx-2.)*(xx-3.)/24.
poids(0)=xx**4/12.+xx**3/3.-13.*xx**2/12.-xx/3.+1.
- poids(1)=-xx*(xx+1.)*(xx+2.)*(xx-2.)/6.
- case(10)
- poids(-3)=0.
- poids(-2)=0.
- poids(-1)=-xx**4/8.+xx**3/12.+5.*xx**2/8.-7.*xx/12.
- poids(0)=(xx**2-4.)*(xx**2-1.)/4.
- case(11)
+ case(5)
poids(-3)=0.
poids(-2)=0.
poids(-1)=0.
poids(0)=1.-xx**4/8.+7.*xx**3/12.-3.*xx**2/8.-13.*xx/12.
poids(1)=xx*(xx+1.)*(xx-3.)*(xx-2.)/4.
+ case(6)
+ poids(-3)=0.
+ poids(-2)=0.
+ poids(-1)=-xx**4/8.+xx**3/12.+5.*xx**2/8.-7.*xx/12.
+ poids(0)=(xx**2-4.)*(xx**2-1.)/4.
+ case(7)
+ poids(-3)=0.
+ poids(-2)=xx*(xx-2.)*(xx**2-1.)/24.
+ poids(-1)=-xx*(xx-1.)*(xx**2-4.)/6.
+ poids(0)=0.5*xx**3-xx**2-0.5*xx+1.
+ poids(1)=xx*(xx+1.)*(xx-3.)*(xx-2.)/4.
+ case(8)
+ poids(-3)=(xx+2.)*xx*(xx**2-1.)/24.
+ poids(-2)=-xx*(xx+3.)*(xx**2-1.)/6.
+ poids(-1)=0.5*xx**3+0.5*xx**2-xx
+ poids(0)=(xx**2-4.)*(xx**2-1.)/4.
+ case(9)
+ poids(-3)=0.
+ poids(-2)=xx*(xx-2.)*(xx**2-1.)/24.
+ poids(-1)=-xx**3/6.+0.5*xx**2-xx/3.
+ poids(0)=-(xx-1.)*(xx-2.)*(xx-3.)*(xx+1.)/6.
+ poids(1)=xx*(xx+1.)*(xx-3.)*(xx-2.)/4.
+ case(10)
+ poids(-3)=(xx+2.)*xx*(xx**2-1.)/24.
+ poids(-2)=xx*(1.-xx**2)/6.
+ poids(-1)=-xx*(xx-1.)*(xx**2-4.)/6.
+ poids(0)=(xx**2-4.)*(xx**2-1.)/4.
end select
select case (blocd(i))
@@ -3767,54 +3788,43 @@ subroutine remaill_l8_z (donne,posx,posy,posz,remaille)
poids(4)=0.
end if
case(1)
- poids(0)=(xx**2-4.)*(xx**2-1.)/4. !c
poids(1)=-xx*(xx+1.)*(xx+2.)*(xx-2.)/6. !d
poids(2)=xx*(xx+1.)*(xx+2.)*(xx-1.)/24. !e
poids(3)=0.
poids(4)=0.
case(2)
- poids(1)=xx*(xx+1.)*(xx-3.)*(xx-2.)/4. !c'
- poids(2)=-xx*(xx+1.)*(xx-3.)*(xx-1.)/6.
- poids(3)=xx*(xx-2.)*(xx**2-1.)/24.
+ poids(2)=-xx**4/8.+5.*xx**3/12.+xx**2/8.-5.*xx/12.
+ poids(3)=0.
poids(4)=0.
case(3)
- poids(1)=-0.5*xx**3+0.5*xx**2+xx
- poids(2)=-xx*(xx-3)*(xx**2-1.)/6.
- poids(3)=xx*(xx-2.)*(xx**2-1.)/24.
+ poids(1)=-xx**4/8.-xx**3/12.+5.*xx**2/8.+7.*xx/12.
+ poids(2)=0.
+ poids(3)=0.
poids(4)=0.
case(4)
- poids(1)=-xx*(xx+1.)*(xx+2.)*(xx-2.)/6.
- poids(2)=xx*(xx**2-1.)/6.
- poids(3)=xx*(xx-2.)*(xx**2-1.)/24.
+ poids(2)=xx*(xx+2.)*(xx**2-1.)/24.
+ poids(3)=0.
poids(4)=0.
case(5)
- poids(2)=-xx*(xx+1.)*(xx-3.)*(xx-1.)/6.
- poids(3)=xx**3/6.-0.5*xx*2+xx/3.
- poids(4)=xx*(xx-1)*(xx-2.)*(xx-3.)/24.
- case(6)
- poids(1)=-xx**4/8.-xx**3/12.+5.*xx**2/8.+7.*xx/12.
+ poids(1)=xx*(xx+1.)*(xx+2.)*(xx+3.)/24.
poids(2)=0.
poids(3)=0.
poids(4)=0.
+ case(6)
+ poids(2)=-xx*(xx+1.)*(xx-3.)*(xx-1.)/6.
+ poids(3)=xx**3/6.-0.5*xx*2+xx/3.
+ poids(4)=xx*(xx-1)*(xx-2.)*(xx-3.)/24.
case(7)
- poids(2)=-xx**4/8.+5.*xx**3/12.+xx**2/8.-5.*xx/12.
- poids(3)=0.
+ poids(1)=-xx*(xx+1.)*(xx+2.)*(xx-2.)/6.
+ poids(2)=xx*(xx**2-1.)/6.
+ poids(3)=xx*(xx-2.)*(xx**2-1.)/24.
poids(4)=0.
case(8)
- poids(1)=xx*(xx+1.)*(xx+2.)*(xx+3.)/24.
- poids(2)=0.
- poids(3)=0.
- poids(4)=0.
- case(9)
- poids(2)=xx*(xx+2.)*(xx**2-1.)/24.
- poids(3)=0.
- poids(4)=0.
- case(10)
- poids(2)=xx*(xx+1.)*(xx+2.)*(xx-1.)/24.
- poids(3)=0.
+ poids(1)=-0.5*xx**3+0.5*xx**2+xx
+ poids(2)=-xx*(xx-3)*(xx**2-1.)/6.
+ poids(3)=xx*(xx-2.)*(xx**2-1.)/24.
poids(4)=0.
end select
-
!remaillage à l' interrieur domaine
@@ -3872,9 +3882,7 @@ subroutine remaill_l8_z (donne,posx,posy,posz,remaille)
!calcul des poidss
!----------------
-
- select case (blocg(i))
-
+ select case (blocg(i))
case(0)
if (xx<=0.5) then
poids(-3)=0.
@@ -3896,57 +3904,52 @@ subroutine remaill_l8_z (donne,posx,posy,posz,remaille)
case(2)
poids(-3)=0.
poids(-2)=0.
- poids(-1)=(xx-1.)*(xx-2.)*(xx-3.)*xx/24. !a'
- poids(0)=-(xx-1.)*(xx-2.)*(xx-3.)*(xx+1.)/6. !b'
- poids(1)=xx*(xx+1.)*(xx-3.)*(xx-2.)/4. !c'
- case(3)
- poids(-3)=(xx+2.)*xx*(xx**2-1.)/24.
- poids(-2)=-xx*(xx+3.)*(xx**2-1.)/6.
- poids(-1)=0.5*xx**3+0.5*xx**2-xx
- case(4)
- poids(-3)=0.
- poids(-2)=xx*(xx-2.)*(xx**2-1.)/24.
- poids(-1)=-xx*(xx-1.)*(xx**2-4.)/6.
- poids(0)=0.5*xx**3-xx**2-0.5*xx+1.
- case(5)
- poids(-3)=(xx+2.)*xx*(xx**2-1.)/24.
- poids(-2)=xx*(1.-xx**2)/6.
- poids(-1)=-xx*(xx-1.)*(xx**2-4.)/6.
- poids(0)=(xx**2-4.)*(xx**2-1.)/4.
- case(6)
- poids(-3)=0.
- poids(-2)=xx*(xx-2.)*(xx**2-1.)/24.
- poids(-1)=-xx**3/6.+0.5*xx**2-xx/3.
+ poids(-1)=(xx-1.)*(xx-2.)*(xx-3.)*xx/24.
poids(0)=-(xx-1.)*(xx-2.)*(xx-3.)*(xx+1.)/6.
- poids(1)=xx*(xx+1.)*(xx-3.)*(xx-2.)/4.
- case(7)
+ poids(1)=xx**4/12.-2.*xx**3/3.+5.*xx**2/12.+7.*xx/6.
+ case(3)
poids(-3)=0.
poids(-2)=xx*(xx-1.)*(xx-2.)*(xx+1.)/24.
poids(-1)=-xx*(xx-1.)*(xx**2-4.)/6.
- poids(0)=1.+xx**4/12.-13.*xx**2/12.-xx**3/3.+xx/3.
- case(8)
- poids(-3)=0.
- poids(-2)=0.
- poids(-1)=(xx-1.)*(xx-2.)*(xx-3.)*xx/24.
- poids(0)=-(xx-1.)*(xx-2.)*(xx-3.)*(xx+1.)/6.
- poids(1)=xx**4/12.-2.*xx**3/3.+5.*xx**2/12.+7.*xx/6.
- case(9)
+ poids(0)=1.+xx**4/12.-13.*xx**2/12.-xx**3/3.+xx/3.
+ case(4)
poids(-3)=0.
poids(-2)=0.
poids(-1)=xx*(xx-1.)*(xx-2.)*(xx-3.)/24.
poids(0)=xx**4/12.+xx**3/3.-13.*xx**2/12.-xx/3.+1.
- poids(1)=-xx*(xx+1.)*(xx+2.)*(xx-2.)/6.
- case(10)
- poids(-3)=0.
- poids(-2)=0.
- poids(-1)=-xx**4/8.+xx**3/12.+5.*xx**2/8.-7.*xx/12.
- poids(0)=(xx**2-4.)*(xx**2-1.)/4.
- case(11)
+ case(5)
poids(-3)=0.
poids(-2)=0.
poids(-1)=0.
poids(0)=1.-xx**4/8.+7.*xx**3/12.-3.*xx**2/8.-13.*xx/12.
poids(1)=xx*(xx+1.)*(xx-3.)*(xx-2.)/4.
+ case(6)
+ poids(-3)=0.
+ poids(-2)=0.
+ poids(-1)=-xx**4/8.+xx**3/12.+5.*xx**2/8.-7.*xx/12.
+ poids(0)=(xx**2-4.)*(xx**2-1.)/4.
+ case(7)
+ poids(-3)=0.
+ poids(-2)=xx*(xx-2.)*(xx**2-1.)/24.
+ poids(-1)=-xx*(xx-1.)*(xx**2-4.)/6.
+ poids(0)=0.5*xx**3-xx**2-0.5*xx+1.
+ poids(1)=xx*(xx+1.)*(xx-3.)*(xx-2.)/4.
+ case(8)
+ poids(-3)=(xx+2.)*xx*(xx**2-1.)/24.
+ poids(-2)=-xx*(xx+3.)*(xx**2-1.)/6.
+ poids(-1)=0.5*xx**3+0.5*xx**2-xx
+ poids(0)=(xx**2-4.)*(xx**2-1.)/4.
+ case(9)
+ poids(-3)=0.
+ poids(-2)=xx*(xx-2.)*(xx**2-1.)/24.
+ poids(-1)=-xx**3/6.+0.5*xx**2-xx/3.
+ poids(0)=-(xx-1.)*(xx-2.)*(xx-3.)*(xx+1.)/6.
+ poids(1)=xx*(xx+1.)*(xx-3.)*(xx-2.)/4.
+ case(10)
+ poids(-3)=(xx+2.)*xx*(xx**2-1.)/24.
+ poids(-2)=xx*(1.-xx**2)/6.
+ poids(-1)=-xx*(xx-1.)*(xx**2-4.)/6.
+ poids(0)=(xx**2-4.)*(xx**2-1.)/4.
end select
select case (blocd(i))
@@ -3962,55 +3965,44 @@ subroutine remaill_l8_z (donne,posx,posy,posz,remaille)
poids(4)=0.
end if
case(1)
- poids(0)=(xx**2-4.)*(xx**2-1.)/4. !c
poids(1)=-xx*(xx+1.)*(xx+2.)*(xx-2.)/6. !d
poids(2)=xx*(xx+1.)*(xx+2.)*(xx-1.)/24. !e
poids(3)=0.
poids(4)=0.
case(2)
- poids(1)=xx*(xx+1.)*(xx-3.)*(xx-2.)/4. !c'
- poids(2)=-xx*(xx+1.)*(xx-3.)*(xx-1.)/6.
- poids(3)=xx*(xx-2.)*(xx**2-1.)/24.
+ poids(2)=-xx**4/8.+5.*xx**3/12.+xx**2/8.-5.*xx/12.
+ poids(3)=0.
poids(4)=0.
case(3)
- poids(1)=-0.5*xx**3+0.5*xx**2+xx
- poids(2)=-xx*(xx-3)*(xx**2-1.)/6.
- poids(3)=xx*(xx-2.)*(xx**2-1.)/24.
+ poids(1)=-xx**4/8.-xx**3/12.+5.*xx**2/8.+7.*xx/12.
+ poids(2)=0.
+ poids(3)=0.
poids(4)=0.
case(4)
- poids(1)=-xx*(xx+1.)*(xx+2.)*(xx-2.)/6.
- poids(2)=xx*(xx**2-1.)/6.
- poids(3)=xx*(xx-2.)*(xx**2-1.)/24.
+ poids(2)=xx*(xx+2.)*(xx**2-1.)/24.
+ poids(3)=0.
poids(4)=0.
case(5)
- poids(2)=-xx*(xx+1.)*(xx-3.)*(xx-1.)/6.
- poids(3)=xx**3/6.-0.5*xx*2+xx/3.
- poids(4)=xx*(xx-1)*(xx-2.)*(xx-3.)/24.
- case(6)
- poids(1)=-xx**4/8.-xx**3/12.+5.*xx**2/8.+7.*xx/12.
+ poids(1)=xx*(xx+1.)*(xx+2.)*(xx+3.)/24.
poids(2)=0.
poids(3)=0.
poids(4)=0.
+ case(6)
+ poids(2)=-xx*(xx+1.)*(xx-3.)*(xx-1.)/6.
+ poids(3)=xx**3/6.-0.5*xx*2+xx/3.
+ poids(4)=xx*(xx-1)*(xx-2.)*(xx-3.)/24.
case(7)
- poids(2)=-xx**4/8.+5.*xx**3/12.+xx**2/8.-5.*xx/12.
- poids(3)=0.
+ poids(1)=-xx*(xx+1.)*(xx+2.)*(xx-2.)/6.
+ poids(2)=xx*(xx**2-1.)/6.
+ poids(3)=xx*(xx-2.)*(xx**2-1.)/24.
poids(4)=0.
case(8)
- poids(1)=xx*(xx+1.)*(xx+2.)*(xx+3.)/24.
- poids(2)=0.
- poids(3)=0.
- poids(4)=0.
- case(9)
- poids(2)=xx*(xx+2.)*(xx**2-1.)/24.
- poids(3)=0.
- poids(4)=0.
- case(10)
- poids(2)=xx*(xx+1.)*(xx+2.)*(xx-1.)/24.
- poids(3)=0.
+ poids(1)=-0.5*xx**3+0.5*xx**2+xx
+ poids(2)=-xx*(xx-3)*(xx**2-1.)/6.
+ poids(3)=xx*(xx-2.)*(xx**2-1.)/24.
poids(4)=0.
end select
-
!remaillage à l' interrieur domaine
!---------------------------------
do c=-3,4
@@ -4024,6 +4016,918 @@ subroutine remaill_l8_z (donne,posx,posy,posz,remaille)
+ subroutine limit_l2_bloc_x (donne,posx,posy,posz,remaille)
+ implicit none
+ real(kind=8),dimension(1:npart),intent(in) :: donne
+ real(kind=8),dimension(1:npart),intent(in) :: posx,posy,posz
+ real(kind=8),dimension(1:nx,1:ny,1:nz),intent(out) :: remaille
+ integer :: i,c,d,e,ib,jb,kb,ivar,jvar,kvar,j,per
+ integer,dimension(-2:2) :: ip
+ real(kind=8),dimension(-2:2) :: poidx
+ real(kind=8) :: xx,tmp_pos
+ real(kind=8),dimension(0:2) :: phi_mdemi,phi_pdemi
+ real(kind=8),dimension(0:1) :: phi
+
+
+ remaille=0.
+
+
+ do i=1,npart
+
+ jvar=nint((posy(i)-yb)/dy)+1
+ kvar=nint((posz(i)-zd)/dz)+1
+
+ poidx=0.
+
+ !numero des points sur le maillage
+ !--------------------------------
+ ip(0) = floor((posx(i)-xg)/dx)
+ ip(-1) = ip(0) - 1
+ ip(-2) = ip(0) - 2
+ ip(1) = ip(0) + 1
+ ip(2) = ip(0) + 2
+
+
+
+
+
+ !distance de la particule à remailler au second point (des quatres utilisés pour le remaillage)
+ !----------------------------------------------------
+ xx = (posx(i) - real(ip(0),kind=8)*dx-xg)/dx !relatif
+
+
+
+
+ !conditions au bord
+ !------------------
+ !periodique:
+ do c=-2,2
+ ip(c)=mod(ip(c)+nx,nx)
+ end do
+
+ !calcul du limiteur pour cas de bloc centre avec xx>0.5
+ !--------------------------------------------------------
+ call calcul_phi_tvd(phi,1)
+ phi_mdemi(1)=phi(0)
+ phi_pdemi(1)=phi(1)
+
+ !les autres cas ...
+ !---------------------
+ call calcul_phi_tvd(phi,0)
+ phi_mdemi(0)=phi(0)
+ phi_pdemi(0)=phi(1)
+
+ phi_mdemi(2)=0.
+ phi_pdemi(2)=0.
+
+ select case (blocg(i))
+
+ case(0)
+ if (xx<=0.5) then
+
+ poidx(-1)=ca(xx,0)
+ if (test (i,1) ) then
+ poidx(0)=cb(xx,0,1)
+ poidx(1)=cc(xx,1)
+ else
+ poidx(0)=cb(xx,0,0)
+ poidx(1)=cc(xx,0)
+ end if
+
+ else
+
+ poidx(0)=ca(xx-1.,1)
+ if (test (i,1) ) then
+ poidx(1)=cb(xx-1.,1,1)
+ poidx(2)=cc(xx-1.,1)
+ else
+ poidx(1)=cb(xx-1.,1,0)
+ poidx(2)=cc(xx-1.,0)
+ end if
+
+ end if
+
+ case(1)
+
+ poidx(-1)=ca(xx,0)
+ if (test (i,1) ) then
+ poidx(0)=cb(xx,0,1)
+ poidx(1)=cc(xx,1)
+ else
+ poidx(0)=cb(xx,0,0)
+ poidx(1)=cc(xx,0)
+ end if
+
+
+ case(2)
+ poidx(0)=ca(xx-1.,1)
+ poidx(1)=cb(xx-1.,1,0)+cc(xx-1.,0)
+ case(3)
+ poidx(-1)=ca(xx,0)
+ poidx(0)=cb(xx,0)+cc(xx,0)
+ case(4)
+ poidx(0)=ca(xx,1)+cb(xx,1,0)
+ poidx(1)=cc(xx,0)
+
+
+
+ case(5)
+ poidx(-1)=ca(xx,0)
+ poidx(0)=cb(xx,0)
+ poidx(1)=cc(xx,0)-cc(xx-1.,1)
+ poidx(2)=cc(xx-1.,1)
+
+ case(6)
+ poidx(-1)=ca(xx,0)
+ poidx(0)=ca(xx-1.,1)-ca(xx,0)
+ poidx(1)=cb(xx-1.,1)
+ poidx(2)=cc(xx-1.,1)
+
+
+ case(7)
+ poidx(-2)=ca(xx+1.,0)
+ poidx(-1)=ca(xx,0)-ca(xx+1.,0)
+ poidx(0)=cb(xx,0)
+ poidx(1)=cc(xx,0)
+ end select
+
+ !remaillage à l' interrieur domaine
+ !---------------------------------
+ do c=-2,2
+ remaille(ip(c)+1,jvar,kvar)=remaille(ip(c)+1,jvar,kvar)+donne(i)*poidx(c)
+ end do
+
+ end do
+
+ contains
+
+ function coeff_l2g (num,x) result(c)
+ implicit none
+ integer,intent(in) :: num
+ real(kind=8),intent(in) :: x
+ real(kind=8) :: c
+
+ select case(num)
+ case(-1)
+ c=0.5*x*(x-1.)
+ case(0)
+ c=1.-x**2
+ case(1)
+ c=0.5*x*(1.+x)
+ end select
+
+ end function coeff_l2g
+
+ function coeff_tscg (num,x) result(c)
+ implicit none
+ integer,intent(in) :: num
+ real(kind=8),intent(in) :: x
+ real(kind=8) :: c
+
+ select case(num)
+ case(-1)
+ c=0.5*x*(x-1.) +1./8. ! +1./6.
+ case(0)
+ c=-x**2+1. -2./8. !-2./6.
+ case(1)
+ c=0.5*x*(1.+x) +1./8. ! +1./6.
+ end select
+
+ end function coeff_tscg
+
+ function ca(x,num)
+ integer :: num
+ real(kind=8) :: x
+ real(kind=8) :: ca
+
+ ca=coeff_tscg(-1,x)+phi_mdemi(num)*(coeff_l2g(-1,x)-coeff_tscg(-1,x))
+
+ end function ca
+
+ function cb(x,num,num2)
+ integer :: num
+ real(kind=8) :: x
+ integer,optional :: num2
+ real(kind=8) :: cb
+
+ if (present (num2)) then
+ cb=1.-(coeff_tscg(-1,x)+coeff_tscg(1,x))-phi_pdemi(num2)*(coeff_l2g(1,x)-coeff_tscg(1,x))-phi_mdemi(num)*(coeff_l2g(-1,x)-coeff_tscg(-1,x))
+ else
+ cb=1.-(coeff_tscg(-1,x)+coeff_tscg(1,x))-phi_pdemi(num)*(coeff_l2g(1,x)-coeff_tscg(1,x))-phi_mdemi(num)*(coeff_l2g(-1,x)-coeff_tscg(-1,x))
+ end if
+
+ end function cb
+
+ function cc(x,num)
+ integer :: num
+ real(kind=8) :: x
+ real(kind=8) :: cc
+
+ cc=coeff_tscg(1,x)+phi_pdemi(num)*(coeff_l2g(1,x)-coeff_tscg(1,x))
+
+
+ end function cc
+
+ function test (ipart,vois)
+ implicit none
+ integer :: ipart,vois
+ logical :: test
+ real(kind=8) :: x
+
+ test=.false.
+ x=(posx(mod(ipart-1+vois+npart,npart)+1) - real( floor((posx(mod(ipart-1+vois+npart,npart)+1)-xg)/dx) ,kind=8)*dx-xg)/dx
+ if ( ((blocg(mod(ipart-1+vois+npart,npart)+1)==0).or.(blocg(mod(ipart-1+vois+npart,npart)+1)==2)).and.(x>0.5)) test=.true.
+
+ end function test
+
+
+
+ subroutine calcul_phi_tvd(phi_out,num)
+ implicit none
+ integer :: num
+ real(kind=8),dimension(0:1),intent(out):: phi_out
+ real(kind=8),dimension(-5:5) :: u
+ real(kind=8) :: rp,rm,diff1,diff2,diff3,diff0,rmm,y,r
+ integer :: ib,jb
+ real(kind=8),dimension(-5:5) :: xx_vois,x,tp1,tp2,tp3
+
+ !particules voisines
+ !-------------------
+ do ib=-2,2
+ u(ib)=donne(mod(i-1+ib+npart,npart)+1)
+ end do
+
+ do ib=-1,2
+ x(ib)=(posx(mod(i-1+ib+npart,npart)+1) - real( floor((posx(mod(i-1+ib+npart,npart)+1)-xg)/dx) ,kind=8)*dx-xg)/dx
+ end do
+
+
+ !cas bloc centre xx>0.5 (psi)
+ !-----------------------
+ if (num==1) then
+
+ do ib=0,1
+ x(ib)=x(ib)-1.
+ tp1(ib)=6.-8.*x(ib)*x(ib)
+ tp3(ib)=4.*(x(ib)-0.5)**2
+ end do
+
+ do ib=0,1
+ r=(u(ib+1)-u(ib))/(u(ib)-u(ib-1))
+
+ !anti diffusif
+ !phi_out(ib)=max(0.,min(tp3(ib),tp1(ib)*r))
+
+ !du type Minmod (attention non justifié pour y>=0.8)
+ phi_out(ib)=max(0.,min(1.,4.*r))
+
+ !bornes dépendent de la condition cfl
+ !phi_out(ib)=max(0.,min(1.,tp3(ib),tp1(ib)*r))
+
+ if(abs(u(ib)-u(ib-1))<0.0000001) phi_out(ib)=0.
+ end do
+
+
+
+
+ else
+ !sinon (phi)
+ !-----------------------
+
+
+ do ib=-1,0
+ tp1(ib)=6.-8.*x(ib)*x(ib)
+ tp2(ib)=4.*(x(ib)+0.5)**2
+ end do
+
+
+ do ib=0,1
+ r=(u(ib-1)-u(ib-2))/(u(ib)-u(ib-1))
+ !phi_out(ib)=max(0.,min(tp2(ib-1),tp1(ib-1)*r))
+ phi_out(ib)=max(0.,min(1.,4.*r))
+ !phi_out(ib)=max(0.,min(1.,tp2(ib-1),tp1(ib-1)*r))
+
+ if(abs(u(ib)-u(ib-1))<0.0000001) phi_out(ib)=0.
+ end do
+
+ end if
+
+
+ end subroutine calcul_phi_tvd
+
+
+ end subroutine limit_l2_bloc_x
+
+
+
+
+ subroutine limit_l2_bloc_y (donne,posx,posy,posz,remaille)
+ implicit none
+ real(kind=8),dimension(1:npart),intent(in) :: donne
+ real(kind=8),dimension(1:npart),intent(in) :: posx,posy,posz
+ real(kind=8),dimension(1:nx,1:ny,1:nz),intent(out) :: remaille
+ integer :: i,c,d,e,ib,jb,kb,ivar,jvar,kvar,j,per
+ integer,dimension(-2:2) :: ip
+ real(kind=8),dimension(-2:2) :: poidx
+ real(kind=8) :: xx,tmp_pos
+ real(kind=8),dimension(0:2) :: phi_mdemi,phi_pdemi
+ real(kind=8),dimension(0:1) :: phi
+
+
+ remaille=0.
+
+
+ do i=1,npart
+
+ ivar=nint((posx(i)-xg)/dx)+1
+ kvar=nint((posz(i)-zd)/dz)+1
+
+ poidx=0.
+
+ !numero des points sur le maillage
+ !--------------------------------
+ ip(0) = floor((posy(i)-yb)/dy)
+ ip(-1) = ip(0) - 1
+ ip(-2) = ip(0) - 2
+ ip(1) = ip(0) + 1
+ ip(2) = ip(0) + 2
+
+
+
+
+
+ !distance de la particule à remailler au second point (des quatres utilisés pour le remaillage)
+ !----------------------------------------------------
+ xx = (posy(i) - real(ip(0),kind=8)*dy-yb)/dy !relatif
+
+
+
+
+ !conditions au bord
+ !------------------
+ !periodique:
+ do c=-2,2
+ ip(c)=mod(ip(c)+ny,ny)
+ end do
+
+ !calcul du limiteur pour cas de bloc centre avec xx>0.5
+ !--------------------------------------------------------
+ call calcul_phi_tvd(phi,1)
+ phi_mdemi(1)=phi(0)
+ phi_pdemi(1)=phi(1)
+
+ !les autres cas ...
+ !---------------------
+ call calcul_phi_tvd(phi,0)
+ phi_mdemi(0)=phi(0)
+ phi_pdemi(0)=phi(1)
+
+ phi_mdemi(2)=0.
+ phi_pdemi(2)=0.
+
+ select case (blocg(i))
+
+ case(0)
+ if (xx<=0.5) then
+
+ poidx(-1)=ca(xx,0)
+ if (test (i,1) ) then
+ poidx(0)=cb(xx,0,1)
+ poidx(1)=cc(xx,1)
+ else
+ poidx(0)=cb(xx,0,0)
+ poidx(1)=cc(xx,0)
+ end if
+
+ else
+
+ poidx(0)=ca(xx-1.,1)
+ if (test (i,1) ) then
+ poidx(1)=cb(xx-1.,1,1)
+ poidx(2)=cc(xx-1.,1)
+ else
+ poidx(1)=cb(xx-1.,1,0)
+ poidx(2)=cc(xx-1.,0)
+ end if
+
+ end if
+
+ case(1)
+
+ poidx(-1)=ca(xx,0)
+ if (test (i,1) ) then
+ poidx(0)=cb(xx,0,1)
+ poidx(1)=cc(xx,1)
+ else
+ poidx(0)=cb(xx,0,0)
+ poidx(1)=cc(xx,0)
+ end if
+
+
+ case(2)
+ poidx(0)=ca(xx-1.,1)
+ poidx(1)=cb(xx-1.,1,0)+cc(xx-1.,0)
+ case(3)
+ poidx(-1)=ca(xx,0)
+ poidx(0)=cb(xx,0)+cc(xx,0)
+ case(4)
+ poidx(0)=ca(xx,1)+cb(xx,1,0)
+ poidx(1)=cc(xx,0)
+
+
+
+ case(5)
+ poidx(-1)=ca(xx,0)
+ poidx(0)=cb(xx,0)
+ poidx(1)=cc(xx,0)-cc(xx-1.,1)
+ poidx(2)=cc(xx-1.,1)
+
+ case(6)
+ poidx(-1)=ca(xx,0)
+ poidx(0)=ca(xx-1.,1)-ca(xx,0)
+ poidx(1)=cb(xx-1.,1)
+ poidx(2)=cc(xx-1.,1)
+
+
+ case(7)
+ poidx(-2)=ca(xx+1.,0)
+ poidx(-1)=ca(xx,0)-ca(xx+1.,0)
+ poidx(0)=cb(xx,0)
+ poidx(1)=cc(xx,0)
+ end select
+
+ !remaillage à l' interrieur domaine
+ !---------------------------------
+ do c=-2,2
+ remaille(ivar,ip(c)+1,kvar)=remaille(ivar,ip(c)+1,kvar)+donne(i)*poidx(c)
+ end do
+
+ end do
+
+ contains
+
+ function coeff_l2g (num,x) result(c)
+ implicit none
+ integer,intent(in) :: num
+ real(kind=8),intent(in) :: x
+ real(kind=8) :: c
+
+ select case(num)
+ case(-1)
+ c=0.5*x*(x-1.)
+ case(0)
+ c=1.-x**2
+ case(1)
+ c=0.5*x*(1.+x)
+ end select
+
+ end function coeff_l2g
+
+ function coeff_tscg (num,x) result(c)
+ implicit none
+ integer,intent(in) :: num
+ real(kind=8),intent(in) :: x
+ real(kind=8) :: c
+
+ select case(num)
+ case(-1)
+ c=0.5*x*(x-1.) +1./8.
+ case(0)
+ c=-x**2+1. -2./8.
+ case(1)
+ c=0.5*x*(1.+x) +1./8.
+ end select
+
+ end function coeff_tscg
+
+ function ca(x,num)
+ integer :: num
+ real(kind=8) :: x
+ real(kind=8) :: ca
+
+ ca=coeff_tscg(-1,x)+phi_mdemi(num)*(coeff_l2g(-1,x)-coeff_tscg(-1,x))
+
+ end function ca
+
+ function cb(x,num,num2)
+ integer :: num
+ real(kind=8) :: x
+ integer,optional :: num2
+ real(kind=8) :: cb
+
+ if (present (num2)) then
+ cb=1.-(coeff_tscg(-1,x)+coeff_tscg(1,x))-phi_pdemi(num2)*(coeff_l2g(1,x)-coeff_tscg(1,x))-phi_mdemi(num)*(coeff_l2g(-1,x)-coeff_tscg(-1,x))
+ else
+ cb=1.-(coeff_tscg(-1,x)+coeff_tscg(1,x))-phi_pdemi(num)*(coeff_l2g(1,x)-coeff_tscg(1,x))-phi_mdemi(num)*(coeff_l2g(-1,x)-coeff_tscg(-1,x))
+ end if
+
+ end function cb
+
+ function cc(x,num)
+ integer :: num
+ real(kind=8) :: x
+ real(kind=8) :: cc
+
+ cc=coeff_tscg(1,x)+phi_pdemi(num)*(coeff_l2g(1,x)-coeff_tscg(1,x))
+
+
+ end function cc
+
+ function test (ipart,vois)
+ implicit none
+ integer :: ipart,vois
+ logical :: test
+ real(kind=8) :: x
+
+ test=.false.
+ x=(posy(mod(ipart-1+vois+npart,npart)+1) - real( floor((posy(mod(ipart-1+vois+npart,npart)+1)-yb)/dy) ,kind=8)*dy-yb)/dy
+ if ( ((blocg(mod(ipart-1+vois+npart,npart)+1)==0).or.(blocg(mod(ipart-1+vois+npart,npart)+1)==2)).and.(x>0.5)) test=.true.
+
+ end function test
+
+
+
+ subroutine calcul_phi_tvd(phi_out,num)
+ implicit none
+ integer :: num
+ real(kind=8),dimension(0:1),intent(out):: phi_out
+ real(kind=8),dimension(-5:5) :: u
+ real(kind=8) :: rp,rm,diff1,diff2,diff3,diff0,rmm,y,r
+ integer :: ib,jb
+ real(kind=8),dimension(-5:5) :: xx_vois,x,tp1,tp2,tp3
+
+ !particules voisines
+ !-------------------
+ do ib=-2,2
+ u(ib)=donne(mod(i-1+ib+npart,npart)+1)
+ end do
+
+ do ib=-1,2
+ x(ib)=(posy(mod(i-1+ib+npart,npart)+1) - real( floor((posy(mod(i-1+ib+npart,npart)+1)-yb)/dy) ,kind=8)*dy-yb)/dy
+ end do
+
+
+ !cas bloc centre xx>0.5 (psi)
+ !-----------------------
+ if (num==1) then
+
+ do ib=0,1
+ x(ib)=x(ib)-1.
+ tp1(ib)=6.-8.*x(ib)*x(ib)
+ tp3(ib)=4.*(x(ib)-0.5)**2
+ end do
+
+ do ib=0,1
+ r=(u(ib+1)-u(ib))/(u(ib)-u(ib-1))
+
+ !anti diffusif
+ !phi_out(ib)=max(0.,min(tp3(ib),tp1(ib)*r))
+
+ !du type Minmod (attention non justifié pour y>=0.8)
+ phi_out(ib)=max(0.,min(1.,4.*r))
+
+ !bornes dépendent de la condition cfl
+ !phi_out(ib)=max(0.,min(1.,tp3(ib),tp1(ib)*r))
+
+ if(abs(u(ib)-u(ib-1))<0.0000001) phi_out(ib)=0.
+ end do
+
+
+
+
+ else
+ !sinon (phi)
+ !-----------------------
+
+
+ do ib=-1,0
+ tp1(ib)=6.-8.*x(ib)*x(ib)
+ tp2(ib)=4.*(x(ib)+0.5)**2
+ end do
+
+
+ do ib=0,1
+ r=(u(ib-1)-u(ib-2))/(u(ib)-u(ib-1))
+ !phi_out(ib)=max(0.,min(tp2(ib-1),tp1(ib-1)*r))
+ phi_out(ib)=max(0.,min(1.,4.*r))
+ !phi_out(ib)=max(0.,min(1.,tp2(ib-1),tp1(ib-1)*r))
+
+ if(abs(u(ib)-u(ib-1))<0.0000001) phi_out(ib)=0.
+ end do
+
+ end if
+
+
+ end subroutine calcul_phi_tvd
+
+
+ end subroutine limit_l2_bloc_y
+
+
+
+
+ subroutine limit_l2_bloc_z (donne,posx,posy,posz,remaille)
+ implicit none
+ real(kind=8),dimension(1:npart),intent(in) :: donne
+ real(kind=8),dimension(1:npart),intent(in) :: posx,posy,posz
+ real(kind=8),dimension(1:nx,1:ny,1:nz),intent(out) :: remaille
+ integer :: i,c,d,e,ib,jb,kb,ivar,jvar,kvar,j,per
+ integer,dimension(-2:2) :: ip
+ real(kind=8),dimension(-2:2) :: poidx
+ real(kind=8) :: xx,tmp_pos
+ real(kind=8),dimension(0:2) :: phi_mdemi,phi_pdemi
+ real(kind=8),dimension(0:1) :: phi
+
+
+ remaille=0.
+
+
+ do i=1,npart
+
+ ivar=nint((posx(i)-xg)/dx)+1
+ jvar=nint((posy(i)-yb)/dy)+1
+
+ poidx=0.
+
+ !numero des points sur le maillage
+ !--------------------------------
+ ip(0) = floor((posz(i)-zd)/dz)
+ ip(-1) = ip(0) - 1
+ ip(-2) = ip(0) - 2
+ ip(1) = ip(0) + 1
+ ip(2) = ip(0) + 2
+
+
+
+
+
+ !distance de la particule à remailler au second point (des quatres utilisés pour le remaillage)
+ !----------------------------------------------------
+ xx = (posz(i) - real(ip(0),kind=8)*dz-zd)/dz !relatif
+
+
+
+
+ !conditions au bord
+ !------------------
+ !periodique:
+ do c=-2,2
+ ip(c)=mod(ip(c)+nz,nz)
+ end do
+
+ !calcul du limiteur pour cas de bloc centre avec xx>0.5
+ !--------------------------------------------------------
+ call calcul_phi_tvd(phi,1)
+ phi_mdemi(1)=phi(0)
+ phi_pdemi(1)=phi(1)
+
+ !les autres cas ...
+ !---------------------
+ call calcul_phi_tvd(phi,0)
+ phi_mdemi(0)=phi(0)
+ phi_pdemi(0)=phi(1)
+
+ phi_mdemi(2)=0.
+ phi_pdemi(2)=0.
+
+ select case (blocg(i))
+
+ case(0)
+ if (xx<=0.5) then
+
+ poidx(-1)=ca(xx,0)
+ if (test (i,1) ) then
+ poidx(0)=cb(xx,0,1)
+ poidx(1)=cc(xx,1)
+ else
+ poidx(0)=cb(xx,0,0)
+ poidx(1)=cc(xx,0)
+ end if
+
+ else
+
+ poidx(0)=ca(xx-1.,1)
+ if (test (i,1) ) then
+ poidx(1)=cb(xx-1.,1,1)
+ poidx(2)=cc(xx-1.,1)
+ else
+ poidx(1)=cb(xx-1.,1,0)
+ poidx(2)=cc(xx-1.,0)
+ end if
+
+ end if
+
+ case(1)
+
+ poidx(-1)=ca(xx,0)
+ if (test (i,1) ) then
+ poidx(0)=cb(xx,0,1)
+ poidx(1)=cc(xx,1)
+ else
+ poidx(0)=cb(xx,0,0)
+ poidx(1)=cc(xx,0)
+ end if
+
+
+ case(2)
+ poidx(0)=ca(xx-1.,1)
+ poidx(1)=cb(xx-1.,1,0)+cc(xx-1.,0)
+ case(3)
+ poidx(-1)=ca(xx,0)
+ poidx(0)=cb(xx,0)+cc(xx,0)
+ case(4)
+ poidx(0)=ca(xx,1)+cb(xx,1,0)
+ poidx(1)=cc(xx,0)
+
+
+
+ case(5)
+ poidx(-1)=ca(xx,0)
+ poidx(0)=cb(xx,0)
+ poidx(1)=cc(xx,0)-cc(xx-1.,1)
+ poidx(2)=cc(xx-1.,1)
+
+ case(6)
+ poidx(-1)=ca(xx,0)
+ poidx(0)=ca(xx-1.,1)-ca(xx,0)
+ poidx(1)=cb(xx-1.,1)
+ poidx(2)=cc(xx-1.,1)
+
+
+ case(7)
+ poidx(-2)=ca(xx+1.,0)
+ poidx(-1)=ca(xx,0)-ca(xx+1.,0)
+ poidx(0)=cb(xx,0)
+ poidx(1)=cc(xx,0)
+ end select
+
+ !remaillage à l' interrieur domaine
+ !---------------------------------
+ do c=-2,2
+ remaille(ivar,jvar,ip(c)+1)=remaille(ivar,jvar,ip(c)+1)+donne(i)*poidx(c)
+ end do
+
+ end do
+
+ contains
+
+ function coeff_l2g (num,x) result(c)
+ implicit none
+ integer,intent(in) :: num
+ real(kind=8),intent(in) :: x
+ real(kind=8) :: c
+
+ select case(num)
+ case(-1)
+ c=0.5*x*(x-1.)
+ case(0)
+ c=1.-x**2
+ case(1)
+ c=0.5*x*(1.+x)
+ end select
+
+ end function coeff_l2g
+
+ function coeff_tscg (num,x) result(c)
+ implicit none
+ integer,intent(in) :: num
+ real(kind=8),intent(in) :: x
+ real(kind=8) :: c
+
+ select case(num)
+ case(-1)
+ c=0.5*x*(x-1.) +1./8.
+ case(0)
+ c=-x**2+1. -2./8.
+ case(1)
+ c=0.5*x*(1.+x) +1./8.
+ end select
+
+ end function coeff_tscg
+
+ function ca(x,num)
+ integer :: num
+ real(kind=8) :: x
+ real(kind=8) :: ca
+
+ ca=coeff_tscg(-1,x)+phi_mdemi(num)*(coeff_l2g(-1,x)-coeff_tscg(-1,x))
+
+ end function ca
+
+ function cb(x,num,num2)
+ integer :: num
+ real(kind=8) :: x
+ integer,optional :: num2
+ real(kind=8) :: cb
+
+ if (present (num2)) then
+ cb=1.-(coeff_tscg(-1,x)+coeff_tscg(1,x))-phi_pdemi(num2)*(coeff_l2g(1,x)-coeff_tscg(1,x))-phi_mdemi(num)*(coeff_l2g(-1,x)-coeff_tscg(-1,x))
+ else
+ cb=1.-(coeff_tscg(-1,x)+coeff_tscg(1,x))-phi_pdemi(num)*(coeff_l2g(1,x)-coeff_tscg(1,x))-phi_mdemi(num)*(coeff_l2g(-1,x)-coeff_tscg(-1,x))
+ end if
+
+ end function cb
+
+ function cc(x,num)
+ integer :: num
+ real(kind=8) :: x
+ real(kind=8) :: cc
+
+ cc=coeff_tscg(1,x)+phi_pdemi(num)*(coeff_l2g(1,x)-coeff_tscg(1,x))
+
+
+ end function cc
+
+ function test (ipart,vois)
+ implicit none
+ integer :: ipart,vois
+ logical :: test
+ real(kind=8) :: x
+
+ test=.false.
+ x=(posz(mod(ipart-1+vois+npart,npart)+1) - real( floor((posz(mod(ipart-1+vois+npart,npart)+1)-zd)/dz) ,kind=8)*dz-zd)/dz
+ if ( ((blocg(mod(ipart-1+vois+npart,npart)+1)==0).or.(blocg(mod(ipart-1+vois+npart,npart)+1)==2)).and.(x>0.5)) test=.true.
+
+ end function test
+
+
+
+ subroutine calcul_phi_tvd(phi_out,num)
+ implicit none
+ integer :: num
+ real(kind=8),dimension(0:1),intent(out):: phi_out
+ real(kind=8),dimension(-5:5) :: u
+ real(kind=8) :: rp,rm,diff1,diff2,diff3,diff0,rmm,y,r
+ integer :: ib,jb
+ real(kind=8),dimension(-5:5) :: xx_vois,x,tp1,tp2,tp3
+
+ !particules voisines
+ !-------------------
+ do ib=-2,2
+ u(ib)=donne(mod(i-1+ib+npart,npart)+1)
+ end do
+
+ do ib=-1,2
+ x(ib)=(posz(mod(i-1+ib+npart,npart)+1) - real( floor((posz(mod(i-1+ib+npart,npart)+1)-zd)/dz) ,kind=8)*dz-zd)/dz
+ end do
+
+
+ !cas bloc centre xx>0.5 (psi)
+ !-----------------------
+ if (num==1) then
+
+ do ib=0,1
+ x(ib)=x(ib)-1.
+ tp1(ib)=6.-8.*x(ib)*x(ib)
+ tp3(ib)=4.*(x(ib)-0.5)**2
+ end do
+
+ do ib=0,1
+ r=(u(ib+1)-u(ib))/(u(ib)-u(ib-1))
+
+ !anti diffusif
+ !phi_out(ib)=max(0.,min(tp3(ib),tp1(ib)*r))
+
+ !du type Minmod (attention non justifié pour y>=0.8)
+ phi_out(ib)=max(0.,min(1.,4.*r))
+
+ !bornes dépendent de la condition cfl
+ !phi_out(ib)=max(0.,min(1.,tp3(ib),tp1(ib)*r))
+
+ if(abs(u(ib)-u(ib-1))<0.0000001) phi_out(ib)=0.
+ end do
+
+
+
+
+ else
+ !sinon (phi)
+ !-----------------------
+
+
+ do ib=-1,0
+ tp1(ib)=6.-8.*x(ib)*x(ib)
+ tp2(ib)=4.*(x(ib)+0.5)**2
+ end do
+
+
+ do ib=0,1
+ r=(u(ib-1)-u(ib-2))/(u(ib)-u(ib-1))
+ !phi_out(ib)=max(0.,min(tp2(ib-1),tp1(ib-1)*r))
+ phi_out(ib)=max(0.,min(1.,4.*r))
+ !phi_out(ib)=max(0.,min(1.,tp2(ib-1),tp1(ib-1)*r))
+
+ if(abs(u(ib)-u(ib-1))<0.0000001) phi_out(ib)=0.
+ end do
+
+ end if
+
+
+ end subroutine calcul_phi_tvd
+
+
+ end subroutine limit_l2_bloc_z
+
diff --git a/Remesh/split_3d_rapide/utile_mod.f90 b/Remesh/split_3d_rapide/utile_mod.f90
index 1a0fa5f40f74c5b875b43680f5ae5a7db6b3fa2a..0ad2996d2871a2431e486dbc2180b9e323f42427 100644
--- a/Remesh/split_3d_rapide/utile_mod.f90
+++ b/Remesh/split_3d_rapide/utile_mod.f90
@@ -687,33 +687,33 @@ subroutine make_bloc_v2(esp)
if (blocd(partfm)==0) then
li=dt/pas*vit(partfm)
if ( (li<n) ) then
- blocd(partfm)=7
+ blocd(partfm)=2
else
- blocd(partfm)=6
+ blocd(partfm)=3
end if
else
- blocd(partfm)=6
+ blocd(partfm)=3
end if
if (blocg(partf)==0) then
li=dt/pas*vit(partf)
if ( (li<n) ) then
- blocg(partf)=8
- blocd(partf)=9
+ blocg(partf)=2
+ blocd(partf)=4
else
- blocg(partf)=7
- blocd(partf)=8
+ blocg(partf)=3
+ blocd(partf)=5
end if
else
- blocg(partf)=7
- blocd(partf)=8
+ blocg(partf)=3
+ blocd(partf)=5
end if
- blocg(partfp)=9
- blocd(partfp)=10
+ blocg(partfp)=4
+ blocd(partfp)=1
li=dt/pas*vit(partfpp)
if ( (blocg(partfpp)==0).and.(li<np) ) then
- blocg(partfpp)=11
+ blocg(partfpp)=5
else
- blocg(partfpp)=10
+ blocg(partfpp)=6
end if
end if
!N_{m+1}=N_m +1
@@ -721,34 +721,34 @@ subroutine make_bloc_v2(esp)
if ((np) == (n+1)) then
li=dt/pas*vit(partfm)
if ( (blocd(partfm)==0).and.(li<n) ) then
- blocd(partfm)=5
+ blocd(partfm)=6
else
- blocd(partfm)=4
+ blocd(partfm)=7
end if
blocg(partf)=1
- blocd(partf)=3
+ blocd(partf)=8
if (blocg(partfp)==0) then
li=dt/pas*vit(partfp)
if ( (li<np) ) then
- blocg(partfp)=4
- blocd(partfp)=2
+ blocg(partfp)=7
+ blocd(partfp)=0
else
- blocg(partfp)=3
+ blocg(partfp)=8
blocd(partfp)=1
end if
else
- blocg(partfp)=3
+ blocg(partfp)=8
blocd(partfp)=1
end if
if (blocg(partfpp)==0) then
li=dt/pas*vit(partfpp)
if ( (li<np) ) then
- blocg(partfpp)=6
+ blocg(partfpp)=9
else
- blocg(partfpp)=5
+ blocg(partfpp)=10
end if
else
- blocg(partfpp)=5
+ blocg(partfpp)=10
end if
end if