質問編集履歴
3
画像を追加
title
CHANGED
|
File without changes
|
body
CHANGED
|
@@ -185,4 +185,8 @@
|
|
|
185
185
|
|
|
186
186
|
2020/04/22 追記
|
|
187
187
|
法線を解析的に求めた場合の画像．
|
|
188
|
-
|
|
188
|
+
|
|
189
|
+
|
|
190
|
+
2020/04/23 追記
|
|
191
|
+
正規化方法をコメント頂いたように直した場合
|
|
192
|
+
|
2
プログラムの更新
title
CHANGED
|
File without changes
|
body
CHANGED
|
@@ -14,33 +14,38 @@
|
|
|
14
14
|
|
|
15
15
|
```fortran
|
|
16
16
|
|
|
17
|
+
|
|
17
18
|
program main
|
|
18
19
|
! Determine some arrays
|
|
19
20
|
implicit none
|
|
20
|
-
integer,parameter :: i_max=
|
|
21
|
+
integer,parameter :: i_max=100, j_max=100
|
|
21
22
|
double precision :: x(0:i_max,0:j_max), y(0:i_max,0:j_max), z(0:i_max,0:j_max) ! this is the position
|
|
22
23
|
double precision :: v(3) ! this is the velocity (uniform flow)
|
|
23
|
-
double precision :: nx(0:i_max-1,0:j_max-1,0:k_max-1),ny(0:i_max-1,0:j_max-1,0:k_max-1),nz(0:i_max-1,0:j_max-1,0:k_max-1)
|
|
24
24
|
double precision :: cp(0:i_max,0:j_max)
|
|
25
25
|
integer :: i,j
|
|
26
26
|
|
|
27
27
|
call grid
|
|
28
28
|
call normal_and_Cp
|
|
29
29
|
call dataSave
|
|
30
|
+
print *, pi_generator()
|
|
30
31
|
stop
|
|
31
32
|
contains
|
|
32
33
|
|
|
34
|
+
double precision function pi_generator()
|
|
35
|
+
pi_generator = 4*atan(1.0)
|
|
36
|
+
return
|
|
37
|
+
end function
|
|
38
|
+
|
|
33
39
|
subroutine grid
|
|
34
40
|
! Generating grids
|
|
35
41
|
double precision :: r = 2 ! radius
|
|
36
42
|
double precision :: theta(0:i_max)
|
|
37
43
|
double precision :: phi(0:j_max)
|
|
38
44
|
do 1000 j=0,j_max
|
|
39
|
-
theta(j)=float(j)/float((j_max))*(
|
|
45
|
+
theta(j)=float(j)/float((j_max))*(pi_generator() )/2.0
|
|
40
|
-
print *, theta(j), j
|
|
41
46
|
1000 continue
|
|
42
47
|
do 1001 i=0,i_max
|
|
43
|
-
phi(i) = float(i)/float(i_max)*
|
|
48
|
+
phi(i) = float(i)/float(i_max)*(pi_generator() )*2
|
|
44
49
|
1001 continue
|
|
45
50
|
! 球座標系で，r方向の変化はないモデルなので，2変数で書けそうだ
|
|
46
51
|
do 100 i=0,i_max
|
|
@@ -49,6 +54,8 @@
|
|
|
49
54
|
y(i,j) = r*sin(phi(i))*sin(theta(j))
|
|
50
55
|
z(i,j) = r*cos(theta(j))
|
|
51
56
|
101 continue
|
|
57
|
+
! x(i,0)=x(i,j_max)
|
|
58
|
+
! y(i,0)=y(i,j_max)
|
|
52
59
|
100 continue
|
|
53
60
|
end subroutine
|
|
54
61
|
|
|
@@ -60,13 +67,13 @@
|
|
|
60
67
|
! 面と格子がずれることになるが，そこは一度妥協する
|
|
61
68
|
|
|
62
69
|
!
|
|
63
|
-
! o
|
|
70
|
+
! o(U)
|
|
64
71
|
! |
|
|
72
|
+
! (L) |
|
|
73
|
+
! o-------o-------o(R)
|
|
65
74
|
! |
|
|
66
|
-
!(L)o-------o-------o(R)
|
|
67
75
|
! |
|
|
68
|
-
! |
|
|
69
|
-
! o
|
|
76
|
+
! o(B)
|
|
70
77
|
|
|
71
78
|
double precision :: vec_R(3),vec_L(3), vec_U(3), vec_B(3), &
|
|
72
79
|
& cross_RU(3), cross_UL(3), cross_LB(3), cross_BR(3), &
|
|
@@ -94,7 +101,6 @@
|
|
|
94
101
|
vec_B(1) = x(i,j-1) - x(i,j)
|
|
95
102
|
vec_B(2) = y(i,j-1) - y(i,j)
|
|
96
103
|
vec_B(3) = z(i,j-1) - z(i,j)
|
|
97
|
-
print *, vec_B(2)
|
|
98
104
|
|
|
99
105
|
|
|
100
106
|
! print *, vec_B(3)
|
|
@@ -125,10 +131,16 @@
|
|
|
125
131
|
cross_avg(2) = cross_avg(2)/sqrt(cross_avg(2)*cross_avg(2)+1.0d-5)
|
|
126
132
|
cross_avg(3) = cross_avg(3)/sqrt(cross_avg(3)*cross_avg(3)+1.0d-5)
|
|
127
133
|
|
|
134
|
+
! 球の場合，normal=(x,y,z)でした．x**2+y**2+z**2=r**2より，f=x**2+y**2+z**2-r**2としてgradとればokでした...
|
|
135
|
+
! =====================================================================================
|
|
136
|
+
cross_avg(1)=0.5*x(i,j)/sqrt( (0.5*x(i,j) )**2+ (0.5*y(i,j) )**2+ (0.5*z(i,j) )**2+1.0d-5)
|
|
137
|
+
cross_avg(2)=0.5*y(i,j)/sqrt( (0.5*x(i,j) )**2+ (0.5*y(i,j) )**2+ (0.5*z(i,j) )**2+1.0d-5)
|
|
138
|
+
cross_avg(3)=0.5*z(i,j)/sqrt( (0.5*x(i,j) )**2+ (0.5*y(i,j) )**2+ (0.5*z(i,j) )**2+1.0d-5)
|
|
139
|
+
! ====================================================================================
|
|
128
140
|
v(1)=0
|
|
129
141
|
v(2)=0
|
|
130
142
|
v(3)=-100*7
|
|
131
|
-
if (dot_product(v,cross_avg)
|
|
143
|
+
if (dot_product(v,cross_avg) < 0) then
|
|
132
144
|
cp(i,j) = 2*((dot_product(v,cross_avg)) &
|
|
133
145
|
& /sqrt(v(1)*v(1) + v(2)*v(2) + v(3)*v(3)))**2
|
|
134
146
|
else
|
|
@@ -147,7 +159,7 @@
|
|
|
147
159
|
subroutine dataSave
|
|
148
160
|
character (40):: filename
|
|
149
161
|
integer :: fo2=12
|
|
150
|
-
write(filename, "(a,i5.5,a)") "Newtonian"
|
|
162
|
+
write(filename, "(a,i5.5,a)") "vtkdata/Newtonian"
|
|
151
163
|
open(fo2,file=filename)
|
|
152
164
|
do 700 i=1,i_max
|
|
153
165
|
do 701 j=1,j_max
|
1
頂いた回答をもとに計算しなおした画像を掲載
title
CHANGED
|
File without changes
|
body
CHANGED
|
@@ -169,4 +169,8 @@
|
|
|
169
169
|
plotはgnuplotを使用して
|
|
170
170
|
splot 'Newtonian' u 1:2:3:4 w pm3d
|
|
171
171
|
で書いています．
|
|
172
|
-
fortran90を使用しています
|
|
172
|
+
fortran90を使用しています
|
|
173
|
+
|
|
174
|
+
2020/04/22 追記
|
|
175
|
+
法線を解析的に求めた場合の画像．
|
|
176
|
+
|