SFEMaNS  version 4.1 (work in progress)
Reference documentation for SFEMaNS
 All Classes Files Functions Variables Groups Pages
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
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
real(kind=8) function, dimension(size(rr, 2)), public source_in_ns_momentum(TYPE, rr, mode, i, time, Re, ty, opt_density, opt_tempn)
real(kind=8) function, dimension(size(rr, 2)), public pp_exact(TYPE, rr, m, t)
real(kind=8) function, dimension(size(rr, 2)), public vv_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)
subroutine error_petsc(string)
Definition: my_util.f90:15
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