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