SFEMaNS  version 4.1 (work in progress)
Reference documentation for SFEMaNS
 All Classes Files Functions Variables Groups Pages
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 IF(z .GE. z0) THEN
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 IF(z .GE. z1) THEN
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
t
section doc_intro_frame_work_num_app Numerical approximation subsection doc_intro_fram_work_num_app_Fourier_FEM Fourier Finite element representation The SFEMaNS code uses a hybrid Fourier Finite element formulation The Fourier decomposition allows to approximate the problem’s solutions for each Fourier mode modulo nonlinear terms that are made explicit The variables are then approximated on a meridian section of the domain with a finite element method The numerical approximation of a function f $f f is written in the following generic t
Definition: doc_intro.h:199
user_data
Definition: read_user_data.f90:10
boundary_test_1::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_1.f90:116
boundary_test_28::smooth_bottom_supporting_disk
real(kind=8) function smooth_bottom_supporting_disk(r, z)
Definition: condlim_test_28.f90:1037
boundary_test_1::pp_exact
real(kind=8) function, dimension(size(rr, 2)), public pp_exact(TYPE, rr, m, t)
Definition: condlim_test_1.f90:354
boundary_test_28::smooth_jump_down
real(kind=8) function smooth_jump_down(x, x0, x1)
Definition: condlim_test_28.f90:1106
def_type_mesh::mesh_type
Definition: def_type_mesh.f90:73
theta
section doc_intro_frame_work_num_app Numerical approximation subsection doc_intro_fram_work_num_app_Fourier_FEM Fourier Finite element representation The SFEMaNS code uses a hybrid Fourier Finite element formulation The Fourier decomposition allows to approximate the problem’s solutions for each Fourier mode modulo nonlinear terms that are made explicit The variables are then approximated on a meridian section of the domain with a finite element method The numerical approximation of a function f $f f is written in the following generic theta
Definition: doc_intro.h:193
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:765
def_type_mesh
Definition: def_type_mesh.f90:4
my_util
Definition: my_util.f90:1
input_data
Definition: read_my_data.f90:200
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:805
boundary_test_24::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_24.f90:328
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:889
boundary_test_1::vv_exact
real(kind=8) function, dimension(size(rr, 2)), public vv_exact(TYPE, rr, m, t)
Definition: condlim_test_1.f90:301
boundary_test_24::imposed_velocity_by_penalty
real(kind=8) function, dimension(size(rr, 2), 6), public imposed_velocity_by_penalty(rr, t)
Definition: condlim_test_24.f90:254
boundary_test_28
Definition: condlim_test_28.f90:1
boundary_test_1::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_1.f90:39
my_util::error_petsc
subroutine error_petsc(string)
Definition: my_util.f90:15
boundary_test_28::smooth_jump_up
real(kind=8) function smooth_jump_up(x, x0, x1)
Definition: condlim_test_28.f90:1122
boundary_test_28::smooth_top_supporting_disk
real(kind=8) function smooth_top_supporting_disk(r, z)
Definition: condlim_test_28.f90:969
bessel
Definition: bessel.f90:4
z
section doc_intro_frame_work_num_app Numerical approximation subsection doc_intro_fram_work_num_app_Fourier_FEM Fourier Finite element representation The SFEMaNS code uses a hybrid Fourier Finite element formulation The Fourier decomposition allows to approximate the problem’s solutions for each Fourier mode modulo nonlinear terms that are made explicit The variables are then approximated on a meridian section of the domain with a finite element method The numerical approximation of a function f $f f is written in the following generic z
Definition: doc_intro.h:193