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