SFEMaNS  version 5.3
Reference documentation for SFEMaNS
condlim_test_28.f90
Go to the documentation of this file.
1 MODULE boundary_test_28
2  USE my_util
3  USE def_type_mesh
4  USE input_data
5  USE bessel
6 !!$ATTENTION
7 !!$Some subroutines have been commented to avoid warning messages when compiling executable.
8 !!$It can not be done in the module boundary_generic that expects all subroutines to be present.
9 !!$END ATTENTION
10  PUBLIC :: init_velocity_pressure
11 !!$ PUBLIC :: init_temperature
12 !!$ PUBLIC :: init_level_set
13  PUBLIC :: source_in_ns_momentum
14 !!$ PUBLIC :: source_in_temperature
15 !!$ PUBLIC :: source_in_level_set
16  PUBLIC :: vv_exact
17  PUBLIC :: imposed_velocity_by_penalty
18  PUBLIC :: pp_exact
19 !!$ PUBLIC :: temperature_exact
20 !!$ PUBLIC :: level_set_exact
21  PUBLIC :: penal_in_real_space
22 !!$ PUBLIC :: extension_velocity
23 !!$ PUBLIC :: Vexact
24 !!$ PUBLIC :: H_B_quasi_static
25 !!$ PUBLIC :: Hexact
26 !!$ PUBLIC :: Phiexact
27 !!$ PUBLIC :: Jexact_gauss
28 !!$ PUBLIC :: Eexact_gauss
29 !!$ PUBLIC :: init_maxwell
30 !!$ PUBLIC :: mu_bar_in_fourier_space
31 !!$ PUBLIC :: grad_mu_bar_in_fourier_space
32 !!$ PUBLIC :: mu_in_real_space
33  PRIVATE
34 
35  REAL(KIND=8), PRIVATE,PARAMETER :: pi = 3.14159265358979323846d0
36  REAL(KIND=8), PRIVATE ,PARAMETER :: twopi=2*pi
37  REAL(KIND=8), PARAMETER:: mu_disk= 5.0d1 !mu, Permeability of blades and disk
38 
39 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
40 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
41  !===TM73
42  !Parameters for discontinuous blades
43  !some offset to begin at the same vertical axis as the bottom propeller
44  REAL(KIND=8), PARAMETER:: top_propeller_angle_offset = 0.d0
45  INTEGER :: nblades = 8 !number of blades
46  REAL(KIND=8), PARAMETER:: lw = 0.0125 !lw= is the half thickness of the blades
47  REAL(KIND=8), PARAMETER:: disk_bot=-1.1d0, top_of_disk_bot=-0.9d0, top_of_blade_bot=-0.7d0
48  REAL(KIND=8), PARAMETER:: disk_top= 1.1d0, bot_of_disk_top= 0.9d0, bot_of_blade_top= 0.7d0
49  REAL(KIND=8), PARAMETER, PUBLIC:: disk_radius=0.75d0, hole_radius=0.1d0
50  !For straight blades use two_rp=0.d0
51  REAL(KIND=8), PARAMETER:: two_rp = disk_radius/sin(twopi*24.d0/360.d0)
52  !Volume of the cylinder (all domain)
53  REAL(KIND=8), PUBLIC :: omega_vol= pi*(1.0)*(2.0) !pi*r^2*h
54 
55  !Parameters for smooth_blades
56  REAL(KIND=8), PARAMETER:: wjump_hole = 0.06*(1.0), wjump= 0.04*(1.0)
57  REAL(KIND=8), PARAMETER:: hole_r = hole_radius, hole_rp=hole_r - wjump_hole
58  REAL(KIND=8), PUBLIC, PARAMETER:: disk_r = disk_radius, disk_rp= disk_r- wjump
59 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
60  !Bottom Smooth_propeller
61  REAL(KIND=8), PARAMETER:: cyl_bott = -1.0
62  REAL(KIND=8), PARAMETER:: bdisk_z = cyl_bott + 0.3d0, bdisk_z_p= bdisk_z - wjump
63  REAL(KIND=8), PARAMETER:: zbot = cyl_bott + 0.1d0, zbot_p = zbot - wjump
64  REAL(KIND=8), PARAMETER:: zbot_bar = zbot - 0.04d0, zbot_bar_p = zbot_bar - wjump
65 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
66  !Top Smooth_propeller
67  REAL(KIND=8), PARAMETER:: cyl_top = 1.0
68  REAL(KIND=8), PARAMETER:: tdisk_z = cyl_top - 0.3d0, tdisk_z_p= tdisk_z + wjump
69  REAL(KIND=8), PARAMETER:: ztop = cyl_top - 0.1d0, ztop_p = ztop + wjump
70  REAL(KIND=8), PARAMETER:: ztop_bar = ztop + 0.04d0, ztop_bar_p = ztop_bar + wjump
71 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
72  !Parameters for both smooth_blades
73  REAL(KIND=8), PARAMETER:: alpha=200.d0*(1.0), alpha_th = 80.d0*(1.0)
74 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
75  !Do we want both propellers?
76  LOGICAL,PARAMETER:: if_bottom_prop=.true.
77  LOGICAL,PARAMETER:: if_top_prop=.true.
78  REAL(KIND=8), PARAMETER:: solid_vel=1.0;
79 
80 
81 CONTAINS
82  !===============================================================================
83  ! Boundary conditions for Navier-Stokes
84  !===============================================================================
85 
86  !===Initialize velocity, pressure
87  SUBROUTINE init_velocity_pressure(mesh_f, mesh_c, time, dt, list_mode, &
88  un_m1, un, pn_m1, pn, phin_m1, phin)
89  IMPLICIT NONE
90  TYPE(mesh_type) :: mesh_f, mesh_c
91  REAL(KIND=8), INTENT(OUT):: time
92  REAL(KIND=8), INTENT(IN) :: dt
93  INTEGER, DIMENSION(:), INTENT(IN) :: list_mode
94  REAL(KIND=8), DIMENSION(:,:,:), INTENT(OUT):: un_m1, un
95  REAL(KIND=8), DIMENSION(:,:,:), INTENT(OUT):: pn_m1, pn, phin_m1, phin
96  INTEGER :: mode, i, j
97  REAL(KIND=8), DIMENSION(mesh_c%np) :: pn_m2
98 
99  time = 0.d0
100  DO i= 1, SIZE(list_mode)
101  mode = list_mode(i)
102  DO j = 1, 6
103  !===velocity
104  un_m1(:,j,i) = vv_exact(j,mesh_f%rr,mode,time-dt)
105  un(:,j,i) = vv_exact(j,mesh_f%rr,mode,time)
106  END DO
107  DO j = 1, 2
108  !===pressure
109  pn_m2(:) = pp_exact(j,mesh_c%rr,mode,time-2*dt)
110  pn_m1(:,j,i) = pp_exact(j,mesh_c%rr,mode,time-dt)
111  pn(:,j,i) = pp_exact(j,mesh_c%rr,mode,time)
112  phin_m1(:,j,i) = pn_m1(:,j,i) - pn_m2(:)
113  phin(:,j,i) = pn(:,j,i) - pn_m1(:,j,i)
114  ENDDO
115  ENDDO
116  END SUBROUTINE init_velocity_pressure
117 
118 !!$ !===Initialize temperature
119 !!$ SUBROUTINE init_temperature(mesh, time, dt, list_mode, tempn_m1, tempn)
120 !!$ IMPLICIT NONE
121 !!$ TYPE(mesh_type) :: mesh
122 !!$ REAL(KIND=8), INTENT(OUT):: time
123 !!$ REAL(KIND=8), INTENT(IN) :: dt
124 !!$ INTEGER, DIMENSION(:), INTENT(IN) :: list_mode
125 !!$ REAL(KIND=8), DIMENSION(:,:,:), INTENT(OUT):: tempn_m1, tempn
126 !!$ INTEGER :: mode, i, j
127 !!$
128 !!$ time = 0.d0
129 !!$ DO i= 1, SIZE(list_mode)
130 !!$ mode = list_mode(i)
131 !!$ DO j = 1, 2
132 !!$ tempn_m1(:,j,i) = temperature_exact(j, mesh%rr, mode, time-dt)
133 !!$ tempn (:,j,i) = temperature_exact(j, mesh%rr, mode, time)
134 !!$ ENDDO
135 !!$ ENDDO
136 !!$ END SUBROUTINE init_temperature
137 
138 !!$ !===Initialize level_set
139 !!$ SUBROUTINE init_level_set(vv_mesh, time, &
140 !!$ dt, list_mode, level_set_m1, level_set)
141 !!$ IMPLICIT NONE
142 !!$ TYPE(mesh_type) :: vv_mesh
143 !!$ REAL(KIND=8), INTENT(OUT):: time
144 !!$ REAL(KIND=8), INTENT(IN) :: dt
145 !!$ INTEGER, DIMENSION(:), INTENT(IN) :: list_mode
146 !!$ REAL(KIND=8), DIMENSION(:,:,:,:), INTENT(OUT):: level_set, level_set_m1
147 !!$ INTEGER :: mode, i, j, n
148 !!$
149 !!$ time = 0.d0
150 !!$ DO i= 1, SIZE(list_mode)
151 !!$ mode = list_mode(i)
152 !!$ DO j = 1, 2
153 !!$ !===level_set
154 !!$ DO n = 1, inputs%nb_fluid -1
155 !!$ level_set_m1(n,:,j,i) = level_set_exact(n,j,vv_mesh%rr,mode,time-dt)
156 !!$ level_set (n,:,j,i) = level_set_exact(n,j,vv_mesh%rr,mode,time)
157 !!$ END DO
158 !!$ END DO
159 !!$ END DO
160 !!$
161 !!$ END SUBROUTINE init_level_set
162 
163  !===Source in momemtum equation. Always called.
164  FUNCTION source_in_ns_momentum(TYPE, rr, mode, i, time, Re, ty, opt_density, opt_tempn) RESULT(vv)
165  IMPLICIT NONE
166  INTEGER , INTENT(IN) :: TYPE
167  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: rr
168  INTEGER , INTENT(IN) :: mode, i
169  REAL(KIND=8), INTENT(IN) :: time
170  REAL(KIND=8), INTENT(IN) :: Re
171  CHARACTER(LEN=2), INTENT(IN) :: ty
172  REAL(KIND=8), DIMENSION(:,:,:), OPTIONAL, INTENT(IN) :: opt_density
173  REAL(KIND=8), DIMENSION(:,:,:), OPTIONAL, INTENT(IN) :: opt_tempn
174  REAL(KIND=8), DIMENSION(SIZE(rr,2)) :: vv
175  REAL(KIND=8) :: r
176  INTEGER :: n
177  CHARACTER(LEN=2) :: np
178 
179  IF (PRESENT(opt_density)) CALL error_petsc('density should not be present for test 28')
180  IF (PRESENT(opt_tempn)) CALL error_petsc('temperature should not be present for test 28')
181 
182  vv = 0.d0
183  RETURN
184 
185  !===Dummies variables to avoid warning
186  n=type; n=SIZE(rr,1); n=mode; n=i; r=time; r=re;np=ty
187  !===Dummies variables to avoid warning
188  END FUNCTION source_in_ns_momentum
189 
190 !!$ !===Extra source in temperature equation. Always called.
191 !!$ FUNCTION source_in_temperature(TYPE, rr, m, t)RESULT(vv)
192 !!$ IMPLICIT NONE
193 !!$ INTEGER , INTENT(IN) :: TYPE
194 !!$ REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: rr
195 !!$ INTEGER , INTENT(IN) :: m
196 !!$ REAL(KIND=8), INTENT(IN) :: t
197 !!$ REAL(KIND=8), DIMENSION(SIZE(rr,2)) :: vv
198 !!$
199 !!$ vv = 0.d0
200 !!$ CALL error_petsc('source_in_temperature: should not be called for this test')
201 !!$ RETURN
202 !!$ END FUNCTION source_in_temperature
203 
204 !!$ !===Extra source in level set equation. Always called.
205 !!$ FUNCTION source_in_level_set(interface_nb,TYPE, rr, m, t)RESULT(vv)
206 !!$ IMPLICIT NONE
207 !!$ INTEGER , INTENT(IN) :: TYPE
208 !!$ REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: rr
209 !!$ INTEGER , INTENT(IN) :: m, interface_nb
210 !!$ REAL(KIND=8), INTENT(IN) :: t
211 !!$ REAL(KIND=8), DIMENSION(SIZE(rr,2)) :: vv
212 !!$
213 !!$ vv=0.d0
214 !!$ CALL error_petsc('sourece_in_temperature: should not be called for this test')
215 !!$ END FUNCTION source_in_level_set
216 
217  !===Velocity for boundary conditions in Navier-Stokes.
218  !===Can be used also to initialize velocity in: init_velocity_pressure_temperature
219  FUNCTION vv_exact(TYPE,rr,m,t) RESULT(vv)
220  IMPLICIT NONE
221  INTEGER , INTENT(IN) :: TYPE
222  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: rr
223  INTEGER, INTENT(IN) :: m
224  REAL(KIND=8), INTENT(IN) :: t
225  REAL(KIND=8), DIMENSION(SIZE(rr,2)) :: vv
226  REAL(KIND=8) :: r,z
227  INTEGER :: n
228 
229  vv=0.0
230 
231  IF (type==3 .AND. m==0) THEN
232  DO n = 1, SIZE(rr,2)
233  r= rr(1,n)
234  z= rr(2,n)
235  ! Are we in the Bottom propeller?
236  IF ( if_bottom_prop .AND. r <disk_radius .AND. z < top_of_blade_bot ) then
237  vv(n)=solid_vel*r
238  END IF
239 
240  !are we in the top Propeller?
241  IF ( if_top_prop .AND. r <disk_radius .AND. z > bot_of_blade_top) then
242  vv(n)=-solid_vel*r
243 
244  END IF
245  END DO
246  END IF
247  RETURN
248 
249  !===Dummies variables to avoid warning
250  r=t
251  !===Dummies variables to avoid warning
252  END FUNCTION vv_exact
253 
254  !===Solid velocity imposed when using penalty technique
255  !===Defined in Fourier space on mode 0 only.
256  FUNCTION imposed_velocity_by_penalty(rr,t) RESULT(vv)
257  IMPLICIT NONE
258  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: rr
259  REAL(KIND=8), INTENT(IN) :: t
260  REAL(KIND=8), DIMENSION(SIZE(rr,2),6) :: vv
261  REAL(KIND=8) :: r, z
262  INTEGER :: n
263 
264  vv=0.d0
265  DO n = 1, SIZE(rr,2)
266  r= rr(1,n)
267  z= rr(2,n)
268  IF (z<-0.5d0) THEN
269  vv(n,3) = solid_vel*rr(1,n)
270  ELSE
271  vv(n,3) = -solid_vel*rr(1,n)
272  ENDIF
273  END DO
274  RETURN
275 
276  !===Dummies variables to avoid warning
277  r=t
278  !===Dummies variables to avoid warning
279  END FUNCTION imposed_velocity_by_penalty
280 
281  !===Pressure for boundary conditions in Navier-Stokes.
282  !===Can be used also to initialize pressure in the subroutine init_velocity_pressure.
283  !===Use this routine for outflow BCs only.
284  !===CAUTION: Do not enfore BCs on pressure where normal component
285  ! of velocity is prescribed.
286  FUNCTION pp_exact(TYPE,rr,m,t) RESULT (vv)
287  IMPLICIT NONE
288  INTEGER , INTENT(IN) :: TYPE
289  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: rr
290  INTEGER , INTENT(IN) :: m
291  REAL(KIND=8), INTENT(IN) :: t
292  REAL(KIND=8), DIMENSION(SIZE(rr,2)) :: vv
293  REAL(KIND=8) :: r
294  INTEGER :: n
295 
296  vv=0.d0
297  RETURN
298 
299  !===Dummies variables to avoid warning
300  n=type; n=SIZE(rr,1); n=m; r=t
301  !===Dummies variables to avoid warning
302  END FUNCTION pp_exact
303 
304 !!$ !===Temperature for boundary conditions in temperature equation.
305 !!$ FUNCTION temperature_exact(TYPE,rr,m,t) RESULT (vv)
306 !!$ IMPLICIT NONE
307 !!$ INTEGER , INTENT(IN) :: TYPE
308 !!$ REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: rr
309 !!$ INTEGER , INTENT(IN) :: m
310 !!$ REAL(KIND=8), INTENT(IN) :: t
311 !!$ REAL(KIND=8), DIMENSION(SIZE(rr,2)) :: vv
312 !!$
313 !!$ vv = 0.d0
314 !!$ CALL error_petsc('temperature_exact: should not be called for this test')
315 !!$ RETURN
316 !!$ END FUNCTION temperature_exact
317 
318 !!$ !===Can be used to initialize level set in the subroutine init_level_set.
319 !!$ FUNCTION level_set_exact(interface_nb,TYPE,rr,m,t) RESULT (vv)
320 !!$ IMPLICIT NONE
321 !!$ INTEGER , INTENT(IN) :: TYPE
322 !!$ REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: rr
323 !!$ INTEGER , INTENT(IN) :: m, interface_nb
324 !!$ REAL(KIND=8), INTENT(IN) :: t
325 !!$ REAL(KIND=8), DIMENSION(SIZE(rr,2)) :: vv
326 !!$
327 !!$ vv = 0.d0
328 !!$ CALL error_petsc('level_set_exact: should not be called for this test')
329 !!$ RETURN
330 !!$
331 !!$ END FUNCTION level_set_exact
332 
333  !===Penalty coefficient (if needed)
334  !===This coefficient is equal to zero in subdomain
335  !===where penalty is applied (penalty is zero in solid)
336  FUNCTION penal_in_real_space(mesh,rr_gauss,angles,nb_angles,nb,ne,time) RESULT(vv)
337  IMPLICIT NONE
338  TYPE(mesh_type) :: mesh
339  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: rr_gauss
340  REAL(KIND=8), DIMENSION(:), INTENT(IN) :: angles
341  INTEGER, INTENT(IN) :: nb_angles
342  INTEGER, INTENT(IN) :: nb, ne
343  REAL(KIND=8), INTENT(IN) :: time
344  REAL(KIND=8), DIMENSION(nb_angles,ne-nb+1) :: vv
345 
346 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
347  !1) USE Smooth Blades
348  vv=smooth_penal_in_real_space(mesh,rr_gauss,angles,nb_angles,nb,ne,time)
349 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
350  RETURN
351  END FUNCTION penal_in_real_space
352 
353 !!$ !===Extension of the velocity field in the solid.
354 !!$ !===Used when temperature or Maxwell equations are solved.
355 !!$ !===It extends the velocity field on the Navier-Stokes domain to a
356 !!$ !===velocity field on the temperature and the Maxwell domain.
357 !!$ !===It is also used if problem type=mxw and restart velocity
358 !!$ !===is set to true in data (type problem denoted mxx in the code).
359 !!$ FUNCTION extension_velocity(TYPE, H_mesh, mode, t, n_start) RESULT(vv)
360 !!$ IMPLICIT NONE
361 !!$ TYPE(mesh_type), INTENT(IN) :: H_mesh
362 !!$ INTEGER , INTENT(IN) :: TYPE, n_start
363 !!$ INTEGER, INTENT(IN) :: mode
364 !!$ REAL(KIND=8), INTENT(IN) :: t
365 !!$ REAL(KIND=8), DIMENSION(H_Mesh%np) :: vv
366 !!$
367 !!$ vv = 0.d0
368 !!$ RETURN
369 !!$
370 !!$ END FUNCTION extension_velocity
371 
372  !===============================================================================
373  ! Boundary conditions for Maxwell
374  !===============================================================================
375 !!$ !===Velocity used in the induction equation.
376 !!$ !===Used only if problem type is mxw and restart velocity is false
377 !!$ FUNCTION Vexact(m, H_mesh) RESULT(vv) !Set uniquement a l'induction
378 !!$ IMPLICIT NONE
379 !!$ TYPE(mesh_type), INTENT(IN) :: H_mesh
380 !!$ INTEGER, INTENT(IN) :: m
381 !!$ REAL(KIND=8), DIMENSION(H_mesh%np,6) :: vv
382 !!$
383 !!$ vv = 0.d0
384 !!$ END FUNCTION Vexact
385 
386 !!$ !===Magnetic field and magnetic induction for quasi-static approximation
387 !!$ !===if needed
388 !!$ FUNCTION H_B_quasi_static(char_h_b, rr, m) RESULT(vv)
389 !!$ IMPLICIT NONE
390 !!$ CHARACTER(LEN=1), INTENT(IN) :: char_h_b
391 !!$ REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: rr
392 !!$ INTEGER, INTENT(IN) :: m
393 !!$ REAL(KIND=8), DIMENSION(SIZE(rr,2),6) :: vv
394 !!$
395 !!$ vv = 0.d0
396 !!$ RETURN
397 !!$ END FUNCTION H_B_quasi_static
398 
399 !!$ !===Magnetic field for boundary conditions in the Maxwell equations.
400 !!$ FUNCTION Hexact(H_mesh,TYPE, rr, m, mu_H_field, t) RESULT(vv)
401 !!$ IMPLICIT NONE
402 !!$ TYPE(mesh_type), INTENT(IN) :: H_mesh
403 !!$ INTEGER , INTENT(IN) :: TYPE
404 !!$ REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: rr
405 !!$ INTEGER , INTENT(IN) :: m
406 !!$ REAL(KIND=8), INTENT(IN) :: t
407 !!$ REAL(KIND=8), DIMENSION(:), INTENT(IN) :: mu_H_field
408 !!$ REAL(KIND=8), DIMENSION(SIZE(rr,2)) :: vv
409 !!$
410 !!$ vv = 0.d0
411 !!$ RETURN
412 !!$ END FUNCTION Hexact
413 
414 !!$ !===Scalar potential for boundary conditions in the Maxwell equations.
415 !!$ FUNCTION Phiexact(TYPE, rr, m, mu_phi,t) RESULT(vv)
416 !!$ IMPLICIT NONE
417 !!$ INTEGER , INTENT(IN) :: TYPE
418 !!$ REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: rr
419 !!$ INTEGER , INTENT(IN) :: m
420 !!$ REAL(KIND=8), INTENT(IN) :: mu_phi, t
421 !!$ REAL(KIND=8), DIMENSION(SIZE(rr,2)) :: vv
422 !!$
423 !!$ vv = 0.d0
424 !!$ RETURN
425 !!$ END FUNCTION Phiexact
426 
427 !!$ !===Current in Ohm's law. Curl(H) = sigma(E + uxB) + current
428 !!$ FUNCTION Jexact_gauss(TYPE, rr, m, mu_phi, sigma, mu_H, t, mesh_id, opt_B_ext) RESULT(vv)
429 !!$ IMPLICIT NONE
430 !!$ INTEGER , INTENT(IN) :: TYPE
431 !!$ REAL(KIND=8), DIMENSION(:), INTENT(IN) :: rr
432 !!$ INTEGER , INTENT(IN) :: m
433 !!$ REAL(KIND=8), INTENT(IN) :: mu_phi, sigma, mu_H, t
434 !!$ INTEGER , INTENT(IN) :: mesh_id
435 !!$ REAL(KIND=8), DIMENSION(6), OPTIONAL,INTENT(IN) :: opt_B_ext
436 !!$ REAL(KIND=8) :: vv
437 !!$
438 !!$ vv = 0.d0
439 !!$ RETURN
440 !!$ END FUNCTION Jexact_gauss
441 
442 !!$ !===Electric field for Neumann BC (cf. doc)
443 !!$ FUNCTION Eexact_gauss(TYPE, rr, m, mu_phi, sigma, mu_H, t) RESULT(vv)
444 !!$ IMPLICIT NONE
445 !!$ INTEGER, INTENT(IN) :: TYPE
446 !!$ REAL(KIND=8), DIMENSION(:), INTENT(IN) :: rr
447 !!$ INTEGER, INTENT(IN) :: m
448 !!$ REAL(KIND=8), INTENT(IN) :: mu_phi, sigma, mu_H, t
449 !!$ REAL(KIND=8) :: vv
450 !!$
451 !!$ vv = 0.d0
452 !!$ CALL error_petsc('Eexact: should not be called for this test')
453 !!$ END FUNCTION Eexact_gauss
454 
455 !!$ !===Initialization of magnetic field and scalar potential (if present)
456 !!$ SUBROUTINE init_maxwell(H_mesh, phi_mesh, time, dt, mu_H_field, mu_phi, &
457 !!$ list_mode, Hn1, Hn, phin1, phin)
458 !!$ IMPLICIT NONE
459 !!$ TYPE(mesh_type) :: H_mesh, phi_mesh
460 !!$ REAL(KIND=8), INTENT(OUT):: time
461 !!$ REAL(KIND=8), INTENT(IN) :: dt
462 !!$ REAL(KIND=8), DIMENSION(:), INTENT(IN) :: mu_H_field
463 !!$ REAL(KIND=8), INTENT(IN) :: mu_phi
464 !!$ INTEGER, DIMENSION(:), INTENT(IN) :: list_mode
465 !!$ REAL(KIND=8), DIMENSION(:,:,:), INTENT(OUT):: Hn, Hn1
466 !!$ REAL(KIND=8), DIMENSION(:,:,:), INTENT(OUT):: phin, phin1
467 !!$
468 !!$ Hn1 = 0.d0
469 !!$ Hn = 0.d0
470 !!$ phin1 = 0.d0
471 !!$ phin = 0.d0
472 !!$ time = 0.d0
473 !!$ RETURN
474 !!$ END SUBROUTINE init_maxwell
475 
476 !!$ !===Analytical permeability (if needed)
477 !!$ !===This function is not needed unless the flag
478 !!$ !=== ===Use FEM Interpolation for magnetic permeability (true/false)
479 !!$ !===is activated and set to .FALSE. in the data data file. Default is .TRUE.
480 !!$ FUNCTION mu_bar_in_fourier_space(H_mesh,nb,ne,pts,pts_ids) RESULT(vv)
481 !!$ IMPLICIT NONE
482 !!$ TYPE(mesh_type), INTENT(IN) :: H_mesh
483 !!$ REAL(KIND=8), DIMENSION(ne-nb+1) :: vv
484 !!$ INTEGER, INTENT(IN) :: nb, ne
485 !!$ REAL(KIND=8),DIMENSION(2,ne-nb+1),OPTIONAL :: pts
486 !!$ INTEGER, DIMENSION(ne-nb+1), OPTIONAL :: pts_ids
487 !!$ INTEGER :: n
488 !!$ REAL(KIND=8),DIMENSION(ne-nb+1) :: r,z
489 !!$
490 !!$ IF( PRESENT(pts) .AND. PRESENT(pts_ids) ) THEN !Computing mu at pts
491 !!$ r=pts(1,nb:ne)
492 !!$ z=pts(2,nb:ne)
493 !!$ ELSE
494 !!$ r=H_mesh%rr(1,nb:ne) !Computing mu at nodes
495 !!$ z=H_mesh%rr(2,nb:ne)
496 !!$ END IF
497 !!$
498 !!$ DO n = 1, ne - nb + 1
499 !!$ vv(n)=mu_bar_func(r(n),z(n))
500 !!$ END DO
501 !!$
502 !!$ RETURN
503 !!$ END FUNCTION mu_bar_in_fourier_space
504 
505 !!$ !===Analytical mu_in_fourier_space (if needed)
506 !!$ !===This function is not needed unless the flag
507 !!$ !=== ===Use FEM Interpolation for magnetic permeability (true/false)
508 !!$ !===is activated and set to .FALSE. in the data data file. Default is .TRUE.
509 !!$ FUNCTION grad_mu_bar_in_fourier_space(pt,pt_id) RESULT(vv)
510 !!$ IMPLICIT NONE
511 !!$ REAL(KIND=8),DIMENSION(2), INTENT(in):: pt
512 !!$ INTEGER,DIMENSION(1), INTENT(in) :: pt_id
513 !!$ REAL(KIND=8),DIMENSION(2) :: vv
514 !!$
515 !!$ vv=grad_mu_bar_func(pt(1),pt(2))
516 !!$
517 !!$ RETURN
518 !!$ END FUNCTION grad_mu_bar_in_fourier_space
519 
520 !!$ !===Analytical permeability, mu in real space (if needed)
521 !!$ FUNCTION mu_in_real_space(H_mesh,angles,nb_angles,nb,ne,time) RESULT(vv)
522 !!$ IMPLICIT NONE
523 !!$ TYPE(mesh_type), INTENT(IN) :: H_mesh
524 !!$ REAL(KIND=8), DIMENSION(:), INTENT(IN) :: angles
525 !!$ INTEGER, INTENT(IN) :: nb_angles
526 !!$ INTEGER, INTENT(IN) :: nb, ne
527 !!$ REAL(KIND=8), INTENT(IN) :: time
528 !!$ REAL(KIND=8), DIMENSION(nb_angles,ne-nb+1) :: vv
529 !!$
530 !!$ vv = (1.d0 - penal_in_real_space(H_mesh,H_mesh%rr(:,nb:ne),angles,nb_angles,nb,ne,time))*(mu_disk - 1.d0 ) +1.d0
531 !!$ RETURN
532 !!$ END FUNCTION mu_in_real_space
533 
534 !!$ FUNCTION mu_bar_func(r,z) RESULT(vv)
535 !!$ USE def_type_mesh
536 !!$ USE input_data
537 !!$ USE my_util
538 !!$ IMPLICIT NONE
539 !!$ REAL(KIND=8) :: r,z,vv
540 !!$
541 !!$!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
542 !!$ !mu bar for Blades , a disks with no hole
543 !!$
544 !!$ IF ( if_bottom_prop .AND. if_top_prop) THEN
545 !!$ vv=(top_mu_bar_func(r,z) + bottom_mu_bar_func(r,z))*(mu_disk-1.0) + 1.0
546 !!$ ELSE IF (if_bottom_prop) THEN
547 !!$ vv= bottom_mu_bar_func(r,z)*(mu_disk-1.0) + 1.0
548 !!$ ELSE
549 !!$ vv= top_mu_bar_func(r,z)*(mu_disk-1.0) + 1.0
550 !!$ END IF
551 !!$ RETURN
552 !!$ END FUNCTION mu_bar_func
553 
554 !!$ FUNCTION bottom_mu_bar_func(r,z) RESULT(vv)
555 !!$ USE def_type_mesh
556 !!$ USE input_data
557 !!$ USE my_util
558 !!$ IMPLICIT NONE
559 !!$ REAL(KIND=8) :: r,z,vv
560 !!$ REAL(KIND=8) :: r2,r3,z0,z1
561 !!$ REAL(KIND=8) :: psi
562 !!$
563 !!$!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
564 !!$ !mu bar for Blades , a disks with no hole
565 !!$ r2=disk_rp
566 !!$ r3=disk_r
567 !!$
568 !!$ z0=zbot_bar_p
569 !!$ z1=zbot_bar
570 !!$
571 !!$ psi=0.d0
572 !!$ IF ( z .GE. z1 .AND. r .GE. r3 ) THEN
573 !!$ psi = 0.d0 ;
574 !!$ ELSE IF (r.LE.r2) THEN
575 !!$ IF(z.LE.z0) THEN
576 !!$ psi=1.0;
577 !!$ ELSE IF (z.GE.z1) THEN
578 !!$ psi=0.0;
579 !!$ ELSE
580 !!$ psi=smooth_jump_down(z,z0,z1);
581 !!$ END IF
582 !!$ ELSE IF (r.GE.r2 .AND. r.LE.r3) THEN
583 !!$ IF(z.GE.z1) THEN
584 !!$ psi=0.0;
585 !!$ ELSE IF(z.LE.z0) THEN
586 !!$ psi=smooth_jump_down(r,r2,r3);
587 !!$ ELSE
588 !!$ psi=smooth_jump_down(r,r2,r3)*smooth_jump_down(z,z0,z1);
589 !!$ END IF
590 !!$ END IF
591 !!$
592 !!$ vv = psi
593 !!$ RETURN
594 !!$ END FUNCTION bottom_mu_bar_func
595 
596 !!$ FUNCTION top_mu_bar_func(r,z) RESULT(vv)
597 !!$ USE def_type_mesh
598 !!$ USE input_data
599 !!$ USE my_util
600 !!$ IMPLICIT NONE
601 !!$ REAL(KIND=8) :: r,z,vv,psi
602 !!$ REAL(KIND=8) :: r2,r3,z2,z3
603 !!$
604 !!$!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
605 !!$ !mu bar for Blades , a disks with no hole
606 !!$ r2=disk_rp
607 !!$ r3=disk_r
608 !!$
609 !!$ z2=ztop_bar_p
610 !!$ z3=ztop_bar
611 !!$
612 !!$ psi=0.d0
613 !!$ IF ( z .LE. z3 .AND. r .GE. r3) THEN
614 !!$ psi = 0.d0 ;
615 !!$ ELSE IF(r.LE.r2) THEN
616 !!$ IF(z.GE.z2) THEN
617 !!$ psi=1.0;
618 !!$ ELSE IF (z.LE.z3) THEN
619 !!$ psi=0.0;
620 !!$ ELSE
621 !!$ psi=smooth_jump_up(z,z3,z2);
622 !!$ END IF
623 !!$ ELSE IF (r.GE.r2 .AND. r.LE.r3) THEN
624 !!$ IF(z.LE.z3) THEN
625 !!$ psi=0.0;
626 !!$ ELSE IF(z.GE.z2) THEN
627 !!$ psi=smooth_jump_down(r,r2,r3);
628 !!$ ELSE
629 !!$ psi=smooth_jump_down(r,r2,r3)*smooth_jump_up(z,z3,z2);
630 !!$ END IF
631 !!$ END IF
632 !!$
633 !!$ vv=psi
634 !!$ RETURN
635 !!$ END FUNCTION top_mu_bar_func
636 
637 !!$ FUNCTION grad_mu_bar_func(r,z) RESULT(vv)
638 !!$ USE def_type_mesh
639 !!$ USE input_data
640 !!$ USE my_util
641 !!$ IMPLICIT NONE
642 !!$ REAL(KIND=8) :: r,z
643 !!$ REAL(KIND=8),DIMENSION(2) :: vv
644 !!$
645 !!$!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
646 !!$ !mu bar for Blades , a disks with no hole
647 !!$
648 !!$ IF ( if_bottom_prop .AND. if_top_prop) THEN
649 !!$ vv=( grad_top_mu_bar_func(r,z) + grad_bottom_mu_bar_func(r,z) )*(mu_disk-1.0)
650 !!$ ELSE IF (if_bottom_prop) THEN
651 !!$ vv= grad_bottom_mu_bar_func(r,z)*(mu_disk-1.0)
652 !!$ ELSE
653 !!$ vv= grad_top_mu_bar_func(r,z)*(mu_disk-1.0)
654 !!$ END IF
655 !!$ RETURN
656 !!$ END FUNCTION grad_mu_bar_func
657 
658 !!$ FUNCTION grad_bottom_mu_bar_func(r,z) RESULT(vv)
659 !!$ USE def_type_mesh
660 !!$ USE input_data
661 !!$ USE my_util
662 !!$ IMPLICIT NONE
663 !!$ REAL(KIND=8) :: r,z
664 !!$ REAL(KIND=8),DIMENSION(2) :: vv
665 !!$ REAL(KIND=8) :: r2,r3,z0,z1
666 !!$ REAL(KIND=8) :: DFr,DFz
667 !!$
668 !!$!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
669 !!$ !mu bar for Blades , a disks with no hole
670 !!$ r2=disk_rp
671 !!$ r3=disk_r
672 !!$
673 !!$ z0=zbot_bar_p
674 !!$ z1=zbot_bar
675 !!$
676 !!$ DFr=0.d0
677 !!$ DFz=0.d0
678 !!$
679 !!$ IF ( z .GE. z1 .AND. r .GE. r3 ) THEN
680 !!$ DFr=0.d0
681 !!$ DFz=0.d0
682 !!$ ELSE IF (r.LE.r2) THEN
683 !!$ IF(z.LE.z0) THEN
684 !!$ DFr=0.d0
685 !!$ DFz=0.d0
686 !!$ ELSE IF (z.GE.z1) THEN
687 !!$ DFr=0.d0
688 !!$ DFz=0.d0
689 !!$ ELSE
690 !!$ DFr=0.d0
691 !!$ DFz=Dsmooth_jump_down(z,z0,z1);
692 !!$ END IF
693 !!$ ELSE IF (r.GE.r2 .AND. r.LE.r3) THEN
694 !!$ IF(z.GE.z1) THEN
695 !!$ DFr=0.d0
696 !!$ DFz=0.d0
697 !!$ ELSE IF(z.LE.z0) THEN
698 !!$ DFr=Dsmooth_jump_down(r,r2,r3);
699 !!$ DFz=0.d0
700 !!$ ELSE
701 !!$ DFr=Dsmooth_jump_down(r,r2,r3)*smooth_jump_down(z,z0,z1);
702 !!$ DFz=smooth_jump_down(r,r2,r3)*Dsmooth_jump_down(z,z0,z1);
703 !!$ END IF
704 !!$ END IF
705 !!$
706 !!$ vv(1)=DFr
707 !!$ vv(2)=DFz
708 !!$ RETURN
709 !!$ END FUNCTION grad_bottom_mu_bar_func
710 
711 !!$ FUNCTION grad_top_mu_bar_func(r,z) RESULT(vv)
712 !!$ USE def_type_mesh
713 !!$ USE input_data
714 !!$ USE my_util
715 !!$ IMPLICIT NONE
716 !!$ REAL(KIND=8) :: r,z
717 !!$ REAL(KIND=8),DIMENSION(2) :: vv
718 !!$ REAL(KIND=8) :: r2,r3,z2,z3
719 !!$ REAL(KIND=8) :: DFr,DFz
720 !!$
721 !!$!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
722 !!$ !mu bar for Blades , a disks with no hole
723 !!$ r2=disk_rp
724 !!$ r3=disk_r
725 !!$
726 !!$ z2=ztop_bar_p
727 !!$ z3=ztop_bar
728 !!$
729 !!$ DFr=0.d0
730 !!$ DFz=0.d0
731 !!$ IF ( z .LE. z3 .AND. r .GE. r3) THEN
732 !!$ DFr=0.d0
733 !!$ DFz=0.d0
734 !!$ ELSE IF(r.LE.r2) THEN
735 !!$ IF(z.GE.z2) THEN
736 !!$ DFr=0.d0
737 !!$ DFz=0.d0
738 !!$ ELSE IF (z.LE.z3) THEN
739 !!$ DFr=0.d0
740 !!$ DFz=0.d0
741 !!$ ELSE
742 !!$ DFr=0.d0
743 !!$ DFz=Dsmooth_jump_up(z,z3,z2);
744 !!$ END IF
745 !!$ ELSE IF (r.GE.r2 .AND. r.LE.r3) THEN
746 !!$ IF(z.LE.z3) THEN
747 !!$ DFr=0.d0
748 !!$ DFz=0.d0
749 !!$ ELSE IF(z.GE.z2) THEN
750 !!$ DFr=Dsmooth_jump_down(r,r2,r3);
751 !!$ DFz=0.d0
752 !!$ ELSE
753 !!$ DFr=Dsmooth_jump_down(r,r2,r3)*smooth_jump_up(z,z3,z2);
754 !!$ DFz=smooth_jump_down(r,r2,r3)*Dsmooth_jump_up(z,z3,z2);
755 !!$ END IF
756 !!$ END IF
757 !!$
758 !!$ vv(1)=DFr
759 !!$ vv(2)=DFz
760 !!$ RETURN
761 !!$ END FUNCTION grad_top_mu_bar_func
762 
763  !This is 1 in the fluid, and 0 in the solid
764 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
765  FUNCTION smooth_penal_in_real_space(mesh,rr_gauss,angles,nb_angles,nb,ne,time) RESULT(vv)
766  USE def_type_mesh
767  USE input_data
768  USE user_data
769  USE my_util
770  IMPLICIT NONE
771  TYPE(mesh_type) :: mesh
772  REAL(KIND=8), DIMENSION(:,:) :: rr_gauss
773  REAL(KIND=8), DIMENSION(:) :: angles
774  INTEGER :: nb_angles
775  INTEGER :: nb, ne
776  REAL(KIND=8), DIMENSION(nb_angles,ne-nb+1) :: vv
777  INTEGER :: n, n_loc
778  REAL(KIND=8) :: r,z,time
779 
780  DO n = nb, ne
781  n_loc = n - nb + 1
782 
783  r = rr_gauss(1,n_loc)
784  z = rr_gauss(2,n_loc)
785 
786  ! DCQ 14/08/2014
787  IF (if_bottom_prop .AND. if_top_prop) THEN
788  vv(:,n_loc) = smooth_top_propeller(r,z,angles, nb_angles,time)&
789  *smooth_bottom_propeller(r,z,angles, nb_angles,time)
790  ELSE IF (if_bottom_prop) THEN
791  vv(:,n_loc) = smooth_bottom_propeller(r,z,angles, nb_angles,time)
792  ELSE
793  vv(:,n_loc) = smooth_top_propeller(r,z,angles, nb_angles,time)
794  END IF
795  END DO
796  RETURN
797 
798  !===Dummies variables to avoid warning
799  n=mesh%np
800  !===Dummies variables to avoid warning
801  END FUNCTION smooth_penal_in_real_space
802 
803  ! DCQ 14/08/2014
804  ! 1 - Characteristic Function of top_propeller
805  FUNCTION smooth_top_propeller(r,z,angles,nb_angles,time) RESULT(vv)
806  USE def_type_mesh
807  USE input_data
808  USE user_data
809  USE my_util
810  IMPLICIT NONE
811  REAL(KIND=8) :: r,theta,z,time
812  REAL(KIND=8), DIMENSION(:) :: angles
813  INTEGER :: nb_angles
814  REAL(KIND=8), DIMENSION(nb_angles) :: vv
815  INTEGER :: na
816  REAL(KIND=8) :: g, a,alphaz,alphar
817  REAL(KIND=8) :: r0,r1,r2,r3,z0,z1,z2,z3
818  REAL(KIND=8) :: psi, tanhp,tanhm, tanhd, r_theta
819 
820  !Get characteristic function of the supporting disk-cylinder
821  psi=smooth_top_supporting_disk(r,z)
822 
823  !If we are outside of the supporting disk (respect to r)
824  IF(abs(psi) .LE. 1.d-8) THEN
825  DO na = 1, nb_angles
826  vv(na) = 1-psi
827  END DO
828  RETURN
829  END IF
830 
831  !Do Blades stuff
832  r0=hole_rp
833  r1=hole_r
834  r2=disk_rp
835  r3=disk_r
836 
837  z0=ztop_p - wjump
838  z1=ztop - wjump
839  z2=tdisk_z_p
840  z3=tdisk_z
841 
842  ! Parabolic jump
843  IF (z .LE. z1 ) THEN
844  alphaz = 1.d0;
845  ELSE IF(z .LE. z0 .AND. z .GE. z1) THEN
846  alphaz= smooth_jump_down(z,z1,z0);
847  ELSE
848  alphaz=0.0;
849  END IF
850 
851  If ( r .LE. r0 ) THEN
852  alphar = 0.d0;
853  ELSE IF(r .GE. r0 .AND. r .LE. r1) THEN
854  alphar= smooth_jump_up(r,r0,r1);
855  ELSE
856  alphar=1.0;
857  END IF
858  alphaz= alpha_th*alphaz*alphar
859 
860  r_theta = asin(r/two_rp)
861 
862  DO na = 1, nb_angles
863 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
864  !DCQ go backwards and do the test
865  !These blades rotate the other way (-lbd), (+user%solid_vel*time) and they begin at the same angle,
866 
867  theta= angles(na) + solid_vel*time
868  theta= theta - top_propeller_angle_offset !some offset to begin at the same vertical axis as the bottom propeller
869 
870  !a = theta + (-lbd*r) - FLOOR(( theta + (-lbd)*r)/(twopi/nblades))*twopi/nblades &
871  ! - twopi/(2*nblades)
872 
873  !JL-CN 01/2015
874  a = theta - r_theta - floor(( theta - r_theta)/(twopi/nblades))*twopi/nblades &
875  - twopi/(2*nblades)
876 
877  tanhp = tanh(alphaz*(r*a+lw+lw*r))
878  tanhm = tanh(alphaz*(r*a-lw-lw*r))
879  tanhd = tanh(alphaz*(lw+lw*r))
880  g=(1+tanhp)*(1-tanhm)/(1+tanhd)**2
881  vv(na) = 1-g*psi
882 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
883  END DO
884  RETURN
885  END FUNCTION smooth_top_propeller
886 
887  ! DCQ 14/08/2014
888  ! 1 - Characteristic Function of bottom_propeller
889  FUNCTION smooth_bottom_propeller(r,z,angles,nb_angles,time) RESULT(vv)
890  USE def_type_mesh
891  USE input_data
892  USE user_data
893  USE my_util
894  IMPLICIT NONE
895  REAL(KIND=8) :: r,theta,z,time
896  REAL(KIND=8), DIMENSION(:) :: angles
897  INTEGER :: nb_angles
898  REAL(KIND=8), DIMENSION(nb_angles) :: vv
899  INTEGER :: na
900  REAL(KIND=8) :: g, a,alphaz,alphar
901  REAL(KIND=8) :: r0,r1,r2,r3,z0,z1,z2,z3
902  REAL(KIND=8) :: psi, tanhp,tanhm, tanhd, r_theta
903 
904  !Supporting disk stuff
905  !Get characteristic function of the supporting disk-cylinder
906  psi=smooth_bottom_supporting_disk(r,z)
907 
908  !If we are outside of the supporting disk (respect to r)
909  IF(abs(psi) .LE. 1.d-8) THEN
910  DO na = 1, nb_angles
911  vv(na) = 1-psi
912  END DO
913  RETURN
914  END IF
915 
916  ! Do blades stuff
917  !! Blades with no hole in the disk
918  r0=hole_rp
919  r1=hole_r
920  r2=disk_rp
921  r3=disk_r
922 
923  z0=zbot_p + wjump
924  z1=zbot + wjump
925  z2=bdisk_z_p
926  z3=bdisk_z
927 
928  ! Parabolic jump
929  IF (z .LE. z0 ) THEN
930  alphaz = 0.d0;
931  ELSE IF(z .GE. z0 .AND. z .LE. z1) THEN
932  alphaz= smooth_jump_up(z,z0,z1);
933  ELSE
934  alphaz=1.0;
935  END IF
936 
937  IF ( r .LE. r0 ) THEN
938  alphar = 0.d0;
939  ELSE IF(r .GE. r0 .AND. r .LE. r1) THEN
940  alphar= smooth_jump_up(r,r0,r1);
941  ELSE
942  alphar=1.0;
943  END IF
944  alphaz= alpha_th*alphaz*alphar
945 
946  r_theta = asin(r/two_rp)
947 
948  DO na = 1, nb_angles
949 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
950  !DCQ go backwards and do the test
951  theta= angles(na) - solid_vel*time
952  !a = theta + lbd*r - FLOOR(( theta + lbd*r)/(twopi/nblades))*twopi/nblades &
953  ! - twopi/(2*nblades)
954 
955  !JL-CN 01/2015
956  a = theta + r_theta - floor(( theta + r_theta)/(twopi/nblades))*twopi/nblades &
957  - twopi/(2*nblades)
958 
959  tanhp = tanh(alphaz*(r*a+lw+lw*r))
960  tanhm = tanh(alphaz*(r*a-lw-lw*r))
961  tanhd = tanh(alphaz*(lw+lw*r))
962  g=(1+tanhp)*(1-tanhm)/(1+tanhd)**2
963  vv(na) = 1-g*psi
964  END DO
965  RETURN
966  END FUNCTION smooth_bottom_propeller
967 
968  !Characteristic Function of top_supporting_disk
969  FUNCTION smooth_top_supporting_disk(r,z) RESULT(vv)
970  USE def_type_mesh
971  USE input_data
972  USE user_data
973  USE my_util
974  IMPLICIT NONE
975  REAL(KIND=8) :: r,z
976  REAL(KIND=8) :: vv
977  REAL(KIND=8) :: r0,r1,r2,r3,z0,z1,z2,z3
978  REAL(KIND=8) :: curve_1,curve_2,psi
979 
980 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
981  !! Blades with no hole in the disk
982  !Supporting disk stuff
983  r0=hole_rp
984  r1=hole_r
985  r2=disk_rp
986  r3=disk_r
987 
988  z0=ztop_p - wjump
989  z1=ztop - wjump
990  z2=tdisk_z_p
991  z3=tdisk_z
992 
993  psi=0.d0
994  IF ( z .LE. z3 .AND. r .GE. r3) THEN
995  psi = 0.d0 ;
996  ELSE IF (r.LE.r0) THEN
997  IF(z.GE.z0) THEN
998  psi=1.0;
999  ELSE IF (z.LE.z1) THEN
1000  psi=0.0;
1001  ELSE
1002  psi=smooth_jump_up(z,z1,z0);
1003  END IF
1004  ELSE IF(r.GE.r0 .AND. r.LE.r1) THEN
1005  curve_2= smooth_jump_up(r,r0,r1)*(z3-z1)+z1;
1006  curve_1= smooth_jump_up(r,r0,r1)*(z2-z0)+z0;
1007  IF(z.LE.curve_2) THEN
1008  psi=0.0;
1009  ELSE IF(z.GE.curve_1) THEN
1010  psi=1.0;
1011  ELSE
1012  psi = 1 - smooth_jump_up(z ,curve_1,curve_2);
1013  END IF
1014  ELSE IF(r.GE.r1 .AND. r.LE.r2) THEN
1015  IF(z.GE.z2) THEN
1016  psi=1.0;
1017  ELSE IF (z.LE.z3) THEN
1018  psi=0.0;
1019  ELSE
1020  psi=smooth_jump_up(z,z3,z2);
1021  END IF
1022  ELSE IF (r.GE.r2 .AND. r.LE.r3) THEN
1023  IF(z.LE.z3) THEN
1024  psi=0.0;
1025  ELSE IF(z.GE.z2) THEN
1026  psi=smooth_jump_down(r,r2,r3);
1027  ELSE
1028  psi=smooth_jump_down(r,r2,r3)*smooth_jump_up(z,z3,z2);
1029  END IF
1030  END IF
1031 
1032  vv=psi
1033 
1034  END FUNCTION smooth_top_supporting_disk
1035 
1036  !Characteristic Function of bot_supporting_disk
1037  FUNCTION smooth_bottom_supporting_disk(r,z) RESULT(vv)
1038  USE def_type_mesh
1039  USE input_data
1040  USE user_data
1041  USE my_util
1042  IMPLICIT NONE
1043  REAL(KIND=8) :: r,z
1044  REAL(KIND=8) :: vv
1045  REAL(KIND=8) :: r0,r1,r2,r3,z0,z1,z2,z3
1046  REAL(KIND=8) :: curve_1,curve_2,psi
1047 
1048 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1049  !! Blades with no hole in the disk
1050  !Supporting disk stuff
1051 
1052  r0=hole_rp
1053  r1=hole_r
1054  r2=disk_rp
1055  r3=disk_r
1056 
1057  z0=zbot_p + wjump
1058  z1=zbot + wjump
1059  z2=bdisk_z_p
1060  z3=bdisk_z
1061 
1062  psi=0.d0
1063  IF ( z .GE. z3 .AND. r .GE. r3) THEN
1064  psi = 0.d0 ;
1065  ELSE IF (r.LE.r0) THEN
1066  IF(z.LE.z0) THEN
1067  psi=1.0;
1068  ELSE IF (z.GE.z1) THEN
1069  psi=0.0;
1070  ELSE
1071  psi=smooth_jump_up(z,z1,z0);
1072  END IF
1073  ELSE IF(r.GE.r0 .AND. r.LE.r1) THEN
1074  curve_2= smooth_jump_up(r,r0,r1)*(z3-z1)+z1;
1075  curve_1= smooth_jump_up(r,r0,r1)*(z2-z0)+z0;
1076  IF(z.GE.curve_2) THEN
1077  psi=0.0;
1078  ELSE IF(z.LE.curve_1) THEN
1079  psi=1.0;
1080  ELSE
1081  psi = 1 - smooth_jump_up(z ,curve_1,curve_2);
1082  END IF
1083  ELSE IF(r.GE.r1 .AND. r.LE.r2) THEN
1084  IF(z.LE.z2) THEN
1085  psi=1.0;
1086  ELSE IF (z.GE.z3) THEN
1087  psi=0.0;
1088  ELSE
1089  psi=smooth_jump_up(z,z3,z2);
1090  END IF
1091  ELSE IF (r.GE.r2 .AND. r.LE.r3) THEN
1092  IF(z.GE.z3) THEN
1093  psi=0.0;
1094  ELSE IF(z.LE.z2) THEN
1095  psi=smooth_jump_down(r,r2,r3);
1096  ELSE
1097  psi=smooth_jump_down(r,r2,r3)*smooth_jump_up(z,z3,z2);
1098  END IF
1099  END IF
1100 
1101  vv=psi
1102 
1103  END FUNCTION smooth_bottom_supporting_disk
1104 
1105  !A cubic profile, which is 1 at x0 and 0 at x1
1106  FUNCTION smooth_jump_down(x,x0,x1) RESULT(vv)
1107  IMPLICIT NONE
1108  REAL(KIND=8) :: x,x0,x1
1109  REAL(KIND=8) :: vv
1110  REAL(KIND=8) :: a0,a1,a2,a3
1111 
1112  !Cubic
1113  a0 = x1**2*(3*x0-x1)/(x0-x1)**3;
1114  a1 = -6.0*x0*x1/(x0-x1)**3;
1115  a2 = (3.0*(x0+x1))/(x0-x1)**3;
1116  a3 = -2.0/(x0-x1)**3;
1117  vv = a0+a1*x+a2*x*x + a3*x*x*x
1118  RETURN
1119  END FUNCTION smooth_jump_down
1120 
1121  !A cubic profile, which is 0 at x0 and 1 at x1
1122  FUNCTION smooth_jump_up(x,x0,x1) RESULT(vv)
1123  IMPLICIT NONE
1124  REAL(KIND=8) :: x,x0,x1
1125  REAL(KIND=8) :: vv
1126 
1127  vv = 1.d0 - smooth_jump_down( x,x0,x1 );
1128  RETURN
1129  END FUNCTION smooth_jump_up
1130 
1131 !!$ !derivative with respect to x
1132 !!$ FUNCTION Dsmooth_jump_down(x,x0,x1) RESULT(vv)
1133 !!$ IMPLICIT NONE
1134 !!$ REAL(KIND=8) :: x,x0,x1
1135 !!$ REAL(KIND=8) :: vv
1136 !!$ REAL(KIND=8) :: a0,a1,a2,a3
1137 !!$
1138 !!$ !Cubic Factorized
1139 !!$ a0 = x1**2*(3*x0-x1)/(x0-x1)**3;
1140 !!$ a1 = -6.0*x0*x1/(x0-x1)**3;
1141 !!$ a2 = (3.0*(x0+x1))/(x0-x1)**3;
1142 !!$ a3 = -2.0/(x0-x1)**3;
1143 !!$
1144 !!$ vv = a1+2.d0*a2*x + 3.d0*a3*x*x
1145 !!$ RETURN
1146 !!$ END FUNCTION Dsmooth_jump_down
1147 
1148 !!$ !derivative with respect to x
1149 !!$ FUNCTION Dsmooth_jump_up(x,x0,x1) RESULT(vv)
1150 !!$ IMPLICIT NONE
1151 !!$ REAL(KIND=8) :: x,x0,x1
1152 !!$ REAL(KIND=8) :: vv
1153 !!$
1154 !!$ vv = - Dsmooth_jump_down( x,x0,x1 );
1155 !!$ RETURN
1156 !!$ END FUNCTION Dsmooth_jump_up
1157 
1158 END MODULE boundary_test_28
boundary_test_28::ztop_bar_p
real(kind=8), parameter ztop_bar_p
Definition: condlim_test_28.f90:70
boundary_test_28::vv_exact
real(kind=8) function, dimension(size(rr, 2)), public vv_exact(TYPE, rr, m, t)
Definition: condlim_test_28.f90:220
boundary_test_28::init_velocity_pressure
subroutine, public init_velocity_pressure(mesh_f, mesh_c, time, dt, list_mode, un_m1, un, pn_m1, pn, phin_m1, phin)
Definition: condlim_test_28.f90:89
boundary_test_28::penal_in_real_space
real(kind=8) function, dimension(nb_angles, ne-nb+1), public penal_in_real_space(mesh, rr_gauss, angles, nb_angles, nb, ne, time)
Definition: condlim_test_28.f90:337
boundary_test_28::top_of_blade_bot
real(kind=8), parameter top_of_blade_bot
Definition: condlim_test_28.f90:47
boundary_test_28::disk_r
real(kind=8), parameter, public disk_r
Definition: condlim_test_28.f90:58
boundary_test_28::zbot
real(kind=8), parameter zbot
Definition: condlim_test_28.f90:63
boundary_test_28::zbot_bar
real(kind=8), parameter zbot_bar
Definition: condlim_test_28.f90:64
boundary_test_28::omega_vol
real(kind=8), public omega_vol
Definition: condlim_test_28.f90:53
bessel
Definition: bessel.f90:4
boundary_test_28::tdisk_z
real(kind=8), parameter tdisk_z
Definition: condlim_test_28.f90:68
boundary_test_28::two_rp
real(kind=8), parameter two_rp
Definition: condlim_test_28.f90:51
boundary_test_28::alpha
real(kind=8), parameter alpha
Definition: condlim_test_28.f90:73
boundary_test_28
Definition: condlim_test_28.f90:1
def_type_mesh::mesh_type
Definition: def_type_mesh.f90:74
my_util
Definition: my_util.f90:1
my_util::error_petsc
subroutine error_petsc(string)
Definition: my_util.f90:16
boundary_test_28::if_top_prop
logical, parameter if_top_prop
Definition: condlim_test_28.f90:77
boundary_test_28::hole_r
real(kind=8), parameter hole_r
Definition: condlim_test_28.f90:57
boundary_test_28::smooth_top_propeller
real(kind=8) function, dimension(nb_angles) smooth_top_propeller(r, z, angles, nb_angles, time)
Definition: condlim_test_28.f90:806
boundary_test_28::bot_of_blade_top
real(kind=8), parameter bot_of_blade_top
Definition: condlim_test_28.f90:48
boundary_test_28::disk_bot
real(kind=8), parameter disk_bot
Definition: condlim_test_28.f90:47
boundary_test_28::source_in_ns_momentum
real(kind=8) function, dimension(size(rr, 2)), public source_in_ns_momentum(TYPE, rr, mode, i, time, Re, ty, opt_density, opt_tempn)
Definition: condlim_test_28.f90:165
boundary_test_28::top_propeller_angle_offset
real(kind=8), parameter top_propeller_angle_offset
Definition: condlim_test_28.f90:44
boundary_test_28::mu_disk
real(kind=8), parameter mu_disk
Definition: condlim_test_28.f90:37
boundary_test_28::zbot_p
real(kind=8), parameter zbot_p
Definition: condlim_test_28.f90:63
boundary_test_28::hole_radius
real(kind=8), parameter, public hole_radius
Definition: condlim_test_28.f90:49
boundary_test_28::smooth_jump_up
real(kind=8) function smooth_jump_up(x, x0, x1)
Definition: condlim_test_28.f90:1123
boundary_test_28::smooth_top_supporting_disk
real(kind=8) function smooth_top_supporting_disk(r, z)
Definition: condlim_test_28.f90:970
boundary_test_28::smooth_bottom_supporting_disk
real(kind=8) function smooth_bottom_supporting_disk(r, z)
Definition: condlim_test_28.f90:1038
boundary_test_28::smooth_bottom_propeller
real(kind=8) function, dimension(nb_angles) smooth_bottom_propeller(r, z, angles, nb_angles, time)
Definition: condlim_test_28.f90:890
boundary_test_28::imposed_velocity_by_penalty
real(kind=8) function, dimension(size(rr, 2), 6), public imposed_velocity_by_penalty(rr, t)
Definition: condlim_test_28.f90:257
boundary_test_28::ztop
real(kind=8), parameter ztop
Definition: condlim_test_28.f90:69
boundary_test_28::bot_of_disk_top
real(kind=8), parameter bot_of_disk_top
Definition: condlim_test_28.f90:48
boundary_test_28::pp_exact
real(kind=8) function, dimension(size(rr, 2)), public pp_exact(TYPE, rr, m, t)
Definition: condlim_test_28.f90:287
boundary_test_28::cyl_top
real(kind=8), parameter cyl_top
Definition: condlim_test_28.f90:67
boundary_test_28::zbot_bar_p
real(kind=8), parameter zbot_bar_p
Definition: condlim_test_28.f90:64
user_data
Definition: read_user_data.f90:10
boundary_test_28::alpha_th
real(kind=8), parameter alpha_th
Definition: condlim_test_28.f90:73
boundary_test_28::tdisk_z_p
real(kind=8), parameter tdisk_z_p
Definition: condlim_test_28.f90:68
boundary_test_28::smooth_penal_in_real_space
real(kind=8) function, dimension(nb_angles, ne-nb+1) smooth_penal_in_real_space(mesh, rr_gauss, angles, nb_angles, nb, ne, time)
Definition: condlim_test_28.f90:766
boundary_test_28::bdisk_z
real(kind=8), parameter bdisk_z
Definition: condlim_test_28.f90:62
boundary_test_28::smooth_jump_down
real(kind=8) function smooth_jump_down(x, x0, x1)
Definition: condlim_test_28.f90:1107
boundary_test_28::lw
real(kind=8), parameter lw
Definition: condlim_test_28.f90:46
boundary_test_28::solid_vel
real(kind=8), parameter solid_vel
Definition: condlim_test_28.f90:78
boundary_test_28::bdisk_z_p
real(kind=8), parameter bdisk_z_p
Definition: condlim_test_28.f90:62
boundary_test_28::wjump
real(kind=8), parameter wjump
Definition: condlim_test_28.f90:56
boundary_test_28::hole_rp
real(kind=8), parameter hole_rp
Definition: condlim_test_28.f90:57
boundary_test_28::ztop_bar
real(kind=8), parameter ztop_bar
Definition: condlim_test_28.f90:70
boundary_test_28::top_of_disk_bot
real(kind=8), parameter top_of_disk_bot
Definition: condlim_test_28.f90:47
input_data
Definition: read_my_data.f90:233
boundary_test_28::disk_top
real(kind=8), parameter disk_top
Definition: condlim_test_28.f90:48
boundary_test_28::if_bottom_prop
logical, parameter if_bottom_prop
Definition: condlim_test_28.f90:76
boundary_test_28::ztop_p
real(kind=8), parameter ztop_p
Definition: condlim_test_28.f90:69
boundary_test_28::cyl_bott
real(kind=8), parameter cyl_bott
Definition: condlim_test_28.f90:61
boundary_test_28::disk_radius
real(kind=8), parameter, public disk_radius
Definition: condlim_test_28.f90:49
def_type_mesh
Definition: def_type_mesh.f90:4
boundary_test_28::disk_rp
real(kind=8), parameter, public disk_rp
Definition: condlim_test_28.f90:58
boundary_test_28::twopi
real(kind=8), parameter, private twopi
Definition: condlim_test_28.f90:36
boundary_test_28::wjump_hole
real(kind=8), parameter wjump_hole
Definition: condlim_test_28.f90:56
boundary_test_28::nblades
integer nblades
Definition: condlim_test_28.f90:45