40 un_m1, un, pn_m1, pn, phin_m1, phin)
43 REAL(KIND=8),
INTENT(OUT):: time
44 REAL(KIND=8),
INTENT(IN) :: dt
45 INTEGER,
DIMENSION(:),
INTENT(IN) :: list_mode
46 REAL(KIND=8),
DIMENSION(:,:,:),
INTENT(OUT):: un_m1, un
47 REAL(KIND=8),
DIMENSION(:,:,:),
INTENT(OUT):: pn_m1, pn, phin_m1, phin
49 REAL(KIND=8),
DIMENSION(mesh_c%np) :: pn_m2
52 DO i= 1,
SIZE(list_mode)
56 un_m1(:,j,i) =
vv_exact(j,mesh_f%rr,mode,time-dt)
57 un(:,j,i) =
vv_exact(j,mesh_f%rr,mode,time)
61 pn_m2(:) =
pp_exact(j,mesh_c%rr,mode,time-2*dt)
62 pn_m1(:,j,i) =
pp_exact(j,mesh_c%rr,mode,time-dt)
63 pn(:,j,i) =
pp_exact(j,mesh_c%rr,mode,time)
64 phin_m1(:,j,i) = pn_m1(:,j,i) - pn_m2(:)
65 phin(:,j,i) = pn(:,j,i) - pn_m1(:,j,i)
118 INTEGER ,
INTENT(IN) :: TYPE
119 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: rr
120 INTEGER ,
INTENT(IN) :: mode, i
121 REAL(KIND=8),
INTENT(IN) :: time
122 REAL(KIND=8),
INTENT(IN) :: Re
123 CHARACTER(LEN=2),
INTENT(IN) :: ty
124 REAL(KIND=8),
DIMENSION(:,:,:),
OPTIONAL,
INTENT(IN) :: opt_density
125 REAL(KIND=8),
DIMENSION(:,:,:),
OPTIONAL,
INTENT(IN) :: opt_tempn
126 REAL(KIND=8),
DIMENSION(SIZE(rr,2)) :: vv
127 REAL(KIND=8),
DIMENSION(SIZE(rr,2)) :: LAP, TEMP,GPRESS
128 REAL(KIND=8),
DIMENSION(SIZE(rr,2)) :: r, z
129 REAL(KIND=8),
DIMENSION(SIZE(rr,2),6) :: Rot, V, vt
132 CHARACTER(LEN=2) :: np
134 IF (
PRESENT(opt_density))
CALL error_petsc(
'density should not be present for test 1')
135 IF (
PRESENT(opt_tempn))
CALL error_petsc(
'temperature should not be present for test 1')
141 IF (m==0 .OR. m==2)
THEN 144 rot(:,1) = m1*r*z**3*(4+m1) &
146 rot(:,2) =-m1*r*z**3*(4-m1)&
148 rot(:,3) = 3*m1*(r*z)**2-6*r**3*z &
150 rot(:,4) =-3*m1*(r*z)**2-6*r**3*z &
152 rot(:,5) = 12*(r*z)**2-9*r*z**3 + &
153 m1*(m1*r*z**3+3*(r*z)**2)
154 rot(:,6) = 12*(r*z)**2-9*r*z**3 + &
155 m1*(m1*r*z**3-3*(r*z)**2)
158 vt(:,1) = m1*z**3 - 3*r*z**2
159 vt(:,4) =3*(r-z)*z**2
161 v(:,1) = vt(:,1)*r**2
162 v(:,4) = vt(:,4)*r**2
163 v(:,5) = vt(:,5)*r**2
165 vt(:,2) =-m1*z**3 - 3*r*z**2
166 vt(:,3) =3*(r-z)*z**2
168 v(:,2) = vt(:,2)*r**2
169 v(:,3) = vt(:,3)*r**2
170 v(:,6) = vt(:,6)*r**2
174 IF (mod(
TYPE,2)==0) then
176 ELSEIF (
TYPE == 1) then
177 vv = -0.5d0*(v(:,3)*rot(:,5)-v(:,5)*rot(:,3) &
178 + v(:,4)*rot(:,6)-v(:,6)*rot(:,4))
179 ELSEIF (
TYPE == 3) then
180 vv = -0.5d0*(v(:,5)*rot(:,1)-v(:,1)*rot(:,5) &
181 + v(:,6)*rot(:,2)-v(:,2)*rot(:,6))
182 ELSEIF (
TYPE == 5) then
183 vv = -0.5d0*(v(:,1)*rot(:,3)-v(:,3)*rot(:,1) &
184 + v(:,2)*rot(:,4)-v(:,4)*rot(:,2))
191 vv = -0.5d0*(v(:,3)*rot(:,5)-v(:,5)*rot(:,3) &
192 - (v(:,4)*rot(:,6)-v(:,6)*rot(:,4)))
193 ELSE IF (
TYPE == 2) then
194 vv = -0.5d0*(v(:,3)*rot(:,6)-v(:,5)*rot(:,4) &
195 + (v(:,4)*rot(:,5)-v(:,6)*rot(:,3)))
196 ELSEIF (
TYPE == 3) then
197 vv = -0.5d0*(v(:,5)*rot(:,1)-v(:,1)*rot(:,5) &
198 - (v(:,6)*rot(:,2)-v(:,2)*rot(:,6)))
199 ELSEIF (
TYPE == 4) then
200 vv = -0.5d0*(v(:,5)*rot(:,2)-v(:,1)*rot(:,6) &
201 + (v(:,6)*rot(:,1)-v(:,2)*rot(:,5)))
202 ELSEIF (
TYPE == 5) then
203 vv = -0.5d0*(v(:,1)*rot(:,3)-v(:,3)*rot(:,1) &
204 - (v(:,2)*rot(:,4)-v(:,4)*rot(:,2)))
205 ELSEIF (
TYPE == 6) then
206 vv = -0.5d0*(v(:,1)*rot(:,4)-v(:,3)*rot(:,2) &
207 + (v(:,2)*rot(:,3)-v(:,4)*rot(:,1)))
220 vt(:,1) = m*z**3 - 3*r*z**2
221 vt(:,4) =3*(r-z)*z**2
223 v(:,1) = vt(:,1)*r**2
224 v(:,4) = vt(:,4)*r**2
225 v(:,5) = vt(:,5)*r**2
227 vt(:,2) =-m*z**3 - 3*r*z**2
228 vt(:,3) =3*(r-z)*z**2
230 v(:,2) = vt(:,2)*r**2
231 v(:,3) = vt(:,3)*r**2
232 v(:,6) = vt(:,6)*r**2
235 lap(:) = m*4*z**3 - 27*r*z**2 + m*6*r**2*z - 6*r**3 &
236 -m**2*vt(:,1) - 2*m*vt(:,4) - vt(:,1)
238 gpress(:) = 2*r*z**3*cos(t)
239 ELSEIF (
TYPE == 2) then
240 lap(:) =-m*4*z**3 - 27*r*z**2 - m*6*r**2*z - 6*r**3 &
241 -m**2*vt(:,2) + 2*m*vt(:,3) - vt(:,2)
243 gpress(:) = 2*r*z**3*cos(t)
244 ELSEIF (
TYPE == 3) then
245 lap(:) = 3*(9*r*z**2 - 4*z**3 + 2*r**3- 6*r**2*z) &
246 -m**2*vt(:,3) + 2*m*vt(:,2) -vt(:,3)
248 gpress(:) = m*r*z**3*cos(t)
249 ELSEIF (
TYPE == 4) then
250 lap(:) = 3*(9*r*z**2 - 4*z**3 + 2*r**3- 6*r**2*z) &
251 -m**2*vt(:,4) -2*m*vt(:,1) -vt(:,4)
253 gpress(:) = -m*r*z**3*cos(t)
254 ELSEIF (
TYPE == 5) then
255 lap(:) = (4-m)*(4*z**3 + 6*r**2*z) -m**2*vt(:,5)
257 gpress(:) = 3*r**2*z**2*cos(t)
258 ELSEIF (
TYPE == 6) then
259 lap(:) = (4+m)*(4*z**3 + 6*r**2*z) -m**2*vt(:,6)
261 gpress(:) = 3*r**2*z**2*cos(t)
263 vv(:) = -sin(t)*temp + cos(t)*(-lap/re) + gpress
301 FUNCTION vv_exact(TYPE,rr,m,t) RESULT(vv)
303 INTEGER ,
INTENT(IN) :: type
304 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: rr
305 INTEGER,
INTENT(IN) :: m
306 REAL(KIND=8),
INTENT(IN) :: t
307 REAL(KIND=8),
DIMENSION(SIZE(rr,2)) :: vv
308 REAL(KIND=8),
DIMENSION(SIZE(rr,2)) :: r, z
319 vv(:) = m*r**2*z**3 - 3*r**3*z**2
320 ELSEIF (
TYPE == 2) then
321 vv(:) = -m*r**2*z**3 - 3*r**3*z**2
322 ELSEIF (
TYPE == 3) then
323 vv(:) = 3*r**3*z**2 - 3 *r**2*z**3
324 ELSEIF (
TYPE == 4) then
325 vv(:) = 3*r**3*z**2 - 3 *r**2*z**3
326 ELSEIF (
TYPE == 5) then
327 vv(:) = r**2*z**3*(4-m)
328 ELSEIF (
TYPE == 6) then
329 vv(:) = r**2*z**3*(4+m)
332 vv(:) = vv(:) * cos(t)
354 FUNCTION pp_exact(TYPE,rr,m,t) RESULT (vv)
356 INTEGER ,
INTENT(IN) :: TYPE
357 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: rr
358 INTEGER ,
INTENT(IN) :: m
359 REAL(KIND=8),
INTENT(IN) :: t
360 REAL(KIND=8),
DIMENSION(SIZE(rr,2)) :: vv
366 IF (type==1.OR.type==2)
THEN 367 vv(:) = rr(1,:)**2*rr(2,:)**3*cos(t)
369 CALL error_petsc(
'pp_exact: TYPE should be equal to 1 or 2')
real(kind=8) function, dimension(size(rr, 2)), public pp_exact(TYPE, rr, m, t)
subroutine, public init_velocity_pressure(mesh_f, mesh_c, time, dt, list_mode, un_m1, un, pn_m1, pn, phin_m1, phin)
subroutine error_petsc(string)
real(kind=8) function, dimension(size(rr, 2)), public source_in_ns_momentum(TYPE, rr, mode, i, time, Re, ty, opt_density, opt_tempn)
real(kind=8) function, dimension(size(rr, 2)), public vv_exact(TYPE, rr, m, t)