[Fortran 90] Expert en Fishpack?
Bonjour tout le monde,
Si quelqu'un peut m'indiquer des erreurs dans mon code fortran 90, cela sera bien bcp apprécié.
Mon code utilise Fishpack("sepeli") et Lapack(LU decomposition) pour résoudre une edp d'ordre 2 sur un domaine 2D (z,r) comme (x,y) dans le sens général.
Je suis presque sur que Lapack le fait bien, du coup je compare les resultats Fishpack avec ceux de Lapack.
L'edp exacte est : D²(u,r) - (1/r)*D(u,r) + D²(u,z) =-r*w
avec w=0 partout.
Les conditions limites sont (Dans un rectangle (a,b)*(c,d) :
u=0 si r=0 ;
u=V0*r² /2 si z=0 ;
u=V0*(d-c)² /2 si r=d ;
D(u,z)=0 si z=b.
Je le colle en bas (excuse-moi pour la longueur):
Vous pouvez négliger la partie "Lapack" et chercher d'abord des erreurs dans la bloc Fishpack.....
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
program valdag87
use fish
implicit none
integer, parameter :: pr=selected_real_kind(15,70) !! Double precision
!! Pour Paraview
real, dimension(,allocatable :: x, y
CHARACTER(len=500) :: car, ca, c1,c2,c3 , bla, blab
!! For Fishpack
real, dimension(
real, dimension(:,
TYPE (fishworkspace) :: w1
INTEGER :: m,n,m1,n1,i,j,mbdcnd,nbdcnd,idmn,intl,iorder,ierror
real :: a, b, c, d, alpha,beta, dum, pertrb
real :: af,bf,cf, df,ef,ff
external :: cofx, cofy
!! Lapack
real(pr),dimension(:,
real(pr),dimension(
real(pr),dimension(:,
INTEGER, dimension(
integer :: INFO, ii
!! Quantites physiques
real(pr), dimension(:,
real(pr) :: drr,dz, err2, temps
real(pr) :: V0, R0, rloc,z
V0 = 5.d0
R0 = 6.d0
! SET LIMITS ON REGION
a = 0.
b = real(6*R0)
c = 0.
d = real(3*R0)
!print*, precision(V0) !! = 15
!print*, precision(c) !! = 6
! SET GRID SIZE
n = 40 !! MODIFICATIONS
m = n*2
m1 = m + 1
n1 = n + 1
dz = (b - a)/FLOAT(m)
drr = (d - c)/FLOAT(n)
Allocate( BDC(m1),BDD(m1) , BDA(n1),BDB(n1) )
Allocate( GRHS(m1,n1), USOL(m1,n1), w(m,n) )
Allocate( USOL2(m1*n1))
Allocate( x(m1),y(n1) )
Allocate( Go(m1*n1,m1*n1),Go1(m1*n1,m1*n1) )
Allocate( Dr(m1*n1),Dr1(m1*n1),Dr3(m1*n1), IPIV(m1*n1) )
Allocate( Dr2(m1,n1) )
w = 0.d0
!!!!!!!!!!!! Bloc FISHPACK !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
DO j = 2, n1 ! SET SPECIFIED BOUNDARY CONDITIONS AT X=A
rloc = (j - 1)*drr
USOL(1,j) = real( V0*rloc**2/2.d0 )
END DO
DO i = 1, m1
z = dble(a) + (i - 1)*dz
usol(i,1) = 0. ! SET SPECIFIED BOUNDARY CONDITIONS AT Y=C,D
usol(i,n1) = real( V0*(d-c)**2/2 )
END DO
DO i = 2, m !! Interieur
CALL cofx (z, af, bf, cf)
DO j = 2, n
rloc = dble(c) + (j - 1)*drr
CALL cofy (rloc, df, ef, ff)
grhs(i,j) = -rloc*( w(i-1,j-1)+w(i-1,j)+w(i,j-1)+w(i,j) )/4.
END DO
END DO
! SET MIXED BOUNDARY CONDITIONS AT X=B
beta = 0.
DO j = 1, n1
bdb(j) = 0.
END DO
IDMN = m+1
IORDER = 2
INTL = 0
!mbdcnd = 2
!nbdcnd = 1
! SUBROUTINE SEPELI (INTL,IORDER,A,B,M,MBDCND,BDA,ALPHA,BDB,BETA, &
! C,D,N,NBDCND,BDC,GAMA ,BDD,XNU, &
! COFX,COFY,GRHS, USOL,IDMN,W,PERTRB,IERROR)
call SEPELI (intl,iorder, a,b, M, 2, BDA, dum, BDB, beta, &
c,d, N, 1, BDC, dum, BDD, dum, &
COFX,COFY, GRHS,USOL, IDMN,W1,PERTRB,IERROR )
print*,"Fishpack_IERROR=",IERROR
CALL fishfin (w1)
!stop
!///////// Bloc pour le calcul par Lapack ///////////
do ii = 1,(M+1)*(N+1)
Go(ii,
Dr(ii) = 0.
end do
do j=1,N+1
do i=1,M+1
rloc = (j-1)*drr
ii = (j-1)*(M+1)+i
if(j==1) then
Go(ii,ii) = 1.
Dr(ii) = 0.
else if(j==N+1) then
Go(ii,ii) = 1.
Dr(ii) = V0*(dble(d-c))**2/2.
else if(i==1) then
Go(ii,ii) = 1.
Dr(ii) = V0*rloc**2/2.
else if(i==M+1) then
Go(ii,ii) = 1.
Go(ii,ii-1) = -1.
Dr(ii) = 0.
else
Go(ii,ii) = -2./(drr**2) -2./(dz**2)
Go(ii,ii-1) = 1./(dz**2)
Go(ii,ii+1) = 1./(dz**2)
Go(ii,ii-(M+1)) = 1./(drr**2) +1./(2.*drr*rloc)
Go(ii,ii+(M+1)) = 1./(drr**2) -1./(2.*drr*rloc)
Dr(ii) = -rloc*( w(i-1,j-1)+w(i-1,j)+w(i,j-1)+w(i,j) )/4.
end if
end do
end do
Go1 = Go !! Membre gauche a conserver, matrice non-LU-décomposée
Dr1 = Dr !! Membre droite a conserver
call DGETRF( (N+1)*(M+1),(N+1)*(M+1), Go, (N+1)*(M+1), IPIV, INFO )
call DGETRS( 'N',(N+1)*(M+1), 1, Go, (N+1)*(M+1), IPIV, Dr, (N+1)*(M+1), INFO )
!///////////// Fin bloc Lapack //////////////////////
!//////// Validation 1 ////////////
err2 = sum( dabs(matmul(Go1,Dr) - Dr1) ) /size(Dr1)
print*, "1 - Erreur du Calcul Lapack :",err2
!//////// Validation 2 ////////////
err2 = 0.
do i=1,M+1
do j=1,N+1
ii = (j-1)*(M+1)+i
Dr2(i,j) = Dr(ii)
! err2 = Amax1( err2 , dabs(Dr2(i,j)-USOL(i,j)) )
end do
end do
err2 = sum( dabs(Dr2(:,
print*, "2 - Difference d'ordre 1 entre Lapack-Fishpack :", err2
!//////// Validation 3 ////////////
print*, "3 - Difference d'ordre 2 entre Lapack-Fishpack :", sqrt( sum( (Dr2(:,
!//////// Validation 4 ////////////
58 do i=1,M+1
do j=1,N+1
ii = (j-1)*(M+1)+i
Dr3(ii) = USOL(i,j)
end do
end do
err2 = sum( dabs(matmul(Go1,Dr3) - Dr1) ) /size(Dr1)
print*, "4 - Erreur Fishpack avec la matrice de Lapack :",err2
print*,"Reference:O(drr**2+dz**2)=",drr**2+dz**2
call cpu_time(temps)
print*, ""
print*,"N et M : ",N," et ", M
print*,"Temps de calcul:",temps,"s"
!! Ecriture de fichiers de sortie pour Paraview -----> Au moins, le dernier fichier s'affiche bien.
do i=1,M+1
x(i) = (i-1)*dz + a
end do
do j=1,N+1
y(j) = (j-1)*drr + c
end do
do j=1,N+1
do i=1,M+1
USOL2((j-1)*(M+1)+i) = dble(USOL(i,j))
end do
end do
Deallocate( BDC,BDD , BDA,BDB )
Deallocate( GRHS, USOL, w , USOL2,x,y)
Deallocate( Go,Go1,Dr,Dr1,Dr2,Dr3, IPIV )
end program
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
SUBROUTINE cofx(x, af, bf, cf)
implicit none
integer, parameter :: pr=selected_real_kind(15,70)
!-----------------------------------------------
! D u m m y A r g u m e n t s
!-----------------------------------------------
real(pr), INTENT(IN) :: x
real, INTENT(OUT) :: af, bf, cf
!-----------------------------------------------
! SET COEFFICIENTS IN THE X-DIRECTION.
af = 1.
bf = 0.
cf = 0.
END SUBROUTINE cofx
!
SUBROUTINE cofy(y, df, ef, ff)
implicit none
integer, parameter :: pr=selected_real_kind(15,70)
!-----------------------------------------------
! D u m m y A r g u m e n t s
!-----------------------------------------------
real(pr), INTENT(IN) :: y
real, INTENT(OUT) :: df, ef, ff
!-----------------------------------------------
! SET COEFFICIENTS IN Y DIRECTION
df = 1.
ef = -1./real(y)
ff = 0.
END SUBROUTINE cofy
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Répondre avec citation
Partager