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)
93 dt, list_mode, level_set_m1, level_set)
96 REAL(KIND=8),
INTENT(OUT):: time
97 REAL(KIND=8),
INTENT(IN) :: dt
98 INTEGER,
DIMENSION(:),
INTENT(IN) :: list_mode
99 REAL(KIND=8),
DIMENSION(:,:,:,:),
INTENT(OUT):: level_set, level_set_m1
100 INTEGER :: mode, i, j, n
103 DO i= 1,
SIZE(list_mode)
107 DO n = 1,
inputs%nb_fluid -1
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)) :: r, z
131 CHARACTER(LEN=2) :: np
132 REAL(KIND=8),
DIMENSION(SIZE(rr,2)) :: ft, fd, fnl, fp
133 REAL(KIND=8) :: rho1, rho2, eta1, eta2
135 IF (
PRESENT(opt_tempn))
CALL error_petsc(
'temperature should not be present for test 40')
149 IF (
inputs%if_level_set)
THEN 151 rho1=
inputs%density_fluid(1)
153 eta1=
inputs%dyna_visc_fluid(1)
154 eta2=
inputs%dyna_visc_fluid(2)-
inputs%dyna_visc_fluid(1)
157 IF (m==1 .AND. type==1)
THEN 158 ft = rho2/32.d0 * (sqrt(2.d0)*cos(sqrt(2.d0)*sin(t))*cos(t)**2 - sin(sqrt(2.d0)*sin(t))*sin(t)) * r**2
159 ELSE IF (m==1 .AND. type==2)
THEN 160 ft = -rho2/32.d0*sqrt(2.d0) * (sqrt(2.d0)*sin(sqrt(2.d0)*sin(t))*cos(t)**2 + cos(sqrt(2.d0)*sin(t))*sin(t)) * r**2
161 ELSE IF (m==2 .AND. type==2)
THEN 162 ft = - (rho2/8.d0*sin(t-z)*cos(t) + (rho1 + rho2/2.d0 + rho2/8.d0*cos(t-z))*sin(t)) * r/2.d0
163 ELSE IF (m==3 .AND. type==1)
THEN 164 ft = -rho2/32.d0 * (sqrt(2.d0)*cos(sqrt(2.d0)*sin(t))*cos(t)**2 - sin(sqrt(2.d0)*sin(t))*sin(t)) * r**2
165 ELSE IF (m==3 .AND. type==2)
THEN 166 ft = -rho2/32.d0*sqrt(2.d0) * (sqrt(2.d0)*sin(sqrt(2.d0)*sin(t))*cos(t)**2 + cos(sqrt(2.d0)*sin(t))*sin(t)) * r**2
167 ELSE IF (m==0 .AND. type==3)
THEN 168 ft = -3.d0 * (rho2/8.d0*sin(t-z)*cos(t) + (rho1 + rho2/2.d0 + rho2/8.d0*cos(t-z))*sin(t)) * r/2.d0
169 ELSE IF (m==1 .AND. type==3)
THEN 170 ft = -7.d0/32.d0*rho2*sqrt(2.d0) * (sqrt(2.d0)*sin(sqrt(2.d0)*sin(t))*cos(t)**2 +cos(sqrt(2.d0)*sin(t))*sin(t)) * r**2
171 ELSE IF (m==1 .AND. type==4)
THEN 172 ft = 5/32.d0*rho2 * (sqrt(2.d0)*cos(sqrt(2.d0)*sin(t))*cos(t)**2 - sin(sqrt(2.d0)*sin(t))*sin(t)) * r**2
173 ELSE IF (m==2 .AND. type==3)
THEN 174 ft = -(rho2/8.d0*sin(t-z)*cos(t) + (rho1 + rho2/2.d0 + rho2/8.d0*cos(t-z))*sin(t)) * r/2.d0
175 ELSE IF (m==3 .AND. type==3)
THEN 176 ft = -rho2/32.d0*sqrt(2.d0) * (sqrt(2.d0)*sin(sqrt(2.d0)*sin(t))*cos(t)**2 +cos(sqrt(2.d0)*sin(t))*sin(t)) * r**2
177 ELSE IF (m==3 .AND. type==4)
THEN 178 ft = rho2/32.d0 * (sqrt(2.d0)*cos(sqrt(2.d0)*sin(t))*cos(t)**2 - sin(sqrt(2.d0)*sin(t))*sin(t)) * r**2
179 ELSE IF (m ==0 .AND. type==5)
THEN 180 ft = -rho2/8.d0*sin(t-z)
181 ELSE IF (m==1 .AND. type==5)
THEN 182 ft = -rho2/4.d0*sin(sqrt(2.d0)*sin(t))*cos(t)*r
183 ELSE IF (m==1 .AND. type==6)
THEN 184 ft = rho2/8.d0*sqrt(2.d0)*cos(sqrt(2.d0)*sin(t))*cos(t)*r
190 IF (m==0 .AND. type==1)
THEN 191 fnl = -(2.d0*rho1 + rho2 + rho2/4.d0*cos(t-z)) * cos(t)**2 * r
192 ELSE IF (m==1 .AND. type==1)
THEN 193 fnl = -9.d0/32.d0*rho2*sqrt(2.d0) * cos(sqrt(2.d0)*sin(t)) * cos(t)**2 * r**2
194 ELSE IF (m==1 .AND. type==2)
THEN 195 fnl = -3.d0/16.d0*rho2 * sin(sqrt(2.d0)*sin(t)) * cos(t)**2 * r**2
196 ELSE IF (m==2 .AND. type==2)
THEN 197 fnl = rho2/16.d0 * sin(t-z)*cos(t) * r
198 ELSE IF (m==3 .AND. type==1)
THEN 199 fnl = rho2/32.d0*sqrt(2.d0) * cos(sqrt(2.d0)*sin(t)) * cos(t)**2 * r**2
200 ELSE IF (m==3 .AND. type==2)
THEN 201 fnl = rho2/16.d0 * sin(sqrt(2.d0)*sin(t)) * cos(t)**2 * r**2
202 ELSE IF (m==0 .AND. type==3)
THEN 203 fnl = 3.d0*rho2/16.d0 * sin(t-z) * cos(t) * r
204 ELSE IF (m==1 .AND. type==3)
THEN 205 fnl = 7.d0*rho2/16.d0 * sin(sqrt(2.d0)*sin(t)) * cos(t)**2 * r**2
206 ELSE IF (m==1 .AND. type==4)
THEN 207 fnl = -5.d0*rho2/32.d0*sqrt(2.d0) * cos(sqrt(2.d0)*sin(t)) * cos(t)**2 * r**2
208 ELSE IF (m==2 .AND. type==3)
THEN 209 fnl = rho2/16.d0 * sin(t-z) * cos(t) * r
210 ELSE IF (m==3 .AND. type==3)
THEN 211 fnl = rho2/16.d0 * sin(sqrt(2.d0)*sin(t)) * cos(t)**2 * r**2
212 ELSE IF (m==3 .AND. type==4)
THEN 213 fnl = -rho2/32.d0*sqrt(2.d0) * cos(sqrt(2.d0)*sin(t)) * cos(t)**2 * r**2
214 ELSE IF (m==0 .AND. type==5)
THEN 215 fnl = rho2/8.d0*sin(t-z)
216 ELSE IF (m==1 .AND. type==5)
THEN 217 fnl = rho2/4.0*sin(sqrt(2.d0)*sin(t))*cos(t)*r
218 ELSE IF (m==1 .AND. type==6)
THEN 219 fnl = -rho2/8.d0*sqrt(2.d0)*cos(sqrt(2.d0)*sin(t))*cos(t)*r
225 IF (m==1 .AND. type==1)
THEN 226 fd = -eta2/(8.d0*
inputs%Re)*sin(sqrt(2.d0)*sin(t))*cos(t)
227 ELSE IF (m==1 .AND. type==2)
THEN 228 fd = -eta2/(8.d0*
inputs%Re)*sqrt(2.d0)*cos(sqrt(2.d0)*sin(t))*cos(t)
229 ELSE IF (m==1 .AND. type==3)
THEN 230 fd = -eta2/(8.d0*
inputs%Re)*sqrt(2.d0)*cos(sqrt(2.d0)*sin(t))*cos(t)
231 ELSE IF (m==1 .AND. type==4)
THEN 232 fd = eta2/(8.d0*
inputs%Re)*sin(sqrt(2.d0)*sin(t))*cos(t)
237 IF (m==0 .AND. type==1)
THEN 238 fp = 2.d0*r*z**3*cos(t)
239 ELSE IF (m==1 .AND. type==2)
THEN 241 ELSE IF (m==2 .AND. type==1)
THEN 242 fp = z*sin(t-r) - r*z*cos(t-r)
243 ELSE IF (m==1 .AND. type==3)
THEN 245 ELSE IF (m==2 .AND. type==4)
THEN 246 fp = -2.d0*z*sin(t-r)
247 ELSE IF (m==0 .AND. type==5)
THEN 248 fp = 3.d0*r**2*z**2*cos(t)
249 ELSE IF (m==1 .AND. type==6)
THEN 251 ELSE IF (m==2 .AND. type==5)
THEN 258 vv = ft + fd + fnl+ fp
260 CALL error_petsc(
'Error in condlim: if_level_set should be true')
265 m=i; m=
SIZE(opt_density,1); np=ty
286 INTEGER ,
INTENT(IN) :: TYPE
287 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: rr
288 INTEGER ,
INTENT(IN) :: m, interface_nb
289 REAL(KIND=8),
INTENT(IN) :: t
290 REAL(KIND=8),
DIMENSION(SIZE(rr,2)) :: vv
298 n=type; n=
SIZE(rr,1); n=m; n=interface_nb; r=t
304 FUNCTION vv_exact(TYPE,rr,m,t) RESULT(vv)
306 INTEGER ,
INTENT(IN) :: TYPE
307 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: rr
308 INTEGER,
INTENT(IN) :: m
309 REAL(KIND=8),
INTENT(IN) :: t
310 REAL(KIND=8),
DIMENSION(SIZE(rr,2)) :: vv
311 REAL(KIND=8),
DIMENSION(SIZE(rr,2)) :: r, z
320 ELSE IF (type==5)
THEN 328 ELSE IF (type==3)
THEN 356 FUNCTION pp_exact(TYPE,rr,m,t) RESULT (vv)
358 INTEGER ,
INTENT(IN) :: TYPE
359 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: rr
360 INTEGER ,
INTENT(IN) :: m
361 REAL(KIND=8),
INTENT(IN) :: t
362 REAL(KIND=8),
DIMENSION(SIZE(rr,2)) :: vv
364 IF (m==0.AND.type==1)
THEN 365 vv(:) = rr(1,:)**2*rr(2,:)**3*cos(t)
366 ELSE IF (m==1 .AND. type==2)
THEN 367 vv = rr(1,:)**2*cos(t-rr(2,:))
368 ELSE IF (m==2 .AND. type==1)
THEN 369 vv = rr(1,:)*rr(2,:)*sin(t-rr(1,:))
393 INTEGER ,
INTENT(IN) :: TYPE
394 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: rr
395 INTEGER ,
INTENT(IN) :: m, interface_nb
396 REAL(KIND=8),
INTENT(IN) :: t
397 REAL(KIND=8),
DIMENSION(SIZE(rr,2)) :: vv
398 REAL(KIND=8),
DIMENSION(SIZE(rr,2)) :: r, z
403 IF (interface_nb==1)
THEN 404 IF (m==0 .AND. type==1)
THEN 405 vv = 0.5d0 + 0.125d0*cos(t-z)
406 ELSE IF (m==1 .AND. type==1)
THEN 407 vv = 0.125d0*sqrt(2.d0)*cos(sqrt(2.d0)*sin(t))*r
408 ELSE IF (m==1 .AND. type==2)
THEN 409 vv = 0.125d0*sin(sqrt(2.d0)*sin(t))*r
414 CALL error_petsc(
' BUG in level_set_exact, we should compute only 1 level set')
real(kind=8) function, dimension(size(rr, 2)), public level_set_exact(interface_nb, TYPE, rr, m, t)
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 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 vv_exact(TYPE, rr, m, t)
subroutine, public init_level_set(pp_mesh, time, dt, list_mode, level_set_m1, level_set)
real(kind=8) function, dimension(size(rr, 2)), public source_in_level_set(interface_nb, TYPE, rr, m, t)