SFEMaNS  version 5.3
Reference documentation for SFEMaNS
condlim_test_14.f90
Go to the documentation of this file.
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  REAL(KIND=8) :: pi=acos(-1.d0)
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 !!$
129 !!$ vv = 0.d0
130 !!$ CALL error_petsc('source_in_NS_momentum: should not be called for this test')
131 !!$ RETURN
132 !!$ END FUNCTION source_in_NS_momentum
133 
134 !!$ !===Extra source in temperature equation. Always called.
135 !!$ FUNCTION source_in_temperature(TYPE, rr, m, t)RESULT(vv)
136 !!$ IMPLICIT NONE
137 !!$ INTEGER , INTENT(IN) :: TYPE
138 !!$ REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: rr
139 !!$ INTEGER , INTENT(IN) :: m
140 !!$ REAL(KIND=8), INTENT(IN) :: t
141 !!$ REAL(KIND=8), DIMENSION(SIZE(rr,2)) :: vv
142 !!$
143 !!$ vv = 0.d0
144 !!$ CALL error_petsc('source_in_temperature: should not be called for this test')
145 !!$ RETURN
146 !!$ END FUNCTION source_in_temperature
147 
148 !!$ !===Extra source in level set equation. Always called.
149 !!$ FUNCTION source_in_level_set(interface_nb,TYPE, rr, m, t)RESULT(vv)
150 !!$ IMPLICIT NONE
151 !!$ INTEGER , INTENT(IN) :: TYPE
152 !!$ REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: rr
153 !!$ INTEGER , INTENT(IN) :: m, interface_nb
154 !!$ REAL(KIND=8), INTENT(IN) :: t
155 !!$ REAL(KIND=8), DIMENSION(SIZE(rr,2)) :: vv
156 !!$
157 !!$ vv=0.d0
158 !!$ CALL error_petsc('sourece_in_temperature: should not be called for this test')
159 !!$ END FUNCTION source_in_level_set
160 
161 !!$ !===Velocity for boundary conditions in Navier-Stokes.
162 !!$ !===Can be used also to initialize velocity in: init_velocity_pressure_temperature
163 !!$ FUNCTION vv_exact(TYPE,rr,m,t) RESULT(vv)
164 !!$
165 !!$ IMPLICIT NONE
166 !!$ INTEGER , INTENT(IN) :: TYPE
167 !!$ REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: rr
168 !!$ INTEGER, INTENT(IN) :: m
169 !!$ REAL(KIND=8), INTENT(IN) :: t
170 !!$ REAL(KIND=8), DIMENSION(SIZE(rr,2)) :: vv
171 !!$
172 !!$ vv(:) = 0.d0
173 !!$ CALL error_petsc('vv_exact: should not be called for this test')
174 !!$ RETURN
175 !!$ END FUNCTION vv_exact
176 
177 !!$ !===Solid velocity imposed when using penalty technique
178 !!$ !===Defined in Fourier space on mode 0 only.
179 !!$ FUNCTION imposed_velocity_by_penalty(rr,t) RESULT(vv)
180 !!$ IMPLICIT NONE
181 !!$ REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: rr
182 !!$ REAL(KIND=8), INTENT(IN) :: t
183 !!$ REAL(KIND=8), DIMENSION(SIZE(rr,2),6) :: vv
184 !!$
185 !!$ vv=0.d0
186 !!$ RETURN
187 !!$ END FUNCTION imposed_velocity_by_penalty
188 
189 !!$ !===Pressure for boundary conditions in Navier-Stokes.
190 !!$ !===Can be used also to initialize pressure in the subroutine init_velocity_pressure.
191 !!$ !===Use this routine for outflow BCs only.
192 !!$ !===CAUTION: Do not enfore BCs on pressure where normal component
193 !!$ ! of velocity is prescribed.
194 !!$ FUNCTION pp_exact(TYPE,rr,m,t) RESULT (vv)
195 !!$ IMPLICIT NONE
196 !!$ INTEGER , INTENT(IN) :: TYPE
197 !!$ REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: rr
198 !!$ INTEGER , INTENT(IN) :: m
199 !!$ REAL(KIND=8), INTENT(IN) :: t
200 !!$ REAL(KIND=8), DIMENSION(SIZE(rr,2)) :: vv
201 !!$
202 !!$ vv=0.d0
203 !!$ CALL error_petsc('pp_exact: should not be called for this test')
204 !!$ RETURN
205 !!$ END FUNCTION pp_exact
206 
207 !!$ !===Temperature for boundary conditions in temperature equation.
208 !!$ FUNCTION temperature_exact(TYPE,rr,m,t) RESULT (vv)
209 !!$ IMPLICIT NONE
210 !!$ INTEGER , INTENT(IN) :: TYPE
211 !!$ REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: rr
212 !!$ INTEGER , INTENT(IN) :: m
213 !!$ REAL(KIND=8), INTENT(IN) :: t
214 !!$ REAL(KIND=8), DIMENSION(SIZE(rr,2)) :: vv
215 !!$
216 !!$ vv = 0.d0
217 !!$ CALL error_petsc('temperature_exact: should not be called for this test')
218 !!$ RETURN
219 !!$ END FUNCTION temperature_exact
220 
221 !!$ !===Can be used to initialize level set in the subroutine init_level_set.
222 !!$ FUNCTION level_set_exact(interface_nb,TYPE,rr,m,t) RESULT (vv)
223 !!$ IMPLICIT NONE
224 !!$ INTEGER , INTENT(IN) :: TYPE
225 !!$ REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: rr
226 !!$ INTEGER , INTENT(IN) :: m, interface_nb
227 !!$ REAL(KIND=8), INTENT(IN) :: t
228 !!$ REAL(KIND=8), DIMENSION(SIZE(rr,2)) :: vv
229 !!$
230 !!$ vv = 0.d0
231 !!$ CALL error_petsc('level_set_exact: should not be called for this test')
232 !!$ RETURN
233 !!$
234 !!$ END FUNCTION level_set_exact
235 
236 !!$ !===Penalty coefficient (if needed)
237 !!$ !===This coefficient is equal to zero in subdomain
238 !!$ !===where penalty is applied (penalty is zero in solid)
239 !!$ FUNCTION penal_in_real_space(mesh,rr_gauss,angles,nb_angles,nb,ne,time) RESULT(vv)
240 !!$ IMPLICIT NONE
241 !!$ TYPE(mesh_type) :: mesh
242 !!$ REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: rr_gauss
243 !!$ REAL(KIND=8), DIMENSION(:), INTENT(IN) :: angles
244 !!$ INTEGER, INTENT(IN) :: nb_angles
245 !!$ INTEGER, INTENT(IN) :: nb, ne
246 !!$ REAL(KIND=8), INTENT(IN) :: time
247 !!$ REAL(KIND=8), DIMENSION(nb_angles,ne-nb+1) :: vv
248 !!$
249 !!$ vv = 1.d0
250 !!$ CALL error_petsc('penal_in_real_space: should not be called for this test')
251 !!$ RETURN
252 !!$ END FUNCTION penal_in_real_space
253 
254 !!$ !===Extension of the velocity field in the solid.
255 !!$ !===Used when temperature or Maxwell equations are solved.
256 !!$ !===It extends the velocity field on the Navier-Stokes domain to a
257 !!$ !===velocity field on the temperature and the Maxwell domain.
258 !!$ !===It is also used if problem type=mxw and restart velocity
259 !!$ !===is set to true in data (type problem denoted mxx in the code).
260 !!$ FUNCTION extension_velocity(TYPE, H_mesh, mode, t, n_start) RESULT(vv)
261 !!$ IMPLICIT NONE
262 !!$ TYPE(mesh_type), INTENT(IN) :: H_mesh
263 !!$ INTEGER , INTENT(IN) :: TYPE, n_start
264 !!$ INTEGER, INTENT(IN) :: mode
265 !!$ REAL(KIND=8), INTENT(IN) :: t
266 !!$ REAL(KIND=8), DIMENSION(H_Mesh%np) :: vv
267 !!$
268 !!$ vv = 0.d0
269 !!$ RETURN
270 !!$
271 !!$ END FUNCTION extension_velocity
272 
273  !===============================================================================
274  ! Boundary conditions for Maxwell
275  !===============================================================================
276  !===Velocity used in the induction equation.
277  !===Used only if problem type is mxw and restart velocity is false
278  FUNCTION vexact(m, H_mesh) RESULT(vv) !Set uniquement a l'induction
279  IMPLICIT NONE
280  TYPE(mesh_type), INTENT(IN) :: H_mesh
281  INTEGER, INTENT(IN) :: m
282  REAL(KIND=8), DIMENSION(H_mesh%np,6) :: vv
283  LOGICAL, DIMENSION(H_mesh%np) :: virgin
284  INTEGER :: i, mjl, njl
285  REAL(Kind=8) :: rr, zz
286  REAL(KIND=8) :: eps=1.d-5, height=2.d0, zeta=30.d0, amp_mnd=1.d0, amp_fl=0.d0
287 
288  vv = 0.d0
289  IF (m==0) THEN
290  virgin = .true.
291  DO mjl = 1, h_mesh%me
292  !IF (H_mesh%i_d(mjl)/= 4) CYCLE
293  !We are in the sodium
294  DO njl = 1, h_mesh%gauss%n_w
295  i = h_mesh%jj(njl,mjl)
296  IF (.NOT.virgin(i)) cycle
297  virgin(i) = .false.
298 
299  rr = h_mesh%rr(1,i)
300  zz = h_mesh%rr(2,i)
301 
302  vv(i,1) = amp_mnd*(-(pi/height)*rr*(1.d0-rr)**2*(1.d0+2.d0*rr)*cos(2.d0*pi*zz/height))
303  vv(i,3) = amp_mnd*(4.d0*height/pi)*rr*(1.d0-rr)*asin(zz)
304  IF (zz .LE. (-eps)) THEN ! On est en bas
305  vv(i,3) = vv(i,3)+ amp_fl*rr*sin(pi*rr)*(1.-2.*zeta*(zz+1.)**2)*exp(-zeta*(zz+1.)**2)
306  ELSE IF (zz .GE. eps) THEN ! On est en haut
307  vv(i,3) = vv(i,3)+ amp_fl*rr*sin(pi*rr)*(1.-2.*zeta*(zz-1.)**2)*exp(-zeta*(zz-1.)**2)
308  ENDIF
309  vv(i,5) = amp_mnd*(1.-rr)*(1.+rr-5.*rr**2)*sin(2.*pi*zz/height)
310 !!$ vel_loc(i) = sqrt(vv(i,1)**2 + vv(i,3)**2 + vv(i,5)**2)
311  END DO !njl
312  END DO
313  END IF
314  END FUNCTION vexact
315 
316 !!$ !===Magnetic field and magnetic induction for quasi-static approximation
317 !!$ !===if needed
318 !!$ FUNCTION H_B_quasi_static(char_h_b, rr, m) RESULT(vv)
319 !!$ IMPLICIT NONE
320 !!$ CHARACTER(LEN=1), INTENT(IN) :: char_h_b
321 !!$ REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: rr
322 !!$ INTEGER, INTENT(IN) :: m
323 !!$ REAL(KIND=8), DIMENSION(SIZE(rr,2),6) :: vv
324 !!$
325 !!$ vv = 0.d0
326 !!$ RETURN
327 !!$ END FUNCTION H_B_quasi_static
328 
329  !===Magnetic field for boundary conditions in the Maxwell equations.
330  FUNCTION hexact(H_mesh,TYPE, rr, m, mu_H_field, t) RESULT(vv)
331  IMPLICIT NONE
332  TYPE(mesh_type), INTENT(IN) :: H_mesh
333  INTEGER , INTENT(IN) :: TYPE
334  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: rr
335  INTEGER , INTENT(IN) :: m
336  REAL(KIND=8), INTENT(IN) :: t
337  REAL(KIND=8), DIMENSION(:), INTENT(IN) :: mu_H_field
338  REAL(KIND=8), DIMENSION(SIZE(rr,2)) :: vv
339  REAL(KIND=8) :: r
340  INTEGER :: n
341 
342  vv = 0.d0
343  RETURN
344 
345  !===Dummies variables to avoid warning
346  n=h_mesh%np; n=type; n=SIZE(rr,1); n=m; r=mu_h_field(1); r=t
347  !===Dummies variables to avoid warning
348  END FUNCTION hexact
349 
350  !===Scalar potential for boundary conditions in the Maxwell equations.
351  FUNCTION phiexact(TYPE, rr, m, mu_phi,t) RESULT(vv)
352  IMPLICIT NONE
353  INTEGER , INTENT(IN) :: TYPE
354  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: rr
355  INTEGER , INTENT(IN) :: m
356  REAL(KIND=8), INTENT(IN) :: mu_phi, t
357  REAL(KIND=8), DIMENSION(SIZE(rr,2)) :: vv
358  REAL(KIND=8) :: r
359  INTEGER :: n
360 
361  vv = 0.d0
362  CALL error_petsc('Phiexact: should not be called for this test')
363  RETURN
364 
365  !===Dummies variables to avoid warning
366  n=type; n=SIZE(rr,1); n=m; r=mu_phi; r=t
367  !===Dummies variables to avoid warning
368  END FUNCTION phiexact
369 
370  !===Current in Ohm's law. Curl(H) = sigma(E + uxB) + current
371  FUNCTION jexact_gauss(TYPE, rr, m, mu_phi, sigma, mu_H, t, mesh_id, opt_B_ext) RESULT(vv)
372  IMPLICIT NONE
373  INTEGER , INTENT(IN) :: TYPE
374  REAL(KIND=8), DIMENSION(:), INTENT(IN) :: rr
375  INTEGER , INTENT(IN) :: m
376  REAL(KIND=8), INTENT(IN) :: mu_phi, sigma, mu_H, t
377  INTEGER , INTENT(IN) :: mesh_id
378  REAL(KIND=8), DIMENSION(6), OPTIONAL,INTENT(IN) :: opt_B_ext
379  REAL(KIND=8) :: vv
380  REAL(KIND=8) :: r
381  INTEGER :: n
382 
383  vv = 0.d0
384  RETURN
385 
386  !===Dummies variables to avoid warning
387  r=rr(1); r=mu_phi; r=sigma; r=mu_h; r=t; n=type; n=m; n=mesh_id
388  IF (PRESENT(opt_b_ext)) r=opt_b_ext(1)
389  !===Dummies variables to avoid warning
390  END FUNCTION jexact_gauss
391 
392  !===Electric field for Neumann BC (cf. doc)
393  FUNCTION eexact_gauss(TYPE, rr, m, mu_phi, sigma, mu_H, t) RESULT(vv)
394  IMPLICIT NONE
395  INTEGER, INTENT(IN) :: TYPE
396  REAL(KIND=8), DIMENSION(:), INTENT(IN) :: rr
397  INTEGER, INTENT(IN) :: m
398  REAL(KIND=8), INTENT(IN) :: mu_phi, sigma, mu_H, t
399  REAL(KIND=8) :: vv
400  REAL(KIND=8) :: r
401  INTEGER :: n
402 
403  vv = 0.d0
404  RETURN
405 
406  !===Dummies variables to avoid warning
407  r=rr(1); r=mu_phi; r=sigma; r=mu_h; r=t; n=type; n=m
408  !===Dummies variables to avoid warning
409  END FUNCTION eexact_gauss
410 
411  !===Initialization of magnetic field and scalar potential (if present)
412  SUBROUTINE init_maxwell(H_mesh, phi_mesh, time, dt, mu_H_field, mu_phi, &
413  list_mode, hn1, hn, phin1, phin)
414  IMPLICIT NONE
415  TYPE(mesh_type) :: H_mesh, phi_mesh
416  REAL(KIND=8), INTENT(OUT):: time
417  REAL(KIND=8), INTENT(IN) :: dt
418  REAL(KIND=8), DIMENSION(:), INTENT(IN) :: mu_H_field
419  REAL(KIND=8), INTENT(IN) :: mu_phi
420  INTEGER, DIMENSION(:), INTENT(IN) :: list_mode
421  REAL(KIND=8), DIMENSION(:,:,:), INTENT(OUT):: Hn, Hn1
422  REAL(KIND=8), DIMENSION(:,:,:), INTENT(OUT):: phin, phin1
423  INTEGER :: i, n0
424  REAL(KIND=8) :: Ri, Ro, A
425 
426  time = 0.d0
427  hn = 0.d0
428  hn1 = 0.d0
429  phin = 0.d0
430  phin1 = 0.d0
431  n0 = 1
432  ri = 1.d0
433  ro = 2.d0
434  a = 1d-3
435  IF (h_mesh%me /= 0) THEN
436  DO i = 1, SIZE(list_mode)
437  IF (list_mode(i) == 1) THEN
438  hn1(:,1,i) = a*(sin(n0*h_mesh%rr(2,:))/h_mesh%rr(1,:))*(ri-h_mesh%rr(1,:))**2 &
439  * (ro-h_mesh%rr(1,:))**2
440  hn(:,1,i) = a*(sin(n0*h_mesh%rr(2,:))/h_mesh%rr(1,:))*(ri-h_mesh%rr(1,:))**2 &
441  * (ro-h_mesh%rr(1,:))**2
442  hn1(:,4,i) = -2*a*sin(n0*h_mesh%rr(2,:))*(h_mesh%rr(1,:)-ri)&
443  *(h_mesh%rr(1,:)-ro)*(2.d0*h_mesh%rr(1,:)-ro-ri)
444  hn(:,4,i) = -2*a*sin(n0*h_mesh%rr(2,:))*(h_mesh%rr(1,:)-ri)&
445  *(h_mesh%rr(1,:)-ro)*(2.d0*h_mesh%rr(1,:)-ro-ri)
446  END IF
447  END DO
448  END IF
449  RETURN
450 
451  !===Dummies variables to avoid warning
452  i=phi_mesh%np; a=time; a=dt; a=mu_h_field(1); a=mu_phi
453  !===Dummies variables to avoid warning
454  END SUBROUTINE init_maxwell
455 
456 !!$ !===Analytical permeability (if needed)
457 !!$ !===This function is not needed unless the flag
458 !!$ !=== ===Use FEM Interpolation for magnetic permeability (true/false)
459 !!$ !===is activated and set to .FALSE. in the data data file. Default is .TRUE.
460 !!$ FUNCTION mu_bar_in_fourier_space(H_mesh,nb,ne,pts,pts_ids) RESULT(vv)
461 !!$ IMPLICIT NONE
462 !!$ TYPE(mesh_type), INTENT(IN) :: H_mesh
463 !!$ REAL(KIND=8), DIMENSION(ne-nb+1) :: vv
464 !!$ INTEGER, INTENT(IN) :: nb, ne
465 !!$ REAL(KIND=8),DIMENSION(2,ne-nb+1),OPTIONAL :: pts
466 !!$ INTEGER, DIMENSION(ne-nb+1), OPTIONAL :: pts_ids
467 !!$
468 !!$ vv = 1.d0
469 !!$ CALL error_petsc('mu_bar_in_fourier_space: should not be called for this test')
470 !!$ RETURN
471 !!$ END FUNCTION mu_bar_in_fourier_space
472 
473 !!$ !===Analytical mu_in_fourier_space (if needed)
474 !!$ !===This function is not needed unless the flag
475 !!$ !=== ===Use FEM Interpolation for magnetic permeability (true/false)
476 !!$ !===is activated and set to .FALSE. in the data data file. Default is .TRUE.
477 !!$ FUNCTION grad_mu_bar_in_fourier_space(pt,pt_id) RESULT(vv)
478 !!$ IMPLICIT NONE
479 !!$ REAL(KIND=8),DIMENSION(2), INTENT(in):: pt
480 !!$ INTEGER,DIMENSION(1), INTENT(in) :: pt_id
481 !!$ REAL(KIND=8),DIMENSION(2) :: vv
482 !!$
483 !!$ vv = 0.d0
484 !!$ CALL error_petsc('grad_mu_bar_in_fourier_space: should not be called for this test')
485 !!$ RETURN
486 !!$ END FUNCTION grad_mu_bar_in_fourier_space
487 
488 !!$ !===Analytical permeability, mu in real space (if needed)
489 !!$ FUNCTION mu_in_real_space(H_mesh,angles,nb_angles,nb,ne,time) RESULT(vv)
490 !!$ IMPLICIT NONE
491 !!$ TYPE(mesh_type), INTENT(IN) :: H_mesh
492 !!$ REAL(KIND=8), DIMENSION(:), INTENT(IN) :: angles
493 !!$ INTEGER, INTENT(IN) :: nb_angles
494 !!$ INTEGER, INTENT(IN) :: nb, ne
495 !!$ REAL(KIND=8), INTENT(IN) :: time
496 !!$ REAL(KIND=8), DIMENSION(nb_angles,ne-nb+1) :: vv
497 !!$
498 !!$ vv = 1.d0
499 !!$ CALL error_petsc('mu_in_real_space: should not be called for this test')
500 !!$ RETURN
501 !!$ END FUNCTION mu_in_real_space
502 
503 END MODULE boundary_test_14
real(kind=8) function, public eexact_gauss(TYPE, rr, m, mu_phi, sigma, mu_H, t)
subroutine error_petsc(string)
Definition: my_util.f90:16
real(kind=8) function, dimension(size(rr, 2)), public phiexact(TYPE, rr, m, mu_phi, t)
real(kind=8) function, public jexact_gauss(TYPE, rr, m, mu_phi, sigma, mu_H, t, mesh_id, opt_B_ext)
subroutine, public init_maxwell(H_mesh, phi_mesh, time, dt, mu_H_field, mu_phi, list_mode, Hn1, Hn, phin1, phin)
real(kind=8) function, dimension(size(rr, 2)), public hexact(H_mesh, TYPE, rr, m, mu_H_field, t)
real(kind=8) function, dimension(h_mesh%np, 6), public vexact(m, H_mesh)