前提・実現したいこと
fortran初学者です.
表題の通り,2次元ラプラス方程式を有限要素法によって近似し,2次元ポテンシャル流れを解析(速度ポテンシャルを求める)しようと考えています.質問としましては,作成したプログラムのエラーを自力で解決できないのでその点に関してご助言頂きたい所存で御座います.特に配列の宣言のルールなどは参考書をみながら手を進めてまいりましたがいささか怪しいです.
有限要素法の理論的な部分は理解し,計算の方法も分かるのですが,それをfortranプログラムで表現するところで躓いております.
そもそもfortranを学び始めたばかりで初歩的な部分かもしれませんが,次の様なプログラムを作成したところプログラム全体にわたって似たようなエラーが複数発生いたしました.
発生している問題・エラーメッセージ(一部)
189 | g(ik,j) = 0.0 | 1 Error: Unclassifiable statement at (1) 213 | a(k,k+1:n) = ar * a(k,k+1:n) | 1 217 | a(i, k+1 :n) = a(i, k+1 : n) - a(i, k) * a(k,k+1:n) | 1 Error: 'a' at (1) is not a variable 158 | g(n3,n3) = (b3*b3 + c3*c3) * S * g(n3,n3) | 1 Error: The function result on the lhs of the assignment at (1) must have the pointer attribute. 13 | call datain(nn,ne,ne1,nn2,x,y,np1,np2,np3,na,nb,qhata,qhatb,ifix,phat) !・ 1 Error: Type mismatch in argument 'x' at (1); passed REAL(8) to REAL(4)
該当のソースコード
fortran
1program FEM3 2 3implicit real(8)(a-h, o-z) 4real(8), dimension(1000,1000)::m1, m2 5dimension a(1000,1000) 6dimension g(1000,1000) 7dimension r(1000), x(1000), y(1000), qhata(1000), qhatb(1000), phat(1000) 8dimension b(1000), w(1000), ip(1000) 9dimension np1(1000), np2(1000), np3(1000) 10dimension na(1000), nb(1000), ifix(1000) 11dimension ug(1000), vg(1000), ru(1000), rv(1000) 12 13call datain(nn,ne,ne1,nn2,x,y,np1,np2,np3,na,nb,qhata,qhatb,ifix,phat) !入力 14 15call asmbly(nn,ne,ne1,ne2,nn2,x,y,np1,np2,np3,na,nb,qhata,qhatb) !行列GとベクトルRの組み立て 16 17call bound(ifix,g,r,phat) !ディリクレ型境界条件の導入 18 19call lineq(g, phat, r, n) !連立1次代数方程式の求解 20 21call output(phat,ug,vg) !出力 22 23 24end program FEM3 25 26!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 27 28subroutine datain(nn,ne,ne1,nn2,x,y,np1,np2,np3,na,nb,qhata,qhatb,ifix,phat) !入力 29 30open (11, file='NP.dat', status='old') 31!implicit real(8)(a-h, o-z) 32!dimension g(1000,1000) 33!dimension r(1000), x(1000), y(1000), qhata(1000), qhatb(1000), phat(1000) 34!dimension b(1000) 35!dimension np1(1000), np2(1000), np3(1000) 36!dimension na(1000), nb(1000), ifix(1000) 37 38 k = 0 39 read (11, '()') 40 do 41 read (11, *) np1, np2, np3 42 k = k + 1 43 end do 44 45 rewind (11) 46 read (11, '()') 47 do n = 1, k 48 read (11, *) k, np1, np2, np3 49 end do 50 51 close (11) 52 53 54 open(21,file='xy.dat',status='old') 55 56 i = 0 57 read (21, '()') 58 do 59 read (21, *) x, y 60 i = i + 1 61 end do 62 63 rewind (21) 64 read (21, '()') 65 do j = 1, i 66 read (21, *) i, x, y 67 end do 68 69 close(21) 70 71 72 open(31,file='k.dat',status='old') 73 74 k = 0 75 read (31, '()') 76 do 77 read (31, *) na,nb,qhata,qhatb 78 k = k + 1 79 end do 80 81 rewind (31) 82 read (31, '()') 83 do l = 1, k 84 read (31, *) k,na,nb,qhata,qhatb 85 end do 86 87 close(31) 88 89 90 open(41,file='I.dat',status='old') 91 92 i = 0 93 read (41, '()') 94 do 95 read (41, *) ifix,phat 96 i = i + 1 97 end do 98 99 rewind (41) 100 read (41, '()') 101 do m = 1, i 102 read (41, *) i,ifix,phat 103 end do 104 105 close(41) 106 107end subroutine datain 108 109 110subroutine asmbly(nn,ne,ne1,ne2,nn2,x,y,np1,np2,np3,na,nb,qhata,qhatb) !行列GとベクトルRの組み立て 111 112!implicit real(8)(a-h, o-z) 113!dimension g(1000,1000) 114!dimension r(1000), x(1000), y(1000), qhata(1000), qhatb(1000), phat(1000) 115!dimension b(1000) 116!dimension np1(1000), np2(1000), np3(1000) 117!dimension na(1000), nb(1000), ifix(1000) 118 119 do j=1,nn 120 do i=1,nn 121 g(i,j) = 0.0 122 end do 123 end do 124 125 do je=1, ne 126 127 n1 = np1(je) 128 n2 = np2(je) 129 n3 = np3(je) 130 131 x1 = x(n1) 132 x2 = x(n2) 133 x3 = x(n3) 134 135 y1 = y(n1) 136 y2 = y(n2) 137 y3 = y(n3) 138 139 S2 = x1*(y2-y3) + x2*(y3-y1) + x3*(y1-y2) 140 S = S2 * 0.5 141 S2 = 1.0 / S2 142 143 b1 = (y2-y3) * S2 !Bを作成 144 b2 = (y3-y1) * S2 145 b3 = (y1-y2) * S2 146 c1 = (x3-x2) * S2 !Cを作成 147 c2 = (x1-x3) * S2 148 c3 = (x2-x1) * S2 149 150 g(n1,n1) = (b1*b1 + c1*c1) * S * g(n1,n1) !Gを作成 151 g(n2,n1) = (b2*b1 + c2*c1) * S * g(n2,n1) 152 g(n3,n1) = (b3*b1 + c3*c1) * S * g(n3,n1) 153 g(n1,n2) = (b1*b2 + c1*c2) * S * g(n1,n2) 154 g(n2,n2) = (b2*b2 + c2*c2) * S * g(n2,n2) 155 g(n3,n2) = (b3*b2 + c3*c2) * S * g(n3,n2) 156 g(n1,n3) = (b1*b3 + c1*c3) * S * g(n1,n3) 157 g(n2,n3) = (b2*b3 + c2*c3) * S * g(n2,n3) 158 g(n3,n3) = (b3*b3 + c3*c3) * S * g(n3,n3) 159 160 end do 161 162 163 do i=1, nn 164 R(i) = 0.0 165 end do 166 167 do k=1, ne1 168 169 n1 = na(k) 170 n2 = nb(k) 171 dx = x(n1) - x(n2) 172 dy = y(n1) - y(n2) 173 dlk = sqrt(dx*dx + dy*dy) 174 qa = qhata(k) 175 qb = qhatb(k) 176 177 r(n1) = (2.0*qa + qb)*dlk / 6.0 + r(n1) 178 r(n2) = (qa + 2.0*qb)*dlk / 6.0 + r(n2) 179 180 end do 181 182end subroutine asmbly 183 184subroutine bound(ifix,g,r,phat) !ディリクレ型境界条件の導入 185 186 do k=1 , nn2 187 ik = ifix(k) 188 do j=1 , nn 189 g(ik,j) = 0.0 190 end do 191 g(ik,ik) = 1.0 192 r(ik) = phat(k) 193 end do 194 195end subroutine bound 196 197subroutine lineq(g, phat, r, n) !連立1次代数方程式の求解(ガウスージョルダン法) 198!implicit real(8)(a-h, o-z) 199!dimension g(1000,1000) 200!dimension r(1000), x(1000), y(1000), qhata(1000), qhatb(1000), phat(1000) 201!dimension b(1000) 202!dimension np1(1000), np2(1000), np3(1000) 203!dimension na(1000), nb(1000), ifix(1000) 204real(8) ar, a 205 206 a(:,:) = g(:,:) 207 x(:) = r(:) 208 209 do k=1, n 210 if(a(k,k) == 0.0d0) stop 211 ar = 1.0d0 / a(k,k) 212 a(k,k) = 1.0d0 213 a(k,k+1:n) = ar * a(k,k+1:n) 214 x(k) = ar * x(k) 215 do i = 1,n 216 if(i /= k) then 217 a(i, k+1 :n) = a(i, k+1 : n) - a(i, k) * a(k,k+1:n) 218 x(i) = x(i) 219 a(i,k) = 0.0d0 220 end if 221 end do 222 end do 223 224end subroutine lineq 225 226subroutine output(phat,k,nn) !出力 227 228!implicit real(8)(a-h, o-z) 229!dimention phat(1000) 230 231 open(51,file='output.dat',status='replace') 232 233 do k=1 , nn 234 write(51,*)'phat(k)=',phat(k) 235 end do 236 237 close(51) 238 239end subroutine output
試したこと
*サブルーチン内での宣言がよく分からなかったことと,他のエラーを解決するためにサブルーチン内での再宣言をコメント化しました.
*サブルーチンを使わないで考えてみました.
補足情報(FW/ツールのバージョンなど)
サブルーチン側での変数宣言は必要です。
"Error: Unclassifiable statement at (1)" というのもそのエラーです。
Fortran の学習はどのように進めているのでしょうか。
参考書を読んでいるのであれば、その内容をまず把握してみましょう。参考書の記述の意味が解らなければ内容について質問することもよいかと思います。
講義を受けているのであれば、指導教官に質問するとか、クラスメイトと相談するとかしてみてはいかがでしょうか。
あなたの回答
tips
プレビュー