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