SFEMaNS  version 5.3
Reference documentation for SFEMaNS
condlim_test_2.f90
Go to the documentation of this file.
2  USE my_util
3  USE def_type_mesh
4  USE input_data
5 !!$Some subroutines have been commented to avoid warning messages when compiling executable.
6 !!$It can not be done in the module boundary_generic that expects all subroutines to be present.
7 !!$END ATTENTION
8  PUBLIC :: init_velocity_pressure
9 !!$ PUBLIC :: init_temperature
10 !!$ PUBLIC :: init_level_set
11  PUBLIC :: source_in_ns_momentum
12 !!$ PUBLIC :: source_in_temperature
13 !!$ PUBLIC :: source_in_level_set
14  PUBLIC :: vv_exact
15 !!$ PUBLIC :: imposed_velocity_by_penalty
16  PUBLIC :: pp_exact
17 !!$ PUBLIC :: temperature_exact
18 !!$ PUBLIC :: level_set_exact
19 !!$ PUBLIC :: penal_in_real_space
20 !!$ PUBLIC :: extension_velocity
21 !!$ PUBLIC :: Vexact
22 !!$ PUBLIC :: H_B_quasi_static
23 !!$ PUBLIC :: Hexact
24 !!$ PUBLIC :: Phiexact
25 !!$ PUBLIC :: Jexact_gauss
26 !!$ PUBLIC :: Eexact_gauss
27 !!$ PUBLIC :: init_maxwell
28 !!$ PUBLIC :: mu_bar_in_fourier_space
29 !!$ PUBLIC :: grad_mu_bar_in_fourier_space
30 !!$ PUBLIC :: mu_in_real_space
31 !!$ATTENTION
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 !!$
102 !!$ time = 0.d0
103 !!$ DO i= 1, SIZE(list_mode)
104 !!$ mode = list_mode(i)
105 !!$ DO j = 1, 2
106 !!$ !===level_set
107 !!$ DO n = 1, inputs%nb_fluid -1
108 !!$ level_set_m1(n,:,j,i) = level_set_exact(n,j,vv_mesh%rr,mode,time-dt)
109 !!$ level_set (n,:,j,i) = level_set_exact(n,j,vv_mesh%rr,mode,time)
110 !!$ END DO
111 !!$ END DO
112 !!$ END DO
113 !!$
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
128  REAL(KIND=8), DIMENSION(SIZE(rr,2)) :: LAP, TEMP,GPRESS
129  REAL(KIND=8), DIMENSION(SIZE(rr,2)) :: r, z
130  REAL(KIND=8), DIMENSION(SIZE(rr,2),6) :: Rot, V, vt
131  INTEGER :: m
132  REAL(KIND=8) :: t, k
133  CHARACTER(LEN=2) :: np
134 
135  IF (PRESENT(opt_density)) CALL error_petsc('density should not be present for test 2')
136  IF (PRESENT(opt_tempn)) CALL error_petsc('temperature should not be present for test 2')
137 
138  m = mode
139  t = time
140  k = 2*acos(-1.d0)
141  r = rr(1,:)
142  z = rr(2,:)
143 
144  IF (m==0 .OR. m==2) THEN
145  !Nonlinear contribution
146  rot(:,1) = r*(4.d0*cos(k*z)+1)
147  rot(:,2) = 0.d0
148  rot(:,3) = 0.d0
149  rot(:,4) = r*(k**2*r**2*cos(k*z)-8.d0*cos(k*z)-2.d0)
150  rot(:,5) = r*(-8.d0-k*r*sin(k*z))
151  rot(:,6) = 0.d0
152  rot = rot*cos(t)
153 
154  vt(:,1) = 0.d0
155  vt(:,4) = 0.d0
156  vt(:,5) = 0.d0
157  vt(:,2) = -r**2*(1.d0-k*r*sin(k*z))
158  vt(:,3) = -3.d0*r**2
159  vt(:,6) = r**2*(4.d0*cos(k*z)+1.d0)
160  v = vt*cos(t)
161 
162  IF (m==0) THEN ! vv = -VxRot(v) pour mode 0
163  IF (mod(TYPE,2)==0) then
164  vv = 0.d0
165  ELSEIF (TYPE == 1) then
166  vv = -0.5d0*(v(:,3)*rot(:,5)-v(:,5)*rot(:,3) &
167  + v(:,4)*rot(:,6)-v(:,6)*rot(:,4))
168  ELSEIF (TYPE == 3) then
169  vv = -0.5d0*(v(:,5)*rot(:,1)-v(:,1)*rot(:,5) &
170  + v(:,6)*rot(:,2)-v(:,2)*rot(:,6))
171  ELSEIF (TYPE == 5) then
172  vv = -0.5d0*(v(:,1)*rot(:,3)-v(:,3)*rot(:,1) &
173  + v(:,2)*rot(:,4)-v(:,4)*rot(:,2))
174  ENDIF
175  RETURN
176  END IF
177 
178  IF (m==2) THEN ! vv = -VxRot(v) pour mode 2
179  IF (TYPE == 1) then
180  vv = -0.5d0*(v(:,3)*rot(:,5)-v(:,5)*rot(:,3) &
181  - (v(:,4)*rot(:,6)-v(:,6)*rot(:,4)))
182  ELSE IF (TYPE == 2) then
183  vv = -0.5d0*(v(:,3)*rot(:,6)-v(:,5)*rot(:,4) &
184  + (v(:,4)*rot(:,5)-v(:,6)*rot(:,3)))
185  ELSEIF (TYPE == 3) then
186  vv = -0.5d0*(v(:,5)*rot(:,1)-v(:,1)*rot(:,5) &
187  - (v(:,6)*rot(:,2)-v(:,2)*rot(:,6)))
188  ELSEIF (TYPE == 4) then
189  vv = -0.5d0*(v(:,5)*rot(:,2)-v(:,1)*rot(:,6) &
190  + (v(:,6)*rot(:,1)-v(:,2)*rot(:,5)))
191  ELSEIF (TYPE == 5) then
192  vv = -0.5d0*(v(:,1)*rot(:,3)-v(:,3)*rot(:,1) &
193  - (v(:,2)*rot(:,4)-v(:,4)*rot(:,2)))
194  ELSEIF (TYPE == 6) then
195  vv = -0.5d0*(v(:,1)*rot(:,4)-v(:,3)*rot(:,2) &
196  + (v(:,2)*rot(:,3)-v(:,4)*rot(:,1)))
197  END IF
198 
199  RETURN
200  END IF
201 
202  END IF
203  IF (m/=1) THEN
204  vv = 0.d0
205  RETURN
206  END IF
207 
208  vt(:,1) = 0.d0
209  vt(:,4) = 0.d0
210  vt(:,5) = 0.d0
211  vt(:,2) = -r**2*(1.d0-k*r*sin(k*z))
212  vt(:,3) = -3.d0*r**2
213  vt(:,6) = r**2*(4.d0*cos(k*z)+1.d0)
214 
215  IF (TYPE == 1) then
216  lap(:) = 0.d0
217  temp(:) = vt(:,1)
218  gpress(:) = 2*r*cos(k*z)
219  ELSEIF (TYPE == 2) then
220  lap(:) = (-8+7*k*r*sin(k*z)-k**3*r**3*sin(k*z))
221  temp(:) = vt(:,2)
222  gpress(:) = 0.d0
223  ELSEIF (TYPE == 3) then
224  lap(:) = (-8+2*k*r*sin(k*z))
225  temp(:) = vt(:,3)
226  gpress(:) = 0.d0
227  ELSEIF (TYPE == 4) then
228  lap(:) = 0.d0
229  temp(:) = vt(:,4)
230  gpress(:) = -r*cos(k*z)
231  ELSEIF (TYPE == 5) then
232  lap(:) = 0.d0
233  temp(:) = vt(:,5)
234  gpress(:) = -r**2*k*sin(k*z)
235  ELSEIF (TYPE == 6) then
236  lap(:) = (3+12*cos(k*z)-4*k**2*r**2*cos(k*z))
237  temp(:) = vt(:,6)
238  gpress(:) = 0.d0
239  ENDIF
240 
241  vv(:) = -sin(t)*temp + cos(t)*(-lap/re) + gpress*cos(t)
242  RETURN
243 
244  !===Dummies variables to avoid warning
245  m=i; np=ty
246  !===Dummies variables to avoid warning
247  END FUNCTION source_in_ns_momentum
248 
249 !!$ !===Extra source in temperature equation. Always called.
250 !!$ FUNCTION source_in_temperature(TYPE, rr, m, t)RESULT(vv)
251 !!$ IMPLICIT NONE
252 !!$ INTEGER , INTENT(IN) :: TYPE
253 !!$ REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: rr
254 !!$ INTEGER , INTENT(IN) :: m
255 !!$ REAL(KIND=8), INTENT(IN) :: t
256 !!$ REAL(KIND=8), DIMENSION(SIZE(rr,2)) :: vv
257 !!$
258 !!$ vv = 0.d0
259 !!$ CALL error_petsc('source_in_temperature: should not be called for this test')
260 !!$ RETURN
261 !!$ END FUNCTION source_in_temperature
262 
263 !!$ !===Extra source in level set equation. Always called.
264 !!$ FUNCTION source_in_level_set(interface_nb,TYPE, rr, m, t)RESULT(vv)
265 !!$ IMPLICIT NONE
266 !!$ INTEGER , INTENT(IN) :: TYPE
267 !!$ REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: rr
268 !!$ INTEGER , INTENT(IN) :: m, interface_nb
269 !!$ REAL(KIND=8), INTENT(IN) :: t
270 !!$ REAL(KIND=8), DIMENSION(SIZE(rr,2)) :: vv
271 !!$
272 !!$ vv=0.d0
273 !!$ CALL error_petsc('source_in_level_set: should not be called for this test')
274 !!$
275 !!$ END FUNCTION source_in_level_set
276 
277  !===Velocity for boundary conditions in Navier-Stokes.
278  !===Can be used also to initialize velocity in: init_velocity_pressure_temperature
279  FUNCTION vv_exact(TYPE,rr,m,t) RESULT(vv)
280  IMPLICIT NONE
281  INTEGER , INTENT(IN) :: type
282  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: rr
283  INTEGER, INTENT(IN) :: m
284  REAL(KIND=8), INTENT(IN) :: t
285  REAL(KIND=8), DIMENSION(SIZE(rr,2)) :: vv
286  REAL(KIND=8) :: k
287  REAL(KIND=8), DIMENSION(SIZE(rr,2)) :: r, z
288 
289  k = 2*acos(-1.d0)
290  r = rr(1,:)
291  z = rr(2,:)
292 
293  IF (m/=1) THEN
294  vv = 0.d0
295  RETURN
296  END IF
297 
298  IF (TYPE == 1) then
299  vv(:) = 0.d0
300  ELSEIF (TYPE == 2 .AND. m /= 0) THEN
301  vv(:) = -r**2*(1-k*r*sin(k*z))
302  ELSEIF (TYPE == 3) then
303  vv(:) = -3*r**2
304  ELSEIF (TYPE == 4 .AND. m /= 0) THEN
305  vv(:) = 0.d0
306  ELSEIF (TYPE == 5) then
307  vv(:) = 0.d0
308  ELSEIF (TYPE == 6 .AND. m /= 0) THEN
309  vv(:) = r**2*(4*cos(k*z)+1)
310  ENDIF
311 
312  vv(:) = vv(:) * cos(t)
313  RETURN
314  END FUNCTION vv_exact
315 
316 !!$ !===Solid velocity imposed when using penalty technique
317 !!$ !===Defined in Fourier space on mode 0 only.
318 !!$ FUNCTION imposed_velocity_by_penalty(rr,t) RESULT(vv)
319 !!$ IMPLICIT NONE
320 !!$ REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: rr
321 !!$ REAL(KIND=8), INTENT(IN) :: t
322 !!$ REAL(KIND=8), DIMENSION(SIZE(rr,2),6) :: vv
323 !!$
324 !!$ vv=0.d0
325 !!$ RETURN
326 !!$ END FUNCTION imposed_velocity_by_penalty
327 
328  !===Pressure for boundary conditions in Navier-Stokes.
329  !===Can be used also to initialize pressure in the subroutine init_velocity_pressure.
330  !===Use this routine for outflow BCs only.
331  !===CAUTION: Do not enfore BCs on pressure where normal component
332  ! of velocity is prescribed.
333  FUNCTION pp_exact(TYPE,rr,m,t) RESULT (vv)
334  IMPLICIT NONE
335  INTEGER , INTENT(IN) :: TYPE
336  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: rr
337  INTEGER , INTENT(IN) :: m
338  REAL(KIND=8), INTENT(IN) :: t
339  REAL(KIND=8), DIMENSION(SIZE(rr,2)) :: vv
340  REAL(KIND=8) :: k
341 
342  k = 2*acos(-1.d0)
343  IF (m/=1) THEN
344  vv = 0.d0
345  RETURN
346  END IF
347 
348  IF (type==1) THEN
349  vv(:)= rr(1,:)**2*cos(k*rr(2,:))*cos(t)
350  ELSE
351  vv(:) = 0.d0
352  END IF
353  RETURN
354  END FUNCTION pp_exact
355 
356 !!$ !===Temperature for boundary conditions in temperature equation.
357 !!$ FUNCTION temperature_exact(TYPE,rr,m,t) RESULT (vv)
358 !!$ IMPLICIT NONE
359 !!$ INTEGER , INTENT(IN) :: TYPE
360 !!$ REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: rr
361 !!$ INTEGER , INTENT(IN) :: m
362 !!$ REAL(KIND=8), INTENT(IN) :: t
363 !!$ REAL(KIND=8), DIMENSION(SIZE(rr,2)) :: vv
364 !!$
365 !!$ vv = 0.d0
366 !!$ CALL error_petsc('temperature_exact: should not be called for this test')
367 !!$ RETURN
368 !!$ END FUNCTION temperature_exact
369 
370 !!$ !===Can be used to initialize level set in the subroutine init_level_set.
371 !!$ FUNCTION level_set_exact(interface_nb,TYPE,rr,m,t) RESULT (vv)
372 !!$ IMPLICIT NONE
373 !!$ INTEGER , INTENT(IN) :: TYPE
374 !!$ REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: rr
375 !!$ INTEGER , INTENT(IN) :: m, interface_nb
376 !!$ REAL(KIND=8), INTENT(IN) :: t
377 !!$ REAL(KIND=8), DIMENSION(SIZE(rr,2)) :: vv
378 !!$
379 !!$ vv = 0.d0
380 !!$ CALL error_petsc('level_set_exact: should not be called for this test')
381 !!$ RETURN
382 !!$
383 !!$ END FUNCTION level_set_exact
384 
385 !!$ !===Penalty coefficient (if needed)
386 !!$ !===This coefficient is equal to zero in subdomain
387 !!$ !===where penalty is applied (penalty is zero in solid)
388 !!$ FUNCTION penal_in_real_space(mesh,rr_gauss,angles,nb_angles,nb,ne,time) RESULT(vv)
389 !!$ IMPLICIT NONE
390 !!$ TYPE(mesh_type) :: mesh
391 !!$ REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: rr_gauss
392 !!$ REAL(KIND=8), DIMENSION(:), INTENT(IN) :: angles
393 !!$ INTEGER, INTENT(IN) :: nb_angles
394 !!$ INTEGER, INTENT(IN) :: nb, ne
395 !!$ REAL(KIND=8), INTENT(IN) :: time
396 !!$ REAL(KIND=8), DIMENSION(nb_angles,ne-nb+1) :: vv
397 !!$
398 !!$ vv = 1.d0
399 !!$ CALL error_petsc('penal_in_real_space: should not be called for this test')
400 !!$ RETURN
401 !!$ END FUNCTION penal_in_real_space
402 
403 !!$ !===Extension of the velocity field in the solid.
404 !!$ !===Used when temperature or Maxwell equations are solved.
405 !!$ !===It extends the velocity field on the Navier-Stokes domain to a
406 !!$ !===velocity field on the temperature and the Maxwell domain.
407 !!$ !===It is also used if problem type=mxw and restart velocity
408 !!$ !===is set to true in data (type problem denoted mxx in the code).
409 !!$ FUNCTION extension_velocity(TYPE, H_mesh, mode, t, n_start) RESULT(vv)
410 !!$ IMPLICIT NONE
411 !!$ TYPE(mesh_type), INTENT(IN) :: H_mesh
412 !!$ INTEGER , INTENT(IN) :: TYPE, n_start
413 !!$ INTEGER, INTENT(IN) :: mode
414 !!$ REAL(KIND=8), INTENT(IN) :: t
415 !!$ REAL(KIND=8), DIMENSION(H_Mesh%np) :: vv
416 !!$
417 !!$ vv = 0.d0
418 !!$ RETURN
419 !!$
420 !!$ END FUNCTION extension_velocity
421 
422  !===============================================================================
423  ! Boundary conditions for Maxwell
424  !===============================================================================
425 !!$ !===Velocity used in the induction equation.
426 !!$ !===Used only if problem type is mxw and restart velocity is false
427 !!$ FUNCTION Vexact(m, H_mesh) RESULT(vv) !Set uniquement a l'induction
428 !!$ IMPLICIT NONE
429 !!$ TYPE(mesh_type), INTENT(IN) :: H_mesh
430 !!$ INTEGER, INTENT(IN) :: m
431 !!$ REAL(KIND=8), DIMENSION(H_mesh%np,6) :: vv
432 !!$
433 !!$ vv = 0.d0
434 !!$ CALL error_petsc('Vexact: should not be called for this test')
435 !!$ END FUNCTION Vexact
436 
437 !!$ !===Magnetic field and magnetic induction for quasi-static approximation
438 !!$ !===if needed
439 !!$ FUNCTION H_B_quasi_static(char_h_b, rr, m) RESULT(vv)
440 !!$ IMPLICIT NONE
441 !!$ CHARACTER(LEN=1), INTENT(IN) :: char_h_b
442 !!$ REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: rr
443 !!$ INTEGER, INTENT(IN) :: m
444 !!$ REAL(KIND=8), DIMENSION(SIZE(rr,2),6) :: vv
445 !!$
446 !!$ vv = 0.d0
447 !!$ RETURN
448 !!$ END FUNCTION H_B_quasi_static
449 
450 !!$ !===Magnetic field for boundary conditions in the Maxwell equations.
451 !!$ FUNCTION Hexact(H_mesh, TYPE, rr, m, mu_H_field, t) RESULT(vv)
452 !!$ IMPLICIT NONE
453 !!$ TYPE(mesh_type), INTENT(IN) :: H_mesh
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(:), INTENT(IN) :: mu_H_field
459 !!$ REAL(KIND=8), DIMENSION(SIZE(rr,2)) :: vv
460 !!$ REAL(KIND=8), DIMENSION(SIZE(rr,2)) :: r
461 !!$ INTEGER :: n
462 !!$
463 !!$ vv=0.d0
464 !!$ CALL error_petsc('Hexact: should not be called for this test')
465 !!$ RETURN
466 !!$ END FUNCTION Hexact
467 
468 !!$ !===Scalar potential for boundary conditions in the Maxwell equations.
469 !!$ FUNCTION Phiexact(TYPE, rr, m, mu_phi,t) RESULT(vv)
470 !!$ IMPLICIT NONE
471 !!$ INTEGER , INTENT(IN) :: TYPE
472 !!$ REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: rr
473 !!$ INTEGER , INTENT(IN) :: m
474 !!$ REAL(KIND=8), INTENT(IN) :: mu_phi, t
475 !!$ REAL(KIND=8), DIMENSION(SIZE(rr,2)) :: vv
476 !!$
477 !!$ vv=0.d0
478 !!$ CALL error_petsc('Phiexact: should not be called for this test')
479 !!$ RETURN
480 !!$ END FUNCTION Phiexact
481 
482 !!$ !===Current in Ohm's law. Curl(H) = sigma(E + uxB) + current
483 !!$ FUNCTION Jexact_gauss(TYPE, rr, m, mu_phi, sigma, mu_H, t, mesh_id, opt_B_ext) 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) :: mu_phi, sigma, mu_H, t
489 !!$ INTEGER , INTENT(IN) :: mesh_id
490 !!$ REAL(KIND=8), DIMENSION(6), OPTIONAL,INTENT(IN) :: opt_B_ext
491 !!$ REAL(KIND=8) :: vv
492 !!$
493 !!$ vv=0.d0
494 !!$ CALL error_petsc('Jexact_gauss: should not be called for this test')
495 !!$ RETURN
496 !!$ END FUNCTION Jexact_gauss
497 
498 !!$ !===Electric field for Neumann BC (cf. doc)
499 !!$ FUNCTION Eexact_gauss(TYPE, rr, m, mu_phi, sigma, mu_H, t) RESULT(vv)
500 !!$ IMPLICIT NONE
501 !!$ INTEGER, INTENT(IN) :: TYPE
502 !!$ REAL(KIND=8), DIMENSION(:), INTENT(IN) :: rr
503 !!$ INTEGER, INTENT(IN) :: m
504 !!$ REAL(KIND=8), INTENT(IN) :: mu_phi, sigma, mu_H, t
505 !!$ REAL(KIND=8) :: vv
506 !!$
507 !!$ vv = 0.d0
508 !!$ CALL error_petsc('Eexact_gauss: should not be called for this test')
509 !!$ END FUNCTION Eexact_gauss
510 
511 !!$ !===Initialization of magnetic field and scalar potential (if present)
512 !!$ SUBROUTINE init_maxwell(H_mesh, phi_mesh, time, dt, mu_H_field, mu_phi, &
513 !!$ list_mode, Hn1, Hn, phin1, phin)
514 !!$ IMPLICIT NONE
515 !!$ TYPE(mesh_type) :: H_mesh, phi_mesh
516 !!$ REAL(KIND=8), INTENT(OUT):: time
517 !!$ REAL(KIND=8), INTENT(IN) :: dt
518 !!$ REAL(KIND=8), DIMENSION(:), INTENT(IN) :: mu_H_field
519 !!$ REAL(KIND=8), INTENT(IN) :: mu_phi
520 !!$ INTEGER, DIMENSION(:), INTENT(IN) :: list_mode
521 !!$ REAL(KIND=8), DIMENSION(:,:,:), INTENT(OUT):: Hn, Hn1
522 !!$ REAL(KIND=8), DIMENSION(:,:,:), INTENT(OUT):: phin, phin1
523 !!$ INTEGER :: i, k
524 !!$
525 !!$ time = -dt
526 !!$ DO k=1,6
527 !!$ DO i=1, SIZE(list_mode)
528 !!$ Hn1(:,k,i) = Hexact(H_mesh,k, H_mesh%rr, list_mode(i), mu_H_field, time)
529 !!$ IF (inputs%nb_dom_phi>0) THEN
530 !!$ IF (k<3) THEN
531 !!$ phin1(:,k,i) = Phiexact(k, phi_mesh%rr, list_mode(i) , mu_phi, time)
532 !!$ ENDIF
533 !!$ END IF
534 !!$ ENDDO
535 !!$ ENDDO
536 !!$
537 !!$ time = time + dt
538 !!$ DO k=1,6
539 !!$ DO i=1, SIZE(list_mode)
540 !!$ Hn(:,k,i) = Hexact(H_mesh,k, H_mesh%rr, list_mode(i), mu_H_field, time)
541 !!$ IF (inputs%nb_dom_phi>0) THEN
542 !!$ IF (k<3) THEN
543 !!$ phin(:,k,i) = Phiexact(k, phi_mesh%rr, list_mode(i), mu_phi, time)
544 !!$ ENDIF
545 !!$ END IF
546 !!$ ENDDO
547 !!$ ENDDO
548 !!$
549 !!$ END SUBROUTINE init_maxwell
550 
551 !!$ !===Analytical permeability (if needed)
552 !!$ !===This function is not needed unless the flag
553 !!$ !=== ===Use FEM Interpolation for magnetic permeability (true/false)
554 !!$ !===is activated and set to .FALSE. in the data data file. Default is .TRUE.
555 !!$ FUNCTION mu_bar_in_fourier_space(H_mesh,nb,ne,pts,pts_ids) RESULT(vv)
556 !!$ IMPLICIT NONE
557 !!$ TYPE(mesh_type), INTENT(IN) :: H_mesh
558 !!$ REAL(KIND=8), DIMENSION(ne-nb+1) :: vv
559 !!$ INTEGER, INTENT(IN) :: nb, ne
560 !!$ REAL(KIND=8),DIMENSION(2,ne-nb+1),OPTIONAL :: pts
561 !!$ INTEGER, DIMENSION(ne-nb+1), OPTIONAL :: pts_ids
562 !!$
563 !!$ vv = 1.d0
564 !!$ CALL error_petsc('mu_bar_in_fourier_space: should not be called for this test')
565 !!$ RETURN
566 !!$ END FUNCTION mu_bar_in_fourier_space
567 
568 !!$ !===Analytical mu_in_fourier_space (if needed)
569 !!$ !===This function is not needed unless the flag
570 !!$ !=== ===Use FEM Interpolation for magnetic permeability (true/false)
571 !!$ !===is activated and set to .FALSE. in the data data file. Default is .TRUE.
572 !!$ FUNCTION grad_mu_bar_in_fourier_space(pt,pt_id) RESULT(vv)
573 !!$ IMPLICIT NONE
574 !!$ REAL(KIND=8),DIMENSION(2), INTENT(in):: pt
575 !!$ INTEGER,DIMENSION(1), INTENT(in) :: pt_id
576 !!$ REAL(KIND=8),DIMENSION(2) :: vv
577 !!$
578 !!$ vv=0.d0
579 !!$ CALL error_petsc('grad_mu_bar_in_fourier_space: should not be called for this test')
580 !!$ RETURN
581 !!$ END FUNCTION grad_mu_bar_in_fourier_space
582 
583 !!$ !===Analytical permeability, mu in real space (if needed)
584 !!$ FUNCTION mu_in_real_space(H_mesh,angles,nb_angles,nb,ne,time) RESULT(vv)
585 !!$ IMPLICIT NONE
586 !!$ TYPE(mesh_type), INTENT(IN) :: H_mesh
587 !!$ REAL(KIND=8), DIMENSION(:), INTENT(IN) :: angles
588 !!$ INTEGER, INTENT(IN) :: nb_angles
589 !!$ INTEGER, INTENT(IN) :: nb, ne
590 !!$ REAL(KIND=8), INTENT(IN) :: time
591 !!$ REAL(KIND=8), DIMENSION(nb_angles,ne-nb+1) :: vv
592 !!$
593 !!$ vv = 1.d0
594 !!$ CALL error_petsc('mu_in_real_space: should not be called for this test')
595 !!$ RETURN
596 !!$ END FUNCTION mu_in_real_space
597 
598 END MODULE boundary_test_2
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 error_petsc(string)
Definition: my_util.f90:16
real(kind=8) function, dimension(size(rr, 2)), public pp_exact(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)
real(kind=8) function, dimension(size(rr, 2)), public vv_exact(TYPE, rr, m, t)