SFEMaNS  version 5.3
Reference documentation for SFEMaNS
condlim_test_31.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 31')
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  RETURN
288 
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  RETURN
347  END FUNCTION vv_exact
348 
349 !!$ !===Solid velocity imposed when using penalty technique
350 !!$ !===Defined in Fourier space on mode 0 only.
351 !!$ FUNCTION imposed_velocity_by_penalty(rr,t) RESULT(vv)
352 !!$ IMPLICIT NONE
353 !!$ REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: rr
354 !!$ REAL(KIND=8), INTENT(IN) :: t
355 !!$ REAL(KIND=8), DIMENSION(SIZE(rr,2),6) :: vv
356 !!$
357 !!$ vv=0.d0
358 !!$ RETURN
359 !!$ END FUNCTION imposed_velocity_by_penalty
360 
361  !===Pressure for boundary conditions in Navier-Stokes.
362  !===Can be used also to initialize pressure in the subroutine init_velocity_pressure.
363  !===Use this routine for outflow BCs only.
364  !===CAUTION: Do not enfore BCs on pressure where normal component
365  ! of velocity is prescribed.
366  FUNCTION pp_exact(TYPE,rr,m,t) RESULT (vv)
367  IMPLICIT NONE
368  INTEGER , INTENT(IN) :: TYPE
369  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: rr
370  INTEGER , INTENT(IN) :: m
371  REAL(KIND=8), INTENT(IN) :: t
372  REAL(KIND=8), DIMENSION(SIZE(rr,2)) :: vv, r, z
373  REAL(KIND=8) :: pi = 3.1415926535897932d0
374 
375  r = rr(1,:)
376  z = rr(2,:)
377 
378  IF ((type==1).AND.(m==1)) THEN
379  vv = r**3*sin(2*pi*z)*cos(t)
380  ELSE
381  vv = 0.d0
382  END IF
383  RETURN
384  END FUNCTION pp_exact
385 
386  !===Temperature for boundary conditions in temperature equation.
387  FUNCTION temperature_exact(TYPE,rr,m,t) RESULT (vv)
388  IMPLICIT NONE
389  INTEGER , INTENT(IN) :: TYPE
390  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: rr
391  INTEGER , INTENT(IN) :: m
392  REAL(KIND=8), INTENT(IN) :: t
393  REAL(KIND=8), DIMENSION(SIZE(rr,2)) :: vv, r, z
394  REAL(KIND=8) :: r0=0.5d0
395  REAL(KIND=8) :: pi = 3.1415926535897932d0
396 
397  r = rr(1,:)
398  z = rr(2,:)
399 
400  IF ((type==1).AND.((m==0).OR.(m==1))) THEN
401  vv = r**2*(r-r0)**2*sin(2*pi*z)*cos(t)
402  ELSE
403  vv = 0.d0
404  END IF
405 
406  RETURN
407  END FUNCTION temperature_exact
408 
409 !!$ !===Can be used to initialize level set in the subroutine init_level_set.
410 !!$ FUNCTION level_set_exact(interface_nb,TYPE,rr,m,t) RESULT (vv)
411 !!$ IMPLICIT NONE
412 !!$ INTEGER , INTENT(IN) :: TYPE
413 !!$ REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: rr
414 !!$ INTEGER , INTENT(IN) :: m, interface_nb
415 !!$ REAL(KIND=8), INTENT(IN) :: t
416 !!$ REAL(KIND=8), DIMENSION(SIZE(rr,2)) :: vv
417 !!$
418 !!$ vv = 0.d0
419 !!$ RETURN
420 !!$ END FUNCTION level_set_exact
421 
422 !!$ !===Penalty coefficient (if needed)
423 !!$ !===This coefficient is equal to zero in subdomain
424 !!$ !===where penalty is applied
425 !!$ FUNCTION penal_in_real_space(mesh,rr_gauss,angles,nb_angles,nb,ne,time) RESULT(vv)
426 !!$ IMPLICIT NONE
427 !!$ TYPE(mesh_type) :: mesh
428 !!$ REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: rr_gauss
429 !!$ REAL(KIND=8), DIMENSION(:), INTENT(IN) :: angles
430 !!$ INTEGER, INTENT(IN) :: nb_angles
431 !!$ INTEGER, INTENT(IN) :: nb, ne
432 !!$ REAL(KIND=8), INTENT(IN) :: time
433 !!$ REAL(KIND=8), DIMENSION(nb_angles,ne-nb+1) :: vv
434 !!$
435 !!$ vv = 1.d0
436 !!$ RETURN
437 !!$ END FUNCTION penal_in_real_space
438 
439  !===Extension of the velocity field in the solid.
440  !===Used when temperature or Maxwell equations are solved.
441  !===It extends the velocity field on the Navier-Stokes domain to a
442  !===velocity field on the temperature and the Maxwell domain.
443  !===It is also used if problem type=mxw and restart velocity
444  !===is set to true in data (type problem denoted mxx in the code).
445  FUNCTION extension_velocity(TYPE, H_mesh, mode, t, n_start) RESULT(vv)
446  IMPLICIT NONE
447  TYPE(mesh_type), INTENT(IN) :: H_mesh
448  INTEGER , INTENT(IN) :: TYPE, n_start
449  INTEGER, INTENT(IN) :: mode
450  REAL(KIND=8), INTENT(IN) :: t
451  REAL(KIND=8), DIMENSION(H_Mesh%np) :: vv
452  REAL(KIND=8) :: r
453  INTEGER :: n
454 
455  vv = 0.d0
456  RETURN
457 
458  !===Dummies variables to avoid warning
459  n=h_mesh%np; r=t; n=type; n=mode; n=n_start
460  !===Dummies variables to avoid warning
461  END FUNCTION extension_velocity
462 
463  !===============================================================================
464  ! Boundary conditions for Maxwell
465  !===============================================================================
466 !!$ !===Velocity used in the induction equation.
467 !!$ !===Used only if problem type is mxw and restart velocity is false
468 !!$ FUNCTION Vexact(m, H_mesh) RESULT(vv) !Set uniquement a l'induction
469 !!$ IMPLICIT NONE
470 !!$ TYPE(mesh_type), INTENT(IN) :: H_mesh
471 !!$ INTEGER, INTENT(IN) :: m
472 !!$ REAL(KIND=8), DIMENSION(H_mesh%np,6) :: vv
473 !!$
474 !!$ vv = 0.d0
475 !!$ RETURN
476 !!$ END FUNCTION Vexact
477 
478 !!$ FUNCTION H_B_quasi_static(char_h_b, rr, m) RESULT(vv)
479 !!$ IMPLICIT NONE
480 !!$ CHARACTER(LEN=1), INTENT(IN) :: char_h_b
481 !!$ REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: rr
482 !!$ INTEGER, INTENT(IN) :: m
483 !!$ REAL(KIND=8), DIMENSION(SIZE(rr,2),6) :: vv
484 !!$
485 !!$ IF (inputs%if_quasi_static_approx) THEN
486 !!$ vv = 0.d0
487 !!$ ELSE
488 !!$ CALL error_petsc('H_B_quasi_static should not be called')
489 !!$ END IF
490 !!$ RETURN
491 !!$ END FUNCTION H_B_quasi_static
492 
493 !!$ !===Magnetic field for boundary conditions in the Maxwell equations.
494 !!$ FUNCTION Hexact(H_mesh, TYPE, rr, m, mu_H_field, t) RESULT(vv)
495 !!$ IMPLICIT NONE
496 !!$ TYPE(mesh_type), INTENT(IN) :: H_mesh
497 !!$ INTEGER , INTENT(IN) :: TYPE
498 !!$ REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: rr
499 !!$ INTEGER , INTENT(IN) :: m
500 !!$ REAL(KIND=8), INTENT(IN) :: t
501 !!$ REAL(KIND=8), DIMENSION(:), INTENT(IN) :: mu_H_field
502 !!$ REAL(KIND=8), DIMENSION(SIZE(rr,2)) :: vv
503 !!$
504 !!$ vv = 0.d0
505 !!$ RETURN
506 !!$ END FUNCTION Hexact
507 
508 !!$ !===Scalar potential for boundary conditions in the Maxwell equations.
509 !!$ FUNCTION Phiexact(TYPE, rr, m, mu_phi,t) RESULT(vv)
510 !!$ IMPLICIT NONE
511 !!$ INTEGER , INTENT(IN) :: TYPE
512 !!$ REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: rr
513 !!$ INTEGER , INTENT(IN) :: m
514 !!$ REAL(KIND=8), INTENT(IN) :: mu_phi, t
515 !!$ REAL(KIND=8), DIMENSION(SIZE(rr,2)) :: vv
516 !!$
517 !!$ vv=0.d0
518 !!$ RETURN
519 !!$ END FUNCTION Phiexact
520 
521 !!$ !===Current in Ohm's law. Curl(H) = sigma(E + uxB) + current
522 !!$ FUNCTION Jexact_gauss(TYPE, rr, m, mu_phi, sigma, mu_H, t, mesh_id, opt_B_ext) RESULT(vv)
523 !!$ IMPLICIT NONE
524 !!$ INTEGER , INTENT(IN) :: TYPE
525 !!$ REAL(KIND=8), DIMENSION(:), INTENT(IN) :: rr
526 !!$ INTEGER , INTENT(IN) :: m
527 !!$ REAL(KIND=8), INTENT(IN) :: mu_phi, sigma, mu_H, t
528 !!$ INTEGER , INTENT(IN) :: mesh_id
529 !!$ REAL(KIND=8), DIMENSION(6), OPTIONAL,INTENT(IN) :: opt_B_ext
530 !!$ REAL(KIND=8) :: vv
531 !!$
532 !!$ vv = 0.d0
533 !!$ RETURN
534 !!$ END FUNCTION Jexact_gauss
535 
536 !!$ !===Electric field for Neumann BC (cf. doc)
537 !!$ FUNCTION Eexact_gauss(TYPE, rr, m, mu_phi, sigma, mu_H, t) RESULT(vv)
538 !!$ IMPLICIT NONE
539 !!$ INTEGER, INTENT(IN) :: TYPE
540 !!$ REAL(KIND=8), DIMENSION(:), INTENT(IN) :: rr
541 !!$ INTEGER, INTENT(IN) :: m
542 !!$ REAL(KIND=8), INTENT(IN) :: mu_phi, sigma, mu_H, t
543 !!$ REAL(KIND=8) :: vv
544 !!$
545 !!$ vv = 0.d0
546 !!$ RETURN
547 !!$ END FUNCTION Eexact_gauss
548 
549 !!$ !===Initialization of magnetic field and scalar potential (if present)
550 !!$ SUBROUTINE init_maxwell(H_mesh, phi_mesh, time, dt, mu_H_field, mu_phi, &
551 !!$ list_mode, Hn1, Hn, phin1, phin)
552 !!$ IMPLICIT NONE
553 !!$ TYPE(mesh_type) :: H_mesh, phi_mesh
554 !!$ REAL(KIND=8), INTENT(OUT):: time
555 !!$ REAL(KIND=8), INTENT(IN) :: dt
556 !!$ REAL(KIND=8), DIMENSION(:), INTENT(IN) :: mu_H_field
557 !!$ REAL(KIND=8), INTENT(IN) :: mu_phi
558 !!$ INTEGER, DIMENSION(:), INTENT(IN) :: list_mode
559 !!$ REAL(KIND=8), DIMENSION(:,:,:), INTENT(OUT):: Hn, Hn1
560 !!$ REAL(KIND=8), DIMENSION(:,:,:), INTENT(OUT):: phin, phin1
561 !!$
562 !!$ time = 0.d0
563 !!$ Hn1=0.d0
564 !!$ HN=0.d0
565 !!$ phin = 0.d0
566 !!$ phin1 =0.d0
567 !!$ RETURN
568 !!$ END SUBROUTINE init_maxwell
569 
570 !!$ !===Analytical permeability (if needed)
571 !!$ !===This function is not needed unless the flag
572 !!$ !=== ===Use FEM Interpolation for magnetic permeability (true/false)
573 !!$ !===is activated and set to .FALSE. in the data data file. Default is .TRUE.
574 !!$ FUNCTION mu_bar_in_fourier_space(H_mesh,nb,ne,pts,pts_ids) RESULT(vv)
575 !!$ IMPLICIT NONE
576 !!$ TYPE(mesh_type), INTENT(IN) :: H_mesh
577 !!$ REAL(KIND=8), DIMENSION(ne-nb+1) :: vv
578 !!$ INTEGER, INTENT(IN) :: nb, ne
579 !!$ REAL(KIND=8),DIMENSION(2,ne-nb+1),OPTIONAL :: pts
580 !!$ INTEGER, DIMENSION(ne-nb+1), OPTIONAL :: pts_ids
581 !!$
582 !!$ vv = 1.d0
583 !!$ RETURN
584 !!$ END FUNCTION mu_bar_in_fourier_space
585 
586 !!$ !===Analytical mu_in_fourier_space (if needed)
587 !!$ !===This function is not needed unless the flag
588 !!$ !=== ===Use FEM Interpolation for magnetic permeability (true/false)
589 !!$ !===is activated and set to .FALSE. in the data data file. Default is .TRUE.
590 !!$ FUNCTION grad_mu_bar_in_fourier_space(pt,pt_id) RESULT(vv)
591 !!$ IMPLICIT NONE
592 !!$ REAL(KIND=8),DIMENSION(2), INTENT(in):: pt
593 !!$ INTEGER,DIMENSION(1), INTENT(in) :: pt_id
594 !!$ REAL(KIND=8),DIMENSION(2) :: vv
595 !!$
596 !!$ vv=0.d0
597 !!$ RETURN
598 !!$ END FUNCTION grad_mu_bar_in_fourier_space
599 
600 !!$ !===Analytical permeability, mu in real space (if needed)
601 !!$ FUNCTION mu_in_real_space(H_mesh,angles,nb_angles,nb,ne,time) RESULT(vv)
602 !!$ IMPLICIT NONE
603 !!$ TYPE(mesh_type), INTENT(IN) :: H_mesh
604 !!$ REAL(KIND=8), DIMENSION(:), INTENT(IN) :: angles
605 !!$ INTEGER, INTENT(IN) :: nb_angles
606 !!$ INTEGER, INTENT(IN) :: nb, ne
607 !!$ REAL(KIND=8), INTENT(IN) :: time
608 !!$ REAL(KIND=8), DIMENSION(nb_angles,ne-nb+1) :: vv
609 !!$
610 !!$ vv = 1.d0
611 !!$ RETURN
612 !!$ END FUNCTION mu_in_real_space
613 
614 END MODULE boundary_test_31
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)
real(kind=8) function, dimension(size(rr, 2)), public source_in_temperature(TYPE, rr, m, t)
subroutine, public init_velocity_pressure(mesh_f, mesh_c, time, dt, list_mode, un_m1, un, pn_m1, pn, phin_m1, phin)
subroutine error_petsc(string)
Definition: my_util.f90:16
type(my_data), public inputs
real(kind=8) function, dimension(size(rr, 2)), public pp_exact(TYPE, rr, m, t)
real(kind=8) function, dimension(size(rr, 2)), public temperature_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)