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