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)
119 INTEGER ,
INTENT(IN) :: TYPE
120 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: rr
121 INTEGER ,
INTENT(IN) :: mode, i
122 REAL(KIND=8),
INTENT(IN) :: time
123 REAL(KIND=8),
INTENT(IN) :: Re
124 CHARACTER(LEN=2),
INTENT(IN) :: ty
125 REAL(KIND=8),
DIMENSION(:,:,:),
OPTIONAL,
INTENT(IN) :: opt_density
126 REAL(KIND=8),
DIMENSION(:,:,:),
OPTIONAL,
INTENT(IN) :: opt_tempn
127 REAL(KIND=8),
DIMENSION(SIZE(rr,2)) :: vv
128 REAL(KIND=8),
DIMENSION(SIZE(rr,2)) :: LAP, TEMP,GPRESS
129 REAL(KIND=8),
DIMENSION(SIZE(rr,2)) :: r, z
130 REAL(KIND=8),
DIMENSION(SIZE(rr,2),6) :: Rot, V, vt
133 CHARACTER(LEN=2) :: np
135 IF (
PRESENT(opt_density))
CALL error_petsc(
'density should not be present for test 2')
136 IF (
PRESENT(opt_tempn))
CALL error_petsc(
'temperature should not be present for test 2')
144 IF (m==0 .OR. m==2)
THEN 146 rot(:,1) = r*(4.d0*cos(k*z)+1)
149 rot(:,4) = r*(k**2*r**2*cos(k*z)-8.d0*cos(k*z)-2.d0)
150 rot(:,5) = r*(-8.d0-k*r*sin(k*z))
157 vt(:,2) = -r**2*(1.d0-k*r*sin(k*z))
159 vt(:,6) = r**2*(4.d0*cos(k*z)+1.d0)
163 IF (mod(
TYPE,2)==0) then
165 ELSEIF (
TYPE == 1) then
166 vv = -0.5d0*(v(:,3)*rot(:,5)-v(:,5)*rot(:,3) &
167 + v(:,4)*rot(:,6)-v(:,6)*rot(:,4))
168 ELSEIF (
TYPE == 3) then
169 vv = -0.5d0*(v(:,5)*rot(:,1)-v(:,1)*rot(:,5) &
170 + v(:,6)*rot(:,2)-v(:,2)*rot(:,6))
171 ELSEIF (
TYPE == 5) then
172 vv = -0.5d0*(v(:,1)*rot(:,3)-v(:,3)*rot(:,1) &
173 + v(:,2)*rot(:,4)-v(:,4)*rot(:,2))
180 vv = -0.5d0*(v(:,3)*rot(:,5)-v(:,5)*rot(:,3) &
181 - (v(:,4)*rot(:,6)-v(:,6)*rot(:,4)))
182 ELSE IF (
TYPE == 2) then
183 vv = -0.5d0*(v(:,3)*rot(:,6)-v(:,5)*rot(:,4) &
184 + (v(:,4)*rot(:,5)-v(:,6)*rot(:,3)))
185 ELSEIF (
TYPE == 3) then
186 vv = -0.5d0*(v(:,5)*rot(:,1)-v(:,1)*rot(:,5) &
187 - (v(:,6)*rot(:,2)-v(:,2)*rot(:,6)))
188 ELSEIF (
TYPE == 4) then
189 vv = -0.5d0*(v(:,5)*rot(:,2)-v(:,1)*rot(:,6) &
190 + (v(:,6)*rot(:,1)-v(:,2)*rot(:,5)))
191 ELSEIF (
TYPE == 5) then
192 vv = -0.5d0*(v(:,1)*rot(:,3)-v(:,3)*rot(:,1) &
193 - (v(:,2)*rot(:,4)-v(:,4)*rot(:,2)))
194 ELSEIF (
TYPE == 6) then
195 vv = -0.5d0*(v(:,1)*rot(:,4)-v(:,3)*rot(:,2) &
196 + (v(:,2)*rot(:,3)-v(:,4)*rot(:,1)))
211 vt(:,2) = -r**2*(1.d0-k*r*sin(k*z))
213 vt(:,6) = r**2*(4.d0*cos(k*z)+1.d0)
218 gpress(:) = 2*r*cos(k*z)
219 ELSEIF (
TYPE == 2) then
220 lap(:) = (-8+7*k*r*sin(k*z)-k**3*r**3*sin(k*z))
223 ELSEIF (
TYPE == 3) then
224 lap(:) = (-8+2*k*r*sin(k*z))
227 ELSEIF (
TYPE == 4) then
230 gpress(:) = -r*cos(k*z)
231 ELSEIF (
TYPE == 5) then
234 gpress(:) = -r**2*k*sin(k*z)
235 ELSEIF (
TYPE == 6) then
236 lap(:) = (3+12*cos(k*z)-4*k**2*r**2*cos(k*z))
241 vv(:) = -sin(t)*temp + cos(t)*(-lap/re) + gpress*cos(t)
279 FUNCTION vv_exact(TYPE,rr,m,t) RESULT(vv)
281 INTEGER ,
INTENT(IN) :: type
282 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: rr
283 INTEGER,
INTENT(IN) :: m
284 REAL(KIND=8),
INTENT(IN) :: t
285 REAL(KIND=8),
DIMENSION(SIZE(rr,2)) :: vv
287 REAL(KIND=8),
DIMENSION(SIZE(rr,2)) :: r, z
300 ELSEIF (
TYPE == 2 .AND. m /= 0)
THEN 301 vv(:) = -r**2*(1-k*r*sin(k*z))
302 ELSEIF (
TYPE == 3) then
304 ELSEIF (
TYPE == 4 .AND. m /= 0)
THEN 306 ELSEIF (
TYPE == 5) then
308 ELSEIF (
TYPE == 6 .AND. m /= 0)
THEN 309 vv(:) = r**2*(4*cos(k*z)+1)
312 vv(:) = vv(:) * cos(t)
333 FUNCTION pp_exact(TYPE,rr,m,t) RESULT (vv)
335 INTEGER ,
INTENT(IN) :: TYPE
336 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: rr
337 INTEGER ,
INTENT(IN) :: m
338 REAL(KIND=8),
INTENT(IN) :: t
339 REAL(KIND=8),
DIMENSION(SIZE(rr,2)) :: vv
349 vv(:)= rr(1,:)**2*cos(k*rr(2,:))*cos(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 error_petsc(string)
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)
real(kind=8) function, dimension(size(rr, 2)), public vv_exact(TYPE, rr, m, t)