SFEMaNS  version 5.3
Reference documentation for SFEMaNS
condlim_test_32.f90
Go to the documentation of this file.
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  PRIVATE
33 
34 CONTAINS
35  !===============================================================================
36  ! Boundary conditions for Navier-Stokes
37  !===============================================================================
38 
39  !===Initialize velocity, pressure
40  SUBROUTINE init_velocity_pressure(mesh_f, mesh_c, time, dt, list_mode, &
41  un_m1, un, pn_m1, pn, phin_m1, phin)
42  IMPLICIT NONE
43  TYPE(mesh_type) :: mesh_f, mesh_c
44  REAL(KIND=8), INTENT(OUT):: time
45  REAL(KIND=8), INTENT(IN) :: dt
46  INTEGER, DIMENSION(:), INTENT(IN) :: list_mode
47  REAL(KIND=8), DIMENSION(:,:,:), INTENT(OUT):: un_m1, un
48  REAL(KIND=8), DIMENSION(:,:,:), INTENT(OUT):: pn_m1, pn, phin_m1, phin
49  INTEGER :: mode, i, j
50  REAL(KIND=8), DIMENSION(mesh_c%np) :: pn_m2
51 
52  time = 0.d0
53  DO i= 1, SIZE(list_mode)
54  mode = list_mode(i)
55  DO j = 1, 6
56  !===velocity
57  un_m1(:,j,i) = vv_exact(j,mesh_f%rr,mode,time-dt)
58  un(:,j,i) = vv_exact(j,mesh_f%rr,mode,time)
59  END DO
60  DO j = 1, 2
61  !===pressure
62  pn_m2(:) = pp_exact(j,mesh_c%rr,mode,time-2*dt)
63  pn_m1(:,j,i) = pp_exact(j,mesh_c%rr,mode,time-dt)
64  pn(:,j,i) = pp_exact(j,mesh_c%rr,mode,time)
65  phin_m1(:,j,i) = pn_m1(:,j,i) - pn_m2(:)
66  phin(:,j,i) = pn(:,j,i) - pn_m1(:,j,i)
67  ENDDO
68  ENDDO
69  END SUBROUTINE init_velocity_pressure
70 
71  !===Initialize temperature
72  SUBROUTINE init_temperature(mesh, time, dt, list_mode, tempn_m1, tempn)
73  IMPLICIT NONE
74  TYPE(mesh_type) :: mesh
75  REAL(KIND=8), INTENT(OUT):: time
76  REAL(KIND=8), INTENT(IN) :: dt
77  INTEGER, DIMENSION(:), INTENT(IN) :: list_mode
78  REAL(KIND=8), DIMENSION(:,:,:), INTENT(OUT):: tempn_m1, tempn
79  INTEGER :: mode, i, j
80 
81  time = 0.d0
82  DO i= 1, SIZE(list_mode)
83  mode = list_mode(i)
84  DO j = 1, 2
85  tempn_m1(:,j,i) = temperature_exact(j, mesh%rr, mode, time-dt)
86  tempn(:,j,i) = temperature_exact(j, mesh%rr, mode, time)
87  ENDDO
88  ENDDO
89  END SUBROUTINE init_temperature
90 
91 !!$ !===Initialize level_set
92 !!$ SUBROUTINE init_level_set(vv_mesh, time, &
93 !!$ dt, list_mode, level_set_m1, level_set)
94 !!$ IMPLICIT NONE
95 !!$ TYPE(mesh_type) :: vv_mesh
96 !!$ REAL(KIND=8), INTENT(OUT):: time
97 !!$ REAL(KIND=8), INTENT(IN) :: dt
98 !!$ INTEGER, DIMENSION(:), INTENT(IN) :: list_mode
99 !!$ REAL(KIND=8), DIMENSION(:,:,:,:), INTENT(OUT):: level_set, level_set_m1
100 !!$ INTEGER :: mode, i, j, n
101 !!$ time = 0.d0
102 !!$ DO i= 1, SIZE(list_mode)
103 !!$ mode = list_mode(i)
104 !!$ DO j = 1, 2
105 !!$ !===level_set
106 !!$ DO n = 1, inputs%nb_fluid -1
107 !!$ level_set_m1(n,:,j,i) = level_set_exact(n,j,vv_mesh%rr,mode,time-dt)
108 !!$ level_set (n,:,j,i) = level_set_exact(n,j,vv_mesh%rr,mode,time)
109 !!$ END DO
110 !!$ END DO
111 !!$ END DO
112 !!$ END SUBROUTINE init_level_set
113 
114  !===Source in momemtum equation. Always called.
115  FUNCTION source_in_ns_momentum(TYPE, rr, mode, i, time, Re, ty, opt_density, opt_tempn) RESULT(vv)
116  IMPLICIT NONE
117  INTEGER , INTENT(IN) :: TYPE
118  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: rr
119  INTEGER , INTENT(IN) :: mode, i
120  REAL(KIND=8), INTENT(IN) :: time
121  REAL(KIND=8), INTENT(IN) :: Re
122  CHARACTER(LEN=2), INTENT(IN) :: ty
123  REAL(KIND=8), DIMENSION(:,:,:), OPTIONAL, INTENT(IN) :: opt_density
124  REAL(KIND=8), DIMENSION(:,:,:), OPTIONAL, INTENT(IN) :: opt_tempn
125  REAL(KIND=8), DIMENSION(SIZE(rr,2)) :: vv, r, z
126  REAL(KIND=8) :: alpha, r0 = 0.5d0
127  REAL(KIND=8) :: pi = 3.1415926535897932d0
128  CHARACTER(LEN=2) :: np
129 
130  IF (PRESENT(opt_density)) CALL error_petsc('density should not be present for test 32')
131 
132  alpha = inputs%gravity_coefficient
133  r = rr(1,:)
134  z = rr(2,:)
135 
136  IF (type==5) THEN
137  vv = alpha*(opt_tempn(:,1,i) - temperature_exact(1,rr,mode,time))
138  ELSE IF (type==6) THEN
139  vv = alpha*(opt_tempn(:,2,i) - temperature_exact(2,rr,mode,time))
140  ELSE
141  vv = 0.d0
142  END IF
143 
144  IF (type==1) THEN
145  IF (mode==0) THEN
146  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) + &
147  (r - r0)*re*cos(time)**2*(14*r**3 - 5*r**2*r0 - 5*r*r0**2 + 2*r0**3 + &
148  (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) + &
149  r**2*(5*r0 - 12*pi**2*r0**3))*cos(4*pi*z)) - &
150  4*pi*r**3*(r - r0)**2*re*cos(2*pi*z)*sin(time))/(2.*r**3*re)
151  ELSE IF (mode==1) THEN
152  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) + &
153  ((3*r**2 - 4*r*r0 + r0**2)*re*cos(time)**2*(3*r**2 - r0**2 + &
154  (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. - &
155  pi*r**3*(r - r0)**2*re*cos(2*pi*z)*sin(time))/(r**3*re)
156  ELSE IF (mode==2) THEN
157  vv = vv -((r - r0)*cos(time)**2*(4*r**2 - r*r0 - r0**2 + (12*pi**2*r**4 - &
158  28*pi**2*r**3*r0 + r0**2 + 4*r**2*(-1.d0 + 5*pi**2*r0**2) + &
159  r*(r0 - 4*pi**2*r0**3))*cos(4*pi*z)))/(2.*r**2)
160  END IF
161  ELSE IF (type==2) THEN
162  IF (mode==1) THEN
163  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 + &
164  r**2*(-3.d0 + 4*pi**2*r0**2) + (3*r**2 + r*r0 - r0**2)*cos(4*pi*z))))/(r**3*re)
165  ELSE IF (mode==2) THEN
166  vv = vv + ((r - r0)**2*cos(time)**2*(4*pi**2*r**4 - r*r0 - 8*pi**2*r**3*r0 + r0**2 + &
167  r**2*(-3.d0 + 4*pi**2*r0**2) + (3*r**2 + r*r0 - r0**2)*cos(4*pi*z)))/(2.*r**3)
168  END IF
169  ELSE IF (type==3) THEN
170  IF (mode==0) THEN
171  vv = vv + (2*pi*(-3*pi*r*(r - r0)**3*(3*r - r0)*re*cos(time)**2 + &
172  (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) - &
173  r**2*(r - r0)**2*re*cos(2*pi*z)*sin(time)))/(r**2*re)
174  ELSE IF (mode==1) THEN
175  vv = vv + (4*pi*r*(2*pi**2*r**4 - r*r0 - 4*pi**2*r**3*r0 + r0**2 + &
176  r**2*(-1 + 2*pi**2*r0**2))*cos(time)*cos(2*pi*z) + &
177  ((r - r0)**3*(3*r - r0)*re*cos(time)**2*(-1.d0 - 16*pi**2*r**2 + cos(4*pi*z)))/2. - &
178  2*pi*r**3*(r - r0)**2*re*cos(2*pi*z)*sin(time))/(r**3*re)
179  ELSE IF (mode==2) THEN
180  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)
181  END IF
182  ELSE IF (type==4) THEN
183  IF (mode==1) THEN
184  vv = vv + ((r - r0)**2*cos(time)*(-8*pi*r*cos(2*pi*z) + re*cos(time)*((-3*r + r0)**2 + &
185  (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)
186  ELSE IF (mode==2) THEN
187  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
188  END IF
189  ELSE IF (type==5) THEN
190  IF (mode==0) THEN
191  vv = vv + ((4*(12*pi**2*r**4 - 16*pi**2*r**3*r0 - r0**2 + &
192  r**2*(-3.d0 + 4*pi**2*r0**2))*cos(time)*sin(2*pi*z) - &
193  4*r**2*(3*r**2 - 4*r*r0 + r0**2)*re*sin(time)*sin(2*pi*z) + &
194  pi*r*(r - r0)**2*(12*pi**2*r**4 - r*r0 - 24*pi**2*r**3*r0 + 2*r0**2 + &
195  4*r**2*(-1.d0 + 3*pi**2*r0**2))*re*4*cos(time)**2*sin(4*pi*z)))/(4.*r**3*re)
196  ELSE IF (mode==1) THEN
197  vv = vv + (4*(-r0 + pi**2*r*(3*r**2 - 4*r*r0 + r0**2))*cos(time)*sin(2*pi*z) - &
198  r*(3*r**2 - 4*r*r0 + r0**2)*re*sin(time)*sin(2*pi*z) + &
199  pi*(r - r0)**2*(16*pi**2*r**4 - 2*r*r0 - 32*pi**2*r**3*r0 + 3*r0**2 + &
200  r**2*(-5.d0 + 16*pi**2*r0**2))*re*cos(time)**2*sin(4*pi*z))/(r**2*re)
201  ELSE IF (mode==2) THEN
202  vv = vv + (pi*(r - r0)**2*(4*pi**2*r**4 - r*r0 - 8*pi**2*r**3*r0 + r0**2 + &
203  r**2*(-1.d0 + 4*pi**2*r0**2))*cos(time)**2*sin(4*pi*z))/r**2
204  END IF
205  ELSE IF (type==6) THEN
206  IF (mode==1) THEN
207  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) -&
208  4*pi*r*(r - r0)**3*re*cos(time)**2*sin(4*pi*z))/(r**2*re)
209  ELSE IF (mode==2) THEN
210  vv = vv -2*pi*(r - r0)**3*cos(time)**2*sin(4*pi*z)/r
211  END IF
212  END IF
213 
214  IF ((type==1).AND.(mode==1)) THEN
215  vv = vv + 3*r**2*cos(time)*sin(2*pi*z)
216  ELSE IF ((type==4).AND.(mode==1)) THEN
217  vv = vv - r**2*cos(time)*sin(2*pi*z)
218  ELSE IF ((type==5).AND.(mode==1)) THEN
219  vv = vv + 2*pi*r**3*cos(time)*cos(2*pi*z)
220  END IF
221 
222  RETURN
223 
224  !===Dummies variables to avoid warning
225  np=ty
226  !===Dummies variables to avoid warning
227  END FUNCTION source_in_ns_momentum
228 
229  !===Extra source in temperature equation. Always called.
230  FUNCTION source_in_temperature(TYPE, rr, m, t)RESULT(vv)
231  IMPLICIT NONE
232  INTEGER , INTENT(IN) :: TYPE
233  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: rr
234  INTEGER , INTENT(IN) :: m
235  REAL(KIND=8), INTENT(IN) :: t
236  REAL(KIND=8), DIMENSION(SIZE(rr,2)) :: vv, r, z, kappa
237  INTEGER :: i
238  REAL(KIND=8) :: r0 = 0.5d0
239  REAL(KIND=8) :: pi = 3.1415926535897932d0
240 
241  r = rr(1,:)
242  z = rr(2,:)
243 
244  DO i=1,SIZE(rr,2)
245  IF (r(i).LE.r0) THEN
246  kappa(i) = inputs%temperature_diffusivity(1)
247  ELSE
248  kappa(i) = inputs%temperature_diffusivity(2)
249  END IF
250  END DO
251 
252  IF (type==1) THEN
253  IF (m==0) THEN
254  vv = - (-2*((2*pi**2*r**4 + 9*r*r0 - 4*pi**2*r**3*r0 - 2*r0**2 + 2*r**2*(-4.d0 + pi**2*r0**2))*kappa*cos(t)) &
255  + r**2*(r - r0)**2*sin(t))*sin(2*pi*z)
256  ELSE IF (m==1) THEN
257  vv = -(-((4*pi**2*r**4 + 16*r*r0 - 8*pi**2*r**3*r0 - 3*r0**2 + r**2*(-15.d0 + 4*pi**2*r0**2))*kappa*cos(t)) &
258  + r**2*(r - r0)**2*sin(t))*sin(2*pi*z)
259  ELSE
260  vv = 0.d0
261  END IF
262  ELSE
263  vv = 0.d0
264  END IF
265 
266  IF (type==1) THEN
267  IF (m==0) THEN
268  DO i=1,size(rr,2)
269  IF (r(i)>r0) THEN
270  vv(i) = vv(i) + (-3*pi*r(i)*(r(i) - r0)**4*cos(t)**2*sin(4*pi*z(i)))/2.d0
271  END IF
272  END DO
273  ELSE IF (m==1) THEN
274  DO i=1,size(rr,2)
275  IF (r(i)>r0) THEN
276  vv(i) = vv(i) - 2*pi*r(i)*(r(i) - r0)**4*cos(t)**2*sin(4*pi*z(i))
277  END IF
278  END DO
279  ELSE IF (m==2) THEN
280  DO i=1,size(rr,2)
281  IF (r(i)>r0) THEN
282  vv(i) = vv(i) - 0.5*pi*(r(i)*(r(i) - r0)**4*cos(t)**2*sin(4*pi*z(i)))
283  END IF
284  END DO
285  END IF
286  END IF
287 
288  RETURN
289  END FUNCTION source_in_temperature
290 
291 !!$ !===Extra source in level set equation. Always called.
292 !!$ FUNCTION source_in_level_set(interface_nb,TYPE, rr, m, t)RESULT(vv)
293 !!$ IMPLICIT NONE
294 !!$ INTEGER , INTENT(IN) :: TYPE
295 !!$ REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: rr
296 !!$ INTEGER , INTENT(IN) :: m, interface_nb
297 !!$ REAL(KIND=8), INTENT(IN) :: t
298 !!$ REAL(KIND=8), DIMENSION(SIZE(rr,2)) :: vv
299 !!$
300 !!$ vv=0.d0
301 !!$ RETURN
302 !!$ END FUNCTION source_in_level_set
303 
304  !===Velocity for boundary conditions in Navier-Stokes.
305  !===Can be used also to initialize velocity in: init_velocity_pressure_temperature
306  FUNCTION vv_exact(TYPE,rr,m,t) RESULT(vv)
307  IMPLICIT NONE
308  INTEGER , INTENT(IN) :: TYPE
309  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: rr
310  INTEGER, INTENT(IN) :: m
311  REAL(KIND=8), INTENT(IN) :: t
312  REAL(KIND=8), DIMENSION(SIZE(rr,2)) :: vv, r, z
313  REAL(KIND=8) :: r0 = 0.5d0
314  REAL(KIND=8) :: pi = 3.1415926535897932d0
315 
316  r = rr(1,:)
317  z = rr(2,:)
318 
319  IF (type==1) THEN
320  IF ((m==0).OR.(m==1)) THEN
321  vv = -2*pi*(r-r0)**2*cos(t)*cos(2*pi*z)
322  ELSE
323  vv = 0.d0
324  END IF
325  ELSE IF (type==3) THEN
326  IF ((m==0).OR.(m==1)) THEN
327  vv = 2*pi*(r-r0)**2*cos(t)*cos(2*pi*z)
328  ELSE
329  vv = 0.d0
330  END IF
331  ELSE IF (type==5) THEN
332  IF ((m==0).OR.(m==1)) THEN
333  vv = (r-r0)*cos(t)*sin(2*pi*z)/r * (3*r-r0)
334  ELSE
335  vv = 0.d0
336  END IF
337  ELSE IF (type==6) THEN
338  IF (m==1) THEN
339  vv = (r-r0)*cos(t)*sin(2*pi*z)/r * (r-r0)
340  ELSE
341  vv = 0.d0
342  END IF
343  ELSE
344  vv = 0.d0
345  END IF
346 
347  RETURN
348  END FUNCTION vv_exact
349 
350 !!$ !===Solid velocity imposed when using penalty technique
351 !!$ !===Defined in Fourier space on mode 0 only.
352 !!$ FUNCTION imposed_velocity_by_penalty(rr,t) RESULT(vv)
353 !!$ IMPLICIT NONE
354 !!$ REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: rr
355 !!$ REAL(KIND=8), INTENT(IN) :: t
356 !!$ REAL(KIND=8), DIMENSION(SIZE(rr,2),6) :: vv
357 !!$
358 !!$ vv=0.d0
359 !!$ RETURN
360 !!$ END FUNCTION imposed_velocity_by_penalty
361 
362  !===Pressure for boundary conditions in Navier-Stokes.
363  !===Can be used also to initialize pressure in the subroutine init_velocity_pressure.
364  !===Use this routine for outflow BCs only.
365  !===CAUTION: Do not enfore BCs on pressure where normal component
366  ! of velocity is prescribed.
367  FUNCTION pp_exact(TYPE,rr,m,t) RESULT (vv)
368  IMPLICIT NONE
369  INTEGER , INTENT(IN) :: TYPE
370  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: rr
371  INTEGER , INTENT(IN) :: m
372  REAL(KIND=8), INTENT(IN) :: t
373  REAL(KIND=8), DIMENSION(SIZE(rr,2)) :: vv, r, z
374  REAL(KIND=8) :: pi = 3.1415926535897932d0
375 
376  r = rr(1,:)
377  z = rr(2,:)
378 
379  IF ((type==1).AND.(m==1)) THEN
380  vv = r**3*sin(2*pi*z)*cos(t)
381  ELSE
382  vv = 0.d0
383  END IF
384 
385  RETURN
386  END FUNCTION pp_exact
387 
388  !===Temperature for boundary conditions in temperature equation.
389  FUNCTION temperature_exact(TYPE,rr,m,t) RESULT (vv)
390  IMPLICIT NONE
391  INTEGER , INTENT(IN) :: TYPE
392  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: rr
393  INTEGER , INTENT(IN) :: m
394  REAL(KIND=8), INTENT(IN) :: t
395  REAL(KIND=8), DIMENSION(SIZE(rr,2)) :: vv, r, z
396  REAL(KIND=8) :: r0=0.5d0
397  REAL(KIND=8) :: pi = 3.1415926535897932d0
398 
399  r = rr(1,:)
400  z = rr(2,:)
401 
402  IF ((type==1).AND.((m==0).OR.(m==1))) THEN
403  vv = r**2*(r-r0)**2*sin(2*pi*z)*cos(t)
404  ELSE
405  vv = 0.d0
406  END IF
407 
408  RETURN
409  END FUNCTION temperature_exact
410 
411 !!$ !===Can be used to initialize level set in the subroutine init_level_set.
412 !!$ FUNCTION level_set_exact(interface_nb,TYPE,rr,m,t) RESULT (vv)
413 !!$ IMPLICIT NONE
414 !!$ INTEGER , INTENT(IN) :: TYPE
415 !!$ REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: rr
416 !!$ INTEGER , INTENT(IN) :: m, interface_nb
417 !!$ REAL(KIND=8), INTENT(IN) :: t
418 !!$ REAL(KIND=8), DIMENSION(SIZE(rr,2)) :: vv
419 !!$
420 !!$ vv = 0.d0
421 !!$ RETURN
422 !!$ END FUNCTION level_set_exact
423 
424 !!$ !===Penalty coefficient (if needed)
425 !!$ !===This coefficient is equal to zero in subdomain
426 !!$ !===where penalty is applied
427 !!$ FUNCTION penal_in_real_space(mesh,rr_gauss,angles,nb_angles,nb,ne,time) RESULT(vv)
428 !!$ IMPLICIT NONE
429 !!$ TYPE(mesh_type) :: mesh
430 !!$ REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: rr_gauss
431 !!$ REAL(KIND=8), DIMENSION(:), INTENT(IN) :: angles
432 !!$ INTEGER, INTENT(IN) :: nb_angles
433 !!$ INTEGER, INTENT(IN) :: nb, ne
434 !!$ REAL(KIND=8), INTENT(IN) :: time
435 !!$ REAL(KIND=8), DIMENSION(nb_angles,ne-nb+1) :: vv
436 !!$
437 !!$ vv = 1.d0
438 !!$ RETURN
439 !!$ END FUNCTION penal_in_real_space
440 
441  !===Extension of the velocity field in the solid.
442  !===Used when temperature or Maxwell equations are solved.
443  !===It extends the velocity field on the Navier-Stokes domain to a
444  !===velocity field on the temperature and the Maxwell domain.
445  !===It is also used if problem type=mxw and restart velocity
446  !===is set to true in data (type problem denoted mxx in the code).
447  FUNCTION extension_velocity(TYPE, H_mesh, mode, t, n_start) RESULT(vv)
448  IMPLICIT NONE
449  TYPE(mesh_type), INTENT(IN) :: H_mesh
450  INTEGER , INTENT(IN) :: TYPE, n_start
451  INTEGER, INTENT(IN) :: mode
452  REAL(KIND=8), INTENT(IN) :: t
453  REAL(KIND=8), DIMENSION(H_Mesh%np) :: vv
454  REAL(KIND=8) :: r
455  INTEGER :: n
456 
457  vv = 0.d0
458  RETURN
459 
460  !===Dummies variables to avoid warning
461  n=h_mesh%np; r=t; n=type; n=mode; n=n_start
462  !===Dummies variables to avoid warning
463  END FUNCTION extension_velocity
464 
465  !===============================================================================
466  ! Boundary conditions for Maxwell
467  !===============================================================================
468 !!$ !===Velocity used in the induction equation.
469 !!$ !===Used only if problem type is mxw and restart velocity is false
470 !!$ FUNCTION Vexact(m, H_mesh) RESULT(vv) !Set uniquement a l'induction
471 !!$ IMPLICIT NONE
472 !!$ TYPE(mesh_type), INTENT(IN) :: H_mesh
473 !!$ INTEGER, INTENT(IN) :: m
474 !!$ REAL(KIND=8), DIMENSION(H_mesh%np,6) :: vv
475 !!$
476 !!$ vv = 0.d0
477 !!$ RETURN
478 !!$ END FUNCTION Vexact
479 
480 !!$ FUNCTION H_B_quasi_static(char_h_b, rr, m) RESULT(vv)
481 !!$ IMPLICIT NONE
482 !!$ CHARACTER(LEN=1), INTENT(IN) :: char_h_b
483 !!$ REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: rr
484 !!$ INTEGER, INTENT(IN) :: m
485 !!$ REAL(KIND=8), DIMENSION(SIZE(rr,2),6) :: vv
486 !!$
487 !!$ IF (inputs%if_quasi_static_approx) THEN
488 !!$ vv = 0.d0
489 !!$ ELSE
490 !!$ CALL error_petsc('H_B_quasi_static should not be called')
491 !!$ END IF
492 !!$ RETURN
493 !!$ END FUNCTION H_B_quasi_static
494 
495 !!$ !===Magnetic field for boundary conditions in the Maxwell equations.
496 !!$ FUNCTION Hexact(H_mesh, TYPE, rr, m, mu_H_field, t) RESULT(vv)
497 !!$ IMPLICIT NONE
498 !!$ TYPE(mesh_type), INTENT(IN) :: H_mesh
499 !!$ INTEGER , INTENT(IN) :: TYPE
500 !!$ REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: rr
501 !!$ INTEGER , INTENT(IN) :: m
502 !!$ REAL(KIND=8), INTENT(IN) :: t
503 !!$ REAL(KIND=8), DIMENSION(:), INTENT(IN) :: mu_H_field
504 !!$ REAL(KIND=8), DIMENSION(SIZE(rr,2)) :: vv
505 !!$
506 !!$ vv = 0.d0
507 !!$ RETURN
508 !!$ END FUNCTION Hexact
509 
510 !!$ !===Scalar potential for boundary conditions in the Maxwell equations.
511 !!$ FUNCTION Phiexact(TYPE, rr, m, mu_phi,t) RESULT(vv)
512 !!$ IMPLICIT NONE
513 !!$ INTEGER , INTENT(IN) :: TYPE
514 !!$ REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: rr
515 !!$ INTEGER , INTENT(IN) :: m
516 !!$ REAL(KIND=8), INTENT(IN) :: mu_phi, t
517 !!$ REAL(KIND=8), DIMENSION(SIZE(rr,2)) :: vv
518 !!$
519 !!$ vv=0.d0
520 !!$ RETURN
521 !!$ END FUNCTION Phiexact
522 
523 !!$ !===Current in Ohm's law. Curl(H) = sigma(E + uxB) + current
524 !!$ FUNCTION Jexact_gauss(TYPE, rr, m, mu_phi, sigma, mu_H, t, mesh_id, opt_B_ext) RESULT(vv)
525 !!$ IMPLICIT NONE
526 !!$ INTEGER , INTENT(IN) :: TYPE
527 !!$ REAL(KIND=8), DIMENSION(:), INTENT(IN) :: rr
528 !!$ INTEGER , INTENT(IN) :: m
529 !!$ REAL(KIND=8), INTENT(IN) :: mu_phi, sigma, mu_H, t
530 !!$ INTEGER , INTENT(IN) :: mesh_id
531 !!$ REAL(KIND=8), DIMENSION(6), OPTIONAL,INTENT(IN) :: opt_B_ext
532 !!$ REAL(KIND=8) :: vv
533 !!$
534 !!$ vv = 0.d0
535 !!$ RETURN
536 !!$ END FUNCTION Jexact_gauss
537 
538 !!$ !===Electric field for Neumann BC (cf. doc)
539 !!$ FUNCTION Eexact_gauss(TYPE, rr, m, mu_phi, sigma, mu_H, t) RESULT(vv)
540 !!$ IMPLICIT NONE
541 !!$ INTEGER, INTENT(IN) :: TYPE
542 !!$ REAL(KIND=8), DIMENSION(:), INTENT(IN) :: rr
543 !!$ INTEGER, INTENT(IN) :: m
544 !!$ REAL(KIND=8), INTENT(IN) :: mu_phi, sigma, mu_H, t
545 !!$ REAL(KIND=8) :: vv
546 !!$
547 !!$ vv = 0.d0
548 !!$ RETURN
549 !!$ END FUNCTION Eexact_gauss
550 
551 !!$ !===Initialization of magnetic field and scalar potential (if present)
552 !!$ SUBROUTINE init_maxwell(H_mesh, phi_mesh, time, dt, mu_H_field, mu_phi, &
553 !!$ list_mode, Hn1, Hn, phin1, phin)
554 !!$ IMPLICIT NONE
555 !!$ TYPE(mesh_type) :: H_mesh, phi_mesh
556 !!$ REAL(KIND=8), INTENT(OUT):: time
557 !!$ REAL(KIND=8), INTENT(IN) :: dt
558 !!$ REAL(KIND=8), DIMENSION(:), INTENT(IN) :: mu_H_field
559 !!$ REAL(KIND=8), INTENT(IN) :: mu_phi
560 !!$ INTEGER, DIMENSION(:), INTENT(IN) :: list_mode
561 !!$ REAL(KIND=8), DIMENSION(:,:,:), INTENT(OUT):: Hn, Hn1
562 !!$ REAL(KIND=8), DIMENSION(:,:,:), INTENT(OUT):: phin, phin1
563 !!$
564 !!$ time = 0.d0
565 !!$ Hn1=0.d0
566 !!$ HN=0.d0
567 !!$ phin = 0.d0
568 !!$ phin1 =0.d0
569 !!$ RETURN
570 !!$ END SUBROUTINE init_maxwell
571 
572 !!$ !===Analytical permeability (if needed)
573 !!$ !===This function is not needed unless the flag
574 !!$ !=== ===Use FEM Interpolation for magnetic permeability (true/false)
575 !!$ !===is activated and set to .FALSE. in the data data file. Default is .TRUE.
576 !!$ FUNCTION mu_bar_in_fourier_space(H_mesh,nb,ne,pts,pts_ids) RESULT(vv)
577 !!$ IMPLICIT NONE
578 !!$ TYPE(mesh_type), INTENT(IN) :: H_mesh
579 !!$ REAL(KIND=8), DIMENSION(ne-nb+1) :: vv
580 !!$ INTEGER, INTENT(IN) :: nb, ne
581 !!$ REAL(KIND=8),DIMENSION(2,ne-nb+1),OPTIONAL :: pts
582 !!$ INTEGER, DIMENSION(ne-nb+1), OPTIONAL :: pts_ids
583 !!$
584 !!$ vv = 1.d0
585 !!$ RETURN
586 !!$ END FUNCTION mu_bar_in_fourier_space
587 
588 !!$ !===Analytical mu_in_fourier_space (if needed)
589 !!$ !===This function is not needed unless the flag
590 !!$ !=== ===Use FEM Interpolation for magnetic permeability (true/false)
591 !!$ !===is activated and set to .FALSE. in the data data file. Default is .TRUE.
592 !!$ FUNCTION grad_mu_bar_in_fourier_space(pt,pt_id) RESULT(vv)
593 !!$ IMPLICIT NONE
594 !!$ REAL(KIND=8),DIMENSION(2), INTENT(in):: pt
595 !!$ INTEGER,DIMENSION(1), INTENT(in) :: pt_id
596 !!$ REAL(KIND=8),DIMENSION(2) :: vv
597 !!$
598 !!$ vv=0.d0
599 !!$ RETURN
600 !!$ END FUNCTION grad_mu_bar_in_fourier_space
601 
602 !!$ !===Analytical permeability, mu in real space (if needed)
603 !!$ FUNCTION mu_in_real_space(H_mesh,angles,nb_angles,nb,ne,time) RESULT(vv)
604 !!$ IMPLICIT NONE
605 !!$ TYPE(mesh_type), INTENT(IN) :: H_mesh
606 !!$ REAL(KIND=8), DIMENSION(:), INTENT(IN) :: angles
607 !!$ INTEGER, INTENT(IN) :: nb_angles
608 !!$ INTEGER, INTENT(IN) :: nb, ne
609 !!$ REAL(KIND=8), INTENT(IN) :: time
610 !!$ REAL(KIND=8), DIMENSION(nb_angles,ne-nb+1) :: vv
611 !!$
612 !!$ vv = 1.d0
613 !!$ RETURN
614 !!$ END FUNCTION mu_in_real_space
615 
616 END MODULE boundary_test_32
real(kind=8) function, dimension(size(rr, 2)), public source_in_temperature(TYPE, rr, m, t)
real(kind=8) function, dimension(size(rr, 2)), public temperature_exact(TYPE, rr, m, t)
subroutine error_petsc(string)
Definition: my_util.f90:16
type(my_data), public inputs
subroutine, public init_velocity_pressure(mesh_f, mesh_c, time, dt, list_mode, un_m1, un, pn_m1, pn, phin_m1, phin)
real(kind=8) function, dimension(size(rr, 2)), public pp_exact(TYPE, rr, m, t)
real(kind=8) function, dimension(size(rr, 2)), public source_in_ns_momentum(TYPE, rr, mode, i, time, Re, ty, opt_density, opt_tempn)
subroutine, public init_temperature(mesh, time, dt, list_mode, tempn_m1, tempn)
real(kind=8) function, dimension(size(rr, 2)), public vv_exact(TYPE, rr, m, t)
real(kind=8) function, dimension(h_mesh%np), public extension_velocity(TYPE, H_mesh, mode, t, n_start)