Bonsoir,

J'ai débuté en Fortran 90 il y a 2 ou 3 jours et mon objectif est de comparer des temps de calcul entre des codes écrits en Matlab, C et F90. Pour un premier test j'ai réalisé un petit code d'éléments finis en élasticité linéaire qui consiste en gros à construire une matrice K, un vecteur F et de résoudre le système linéaire K*U = F.

Mon code F90 étant beaucoup plus lent que ses copains en C et Matlab, j'aimerais avoir un avis sur ce petit code pour savoir si je réalise des mauvaises opérations, mauvaises initialisations des matrices etc. Voici le code:

Code : Sélectionner tout - Visualiser dans une fenêtre à part
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
program linear_elasticity
 
!---------------------VARIABLES------------------------
implicit none
 
integer :: e, NE, N, ND, i, j, NLIN, NRHS, INFO
integer, parameter :: r=3, s=6
integer, dimension(3) :: nodes
integer, dimension(6) :: connect
integer, dimension(:, :), allocatable :: t
integer, dimension(:), allocatable :: fixeddofs
integer, dimension(:), allocatable :: IPIV
double precision, dimension(:, :), allocatable :: p, K
double precision, dimension(:), allocatable :: F
double precision, dimension(3, 6) :: Be, D
double precision, dimension(6, 3) :: BeT
double precision, dimension(6, 6) :: Ke
double precision, dimension(3, 3) :: L
double precision, dimension(3) :: x, y
double precision, dimension(:, :), allocatable :: IDENTITY
double precision :: A, young, nu, lambda, mu, t1, t2, ALPHA, BETA
 
!--------------------INSTRUCTIONS----------------------
call read_file(10, 'mesh.txt', NE, N, ND, p, t, fixeddofs)
 
call cpu_time (t1)
 
ALPHA = 1.0
BETA = 0.0
NLIN = 2*N
NRHS = 1
 
allocate (K(1:2*N,1:2*N))
allocate (F(1:2*N))
allocate (IDENTITY(1:ND, 1:ND))
allocate (IPIV(1:2*N))
 
IDENTITY = 0.0
DO i = 1,ND
   IDENTITY(i, i) = 1
ENDDO
 
young = 200e-2 !just for testing
nu = 0.3 
lambda = (young*nu)/((nu+1)*(1-2*nu))
mu = young/(2*(nu+1))
lambda = 2*lambda*mu/(lambda+2*mu)
L(1, :) = (/ lambda+2*mu, lambda, 0.0d0 /)
L(2, :) = (/ lambda, lambda+2*mu, 0.0d0 /)
L(3, :) = (/ 0.0d0, 0.0d0, mu /)
 
K = 0.0d0
 
DO e = 1,NE
 
   nodes = t(e, :)
   x = p(nodes, 1)
   y = p(nodes, 2)
 
   A = 0.5*((x(2)-x(1))*(y(3)-y(1)) - (y(2)-y(1))*(x(3)-x(1)))
 
   Be(1, :) = (/ (0.5/A)*(y(2)-y(3)), 0.0d0, (0.5/A)*(y(3)-y(1)), 0.0d0, (0.5/A)*(y(1)-y(2)), 0.0d0 /)
   Be(2, :) = (/ 0.0d0, (0.5/A)*(x(3)-x(2)), 0.0d0, (0.5/A)*(x(1)-x(3)), 0.0d0, (0.5/A)*(x(2)-x(1)) /)
   Be(3, :) = (/ (0.5/A)*(x(3)-x(2)), (0.5/A)*(y(2)-y(3)), (0.5/A)*(x(1)-x(3)), (0.5/A)*(y(3)-y(1)), &
        &        (0.5/A)*(x(2)-x(1)), (0.5/A)*(y(1)-y(2)) /)
 
   D = 0.0d0
   Ke = 0.0d0
 
   DO j = 1,r
      DO i = 1,s
         BeT(i, j) = Be(j, i)
      ENDDO
   ENDDO
 
   CALL DGEMM('N', 'N', r, s, r, ALPHA, L, r, Be, r, BETA, D, r)   
   CALL DGEMM('N', 'N', s, s, r, ALPHA, BeT, s, D, r, BETA, Ke, s) 
 
   connect = (/ 2*nodes(1)-1, 2*nodes(1), 2*nodes(2)-1, 2*nodes(2), 2*nodes(3)-1, 2*nodes(3) /)
 
   K(connect, connect) = K(connect, connect) + A*Ke
 
ENDDO
 
K(fixeddofs, :) = 0.0
K(fixeddofs, fixeddofs) = IDENTITY
F(fixeddofs(1:ND/2)) = p((1/2)*fixeddofs(1:ND/2)+1, 1)
 
CALL DGESV(NLIN, NRHS, K, NLIN, IPIV, F, NLIN, INFO);
 
call cpu_time (t2)
print*, t2-t1
 
if (INFO==0) then
write(*,*)'INFO=0:exit successful'
else
write(*,*)'bad exit, INFO=',INFO,'system not solved'
endif
 
!---------------------SUBROUTINES----------------------
contains
 
!-------------SUBROUTINE FOR READING A FILE------------
subroutine read_file (UnitNum, FileName, NE, N, ND, p, t, fixeddofs)
 
integer, intent (in) :: UnitNum
integer, intent (out) :: Ne, N, ND
integer, dimension(:, :), intent (out),  allocatable :: t
integer, dimension(:), intent (out), allocatable :: fixeddofs
integer :: i
character (len=*), intent (in) :: FileName
double precision, dimension(:, :), allocatable :: p
 
open(unit=UnitNum, file=FileName, status='old', action='read')
 
read(10,*) NE
read(10,*) N
read(10,*) ND
 
allocate (p(1:N,1:2))
allocate (t(1:NE,1:3))
allocate (fixeddofs(1:ND))
 
do i=1,N
read(10,fmt='(2F12.9)')p(i, 1),p(i, 2)
enddo
 
do i=1,NE
read(10,*)t(i, 1),t(i, 2),t(i, 3)
enddo
 
do i=1,ND
read(10,*)fixeddofs(i)
enddo
 
close(10)
 
end subroutine read_file
 
!---------------------END-SUBROUTINES------------------
 
end program linear_elasticity
Merci beaucoup et bonne soirée !