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