41 un_m1, un, pn_m1, pn, phin_m1, phin)
44 REAL(KIND=8),
INTENT(OUT):: time
45 REAL(KIND=8),
INTENT(IN) :: dt
46 INTEGER,
DIMENSION(:),
INTENT(IN) :: list_mode
47 REAL(KIND=8),
DIMENSION(:,:,:),
INTENT(OUT):: un_m1, un
48 REAL(KIND=8),
DIMENSION(:,:,:),
INTENT(OUT):: pn_m1, pn, phin_m1, phin
50 REAL(KIND=8),
DIMENSION(mesh_c%np) :: pn_m2
53 DO i= 1,
SIZE(list_mode)
57 un_m1(:,j,i) =
vv_exact(j,mesh_f%rr,mode,time-dt)
58 un(:,j,i) =
vv_exact(j,mesh_f%rr,mode,time)
62 pn_m2(:) =
pp_exact(j,mesh_c%rr,mode,time-2*dt)
63 pn_m1(:,j,i) =
pp_exact(j,mesh_c%rr,mode,time-dt)
64 pn(:,j,i) =
pp_exact(j,mesh_c%rr,mode,time)
65 phin_m1(:,j,i) = pn_m1(:,j,i) - pn_m2(:)
66 phin(:,j,i) = pn(:,j,i) - pn_m1(:,j,i)
75 REAL(KIND=8),
INTENT(OUT):: time
76 REAL(KIND=8),
INTENT(IN) :: dt
77 INTEGER,
DIMENSION(:),
INTENT(IN) :: list_mode
78 REAL(KIND=8),
DIMENSION(:,:,:),
INTENT(OUT):: tempn_m1, tempn
82 DO i= 1,
SIZE(list_mode)
117 INTEGER ,
INTENT(IN) :: TYPE
118 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: rr
119 INTEGER ,
INTENT(IN) :: mode, i
120 REAL(KIND=8),
INTENT(IN) :: time
121 REAL(KIND=8),
INTENT(IN) :: Re
122 CHARACTER(LEN=2),
INTENT(IN) :: ty
123 REAL(KIND=8),
DIMENSION(:,:,:),
OPTIONAL,
INTENT(IN) :: opt_density
124 REAL(KIND=8),
DIMENSION(:,:,:),
OPTIONAL,
INTENT(IN) :: opt_tempn
125 REAL(KIND=8),
DIMENSION(SIZE(rr,2)) :: vv, r, z
126 REAL(KIND=8) :: alpha, r0 = 0.5d0
127 REAL(KIND=8) :: pi = 3.1415926535897932d0
128 CHARACTER(LEN=2) :: np
130 IF (
PRESENT(opt_density))
CALL error_petsc(
'density should not be present for test 32')
132 alpha =
inputs%gravity_coefficient
138 ELSE IF (type==6)
THEN 146 vv = vv -(4*pi*r*(4*pi**2*r**4 - 8*pi**2*r**3*r0 + r0**2 + r**2*(-3 + 4*pi**2*r0**2))*cos(time)*cos(2*pi*z) + &
147 (r - r0)*re*cos(time)**2*(14*r**3 - 5*r**2*r0 - 5*r*r0**2 + 2*r0**3 + &
148 (36*pi**2*r**5 - 84*pi**2*r**4*r0 + 5*r*r0**2 - 2*r0**3 + 2*r**3*(-7 + 30*pi**2*r0**2) + &
149 r**2*(5*r0 - 12*pi**2*r0**3))*cos(4*pi*z)) - &
150 4*pi*r**3*(r - r0)**2*re*cos(2*pi*z)*sin(time))/(2.*r**3*re)
151 ELSE IF (mode==1)
THEN 152 vv = vv -2*(2*pi*r*(2*pi**2*r**4 - r*r0 - 4*pi**2*r**3*r0 + r0**2 + r**2*(-1 + 2*pi**2*r0**2))*cos(time)*cos(2*pi*z) + &
153 ((3*r**2 - 4*r*r0 + r0**2)*re*cos(time)**2*(3*r**2 - r0**2 + &
154 (8*pi**2*r**4 - 16*pi**2*r**3*r0 + r0**2 + r**2*(-3 + 8*pi**2*r0**2))*cos(4*pi*z)))/2. - &
155 pi*r**3*(r - r0)**2*re*cos(2*pi*z)*sin(time))/(r**3*re)
156 ELSE IF (mode==2)
THEN 157 vv = vv -((r - r0)*cos(time)**2*(4*r**2 - r*r0 - r0**2 + (12*pi**2*r**4 - &
158 28*pi**2*r**3*r0 + r0**2 + 4*r**2*(-1.d0 + 5*pi**2*r0**2) + &
159 r*(r0 - 4*pi**2*r0**3))*cos(4*pi*z)))/(2.*r**2)
161 ELSE IF (type==2)
THEN 163 vv = vv + ((r - r0)**2*cos(time)*(-4*pi*r*cos(2*pi*z) + re*cos(time)*(4*pi**2*r**4 - r*r0 - 8*pi**2*r**3*r0 + r0**2 + &
164 r**2*(-3.d0 + 4*pi**2*r0**2) + (3*r**2 + r*r0 - r0**2)*cos(4*pi*z))))/(r**3*re)
165 ELSE IF (mode==2)
THEN 166 vv = vv + ((r - r0)**2*cos(time)**2*(4*pi**2*r**4 - r*r0 - 8*pi**2*r**3*r0 + r0**2 + &
167 r**2*(-3.d0 + 4*pi**2*r0**2) + (3*r**2 + r*r0 - r0**2)*cos(4*pi*z)))/(2.*r**3)
169 ELSE IF (type==3)
THEN 171 vv = vv + (2*pi*(-3*pi*r*(r - r0)**3*(3*r - r0)*re*cos(time)**2 + &
172 (4*pi**2*r**4 - 8*pi**2*r**3*r0 + r0**2 + r**2*(-3.d0 + 4*pi**2*r0**2))*cos(time)*cos(2*pi*z) - &
173 r**2*(r - r0)**2*re*cos(2*pi*z)*sin(time)))/(r**2*re)
174 ELSE IF (mode==1)
THEN 175 vv = vv + (4*pi*r*(2*pi**2*r**4 - r*r0 - 4*pi**2*r**3*r0 + r0**2 + &
176 r**2*(-1 + 2*pi**2*r0**2))*cos(time)*cos(2*pi*z) + &
177 ((r - r0)**3*(3*r - r0)*re*cos(time)**2*(-1.d0 - 16*pi**2*r**2 + cos(4*pi*z)))/2. - &
178 2*pi*r**3*(r - r0)**2*re*cos(2*pi*z)*sin(time))/(r**3*re)
179 ELSE IF (mode==2)
THEN 180 vv = vv + ((r - r0)**3*(3*r - r0)*cos(time)**2*(-1.d0 - 4*pi**2*r**2 + cos(4*pi*z)))/(2.*r**3)
182 ELSE IF (type==4)
THEN 184 vv = vv + ((r - r0)**2*cos(time)*(-8*pi*r*cos(2*pi*z) + re*cos(time)*((-3*r + r0)**2 + &
185 (8*pi**2*r**4 + 6*r*r0 - 16*pi**2*r**3*r0 - r0**2 + r**2*(-9.d0 + 8*pi**2*r0**2))*cos(4*pi*z))))/(2.*r**3*re)
186 ELSE IF (mode==2)
THEN 187 vv = vv + ((r - r0)**2*cos(time)**2*(2*r - r0 + (2*r*(-1.d0 + pi*(r - r0))*(1 + pi*(r - r0)) + r0)*cos(4*pi*z)))/r**2
189 ELSE IF (type==5)
THEN 191 vv = vv + ((4*(12*pi**2*r**4 - 16*pi**2*r**3*r0 - r0**2 + &
192 r**2*(-3.d0 + 4*pi**2*r0**2))*cos(time)*sin(2*pi*z) - &
193 4*r**2*(3*r**2 - 4*r*r0 + r0**2)*re*sin(time)*sin(2*pi*z) + &
194 pi*r*(r - r0)**2*(12*pi**2*r**4 - r*r0 - 24*pi**2*r**3*r0 + 2*r0**2 + &
195 4*r**2*(-1.d0 + 3*pi**2*r0**2))*re*4*cos(time)**2*sin(4*pi*z)))/(4.*r**3*re)
196 ELSE IF (mode==1)
THEN 197 vv = vv + (4*(-r0 + pi**2*r*(3*r**2 - 4*r*r0 + r0**2))*cos(time)*sin(2*pi*z) - &
198 r*(3*r**2 - 4*r*r0 + r0**2)*re*sin(time)*sin(2*pi*z) + &
199 pi*(r - r0)**2*(16*pi**2*r**4 - 2*r*r0 - 32*pi**2*r**3*r0 + 3*r0**2 + &
200 r**2*(-5.d0 + 16*pi**2*r0**2))*re*cos(time)**2*sin(4*pi*z))/(r**2*re)
201 ELSE IF (mode==2)
THEN 202 vv = vv + (pi*(r - r0)**2*(4*pi**2*r**4 - r*r0 - 8*pi**2*r**3*r0 + r0**2 + &
203 r**2*(-1.d0 + 4*pi**2*r0**2))*cos(time)**2*sin(4*pi*z))/r**2
205 ELSE IF (type==6)
THEN 207 vv = vv + (2*(2*pi**2*r*(r - r0)**2 - r0)*cos(time)*sin(2*pi*z) - r*(r - r0)**2*re*sin(time)*sin(2*pi*z) -&
208 4*pi*r*(r - r0)**3*re*cos(time)**2*sin(4*pi*z))/(r**2*re)
209 ELSE IF (mode==2)
THEN 210 vv = vv -2*pi*(r - r0)**3*cos(time)**2*sin(4*pi*z)/r
214 IF ((type==1).AND.(mode==1))
THEN 215 vv = vv + 3*r**2*cos(time)*sin(2*pi*z)
216 ELSE IF ((type==4).AND.(mode==1))
THEN 217 vv = vv - r**2*cos(time)*sin(2*pi*z)
218 ELSE IF ((type==5).AND.(mode==1))
THEN 219 vv = vv + 2*pi*r**3*cos(time)*cos(2*pi*z)
232 INTEGER ,
INTENT(IN) :: TYPE
233 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: rr
234 INTEGER ,
INTENT(IN) :: m
235 REAL(KIND=8),
INTENT(IN) :: t
236 REAL(KIND=8),
DIMENSION(SIZE(rr,2)) :: vv, r, z, kappa
238 REAL(KIND=8) :: r0 = 0.5d0
239 REAL(KIND=8) :: pi = 3.1415926535897932d0
246 kappa(i) =
inputs%temperature_diffusivity(1)
248 kappa(i) =
inputs%temperature_diffusivity(2)
254 vv = - (-2*((2*pi**2*r**4 + 9*r*r0 - 4*pi**2*r**3*r0 - 2*r0**2 + 2*r**2*(-4.d0 + pi**2*r0**2))*kappa*cos(t)) &
255 + r**2*(r - r0)**2*sin(t))*sin(2*pi*z)
257 vv = -(-((4*pi**2*r**4 + 16*r*r0 - 8*pi**2*r**3*r0 - 3*r0**2 + r**2*(-15.d0 + 4*pi**2*r0**2))*kappa*cos(t)) &
258 + r**2*(r - r0)**2*sin(t))*sin(2*pi*z)
270 vv(i) = vv(i) + (-3*pi*r(i)*(r(i) - r0)**4*cos(t)**2*sin(4*pi*z(i)))/2.d0
276 vv(i) = vv(i) - 2*pi*r(i)*(r(i) - r0)**4*cos(t)**2*sin(4*pi*z(i))
282 vv(i) = vv(i) - 0.5*pi*(r(i)*(r(i) - r0)**4*cos(t)**2*sin(4*pi*z(i)))
306 FUNCTION vv_exact(TYPE,rr,m,t) RESULT(vv)
308 INTEGER ,
INTENT(IN) :: TYPE
309 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: rr
310 INTEGER,
INTENT(IN) :: m
311 REAL(KIND=8),
INTENT(IN) :: t
312 REAL(KIND=8),
DIMENSION(SIZE(rr,2)) :: vv, r, z
313 REAL(KIND=8) :: r0 = 0.5d0
314 REAL(KIND=8) :: pi = 3.1415926535897932d0
320 IF ((m==0).OR.(m==1))
THEN 321 vv = -2*pi*(r-r0)**2*cos(t)*cos(2*pi*z)
325 ELSE IF (type==3)
THEN 326 IF ((m==0).OR.(m==1))
THEN 327 vv = 2*pi*(r-r0)**2*cos(t)*cos(2*pi*z)
331 ELSE IF (type==5)
THEN 332 IF ((m==0).OR.(m==1))
THEN 333 vv = (r-r0)*cos(t)*sin(2*pi*z)/r * (3*r-r0)
337 ELSE IF (type==6)
THEN 339 vv = (r-r0)*cos(t)*sin(2*pi*z)/r * (r-r0)
367 FUNCTION pp_exact(TYPE,rr,m,t) RESULT (vv)
369 INTEGER ,
INTENT(IN) :: TYPE
370 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: rr
371 INTEGER ,
INTENT(IN) :: m
372 REAL(KIND=8),
INTENT(IN) :: t
373 REAL(KIND=8),
DIMENSION(SIZE(rr,2)) :: vv, r, z
374 REAL(KIND=8) :: pi = 3.1415926535897932d0
379 IF ((type==1).AND.(m==1))
THEN 380 vv = r**3*sin(2*pi*z)*cos(t)
391 INTEGER ,
INTENT(IN) :: TYPE
392 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: rr
393 INTEGER ,
INTENT(IN) :: m
394 REAL(KIND=8),
INTENT(IN) :: t
395 REAL(KIND=8),
DIMENSION(SIZE(rr,2)) :: vv, r, z
396 REAL(KIND=8) :: r0=0.5d0
397 REAL(KIND=8) :: pi = 3.1415926535897932d0
402 IF ((type==1).AND.((m==0).OR.(m==1)))
THEN 403 vv = r**2*(r-r0)**2*sin(2*pi*z)*cos(t)
450 INTEGER ,
INTENT(IN) ::
TYPE, n_start
451 INTEGER,
INTENT(IN) :: mode
452 REAL(KIND=8),
INTENT(IN) :: t
453 REAL(KIND=8),
DIMENSION(H_Mesh%np) :: vv
461 n=h_mesh%np; r=t; n=type; n=mode; n=n_start
real(kind=8) function, dimension(size(rr, 2)), public source_in_temperature(TYPE, rr, m, t)
real(kind=8) function, dimension(size(rr, 2)), public temperature_exact(TYPE, rr, m, t)
subroutine error_petsc(string)
subroutine, public init_velocity_pressure(mesh_f, mesh_c, time, dt, list_mode, un_m1, un, pn_m1, pn, phin_m1, phin)
real(kind=8) function, dimension(size(rr, 2)), public pp_exact(TYPE, rr, m, t)
real(kind=8) function, dimension(size(rr, 2)), public source_in_ns_momentum(TYPE, rr, mode, i, time, Re, ty, opt_density, opt_tempn)
subroutine, public init_temperature(mesh, time, dt, list_mode, tempn_m1, tempn)
real(kind=8) function, dimension(size(rr, 2)), public vv_exact(TYPE, rr, m, t)
real(kind=8) function, dimension(h_mesh%np), public extension_velocity(TYPE, H_mesh, mode, t, n_start)