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