SFEMaNS  version 4.1 (work in progress)
Reference documentation for SFEMaNS
 All Classes Files Functions Variables Groups Pages
condlim_test_31.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 31')
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  RETURN
288 
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  RETURN
347  END FUNCTION vv_exact
348 
349 !!$ !===Solid velocity imposed when using penalty technique
350 !!$ !===Defined in Fourier space on mode 0 only.
351 !!$ FUNCTION imposed_velocity_by_penalty(rr,t) RESULT(vv)
352 !!$ IMPLICIT NONE
353 !!$ REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: rr
354 !!$ REAL(KIND=8), INTENT(IN) :: t
355 !!$ REAL(KIND=8), DIMENSION(SIZE(rr,2),6) :: vv
356 !!$
357 !!$ vv=0.d0
358 !!$ RETURN
359 !!$ END FUNCTION imposed_velocity_by_penalty
360 
361  !===Pressure for boundary conditions in Navier-Stokes.
362  !===Can be used also to initialize pressure in the subroutine init_velocity_pressure.
363  !===Use this routine for outflow BCs only.
364  !===CAUTION: Do not enfore BCs on pressure where normal component
365  ! of velocity is prescribed.
366  FUNCTION pp_exact(TYPE,rr,m,t) RESULT (vv)
367  IMPLICIT NONE
368  INTEGER , INTENT(IN) :: type
369  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: rr
370  INTEGER , INTENT(IN) :: m
371  REAL(KIND=8), INTENT(IN) :: t
372  REAL(KIND=8), DIMENSION(SIZE(rr,2)) :: vv, r, z
373  REAL(KIND=8) :: pi = 3.1415926535897932d0
374 
375  r = rr(1,:)
376  z = rr(2,:)
377 
378  IF ((type==1).AND.(m==1)) THEN
379  vv = r**3*sin(2*pi*z)*cos(t)
380  ELSE
381  vv = 0.d0
382  END IF
383  RETURN
384  END FUNCTION pp_exact
385 
386  !===Temperature for boundary conditions in temperature equation.
387  FUNCTION temperature_exact(TYPE,rr,m,t) RESULT (vv)
388  IMPLICIT NONE
389  INTEGER , INTENT(IN) :: type
390  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: rr
391  INTEGER , INTENT(IN) :: m
392  REAL(KIND=8), INTENT(IN) :: t
393  REAL(KIND=8), DIMENSION(SIZE(rr,2)) :: vv, r, z
394  REAL(KIND=8) :: r0=0.5d0
395  REAL(KIND=8) :: pi = 3.1415926535897932d0
396 
397  r = rr(1,:)
398  z = rr(2,:)
399 
400  IF ((type==1).AND.((m==0).OR.(m==1))) THEN
401  vv = r**2*(r-r0)**2*sin(2*pi*z)*cos(t)
402  ELSE
403  vv = 0.d0
404  END IF
405 
406  RETURN
407  END FUNCTION temperature_exact
408 
409 !!$ !===Can be used to initialize level set in the subroutine init_level_set.
410 !!$ FUNCTION level_set_exact(interface_nb,TYPE,rr,m,t) RESULT (vv)
411 !!$ IMPLICIT NONE
412 !!$ INTEGER , INTENT(IN) :: TYPE
413 !!$ REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: rr
414 !!$ INTEGER , INTENT(IN) :: m, interface_nb
415 !!$ REAL(KIND=8), INTENT(IN) :: t
416 !!$ REAL(KIND=8), DIMENSION(SIZE(rr,2)) :: vv
417 !!$
418 !!$ vv = 0.d0
419 !!$ RETURN
420 !!$ END FUNCTION level_set_exact
421 
422 !!$ !===Penalty coefficient (if needed)
423 !!$ !===This coefficient is equal to zero in subdomain
424 !!$ !===where penalty is applied
425 !!$ FUNCTION penal_in_real_space(mesh,rr_gauss,angles,nb_angles,nb,ne,time) RESULT(vv)
426 !!$ IMPLICIT NONE
427 !!$ TYPE(mesh_type) :: mesh
428 !!$ REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: rr_gauss
429 !!$ REAL(KIND=8), DIMENSION(:), INTENT(IN) :: angles
430 !!$ INTEGER, INTENT(IN) :: nb_angles
431 !!$ INTEGER, INTENT(IN) :: nb, ne
432 !!$ REAL(KIND=8), INTENT(IN) :: time
433 !!$ REAL(KIND=8), DIMENSION(nb_angles,ne-nb+1) :: vv
434 !!$
435 !!$ vv = 1.d0
436 !!$ RETURN
437 !!$ END FUNCTION penal_in_real_space
438 
439  !===Extension of the velocity field in the solid.
440  !===Used when temperature or Maxwell equations are solved.
441  !===It extends the velocity field on the Navier-Stokes domain to a
442  !===velocity field on the temperature and the Maxwell domain.
443  !===It is also used if problem type=mxw and restart velocity
444  !===is set to true in data (type problem denoted mxx in the code).
445  FUNCTION extension_velocity(TYPE, H_mesh, mode, t, n_start) RESULT(vv)
446  IMPLICIT NONE
447  TYPE(mesh_type), INTENT(IN) :: h_mesh
448  INTEGER , INTENT(IN) :: type, n_start
449  INTEGER, INTENT(IN) :: mode
450  REAL(KIND=8), INTENT(IN) :: t
451  REAL(KIND=8), DIMENSION(H_Mesh%np) :: vv
452  REAL(KIND=8) :: r
453  INTEGER :: n
454 
455  vv = 0.d0
456  RETURN
457 
458  !===Dummies variables to avoid warning
459  n=h_mesh%np; r=t; n=type; n=mode; n=n_start
460  !===Dummies variables to avoid warning
461  END FUNCTION extension_velocity
462 
463  !===============================================================================
464  ! Boundary conditions for Maxwell
465  !===============================================================================
466 !!$ !===Velocity used in the induction equation.
467 !!$ !===Used only if problem type is mxw and restart velocity is false
468 !!$ FUNCTION Vexact(m, H_mesh) RESULT(vv) !Set uniquement a l'induction
469 !!$ IMPLICIT NONE
470 !!$ TYPE(mesh_type), INTENT(IN) :: H_mesh
471 !!$ INTEGER, INTENT(IN) :: m
472 !!$ REAL(KIND=8), DIMENSION(H_mesh%np,6) :: vv
473 !!$
474 !!$ vv = 0.d0
475 !!$ RETURN
476 !!$ END FUNCTION Vexact
477 
478 !!$ FUNCTION H_B_quasi_static(char_h_b, rr, m) RESULT(vv)
479 !!$ IMPLICIT NONE
480 !!$ CHARACTER(LEN=1), INTENT(IN) :: char_h_b
481 !!$ REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: rr
482 !!$ INTEGER, INTENT(IN) :: m
483 !!$ REAL(KIND=8), DIMENSION(SIZE(rr,2),6) :: vv
484 !!$
485 !!$ IF (inputs%if_quasi_static_approx) THEN
486 !!$ vv = 0.d0
487 !!$ ELSE
488 !!$ CALL error_petsc('H_B_quasi_static should not be called')
489 !!$ END IF
490 !!$ RETURN
491 !!$ END FUNCTION H_B_quasi_static
492 
493 !!$ !===Magnetic field for boundary conditions in the Maxwell equations.
494 !!$ FUNCTION Hexact(H_mesh, TYPE, rr, m, mu_H_field, t) RESULT(vv)
495 !!$ IMPLICIT NONE
496 !!$ TYPE(mesh_type), INTENT(IN) :: H_mesh
497 !!$ INTEGER , INTENT(IN) :: TYPE
498 !!$ REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: rr
499 !!$ INTEGER , INTENT(IN) :: m
500 !!$ REAL(KIND=8), INTENT(IN) :: t
501 !!$ REAL(KIND=8), DIMENSION(:), INTENT(IN) :: mu_H_field
502 !!$ REAL(KIND=8), DIMENSION(SIZE(rr,2)) :: vv
503 !!$
504 !!$ vv = 0.d0
505 !!$ RETURN
506 !!$ END FUNCTION Hexact
507 
508 !!$ !===Scalar potential for boundary conditions in the Maxwell equations.
509 !!$ FUNCTION Phiexact(TYPE, rr, m, mu_phi,t) RESULT(vv)
510 !!$ IMPLICIT NONE
511 !!$ INTEGER , INTENT(IN) :: TYPE
512 !!$ REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: rr
513 !!$ INTEGER , INTENT(IN) :: m
514 !!$ REAL(KIND=8), INTENT(IN) :: mu_phi, t
515 !!$ REAL(KIND=8), DIMENSION(SIZE(rr,2)) :: vv
516 !!$
517 !!$ vv=0.d0
518 !!$ RETURN
519 !!$ END FUNCTION Phiexact
520 
521 !!$ !===Current in Ohm's law. Curl(H) = sigma(E + uxB) + current
522 !!$ FUNCTION Jexact_gauss(TYPE, rr, m, mu_phi, sigma, mu_H, t, mesh_id, opt_B_ext) RESULT(vv)
523 !!$ IMPLICIT NONE
524 !!$ INTEGER , INTENT(IN) :: TYPE
525 !!$ REAL(KIND=8), DIMENSION(:), INTENT(IN) :: rr
526 !!$ INTEGER , INTENT(IN) :: m
527 !!$ REAL(KIND=8), INTENT(IN) :: mu_phi, sigma, mu_H, t
528 !!$ INTEGER , INTENT(IN) :: mesh_id
529 !!$ REAL(KIND=8), DIMENSION(6), OPTIONAL,INTENT(IN) :: opt_B_ext
530 !!$ REAL(KIND=8) :: vv
531 !!$
532 !!$ vv = 0.d0
533 !!$ RETURN
534 !!$ END FUNCTION Jexact_gauss
535 
536 !!$ !===Electric field for Neumann BC (cf. doc)
537 !!$ FUNCTION Eexact_gauss(TYPE, rr, m, mu_phi, sigma, mu_H, t) RESULT(vv)
538 !!$ IMPLICIT NONE
539 !!$ INTEGER, INTENT(IN) :: TYPE
540 !!$ REAL(KIND=8), DIMENSION(:), INTENT(IN) :: rr
541 !!$ INTEGER, INTENT(IN) :: m
542 !!$ REAL(KIND=8), INTENT(IN) :: mu_phi, sigma, mu_H, t
543 !!$ REAL(KIND=8) :: vv
544 !!$
545 !!$ vv = 0.d0
546 !!$ RETURN
547 !!$ END FUNCTION Eexact_gauss
548 
549 !!$ !===Initialization of magnetic field and scalar potential (if present)
550 !!$ SUBROUTINE init_maxwell(H_mesh, phi_mesh, time, dt, mu_H_field, mu_phi, &
551 !!$ list_mode, Hn1, Hn, phin1, phin)
552 !!$ IMPLICIT NONE
553 !!$ TYPE(mesh_type) :: H_mesh, phi_mesh
554 !!$ REAL(KIND=8), INTENT(OUT):: time
555 !!$ REAL(KIND=8), INTENT(IN) :: dt
556 !!$ REAL(KIND=8), DIMENSION(:), INTENT(IN) :: mu_H_field
557 !!$ REAL(KIND=8), INTENT(IN) :: mu_phi
558 !!$ INTEGER, DIMENSION(:), INTENT(IN) :: list_mode
559 !!$ REAL(KIND=8), DIMENSION(:,:,:), INTENT(OUT):: Hn, Hn1
560 !!$ REAL(KIND=8), DIMENSION(:,:,:), INTENT(OUT):: phin, phin1
561 !!$
562 !!$ time = 0.d0
563 !!$ Hn1=0.d0
564 !!$ HN=0.d0
565 !!$ phin = 0.d0
566 !!$ phin1 =0.d0
567 !!$ RETURN
568 !!$ END SUBROUTINE init_maxwell
569 
570 !!$ !===Analytical permeability (if needed)
571 !!$ !===This function is not needed unless the flag
572 !!$ !=== ===Use FEM Interpolation for magnetic permeability (true/false)
573 !!$ !===is activated and set to .FALSE. in the data data file. Default is .TRUE.
574 !!$ FUNCTION mu_bar_in_fourier_space(H_mesh,nb,ne,pts,pts_ids) RESULT(vv)
575 !!$ IMPLICIT NONE
576 !!$ TYPE(mesh_type), INTENT(IN) :: H_mesh
577 !!$ REAL(KIND=8), DIMENSION(ne-nb+1) :: vv
578 !!$ INTEGER, INTENT(IN) :: nb, ne
579 !!$ REAL(KIND=8),DIMENSION(2,ne-nb+1),OPTIONAL :: pts
580 !!$ INTEGER, DIMENSION(ne-nb+1), OPTIONAL :: pts_ids
581 !!$
582 !!$ vv = 1.d0
583 !!$ RETURN
584 !!$ END FUNCTION mu_bar_in_fourier_space
585 
586 !!$ !===Analytical mu_in_fourier_space (if needed)
587 !!$ !===This function is not needed unless the flag
588 !!$ !=== ===Use FEM Interpolation for magnetic permeability (true/false)
589 !!$ !===is activated and set to .FALSE. in the data data file. Default is .TRUE.
590 !!$ FUNCTION grad_mu_bar_in_fourier_space(pt,pt_id) RESULT(vv)
591 !!$ IMPLICIT NONE
592 !!$ REAL(KIND=8),DIMENSION(2), INTENT(in):: pt
593 !!$ INTEGER,DIMENSION(1), INTENT(in) :: pt_id
594 !!$ REAL(KIND=8),DIMENSION(2) :: vv
595 !!$
596 !!$ vv=0.d0
597 !!$ RETURN
598 !!$ END FUNCTION grad_mu_bar_in_fourier_space
599 
600 !!$ !===Analytical permeability, mu in real space (if needed)
601 !!$ FUNCTION mu_in_real_space(H_mesh,angles,nb_angles,nb,ne,time) RESULT(vv)
602 !!$ IMPLICIT NONE
603 !!$ TYPE(mesh_type), INTENT(IN) :: H_mesh
604 !!$ REAL(KIND=8), DIMENSION(:), INTENT(IN) :: angles
605 !!$ INTEGER, INTENT(IN) :: nb_angles
606 !!$ INTEGER, INTENT(IN) :: nb, ne
607 !!$ REAL(KIND=8), INTENT(IN) :: time
608 !!$ REAL(KIND=8), DIMENSION(nb_angles,ne-nb+1) :: vv
609 !!$
610 !!$ vv = 1.d0
611 !!$ RETURN
612 !!$ END FUNCTION mu_in_real_space
613 
614 END MODULE boundary_test_31
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