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