SFEMaNS  version 5.3
Reference documentation for SFEMaNS
condlim_test_34.f90
Go to the documentation of this file.
1 MODULE boundary_test_34
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)
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 
137  r = rr(1,:)
138  z = rr(2,:)
139 
140  DO j=1,SIZE(rr,2)
141  IF (r(j).LE.r0) THEN
142  lambda(j) = inputs%temperature_diffusivity(1)
143  ELSE
144  lambda(j) = inputs%temperature_diffusivity(2)
145  END IF
146  END DO
147 
148  IF (type==5) THEN
149  vv = alpha*(opt_tempn(:,1,i) - temperature_exact(1,rr,mode,time))
150  ELSE IF (type==6) THEN
151  vv = alpha*(opt_tempn(:,2,i) - temperature_exact(2,rr,mode,time))
152  ELSE
153  vv = 0.d0
154  END IF
155 
156  IF (type==1) THEN
157  IF (mode==0) THEN
158  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) + &
159  (r - r0)*re*cos(time)**2*(14*r**3 - 5*r**2*r0 - 5*r*r0**2 + 2*r0**3 + &
160  (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) + &
161  r**2*(5*r0 - 12*pi**2*r0**3))*cos(4*pi*z)) - &
162  4*pi*r**3*(r - r0)**2*re*cos(2*pi*z)*sin(time))/(2.*r**3*re)
163  ELSE IF (mode==1) THEN
164  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) + &
165  ((3*r**2 - 4*r*r0 + r0**2)*re*cos(time)**2*(3*r**2 - r0**2 + &
166  (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. - &
167  pi*r**3*(r - r0)**2*re*cos(2*pi*z)*sin(time))/(r**3*re)
168  ELSE IF (mode==2) THEN
169  vv = vv -((r - r0)*cos(time)**2*(4*r**2 - r*r0 - r0**2 + (12*pi**2*r**4 - &
170  28*pi**2*r**3*r0 + r0**2 + 4*r**2*(-1.d0 + 5*pi**2*r0**2) + &
171  r*(r0 - 4*pi**2*r0**3))*cos(4*pi*z)))/(2.*r**2)
172  END IF
173  ELSE IF (type==2) THEN
174  IF (mode==1) THEN
175  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 + &
176  r**2*(-3.d0 + 4*pi**2*r0**2) + (3*r**2 + r*r0 - r0**2)*cos(4*pi*z))))/(r**3*re)
177  ELSE IF (mode==2) THEN
178  vv = vv + ((r - r0)**2*cos(time)**2*(4*pi**2*r**4 - r*r0 - 8*pi**2*r**3*r0 + r0**2 + &
179  r**2*(-3.d0 + 4*pi**2*r0**2) + (3*r**2 + r*r0 - r0**2)*cos(4*pi*z)))/(2.*r**3)
180  END IF
181  ELSE IF (type==3) THEN
182  IF (mode==0) THEN
183  vv = vv + (2*pi*(-3*pi*r*(r - r0)**3*(3*r - r0)*re*cos(time)**2 + &
184  (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) - &
185  r**2*(r - r0)**2*re*cos(2*pi*z)*sin(time)))/(r**2*re)
186  ELSE IF (mode==1) THEN
187  vv = vv + (4*pi*r*(2*pi**2*r**4 - r*r0 - 4*pi**2*r**3*r0 + r0**2 + &
188  r**2*(-1 + 2*pi**2*r0**2))*cos(time)*cos(2*pi*z) + &
189  ((r - r0)**3*(3*r - r0)*re*cos(time)**2*(-1.d0 - 16*pi**2*r**2 + cos(4*pi*z)))/2. - &
190  2*pi*r**3*(r - r0)**2*re*cos(2*pi*z)*sin(time))/(r**3*re)
191  ELSE IF (mode==2) THEN
192  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)
193  END IF
194  ELSE IF (type==4) THEN
195  IF (mode==1) THEN
196  vv = vv + ((r - r0)**2*cos(time)*(-8*pi*r*cos(2*pi*z) + re*cos(time)*((-3*r + r0)**2 + &
197  (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)
198  ELSE IF (mode==2) THEN
199  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
200  END IF
201  ELSE IF (type==5) THEN
202  IF (mode==0) THEN
203  vv = vv + ((4*(12*pi**2*r**4 - 16*pi**2*r**3*r0 - r0**2 + &
204  r**2*(-3.d0 + 4*pi**2*r0**2))*cos(time)*sin(2*pi*z) - &
205  4*r**2*(3*r**2 - 4*r*r0 + r0**2)*re*sin(time)*sin(2*pi*z) + &
206  pi*r*(r - r0)**2*(12*pi**2*r**4 - r*r0 - 24*pi**2*r**3*r0 + 2*r0**2 + &
207  4*r**2*(-1.d0 + 3*pi**2*r0**2))*re*4*cos(time)**2*sin(4*pi*z)))/(4.*r**3*re)
208  ELSE IF (mode==1) THEN
209  vv = vv + (4*(-r0 + pi**2*r*(3*r**2 - 4*r*r0 + r0**2))*cos(time)*sin(2*pi*z) - &
210  r*(3*r**2 - 4*r*r0 + r0**2)*re*sin(time)*sin(2*pi*z) + &
211  pi*(r - r0)**2*(16*pi**2*r**4 - 2*r*r0 - 32*pi**2*r**3*r0 + 3*r0**2 + &
212  r**2*(-5.d0 + 16*pi**2*r0**2))*re*cos(time)**2*sin(4*pi*z))/(r**2*re)
213  ELSE IF (mode==2) THEN
214  vv = vv + (pi*(r - r0)**2*(4*pi**2*r**4 - r*r0 - 8*pi**2*r**3*r0 + r0**2 + &
215  r**2*(-1.d0 + 4*pi**2*r0**2))*cos(time)**2*sin(4*pi*z))/r**2
216  END IF
217  ELSE IF (type==6) THEN
218  IF (mode==1) THEN
219  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) -&
220  4*pi*r*(r - r0)**3*re*cos(time)**2*sin(4*pi*z))/(r**2*re)
221  ELSE IF (mode==2) THEN
222  vv = vv -2*pi*(r - r0)**3*cos(time)**2*sin(4*pi*z)/r
223  END IF
224  END IF
225 
226  IF ((type==1).AND.(mode==1)) THEN
227  vv = vv + 3*r**2*cos(time)*sin(2*pi*z)
228  ELSE IF ((type==4).AND.(mode==1)) THEN
229  vv = vv - r**2*cos(time)*sin(2*pi*z)
230  ELSE IF ((type==5).AND.(mode==1)) THEN
231  vv = vv + 2*pi*r**3*cos(time)*cos(2*pi*z)
232  END IF
233 
234  IF (TYPE == 1) then
235  IF (mode == 0) THEN
236  vv = vv + r**7 * (r - r0)**2 * cos(time)**4 * (215.d0 + 204 * pi**2 * r**2 + &
237  (215.d0 - 204 * pi**2 * r**2) * cos(4*pi * z)) * sin(2*pi*z)**2 / &
238  (8 * lambda**2)
239  ELSE IF (mode == 1) THEN
240  vv = vv - 5 * r**7 * (r - r0)**2 * cos(time)**4 * (-11.d0 - 12 * pi**2 * r**2 + &
241  (-11.d0 + 12 * pi**2 * r**2) * cos(4*pi*z)) * sin(2*pi*z)**2 / &
242  (2 * lambda**2)
243  ELSE IF (mode == 2) THEN
244  vv = vv - 7 * r**7 * (r - r0)**2 * cos(time)**4 * sin(4*pi*z)**2 / &
245  (2 * lambda**2)
246  ELSE IF (mode == 3) THEN
247  vv = vv + r**7 * (r - r0)**2 * cos(time)**4 * (-19.d0 - 12 * pi**2 *r**2 + &
248  (-19.d0 + 12 * pi**2 *r**2) * cos(4*pi*z)) * sin(2*pi*z)**2 / &
249  (2 * lambda**2)
250  ELSE IF (mode == 4) THEN
251  vv = vv + 3 * r**7 * (r - r0)**2 * cos(time)**4 * (-5.d0 - 4 * pi**2 * r**2 + &
252  (-5.d0 + 4 * pi**2 * r**2) * cos(4*pi*z)) * sin(2*pi*z)**2 / &
253  (8 * lambda**2)
254  END IF
255  ELSE IF (TYPE == 2) then
256  IF (mode == 1) THEN
257  vv = vv - 6 * r**7 * (r - r0)**2 * cos(time)**4 * (-6.d0 - 5 * pi**2 * r**2 + &
258  (-6.d0 + 5.d0 * pi**2 * r**2) * cos(4*pi*z)) * sin(2*pi*z)**2 / &
259  lambda**2
260  ELSE IF (mode == 2) THEN
261  vv = vv + 2 * r**7 * (r - r0)**2 * cos(time)**4 * (13.d0 + 12 * pi**2 * r**2 + &
262  (13.d0 - 12 * pi**2 * r**2) * cos(4*pi*z)) * sin(2*pi*z)**2 / &
263  lambda**2
264  ELSE IF (mode == 3) THEN
265  vv = vv + 2 * r**7 * (r - r0)**2 * cos(time)**4 * (2.d0 + 3 * pi**2 * r**2 + &
266  (2.d0 - 3 * pi**2 * r**2) * cos(4*pi*z)) * sin(2*pi*z)**2 / &
267  lambda**2
268  ELSE IF (mode == 4) THEN
269  vv = vv - r**7 * (r - r0)**2 * cos(time)**4 * sin(4*pi*z)**2 / &
270  (2 * lambda**2)
271  END IF
272  ELSE IF (TYPE == 3) then
273  IF (mode == 0) THEN
274  vv = vv + r**7 * (r - r0)**2 * cos(time)**4 * (15.d0 + 8 * pi**2 * r**2 + &
275  (15.d0 - 8 * pi**2 * r**2) * cos(4*pi*z)) * sin(2*pi*z)**2 / &
276  (2 * lambda**2)
277  ELSE IF (mode == 1) THEN
278  vv = vv + r**7 * (r - r0)**2 * cos(time)**4 * (12.d0 + 7 * pi**2 * r**2 + &
279  (12.d0 - 7 * pi**2 * r**2) * cos(4*pi*z)) * sin(2*pi*z)**2 / &
280  lambda**2
281  ELSE IF (mode == 2) THEN
282  vv = vv + r**7 * (r - r0)**2 * cos(time)**4 * (5.d0 + 4 * pi**2 * r**2 + &
283  (5.d0 - 4 * pi**2 * r**2) * cos(4*pi*z)) * sin(2*pi*z)**2 / &
284  lambda**2
285  ELSE IF (mode == 3) THEN
286  vv = vv + 2 * pi**2 * r**9 * (r - r0)**2 * cos(time)**4 * sin(2*pi*z)**4 / &
287  lambda**2
288  ELSE IF (mode == 4) THEN
289  vv = vv - r**7 * (r - r0)**2 * cos(time)**4 * sin(4*pi*z)**2 / &
290  (4 * lambda**2)
291  END IF
292  ELSE IF (TYPE == 4) then
293  IF (mode == 1) THEN
294  vv = vv + r**7 * (r - r0)**2 * cos(time)**4 * (25.d0 + 8 * pi**2 * r**2 + &
295  (25.d0 - 8 * pi**2 * r**2) * cos(4*pi*z)) * sin(2*pi*z)**2 / &
296  (4 * lambda**2)
297  ELSE IF (mode == 2) THEN
298  vv = vv + r**7 * (r - r0)**2 * cos(time)**4 * (61.d0 + 24 * pi**2 * r**2 + &
299  (61.d0 - 24 * pi**2 * r**2) * cos(4*pi*z)) * sin(2*pi*z)**2 / &
300  (8 * lambda**2)
301  ELSE IF (mode == 3) THEN
302  vv = vv + r**7 * (r - r0)**2 * cos(time)**4 * (17.d0 + 8 * pi**2 * r**2 + &
303  (17.d0 - 8 * pi**2 * r**2) * cos(4*pi*z)) * sin(2*pi*z)**2 / &
304  (4 * lambda**2)
305  ELSE IF (mode == 4) THEN
306  vv = vv + r**7 * (r - r0)**2 * cos(time)**4 * (15.d0 + 8 * pi**2 * r**2 + &
307  (15.d0 - 8 * pi**2 * r**2) * cos(4*pi*z)) * sin(2*pi*z)**2 / &
308  (16 * lambda**2)
309  END IF
310  ELSE IF (TYPE == 5) then
311  IF (mode == 0) THEN
312  vv = vv + pi * r**8 * (-215.d0 + 136 * pi**2 * r**2) * (r - r0)**2 * cos(time)**4 * &
313  cos(2*pi*z) * sin(2*pi*z)**3 / (4 * lambda**2)
314  ELSE IF (mode == 1) THEN
315  vv = vv + 5 * pi * r**8 * (-11.d0 + 8 * pi**2 * r**2) * (r - r0)**2 * cos(time)**4 * &
316  cos(2*pi*z) * sin(2*pi*z)**3 / lambda**2
317  ELSE IF (mode == 2) THEN
318  vv = vv + 14 * pi * r**8 * (r - r0)**2 * cos(time)**4 * &
319  cos(2*pi*z) * sin(2*pi*z)**3 / lambda**2
320  ELSE IF (mode == 3) THEN
321  vv = vv - pi * r**8 * (-19.d0 + 8 * pi**2 * r**2) * (r - r0)**2 * cos(time)**4 * &
322  cos(2*pi*z) * sin(2*pi*z)**3 / lambda**2
323  ELSE IF (mode == 4) THEN
324  vv = vv - pi * r**8 * (-15.d0 + 8 * pi**2 * r**2) * (r - r0)**2 * cos(time)**4 * &
325  cos(2*pi*z) * sin(2*pi*z)**3 / (4 * lambda**2)
326  END IF
327  ELSE IF (TYPE == 6) then
328  IF (mode == 1) THEN
329  vv = vv + 8 * pi * r**8 * (-9.d0 + 5 * pi**2 * r**2) * (r - r0)**2 * cos(time)**4 * &
330  cos(2*pi*z) * sin(2*pi*z)**3 / lambda**2
331  ELSE IF (mode == 2) THEN
332  vv = vv + 4 * pi * r**8 * (-13.d0 + 8 * pi**2 * r**2) * (r - r0)**2 * cos(time)**4 * &
333  cos(2*pi*z) * sin(2*pi*z)**3 / lambda**2
334  ELSE IF (mode == 3) THEN
335  vv = vv + 8 * pi * r**8 * (-1.d0 + pi**2 * r**2) * (r - r0)**2 * cos(time)**4 * &
336  cos(2*pi*z) * sin(2*pi*z)**3 / lambda**2
337  ELSE IF (mode == 4) THEN
338  vv = vv + 2 * pi * r**8 * (r - r0)**2 * cos(time)**4 * &
339  cos(2*pi*z) * sin(2*pi*z)**3 / lambda**2
340  END IF
341  END IF
342 
343  RETURN
344 
345  !===Dummies variables to avoid warning
346  np=ty
347  !===Dummies variables to avoid warning
348  END FUNCTION source_in_ns_momentum
349 
350  !===Extra source in temperature equation. Always called.
351  FUNCTION source_in_temperature(TYPE, rr, m, t)RESULT(vv)
352  IMPLICIT NONE
353  INTEGER , INTENT(IN) :: TYPE
354  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: rr
355  INTEGER , INTENT(IN) :: m
356  REAL(KIND=8), INTENT(IN) :: t
357  REAL(KIND=8), DIMENSION(SIZE(rr,2)) :: vv, r, z, c, lambda
358  INTEGER :: i
359  REAL(KIND=8) :: r0 = 0.5d0, pi = acos(-1.d0)
360 
361  r = rr(1,:)
362  z = rr(2,:)
363 
364  DO i=1,SIZE(rr,2)
365  IF (r(i).LE.r0) THEN
366  c(i) = inputs%vol_heat_capacity(1)
367  ELSE
368  c(i) = inputs%vol_heat_capacity(2)
369  END IF
370  END DO
371 
372  DO i=1,SIZE(rr,2)
373  IF (r(i).LE.r0) THEN
374  lambda(i) = inputs%temperature_diffusivity(1)
375  ELSE
376  lambda(i) = inputs%temperature_diffusivity(2)
377  END IF
378  END DO
379 
380  IF (type==1) THEN
381  IF (m==0) THEN
382  vv = ((-9*r + 4*pi**2*r**3 + 4*r0 - 4*pi**2*r**2*r0) * lambda * cos(t) + &
383  c * r**2 * (-r + r0) * sin(t)) * sin(2*pi*z) / lambda
384  ELSE IF (m==1) THEN
385  vv = - ((-3*r0 + 4*r*(2.d0 + pi**2*r*(-r + r0))) * lambda * cos(t) + &
386  c * r**2 * (r - r0) * sin(t)) * sin(2*pi*z) / lambda
387  ELSE
388  vv = 0.d0
389  END IF
390  ELSE
391  vv = 0.d0
392  END IF
393 
394  IF (type==1) THEN
395  IF (m==0) THEN
396  DO i=1,size(rr,2)
397  IF (r(i)>r0) THEN
398  vv(i) = vv(i) + 3 * c(i) * pi * r(i) * (r(i) - r0)**2 * r0 * &
399  cos(t)**2 * sin(4*pi*z(i)) / (2 * lambda(i))
400  END IF
401  END DO
402  ELSE IF (m==1) THEN
403  DO i=1,size(rr,2)
404  IF (r(i)>r0) THEN
405  vv(i) = vv(i) + 2 * c(i) * pi * r(i) * (r(i) - r0)**2 * r0 * &
406  cos(t)**2 * sin(4*pi*z(i)) / lambda(i)
407  END IF
408  END DO
409  ELSE IF (m==2) THEN
410  DO i=1,size(rr,2)
411  IF (r(i)>r0) THEN
412  vv(i) = vv(i) + c(i) * pi * r(i) * (r(i) - r0)**2 * r0 * &
413  cos(t)**2 * sin(4*pi*z(i)) / (2 *lambda(i))
414  END IF
415  END DO
416  END IF
417  END IF
418 
419  RETURN
420  END FUNCTION source_in_temperature
421 
422 !!$ FUNCTION source_in_level_set(interface_nb,TYPE, rr, m, t)RESULT(vv)
423 !!$ IMPLICIT NONE
424 !!$ INTEGER , INTENT(IN) :: TYPE
425 !!$ REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: rr
426 !!$ INTEGER , INTENT(IN) :: m, interface_nb
427 !!$ REAL(KIND=8), INTENT(IN) :: t
428 !!$ REAL(KIND=8), DIMENSION(SIZE(rr,2)) :: vv
429 !!$
430 !!$ vv=0.d0
431 !!$ RETURN
432 !!$ END FUNCTION source_in_level_set
433 
434  !===Velocity for boundary conditions in Navier-Stokes.
435  !===Can be used also to initialize velocity in: init_velocity_pressure_temperature
436  FUNCTION vv_exact(TYPE,rr,m,t) RESULT(vv)
437  IMPLICIT NONE
438  INTEGER , INTENT(IN) :: TYPE
439  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: rr
440  INTEGER, INTENT(IN) :: m
441  REAL(KIND=8), INTENT(IN) :: t
442  REAL(KIND=8), DIMENSION(SIZE(rr,2)) :: vv, r, z
443  REAL(KIND=8) :: r0 = 0.5d0, pi = acos(-1.d0)
444 
445  r = rr(1,:)
446  z = rr(2,:)
447 
448  IF (type==1) THEN
449  IF ((m==0).OR.(m==1)) THEN
450  vv = -2*pi*(r-r0)**2*cos(t)*cos(2*pi*z)
451  ELSE
452  vv = 0.d0
453  END IF
454  ELSE IF (type==3) THEN
455  IF ((m==0).OR.(m==1)) THEN
456  vv = 2*pi*(r-r0)**2*cos(t)*cos(2*pi*z)
457  ELSE
458  vv = 0.d0
459  END IF
460  ELSE IF (type==5) THEN
461  IF ((m==0).OR.(m==1)) THEN
462  vv = (r-r0)*cos(t)*sin(2*pi*z)/r * (3*r-r0)
463  ELSE
464  vv = 0.d0
465  END IF
466  ELSE IF (type==6) THEN
467  IF (m==1) THEN
468  vv = (r-r0)*cos(t)*sin(2*pi*z)/r * (r-r0)
469  ELSE
470  vv = 0.d0
471  END IF
472  ELSE
473  vv = 0.d0
474  END IF
475 
476  RETURN
477  END FUNCTION vv_exact
478 
479 !!$ !===Solid velocity imposed when using penalty technique
480 !!$ !===Defined in Fourier space on mode 0 only.
481 !!$ FUNCTION imposed_velocity_by_penalty(rr,t) RESULT(vv)
482 !!$ IMPLICIT NONE
483 !!$ REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: rr
484 !!$ REAL(KIND=8), INTENT(IN) :: t
485 !!$ REAL(KIND=8), DIMENSION(SIZE(rr,2),6) :: vv
486 !!$
487 !!$ vv=0.d0
488 !!$ RETURN
489 !!$ END FUNCTION imposed_velocity_by_penalty
490 
491  !===Pressure for boundary conditions in Navier-Stokes.
492  !===Can be used also to initialize pressure in the subroutine init_velocity_pressure.
493  !===Use this routine for outflow BCs only.
494  !===CAUTION: Do not enfore BCs on pressure where normal component
495  ! of velocity is prescribed.
496  FUNCTION pp_exact(TYPE,rr,m,t) RESULT (vv)
497  IMPLICIT NONE
498  INTEGER , INTENT(IN) :: TYPE
499  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: rr
500  INTEGER , INTENT(IN) :: m
501  REAL(KIND=8), INTENT(IN) :: t
502  REAL(KIND=8), DIMENSION(SIZE(rr,2)) :: vv, r, z
503  REAL(KIND=8) :: pi = acos(-1.d0)
504 
505  r = rr(1,:)
506  z = rr(2,:)
507 
508  IF ((type==1).AND.(m==1)) THEN
509  vv = r**3*sin(2*pi*z)*cos(t)
510  ELSE
511  vv = 0.d0
512  END IF
513 
514  RETURN
515  END FUNCTION pp_exact
516 
517  !===Temperature for boundary conditions in temperature equation.
518  FUNCTION temperature_exact(TYPE,rr,m,t) RESULT (vv)
519  IMPLICIT NONE
520  INTEGER , INTENT(IN) :: TYPE
521  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: rr
522  INTEGER , INTENT(IN) :: m
523  REAL(KIND=8), INTENT(IN) :: t
524  REAL(KIND=8), DIMENSION(SIZE(rr,2)) :: vv, r, z, lambda
525  REAL(KIND=8) :: r0=0.5d0, pi = acos(-1.d0)
526  INTEGER :: i
527 
528  r = rr(1,:)
529  z = rr(2,:)
530 
531  DO i=1,SIZE(rr,2)
532  IF (r(i).LE.r0) THEN
533  lambda(i) = inputs%temperature_diffusivity(1)
534  ELSE
535  lambda(i) = inputs%temperature_diffusivity(2)
536  END IF
537  END DO
538 
539  IF ((type==1).AND.((m==0).OR.(m==1))) THEN
540  vv = r**2*(r-r0)*sin(2*pi*z)*cos(t) / lambda
541  ELSE
542  vv = 0.d0
543  END IF
544 
545  RETURN
546  END FUNCTION temperature_exact
547 
548 !!$ !===Can be used to initialize level set in the subroutine init_level_set.
549 !!$ FUNCTION level_set_exact(interface_nb,TYPE,rr,m,t) RESULT (vv)
550 !!$ IMPLICIT NONE
551 !!$ INTEGER , INTENT(IN) :: TYPE
552 !!$ REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: rr
553 !!$ INTEGER , INTENT(IN) :: m, interface_nb
554 !!$ REAL(KIND=8), INTENT(IN) :: t
555 !!$ REAL(KIND=8), DIMENSION(SIZE(rr,2)) :: vv
556 !!$
557 !!$ vv = 0.d0
558 !!$ RETURN
559 !!$ END FUNCTION level_set_exact
560 
561 !!$ !===Penalty coefficient (if needed)
562 !!$ !===This coefficient is equal to zero in subdomain
563 !!$ !===where penalty is applied
564 !!$ FUNCTION penal_in_real_space(mesh,rr_gauss,angles,nb_angles,nb,ne,time) RESULT(vv)
565 !!$ IMPLICIT NONE
566 !!$ TYPE(mesh_type) :: mesh
567 !!$ REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: rr_gauss
568 !!$ REAL(KIND=8), DIMENSION(:), INTENT(IN) :: angles
569 !!$ INTEGER, INTENT(IN) :: nb_angles
570 !!$ INTEGER, INTENT(IN) :: nb, ne
571 !!$ REAL(KIND=8), INTENT(IN) :: time
572 !!$ REAL(KIND=8), DIMENSION(nb_angles,ne-nb+1) :: vv
573 !!$
574 !!$ vv = 1.d0
575 !!$ RETURN
576 !!$ END FUNCTION penal_in_real_space
577 
578  !===Extension of the velocity field in the solid.
579  !===Used when temperature or Maxwell equations are solved.
580  !===It extends the velocity field on the Navier-Stokes domain to a
581  !===velocity field on the temperature and the Maxwell domain.
582  !===It is also used if problem type=mxw and restart velocity
583  !===is set to true in data (type problem denoted mxx in the code).
584  FUNCTION extension_velocity(TYPE, H_mesh, mode, t, n_start) RESULT(vv)
585  IMPLICIT NONE
586  TYPE(mesh_type), INTENT(IN) :: H_mesh
587  INTEGER , INTENT(IN) :: TYPE, n_start
588  INTEGER, INTENT(IN) :: mode
589  REAL(KIND=8), INTENT(IN) :: t
590  REAL(KIND=8), DIMENSION(H_Mesh%np) :: vv
591  REAL(KIND=8) :: r
592  INTEGER :: n
593 
594  vv = 0.d0
595  RETURN
596 
597  !===Dummies variables to avoid warning
598  n=h_mesh%np; r=t; n=type; n=mode; n=n_start
599  !===Dummies variables to avoid warning
600  END FUNCTION extension_velocity
601 
602  !===============================================================================
603  ! Boundary conditions for Maxwell
604  !===============================================================================
605 !!$ !===Velocity used in the induction equation.
606 !!$ !===Used only if problem type is mxw and restart velocity is false
607 !!$ FUNCTION Vexact(m, H_mesh) RESULT(vv) !Set uniquement a l'induction
608 !!$ IMPLICIT NONE
609 !!$ TYPE(mesh_type), INTENT(IN) :: H_mesh
610 !!$ INTEGER, INTENT(IN) :: m
611 !!$ REAL(KIND=8), DIMENSION(H_mesh%np,6) :: vv
612 !!$
613 !!$ vv = 0.d0
614 !!$ RETURN
615 !!$ END FUNCTION Vexact
616 
617 !!$ FUNCTION H_B_quasi_static(char_h_b, rr, m) RESULT(vv)
618 !!$ IMPLICIT NONE
619 !!$ CHARACTER(LEN=1), INTENT(IN) :: char_h_b
620 !!$ REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: rr
621 !!$ INTEGER, INTENT(IN) :: m
622 !!$ REAL(KIND=8), DIMENSION(SIZE(rr,2),6) :: vv
623 !!$
624 !!$ IF (inputs%if_quasi_static_approx) THEN
625 !!$ vv = 0.d0
626 !!$ ELSE
627 !!$ CALL error_petsc('H_B_quasi_static should not be called')
628 !!$ END IF
629 !!$ RETURN
630 !!$ END FUNCTION H_B_quasi_static
631 
632  !===Magnetic field for boundary conditions in the Maxwell equations.
633  FUNCTION hexact(H_mesh, TYPE, rr, m, mu_H_field, t) RESULT(vv)
634  IMPLICIT NONE
635  TYPE(mesh_type), INTENT(IN) :: H_mesh
636  INTEGER , INTENT(IN) :: TYPE
637  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: rr
638  INTEGER , INTENT(IN) :: m
639  REAL(KIND=8), INTENT(IN) :: t
640  REAL(KIND=8), DIMENSION(:), INTENT(IN) :: mu_H_field
641  REAL(KIND=8), DIMENSION(SIZE(rr,2)) :: vv, r, z
642  REAL(KIND=8) :: pi = acos(-1.d0)
643  INTEGER :: i
644 
645  r = rr(1,:)
646  z = rr(2,:)
647 
648  IF (TYPE == 1) then
649  IF (m == 0) THEN
650  vv = 2 * pi * r**3 * sin(2*pi*z) * cos(t)
651  ELSE
652  vv = 0.d0
653  END IF
654  ELSE IF (TYPE == 2) then
655  IF (m == 1) THEN
656  vv = 2 * pi * r**3 * sin(2*pi*z) * cos(t)
657  ELSE
658  vv = 0.d0
659  END IF
660  ELSE IF (TYPE == 3) then
661  IF (m == 0) THEN
662  vv = - 2 * pi * r**3 * sin(2*pi*z) * cos(t)
663  ELSE
664  vv = 0.d0
665  END IF
666  ELSE IF (TYPE == 4) then
667  IF (m == 1) THEN
668  vv = - 2 * pi * r**3 * sin(2*pi*z) * cos(t)
669  ELSE
670  vv = 0.d0
671  END IF
672  ELSE IF (TYPE == 5) then
673  IF (m == 0) THEN
674  vv = 4 * r**2 * cos(2*pi*z) * cos(t)
675  ELSE IF (m == 1) THEN
676  vv = - r**2 * cos(2*pi*z) * cos(t)
677  ELSE
678  vv = 0.d0
679  END IF
680  ELSE
681  IF (m == 1) THEN
682  vv = 4 * r**2 * cos(2*pi*z) * cos(t)
683  ELSE
684  vv = 0.d0
685  END IF
686  END IF
687  RETURN
688 
689  !===Dummies variables to avoid warning
690  i=h_mesh%np; i=SIZE(mu_h_field);
691  !===Dummies variables to avoid warning
692  END FUNCTION hexact
693 
694  !===Scalar potential for boundary conditions in the Maxwell equations.
695  FUNCTION phiexact(TYPE, rr, m, mu_phi,t) RESULT(vv)
696  IMPLICIT NONE
697  INTEGER , INTENT(IN) :: TYPE
698  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: rr
699  INTEGER , INTENT(IN) :: m
700  REAL(KIND=8), INTENT(IN) :: mu_phi, t
701  REAL(KIND=8), DIMENSION(SIZE(rr,2)) :: vv
702  REAL(KIND=8) :: r
703  INTEGER :: n
704 
705  vv=0.d0
706  RETURN
707 
708  !===Dummies variables to avoid warning
709  n=type; n=SIZE(rr,1); n=m; r=mu_phi; r=t
710  !===Dummies variables to avoid warning
711  END FUNCTION phiexact
712 
713  !===Current in Ohm's law. Curl(H) = sigma(E + uxB) + current
714  FUNCTION jexact_gauss(TYPE, rr, m, mu_phi, sigma, mu_H, t, mesh_id, opt_B_ext) RESULT(vv)
715  IMPLICIT NONE
716  INTEGER , INTENT(IN) :: TYPE
717  REAL(KIND=8), DIMENSION(:), INTENT(IN) :: rr
718  INTEGER , INTENT(IN) :: m
719  REAL(KIND=8), INTENT(IN) :: mu_phi, sigma, mu_H, t
720  INTEGER , INTENT(IN) :: mesh_id
721  REAL(KIND=8), DIMENSION(6), OPTIONAL,INTENT(IN) :: opt_B_ext
722  REAL(KIND=8) :: vv, r, z
723  REAL(KIND=8) :: pi = acos(-1.d0)
724  INTEGER :: n
725 
726  r = rr(1)
727  z = rr(2)
728 
729  IF (TYPE == 1) then
730  IF (m == 0) THEN
731  vv = 4 * pi**2 * r**3 * cos(2*pi*z) * cos(t)
732  ELSE IF (m == 1) THEN
733  vv = 4 * r * cos(2*pi*z) * cos(t)
734  ELSE
735  vv = 0.d0
736  END IF
737  ELSE IF (TYPE == 2) then
738  IF (m == 1) THEN
739  vv = r * (1.d0 + 4 * pi**2 * r**2) * cos(2*pi*z) * cos(t)
740  ELSE
741  vv = 0.d0
742  END IF
743  ELSE IF (TYPE == 3) then
744  IF (m == 0) THEN
745  vv = 4 * r * (-2.d0 + pi**2 * r**2) * cos(2*pi*z) * cos(t)
746  ELSE IF (m == 1) THEN
747  vv = 2 * r * cos(2*pi*z) * cos(t)
748  ELSE
749  vv = 0.d0
750  END IF
751  ELSE IF (TYPE == 4) then
752  IF (m == 1) THEN
753  vv = 4 * r * (-2.d0 + pi**2 * r**2) * cos(2*pi*z) * cos(t)
754  ELSE
755  vv = 0.d0
756  END IF
757  ELSE IF (TYPE == 5) then
758  IF (m == 0) THEN
759  vv = - 8 * pi * r**2 * sin(2*pi*z) * cos(t)
760  ELSE IF (m == 1) THEN
761  vv = - 2 * pi * r**2 * sin(2*pi*z) * cos(t)
762  ELSE
763  vv = 0.d0
764  END IF
765  ELSE
766  IF (m == 1) THEN
767  vv = - 8 * pi * r**2 * sin(2*pi*z) * cos(t)
768  ELSE
769  vv = 0.d0
770  END IF
771  END IF
772  RETURN
773 
774  !===Dummies variables to avoid warning
775  r=mu_phi; r=sigma; r=mu_h; n=mesh_id
776  IF (PRESENT(opt_b_ext)) r=opt_b_ext(1)
777  !===Dummies variables to avoid warning
778  END FUNCTION jexact_gauss
779 
780  !===Electric field for Neumann BC (cf. doc)
781  FUNCTION eexact_gauss(TYPE, rr, m, mu_phi, sigma, mu_H, t) RESULT(vv)
782  IMPLICIT NONE
783  INTEGER, INTENT(IN) :: TYPE
784  REAL(KIND=8), DIMENSION(:), INTENT(IN) :: rr
785  INTEGER, INTENT(IN) :: m
786  REAL(KIND=8), INTENT(IN) :: mu_phi, sigma, mu_H, t
787  REAL(KIND=8) :: vv
788  REAL(KIND=8) :: r
789  INTEGER :: n
790 
791  vv = 0.d0
792  RETURN
793 
794  !===Dummies variables to avoid warning
795  r=rr(1); r=mu_phi; r=sigma; r=mu_h; r=t; n=type; n=m
796  !===Dummies variables to avoid warning
797  END FUNCTION eexact_gauss
798 
799  !===Initialization of magnetic field and scalar potential (if present)
800  SUBROUTINE init_maxwell(H_mesh, phi_mesh, time, dt, mu_H_field, mu_phi, &
801  list_mode, hn1, hn, phin1, phin)
802  IMPLICIT NONE
803  TYPE(mesh_type) :: H_mesh, phi_mesh
804  REAL(KIND=8), INTENT(OUT):: time
805  REAL(KIND=8), INTENT(IN) :: dt
806  REAL(KIND=8), DIMENSION(:), INTENT(IN) :: mu_H_field
807  REAL(KIND=8), INTENT(IN) :: mu_phi
808  INTEGER, DIMENSION(:), INTENT(IN) :: list_mode
809  REAL(KIND=8), DIMENSION(:,:,:), INTENT(OUT):: Hn, Hn1
810  REAL(KIND=8), DIMENSION(:,:,:), INTENT(OUT):: phin, phin1
811  INTEGER :: i, k
812 
813  time = -dt
814  DO k=1,6
815  DO i=1, SIZE(list_mode)
816  hn1(:,k,i) = hexact(h_mesh,k, h_mesh%rr, list_mode(i), mu_h_field, time)
817  IF (inputs%nb_dom_phi>0) THEN
818  IF (k<3) THEN
819  phin1(:,k,i) = phiexact(k, phi_mesh%rr, list_mode(i) , mu_phi, time)
820  ENDIF
821  END IF
822  ENDDO
823  ENDDO
824 
825  time = time + dt
826  DO k=1,6
827  DO i=1, SIZE(list_mode)
828  hn(:,k,i) = hexact(h_mesh,k, h_mesh%rr, list_mode(i), mu_h_field, time)
829  IF (inputs%nb_dom_phi>0) THEN
830  IF (k<3) THEN
831  phin(:,k,i) = phiexact(k, phi_mesh%rr, list_mode(i), mu_phi, time)
832  ENDIF
833  END IF
834  ENDDO
835  ENDDO
836  RETURN
837  END SUBROUTINE init_maxwell
838 
839 !!$ !===Analytical permeability (if needed)
840 !!$ !===This function is not needed unless the flag
841 !!$ !=== ===Use FEM Interpolation for magnetic permeability (true/false)
842 !!$ !===is activated and set to .FALSE. in the data data file. Default is .TRUE.
843 !!$ FUNCTION mu_bar_in_fourier_space(H_mesh,nb,ne,pts,pts_ids) RESULT(vv)
844 !!$ IMPLICIT NONE
845 !!$ TYPE(mesh_type), INTENT(IN) :: H_mesh
846 !!$ REAL(KIND=8), DIMENSION(ne-nb+1) :: vv
847 !!$ INTEGER, INTENT(IN) :: nb, ne
848 !!$ REAL(KIND=8),DIMENSION(2,ne-nb+1),OPTIONAL :: pts
849 !!$ INTEGER, DIMENSION(ne-nb+1), OPTIONAL :: pts_ids
850 !!$
851 !!$ vv = 1.d0
852 !!$ RETURN
853 !!$ END FUNCTION mu_bar_in_fourier_space
854 
855 !!$ !===Analytical mu_in_fourier_space (if needed)
856 !!$ !===This function is not needed unless the flag
857 !!$ !=== ===Use FEM Interpolation for magnetic permeability (true/false)
858 !!$ !===is activated and set to .FALSE. in the data data file. Default is .TRUE.
859 !!$ FUNCTION grad_mu_bar_in_fourier_space(pt,pt_id) RESULT(vv)
860 !!$ IMPLICIT NONE
861 !!$ REAL(KIND=8),DIMENSION(2), INTENT(in):: pt
862 !!$ INTEGER,DIMENSION(1), INTENT(in) :: pt_id
863 !!$ REAL(KIND=8),DIMENSION(2) :: vv
864 !!$
865 !!$ vv = 0.d0
866 !!$ RETURN
867 !!$ END FUNCTION grad_mu_bar_in_fourier_space
868 
869 !!$ !===Analytical permeability, mu in real space (if needed)
870 !!$ FUNCTION mu_in_real_space(H_mesh,angles,nb_angles,nb,ne,time) RESULT(vv)
871 !!$ IMPLICIT NONE
872 !!$ TYPE(mesh_type), INTENT(IN) :: H_mesh
873 !!$ REAL(KIND=8), DIMENSION(:), INTENT(IN) :: angles
874 !!$ INTEGER, INTENT(IN) :: nb_angles
875 !!$ INTEGER, INTENT(IN) :: nb, ne
876 !!$ REAL(KIND=8), INTENT(IN) :: time
877 !!$ REAL(KIND=8), DIMENSION(nb_angles,ne-nb+1) :: vv
878 !!$
879 !!$ vv = 1.d0
880 !!$ RETURN
881 !!$ END FUNCTION mu_in_real_space
882 
883  !===Coefficient contaning the magnetic susceptibility for magnetic force in ferrofluids:
884  !===F = chi_coeff(T) * grad(H**2/2) (Kelvin form)
885  !===or F = -H**2/2 * grad(chi_coeff(T)) (Helmholtz form)
886  FUNCTION chi_coeff_law(temp) RESULT(vv)
887  IMPLICIT NONE
888  REAL(KIND=8) :: temp
889  REAL(KIND=8) :: vv
890 
891  vv = - temp**2
892  RETURN
893  END FUNCTION chi_coeff_law
894 
895  !===Coefficient contaning the temperature dependant factor in the
896  !===pyromagnetic coefficient term of the temperature equation for ferrofluids:
897  !===T * dchi/dT(T) * D/Dt(H**2/2)
898  FUNCTION t_dchi_dt_coeff_law(temp) RESULT(vv)
899  IMPLICIT NONE
900  REAL(KIND=8) :: temp
901  REAL(KIND=8) :: vv
902 
903  vv = 0.d0
904  RETURN
905 
906  !===Dummies variables to avoid warning
907  vv=temp
908  !===Dummies variables to avoid warning
909  END FUNCTION t_dchi_dt_coeff_law
910 
911 END MODULE boundary_test_34
boundary_test_34
Definition: condlim_test_34.f90:1
boundary_test_34::init_temperature
subroutine, public init_temperature(mesh, time, dt, list_mode, tempn_m1, tempn)
Definition: condlim_test_34.f90:75
boundary_test_34::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_34.f90:119
boundary_test_34::chi_coeff_law
real(kind=8) function, public chi_coeff_law(temp)
Definition: condlim_test_34.f90:887
boundary_test_34::eexact_gauss
real(kind=8) function, public eexact_gauss(TYPE, rr, m, mu_phi, sigma, mu_H, t)
Definition: condlim_test_34.f90:782
def_type_mesh::mesh_type
Definition: def_type_mesh.f90:74
boundary_test_34::extension_velocity
real(kind=8) function, dimension(h_mesh%np), public extension_velocity(TYPE, H_mesh, mode, t, n_start)
Definition: condlim_test_34.f90:585
my_util
Definition: my_util.f90:1
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_34::t_dchi_dt_coeff_law
real(kind=8) function, public t_dchi_dt_coeff_law(temp)
Definition: condlim_test_34.f90:899
boundary_test_34::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_34.f90:715
boundary_test_34::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_34.f90:802
boundary_test_34::vv_exact
real(kind=8) function, dimension(size(rr, 2)), public vv_exact(TYPE, rr, m, t)
Definition: condlim_test_34.f90:437
boundary_test_34::phiexact
real(kind=8) function, dimension(size(rr, 2)), public phiexact(TYPE, rr, m, mu_phi, t)
Definition: condlim_test_34.f90:696
boundary_test_34::hexact
real(kind=8) function, dimension(size(rr, 2)), public hexact(H_mesh, TYPE, rr, m, mu_H_field, t)
Definition: condlim_test_34.f90:634
boundary_test_34::temperature_exact
real(kind=8) function, dimension(size(rr, 2)), public temperature_exact(TYPE, rr, m, t)
Definition: condlim_test_34.f90:519
boundary_test_34::source_in_temperature
real(kind=8) function, dimension(size(rr, 2)), public source_in_temperature(TYPE, rr, m, t)
Definition: condlim_test_34.f90:352
input_data
Definition: read_my_data.f90:233
def_type_mesh
Definition: def_type_mesh.f90:4
boundary_test_34::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_34.f90:44
boundary_test_34::pp_exact
real(kind=8) function, dimension(size(rr, 2)), public pp_exact(TYPE, rr, m, t)
Definition: condlim_test_34.f90:497