SFEMaNS  version 4.1 (work in progress)
Reference documentation for SFEMaNS
 All Classes Files Functions Variables Groups Pages
condlim_test_33.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  PUBLIC :: kelvin_force_coeff
33  PRIVATE
34 
35 CONTAINS
36  !===============================================================================
37  ! Boundary conditions for Navier-Stokes
38  !===============================================================================
39 
40  !===Initialize velocity, pressure
41  SUBROUTINE init_velocity_pressure(mesh_f, mesh_c, time, dt, list_mode, &
42  un_m1, un, pn_m1, pn, phin_m1, phin)
43  IMPLICIT NONE
44  TYPE(mesh_type) :: mesh_f, mesh_c
45  REAL(KIND=8), INTENT(OUT):: time
46  REAL(KIND=8), INTENT(IN) :: dt
47  INTEGER, DIMENSION(:), INTENT(IN) :: list_mode
48  REAL(KIND=8), DIMENSION(:,:,:), INTENT(OUT):: un_m1, un
49  REAL(KIND=8), DIMENSION(:,:,:), INTENT(OUT):: pn_m1, pn, phin_m1, phin
50  INTEGER :: mode, i, j
51  REAL(KIND=8), DIMENSION(mesh_c%np) :: pn_m2
52 
53  time = 0.d0
54  DO i= 1, SIZE(list_mode)
55  mode = list_mode(i)
56  DO j = 1, 6
57  !===velocity
58  un_m1(:,j,i) = vv_exact(j,mesh_f%rr,mode,time-dt)
59  un(:,j,i) = vv_exact(j,mesh_f%rr,mode,time)
60  END DO
61  DO j = 1, 2
62  !===pressure
63  pn_m2(:) = pp_exact(j,mesh_c%rr,mode,time-2*dt)
64  pn_m1(:,j,i) = pp_exact(j,mesh_c%rr,mode,time-dt)
65  pn(:,j,i) = pp_exact(j,mesh_c%rr,mode,time)
66  phin_m1(:,j,i) = pn_m1(:,j,i) - pn_m2(:)
67  phin(:,j,i) = pn(:,j,i) - pn_m1(:,j,i)
68  ENDDO
69  ENDDO
70  END SUBROUTINE init_velocity_pressure
71 
72  !===Initialize temperature
73  SUBROUTINE init_temperature(mesh, time, dt, list_mode, tempn_m1, tempn)
74  IMPLICIT NONE
75  TYPE(mesh_type) :: mesh
76  REAL(KIND=8), INTENT(OUT):: time
77  REAL(KIND=8), INTENT(IN) :: dt
78  INTEGER, DIMENSION(:), INTENT(IN) :: list_mode
79  REAL(KIND=8), DIMENSION(:,:,:), INTENT(OUT):: tempn_m1, tempn
80  INTEGER :: mode, i, j
81 
82  time = 0.d0
83  DO i= 1, SIZE(list_mode)
84  mode = list_mode(i)
85  DO j = 1, 2
86  tempn_m1(:,j,i) = temperature_exact(j, mesh%rr, mode, time-dt)
87  tempn(:,j,i) = temperature_exact(j, mesh%rr, mode, time)
88  ENDDO
89  ENDDO
90  END SUBROUTINE init_temperature
91 
92 !!$ !===Initialize level_set
93 !!$ SUBROUTINE init_level_set(vv_mesh, time, &
94 !!$ dt, list_mode, level_set_m1, level_set)
95 !!$ IMPLICIT NONE
96 !!$ TYPE(mesh_type) :: vv_mesh
97 !!$ REAL(KIND=8), INTENT(OUT):: time
98 !!$ REAL(KIND=8), INTENT(IN) :: dt
99 !!$ INTEGER, DIMENSION(:), INTENT(IN) :: list_mode
100 !!$ REAL(KIND=8), DIMENSION(:,:,:,:), INTENT(OUT):: level_set, level_set_m1
101 !!$ INTEGER :: mode, i, j, n
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 !!$ END SUBROUTINE init_level_set
114 
115  !===Source in momemtum equation. Always called.
116  FUNCTION source_in_ns_momentum(TYPE, rr, mode, i, time, Re, ty, opt_density, opt_tempn) RESULT(vv)
117  IMPLICIT NONE
118  INTEGER , INTENT(IN) :: type
119  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: rr
120  INTEGER , INTENT(IN) :: mode, i
121  REAL(KIND=8), INTENT(IN) :: time
122  REAL(KIND=8), INTENT(IN) :: re
123  CHARACTER(LEN=2), INTENT(IN) :: ty
124  REAL(KIND=8), DIMENSION(:,:,:), OPTIONAL, INTENT(IN) :: opt_density
125  REAL(KIND=8), DIMENSION(:,:,:), OPTIONAL, INTENT(IN) :: opt_tempn
126  REAL(KIND=8), DIMENSION(SIZE(rr,2)) :: vv, r, z
127  REAL(KIND=8) :: alpha, r0 = 0.5d0
128  CHARACTER(LEN=2) :: np
129 
130  IF (present(opt_density)) CALL error_petsc('density should not be present for test 33')
131 
132  alpha = inputs%gravity_coefficient
133  r = rr(1,:)
134  z = rr(2,:)
135 
136  IF (type==5) THEN
137  vv = alpha*(opt_tempn(:,1,i) - temperature_exact(1,rr,mode,time))
138  ELSE IF (type==6) THEN
139  vv = alpha*(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 -(2*r*(r**4 - 2*r**3*r0 + r0**2 + r**2*(-3 + r0**2))*cos(time)*cos(z) + &
147  (r - r0)*re*cos(time)**2*(14*r**3 - 5*r**2*r0 - 5*r*r0**2 + 2*r0**3 + &
148  (9*r**5 - 21*r**4*r0 + 5*r*r0**2 - 2*r0**3 + r**3*(-14 + 15*r0**2) + r**2*(5*r0 - 3*r0**3))*cos(2*z)) - &
149  2*r**3*(r - r0)**2*re*cos(z)*sin(time))/(2.*r**3*re)
150  ELSE IF (mode==1) THEN
151  vv = vv +(-(r*(r**4 - 2*r*r0 - 2*r**3*r0 + 2*r0**2 + r**2*(-2 + r0**2))*cos(time)*cos(z)) - &
152  (3*r**2 - 4*r*r0 + r0**2)*re*cos(time)**2*(3*r**2 - r0**2 + (2*r**4 - 4*r**3*r0 + r0**2 + &
153  r**2*(-3 + 2*r0**2))*cos(2*z)) + r**3*(r - r0)**2*re*cos(z)*sin(time))/ (r**3*re)
154  ELSE IF (mode==2) THEN
155  vv = vv -((r - r0)*cos(time)**2*(4*r**2 - r*r0 - r0**2 + (3*r**4 - 7*r**3*r0 + r0**2 + &
156  r**2*(-4 + 5*r0**2) + r*(r0 - r0**3))*cos(2*z)))/(2.*r**2)
157  END IF
158  ELSE IF (type==2) THEN
159  IF (mode==1) THEN
160  vv = vv + ((r - r0)**2*cos(time)*(-2*r*cos(z) + re*cos(time)*(r**4 - r*r0 - 2*r**3*r0 + &
161  r0**2 + r**2*(-3 + r0**2) + (3*r**2 + r*r0 - r0**2)*cos(2*z))))/(r**3*re)
162  ELSE IF (mode==2) THEN
163  vv = vv + ((r - r0)**2*cos(time)**2*(r**4 - r*r0 - 2*r**3*r0 + r0**2 + r**2*(-3 + r0**2) + &
164  (3*r**2 + r*r0 - r0**2)*cos(2*z)))/(2.*r**3)
165  END IF
166  ELSE IF (type==3) THEN
167  IF (mode==0) THEN
168  vv = vv + (-3*r*(r - r0)**3*(3*r - r0)*re*cos(time)**2 + 2*(r**4 - 2*r**3*r0 + r0**2 + &
169  r**2*(-3 + r0**2))*cos(time)*cos(z) - 2*r**2*(r - r0)**2*re*cos(z)*sin(time))/(2.*r**2*re)
170  ELSE IF (mode==1) THEN
171  vv = vv -(-2*r*(r**4 - 2*r*r0 - 2*r**3*r0 + 2*r0**2 + r**2*(-2 + r0**2))*cos(time)*cos(z) + &
172  (r - r0)**3*(3*r - r0)*re*cos(time)**2*(1 + 4*r**2 - cos(2*z)) + &
173  2*r**3*(r - r0)**2*re*cos(z)*sin(time))/(2.*r**3*re)
174  ELSE IF (mode==2) THEN
175  vv = vv + ((r - r0)**3*(3*r - r0)*cos(time)**2*(-1 - r**2 + cos(2*z)))/(2.*r**3)
176  END IF
177  ELSE IF (type==4) THEN
178  IF (mode==1) THEN
179  vv = vv + ((r - r0)**2*cos(time)*(-4*r*cos(z) + re*cos(time)*((-3*r + r0)**2 + &
180  (2*r**4 + 6*r*r0 - 4*r**3*r0 - r0**2 + r**2*(-9 + 2*r0**2))*cos(2*z))))/(2.*r**3*re)
181  ELSE IF (mode==2) THEN
182  vv = vv + ((r - r0)**2*cos(time)**2*(4*r - 2*r0 + (r*(-4 + (r - r0)**2) + 2*r0)*cos(2*z)))/(2.*r**2)
183  END IF
184  ELSE IF (type==5) THEN
185  IF (mode==0) THEN
186  vv = vv + (((3*r**4 - 4*r**3*r0 - r0**2 + r**2*(-3 + r0**2))*cos(time) + &
187  r*(r - r0)**2*(3*r**4 - r*r0 - 6*r**3*r0 + 2*r0**2 + r**2*(-4 + 3*r0**2))*re*cos(time)**2*cos(z) - &
188  r**2*(3*r**2 - 4*r*r0 + r0**2)*re*sin(time))*sin(z))/(r**3*re)
189  ELSE IF (mode==1) THEN
190  vv = vv + ((3*r**3 - 4*r0 - 4*r**2*r0 + r*r0**2)*cos(time)*sin(z) - r*(3*r**2 - 4*r*r0 + r0**2)*re*sin(time)*sin(z) + &
191  ((r - r0)**2*(4*r**4 - 2*r*r0 - 8*r**3*r0 + 3*r0**2 + r**2*(-5 + 4*r0**2))*re*cos(time)**2*sin(2*z))/2.)/(r**2*re)
192  ELSE IF (mode==2) THEN
193  vv = vv + ((r - r0)**2*(r**4 - r*r0 - 2*r**3*r0 + r0**2 + r**2*(-1 + r0**2))*cos(time)**2*sin(2*z))/(2.*r**2)
194  END IF
195  ELSE IF (type==6) THEN
196  IF (mode==1) THEN
197  vv = vv -((((-r**3 + 2*r0 + 2*r**2*r0 - r*r0**2)*cos(time) + 4*r*(r - r0)**3*re*cos(time)**2*cos(z) + &
198  r*(r - r0)**2*re*sin(time))*sin(z))/(r**2*re))
199  ELSE IF (mode==2) THEN
200  vv = vv -(((r - r0)**3*cos(time)**2*sin(2*z))/r)
201  END IF
202  END IF
203 
204  IF (TYPE == 1) then
205  IF (mode == 0) THEN
206  vv = vv - 0.25 * exp(2*z) * r**7 * (215.d0 + 51*r**2) * (r - r0)**4 * cos(time)**4 * sin(z)**2
207  ELSE IF (mode == 1) THEN
208  vv = vv - 5 * exp(2*z) * r**7 * (11.d0 + 3*r**2) * (r - r0)**4 * cos(time)**4 * sin(z)**2
209  ELSE IF (mode == 2) THEN
210  vv = vv + 14 * exp(2*z) * r**7 * (r - r0)**4 * cos(time)**4 * sin(z)**2
211  ELSE IF (mode == 3) THEN
212  vv = vv + exp(2*z) * r**7 * (19.d0 + 3*r**2) * (r - r0)**4 * cos(time)**4 * sin(z)**2
213  ELSE IF (mode == 4) THEN
214  vv = vv + 0.75 * exp(2*z) * r**7 * (5.d0 + r**2) * (r - r0)**4 * cos(time)**4 * sin(z)**2
215  ELSE
216  vv = 0.d0
217  END IF
218  ELSE IF (TYPE == 2) then
219  IF (mode == 1) THEN
220  vv = vv - 3 * exp(2*z) * r**7 * (24.d0 + 5*r**2) * (r - r0)**4 * cos(time)**4 * sin(z)**2
221  ELSE IF (mode == 2) THEN
222  vv = vv - 4 * exp(2*z) * r**7 * (13.d0 + 3*r**2) * (r - r0)**4 * cos(time)**4 * sin(z)**2
223  ELSE IF (mode == 3) THEN
224  vv = vv - exp(2*z) * r**7 * (8.d0 + 3*r**2) * (r - r0)**4 * cos(time)**4 * sin(z)**2
225  ELSE IF (mode == 4) THEN
226  vv = vv + 2 * exp(2*z) * r**7 * (r - r0)**4 * cos(time)**4 * sin(z)**2
227  ELSE
228  vv = 0.d0
229  END IF
230  ELSE IF (TYPE == 3) then
231  IF (mode == 0) THEN
232  vv = vv - exp(2*z) * r**7 * (15.d0 + 2*r**2) * (r - r0)**4 * cos(time)**4 * sin(z)**2
233  ELSE IF (mode == 1) THEN
234  vv = vv - 0.5 * exp(2*z) * r**7 * (48.d0 + 7*r**2) * (r - r0)**4 * cos(time)**4 * sin(z)**2
235  ELSE IF (mode == 2) THEN
236  vv = vv - 2 * exp(2*z) * r**7 * (5.d0 + r**2) * (r - r0)**4 * cos(time)**4 * sin(z)**2
237  ELSE IF (mode == 3) THEN
238  vv = vv - 0.5 * exp(2*z) * r**9 * (r - r0)**4 * cos(time)**4 * sin(z)**2
239  ELSE IF (mode == 4) THEN
240  vv = vv + exp(2*z) * r**7 * (r - r0)**4 * cos(time)**4 * sin(z)**2
241  ELSE
242  vv = 0.d0
243  END IF
244  ELSE IF (TYPE == 4) then
245  IF (mode == 1) THEN
246  vv = vv - 0.5 * exp(2*z) * r**7 * (25.d0 + 2*r**2) * (r - r0)**4 * cos(time)**4 * sin(z)**2
247  ELSE IF (mode == 2) THEN
248  vv = vv - 0.25 * exp(2*z) * r**7 * (61.d0 + 6*r**2) * (r - r0)**4 * cos(time)**4 * sin(z)**2
249  ELSE IF (mode == 3) THEN
250  vv = vv - 0.5 * exp(2*z) * r**7 * (17.d0 + 2*r**2) * (r - r0)**4 * cos(time)**4 * sin(z)**2
251  ELSE IF (mode == 4) THEN
252  vv = vv - 0.125 * exp(2*z) * r**7 * (15.d0 + 2*r**2) * (r - r0)**4 * cos(time)**4 * sin(z)**2
253  ELSE
254  vv = 0.d0
255  END IF
256  ELSE IF (TYPE == 5) then
257  IF (mode == 0) THEN
258  vv = vv - 0.125 * exp(2*z) * r**8 * (215.d0 + 34*r**2) * (r - r0)**4 * cos(time)**4 * sin(z)**2
259  ELSE IF (mode == 1) THEN
260  vv = vv - 2.5 * exp(2*z) * r**8 * (11.d0 + 2*r**2) * (r - r0)**4 * cos(time)**4 * sin(z)**2
261  ELSE IF (mode == 2) THEN
262  vv = vv + 7 * exp(2*z) * r**8 * (r - r0)**4 * cos(time)**4 * sin(z)**2
263  ELSE IF (mode == 3) THEN
264  vv = vv + 0.5 * exp(2*z) * r**8 * (19.d0 + 2*r**2) * (r - r0)**4 * cos(time)**4 * sin(z)**2
265  ELSE IF (mode == 4) THEN
266  vv = vv + 0.125 * exp(2*z) * r**8 * (15.d0 + 2*r**2) * (r - r0)**4 * cos(time)**4 * sin(z)**2
267  ELSE
268  vv = 0.d0
269  END IF
270  ELSE IF (TYPE == 6) then
271  IF (mode == 1) THEN
272  vv = vv - exp(2*z) * r**8 * (36.d0 + 5*r**2) * (r - r0)**4 * cos(time)**4 * sin(z)**2
273  ELSE IF (mode == 2) THEN
274  vv = vv - 2 * exp(2*z) * r**8 * (13.d0 + 2*r**2) * (r - r0)**4 * cos(time)**4 * sin(z)**2
275  ELSE IF (mode == 3) THEN
276  vv = vv - exp(2*z) * r**8 * (4.d0 + r**2) * (r - r0)**4 * cos(time)**4 * sin(z)**2
277  ELSE IF (mode == 4) THEN
278  vv = vv + exp(2*z) * r**8 * (r - r0)**4 * cos(time)**4 * sin(z)**2
279  ELSE
280  vv = 0.d0
281  END IF
282  END IF
283 
284  RETURN
285 
286  !===Dummies variables to avoid warning
287  np=ty
288  !===Dummies variables to avoid warning
289  END FUNCTION source_in_ns_momentum
290 
291  !===Extra source in temperature equation. Always called.
292  FUNCTION source_in_temperature(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
297  REAL(KIND=8), INTENT(IN) :: t
298  REAL(KIND=8), DIMENSION(SIZE(rr,2)) :: vv, r, z, c, lambda
299  INTEGER :: i
300  REAL(KIND=8) :: r0 = 0.5d0
301 
302  r = rr(1,:)
303  z = rr(2,:)
304 
305  DO i=1,SIZE(rr,2)
306  IF (r(i).LE.r0) THEN
307  c(i) = inputs%vol_heat_capacity(1)
308  ELSE
309  c(i) = inputs%vol_heat_capacity(2)
310  END IF
311  END DO
312 
313  DO i=1,SIZE(rr,2)
314  IF (r(i).LE.r0) THEN
315  lambda(i) = inputs%temperature_diffusivity(1)
316  ELSE
317  lambda(i) = inputs%temperature_diffusivity(2)
318  END IF
319  END DO
320 
321  IF (type==1) THEN
322  IF (m==0) THEN
323  vv = - (-(r**4 + 18*r*r0 - 2*r**3*r0 - 4*r0**2 + r**2*(-16.d0 + r0**2))*lambda*cos(t) &
324  + c*r**2*(r - r0)**2*sin(t))*sin(z)
325  ELSE IF (m==1) THEN
326  vv = -(-(r**4 + 16*r*r0 - 2*r**3*r0 - 3*r0**2 + r**2*(-15.d0 + r0**2))*lambda*cos(t) &
327  + c*r**2*(r - r0)**2*sin(t))*sin(z)
328  ELSE
329  vv = 0.d0
330  END IF
331  ELSE
332  vv = 0.d0
333  END IF
334 
335  IF (type==1) THEN
336  IF (m==0) THEN
337  DO i=1,size(rr,2)
338  IF (r(i)>r0) THEN
339  vv(i) = vv(i) + c(i)*(-3*r(i)*(r(i) - r0)**4*cos(t)**2*sin(2*z(i)))/4.d0
340  END IF
341  END DO
342  ELSE IF (m==1) THEN
343  DO i=1,size(rr,2)
344  IF (r(i)>r0) THEN
345  vv(i) = vv(i) - c(i)*r(i)*(r(i) - r0)**4*cos(t)**2*sin(2*z(i))
346  END IF
347  END DO
348  ELSE IF (m==2) THEN
349  DO i=1,size(rr,2)
350  IF (r(i)>r0) THEN
351  vv(i) = vv(i) - c(i)/lambda(i)*(r(i)*(r(i) - r0)**4*cos(t)**2*sin(2*z(i)))/4.d0
352  END IF
353  END DO
354  END IF
355  END IF
356  RETURN
357  END FUNCTION source_in_temperature
358 
359 !!$ FUNCTION source_in_level_set(interface_nb,TYPE, rr, m, t)RESULT(vv)
360 !!$ IMPLICIT NONE
361 !!$ INTEGER , INTENT(IN) :: TYPE
362 !!$ REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: rr
363 !!$ INTEGER , INTENT(IN) :: m, interface_nb
364 !!$ REAL(KIND=8), INTENT(IN) :: t
365 !!$ REAL(KIND=8), DIMENSION(SIZE(rr,2)) :: vv
366 !!$
367 !!$ vv=0.d0
368 !!$ RETURN
369 !!$ END FUNCTION source_in_level_set
370 
371  !===Velocity for boundary conditions in Navier-Stokes.
372  !===Can be used also to initialize velocity in: init_velocity_pressure_temperature
373  FUNCTION vv_exact(TYPE,rr,m,t) RESULT(vv)
374  IMPLICIT NONE
375  INTEGER , INTENT(IN) :: type
376  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: rr
377  INTEGER, INTENT(IN) :: m
378  REAL(KIND=8), INTENT(IN) :: t
379  REAL(KIND=8), DIMENSION(SIZE(rr,2)) :: vv, r, z
380  REAL(KIND=8) :: r0 = 0.5d0
381 
382  r = rr(1,:)
383  z = rr(2,:)
384 
385  IF (type==1) THEN
386  IF ((m==0).OR.(m==1)) THEN
387  vv = -(r-r0)**2*cos(t)*cos(z)
388  ELSE
389  vv = 0.d0
390  END IF
391  ELSE IF (type==3) THEN
392  IF ((m==0).OR.(m==1)) THEN
393  vv = (r-r0)**2*cos(t)*cos(z)
394  ELSE
395  vv = 0.d0
396  END IF
397  ELSE IF (type==5) THEN
398  IF ((m==0).OR.(m==1)) THEN
399  vv = (r-r0)*cos(t)*sin(z)/r * (3*r-r0)
400  ELSE
401  vv = 0.d0
402  END IF
403  ELSE IF (type==6) THEN
404  IF (m==1) THEN
405  vv = (r-r0)*cos(t)*sin(z)/r * (r-r0)
406  ELSE
407  vv = 0.d0
408  END IF
409  ELSE
410  vv = 0.d0
411  END IF
412  RETURN
413  END FUNCTION vv_exact
414 
415 !!$ !===Solid velocity imposed when using penalty technique
416 !!$ !===Defined in Fourier space on mode 0 only.
417 !!$ FUNCTION imposed_velocity_by_penalty(rr,t) RESULT(vv)
418 !!$ IMPLICIT NONE
419 !!$ REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: rr
420 !!$ REAL(KIND=8), INTENT(IN) :: t
421 !!$ REAL(KIND=8), DIMENSION(SIZE(rr,2),6) :: vv
422 !!$
423 !!$ vv=0.d0
424 !!$ RETURN
425 !!$ END FUNCTION imposed_velocity_by_penalty
426 
427  !===Pressure for boundary conditions in Navier-Stokes.
428  !===Can be used also to initialize pressure in the subroutine init_velocity_pressure.
429  !===Use this routine for outflow BCs only.
430  !===CAUTION: Do not enfore BCs on pressure where normal component
431  ! of velocity is prescribed.
432  FUNCTION pp_exact(TYPE,rr,m,t) RESULT (vv)
433  IMPLICIT NONE
434  INTEGER , INTENT(IN) :: type
435  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: rr
436  INTEGER , INTENT(IN) :: m
437  REAL(KIND=8), INTENT(IN) :: t
438  REAL(KIND=8), DIMENSION(SIZE(rr,2)) :: vv
439  REAL(KIND=8) :: r
440  INTEGER :: n
441 
442  vv=0.d0
443  RETURN
444 
445  !===Dummies variables to avoid warning
446  n=type; n=SIZE(rr,1); n=m; r=t
447  !===Dummies variables to avoid warning
448  END FUNCTION pp_exact
449 
450  !===Temperature for boundary conditions in temperature equation.
451  FUNCTION temperature_exact(TYPE,rr,m,t) RESULT (vv)
452  IMPLICIT NONE
453  INTEGER , INTENT(IN) :: type
454  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: rr
455  INTEGER , INTENT(IN) :: m
456  REAL(KIND=8), INTENT(IN) :: t
457  REAL(KIND=8), DIMENSION(SIZE(rr,2)) :: vv, r, z
458  REAL(KIND=8) :: r0=0.5d0
459 
460  r = rr(1,:)
461  z = rr(2,:)
462 
463  IF ((type==1).AND.((m==0).OR.(m==1))) THEN
464  vv = r**2*(r-r0)**2*sin(z)*cos(t)
465  ELSE
466  vv = 0.d0
467  END IF
468  RETURN
469  END FUNCTION temperature_exact
470 
471 !!$ !===Can be used to initialize level set in the subroutine init_level_set.
472 !!$ FUNCTION level_set_exact(interface_nb,TYPE,rr,m,t) RESULT (vv)
473 !!$ IMPLICIT NONE
474 !!$ INTEGER , INTENT(IN) :: TYPE
475 !!$ REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: rr
476 !!$ INTEGER , INTENT(IN) :: m, interface_nb
477 !!$ REAL(KIND=8), INTENT(IN) :: t
478 !!$ REAL(KIND=8), DIMENSION(SIZE(rr,2)) :: vv
479 !!$
480 !!$ vv = 0.d0
481 !!$ RETURN
482 !!$ END FUNCTION level_set_exact
483 
484 !!$ !===Penalty coefficient (if needed)
485 !!$ !===This coefficient is equal to zero in subdomain
486 !!$ !===where penalty is applied
487 !!$ FUNCTION penal_in_real_space(mesh,rr_gauss,angles,nb_angles,nb,ne,time) RESULT(vv)
488 !!$ IMPLICIT NONE
489 !!$ TYPE(mesh_type) :: mesh
490 !!$ REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: rr_gauss
491 !!$ REAL(KIND=8), DIMENSION(:), INTENT(IN) :: angles
492 !!$ INTEGER, INTENT(IN) :: nb_angles
493 !!$ INTEGER, INTENT(IN) :: nb, ne
494 !!$ REAL(KIND=8), INTENT(IN) :: time
495 !!$ REAL(KIND=8), DIMENSION(nb_angles,ne-nb+1) :: vv
496 !!$
497 !!$ vv = 1.d0
498 !!$ RETURN
499 !!$ END FUNCTION penal_in_real_space
500 
501  !===Extension of the velocity field in the solid.
502  !===Used when temperature or Maxwell equations are solved.
503  !===It extends the velocity field on the Navier-Stokes domain to a
504  !===velocity field on the temperature and the Maxwell domain.
505  !===It is also used if problem type=mxw and restart velocity
506  !===is set to true in data (type problem denoted mxx in the code).
507  FUNCTION extension_velocity(TYPE, H_mesh, mode, t, n_start) RESULT(vv)
508  IMPLICIT NONE
509  TYPE(mesh_type), INTENT(IN) :: h_mesh
510  INTEGER , INTENT(IN) :: type, n_start
511  INTEGER, INTENT(IN) :: mode
512  REAL(KIND=8), INTENT(IN) :: t
513  REAL(KIND=8), DIMENSION(H_Mesh%np) :: vv
514  REAL(KIND=8) :: r
515  INTEGER :: n
516 
517  vv = 0.d0
518  RETURN
519 
520  !===Dummies variables to avoid warning
521  n=h_mesh%np; r=t; n=type; n=mode; n=n_start
522  !===Dummies variables to avoid warning
523  END FUNCTION extension_velocity
524 
525  !===============================================================================
526  ! Boundary conditions for Maxwell
527  !===============================================================================
528 !!$ !===Velocity used in the induction equation.
529 !!$ !===Used only if problem type is mxw and restart velocity is false
530 !!$ FUNCTION Vexact(m, H_mesh) RESULT(vv) !Set uniquement a l'induction
531 !!$ IMPLICIT NONE
532 !!$ TYPE(mesh_type), INTENT(IN) :: H_mesh
533 !!$ INTEGER, INTENT(IN) :: m
534 !!$ REAL(KIND=8), DIMENSION(H_mesh%np,6) :: vv
535 !!$
536 !!$ vv = 0.d0
537 !!$ RETURN
538 !!$ END FUNCTION Vexact
539 
540 !!$ FUNCTION H_B_quasi_static(char_h_b, rr, m) RESULT(vv)
541 !!$ IMPLICIT NONE
542 !!$ CHARACTER(LEN=1), INTENT(IN) :: char_h_b
543 !!$ REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: rr
544 !!$ INTEGER, INTENT(IN) :: m
545 !!$ REAL(KIND=8), DIMENSION(SIZE(rr,2),6) :: vv
546 !!$
547 !!$ IF (inputs%if_quasi_static_approx) THEN
548 !!$ vv = 0.d0
549 !!$ ELSE
550 !!$ CALL error_petsc('H_B_quasi_static should not be called')
551 !!$ END IF
552 !!$ RETURN
553 !!$ END FUNCTION H_B_quasi_static
554 
555  !===Magnetic field for boundary conditions in the Maxwell equations.
556  FUNCTION hexact(H_mesh, TYPE, rr, m, mu_H_field, t) RESULT(vv)
557  IMPLICIT NONE
558  TYPE(mesh_type), INTENT(IN) :: h_mesh
559  INTEGER , INTENT(IN) :: type
560  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: rr
561  INTEGER , INTENT(IN) :: m
562  REAL(KIND=8), INTENT(IN) :: t
563  REAL(KIND=8), DIMENSION(:), INTENT(IN) :: mu_h_field
564  REAL(KIND=8), DIMENSION(SIZE(rr,2)) :: vv, r, z
565  INTEGER :: n
566 
567  r = rr(1,:)
568  z = rr(2,:)
569 
570  IF (TYPE == 1) then
571  IF (m == 0) THEN
572  vv = - exp(z) * r**3 * cos(t)
573  ELSE
574  vv = 0.d0
575  END IF
576  ELSE IF (TYPE == 2) then
577  IF (m == 1) THEN
578  vv = - exp(z) * r**3 * cos(t)
579  ELSE
580  vv = 0.d0
581  END IF
582  ELSE IF (TYPE == 3) then
583  IF (m == 0) THEN
584  vv = exp(z) * r**3 * cos(t)
585  ELSE
586  vv = 0.d0
587  END IF
588  ELSE IF (TYPE == 4) then
589  IF (m == 1) THEN
590  vv = exp(z) * r**3 * cos(t)
591  ELSE
592  vv = 0.d0
593  END IF
594  ELSE IF (TYPE == 5) then
595  IF (m == 0) THEN
596  vv = 4 * exp(z) * r**2 * cos(t)
597  ELSE IF (m == 1) THEN
598  vv = - exp(z) * r**2 * cos(t)
599  ELSE
600  vv = 0.d0
601  END IF
602  ELSE
603  IF (m == 1) THEN
604  vv = 4 * exp(z) * r**2 * cos(t)
605  ELSE
606  vv = 0.d0
607  END IF
608  END IF
609  RETURN
610 
611  !===Dummies variables to avoid warning
612  n=h_mesh%np; n=SIZE(mu_h_field);
613  !===Dummies variables to avoid warning
614  END FUNCTION hexact
615 
616  !===Scalar potential for boundary conditions in the Maxwell equations.
617  FUNCTION phiexact(TYPE, rr, m, mu_phi,t) RESULT(vv)
618  IMPLICIT NONE
619  INTEGER , INTENT(IN) :: type
620  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: rr
621  INTEGER , INTENT(IN) :: m
622  REAL(KIND=8), INTENT(IN) :: mu_phi, t
623  REAL(KIND=8), DIMENSION(SIZE(rr,2)) :: vv
624  REAL(KIND=8) :: r
625  INTEGER :: n
626 
627  vv=0.d0
628  RETURN
629 
630  !===Dummies variables to avoid warning
631  n=type; n=SIZE(rr,1); n=m; r=mu_phi; r=t
632  !===Dummies variables to avoid warning
633  END FUNCTION phiexact
634 
635  !===Current in Ohm's law. Curl(H) = sigma(E + uxB) + current
636  FUNCTION jexact_gauss(TYPE, rr, m, mu_phi, sigma, mu_H, t, mesh_id, opt_B_ext) RESULT(vv)
637  IMPLICIT NONE
638  INTEGER , INTENT(IN) :: type
639  REAL(KIND=8), DIMENSION(:), INTENT(IN) :: rr
640  INTEGER , INTENT(IN) :: m
641  REAL(KIND=8), INTENT(IN) :: mu_phi, sigma, mu_h, t
642  INTEGER , INTENT(IN) :: mesh_id
643  REAL(KIND=8), DIMENSION(6), OPTIONAL,INTENT(IN) :: opt_b_ext
644  REAL(KIND=8) :: vv
645  REAL(KIND=8) :: r, z
646  INTEGER :: n
647 
648  r = rr(1)
649  z = rr(2)
650 
651  IF (TYPE == 1) then
652  IF (m == 0) THEN
653  vv = - exp(z) * r**3 * cos(t)
654  ELSE IF (m == 1) THEN
655  vv = 4 * exp(z) * r * cos(t)
656  ELSE
657  vv = 0.d0
658  END IF
659  ELSE IF (TYPE == 2) then
660  IF (m == 1) THEN
661  vv = - exp(z) * r * (-1.d0 + r**2) * cos(t)
662  ELSE
663  vv = 0.d0
664  END IF
665  ELSE IF (TYPE == 3) then
666  IF (m == 0) THEN
667  vv = - exp(z) * r * (8.d0 + r**2) * cos(t)
668  ELSE IF (m == 1) THEN
669  vv = 2 * exp(z) * r * cos(t)
670  ELSE
671  vv = 0.d0
672  END IF
673  ELSE IF (TYPE == 4) then
674  IF (m == 1) THEN
675  vv = - exp(z) * r * (8.d0 + r**2) * cos(t)
676  ELSE
677  vv = 0.d0
678  END IF
679  ELSE IF (TYPE == 5) then
680  IF (m == 0) THEN
681  vv = 4 * exp(z) *r**2 * cos(t)
682  ELSE IF (m == 1) THEN
683  vv = exp(z) * r**2 * cos(t)
684  ELSE
685  vv = 0.d0
686  END IF
687  ELSE
688  IF (m == 1) THEN
689  vv = 4 * exp(z) * r**2 * cos(t)
690  ELSE
691  vv = 0.d0
692  END IF
693  END IF
694  RETURN
695 
696  !===Dummies variables to avoid warning
697  r=mu_phi; r=sigma; r=mu_h; n=mesh_id
698  IF (present(opt_b_ext)) r=opt_b_ext(1)
699  !===Dummies variables to avoid warning
700  END FUNCTION jexact_gauss
701 
702  !===Electric field for Neumann BC (cf. doc)
703  FUNCTION eexact_gauss(TYPE, rr, m, mu_phi, sigma, mu_H, t) RESULT(vv)
704  IMPLICIT NONE
705  INTEGER, INTENT(IN) :: type
706  REAL(KIND=8), DIMENSION(:), INTENT(IN) :: rr
707  INTEGER, INTENT(IN) :: m
708  REAL(KIND=8), INTENT(IN) :: mu_phi, sigma, mu_h, t
709  REAL(KIND=8) :: vv
710  REAL(KIND=8) :: r
711  INTEGER :: n
712 
713  vv = 0.d0
714  RETURN
715 
716  !===Dummies variables to avoid warning
717  r=rr(1); r=mu_phi; r=sigma; r=mu_h; r=t; n=type; n=m
718  !===Dummies variables to avoid warning
719  END FUNCTION eexact_gauss
720 
721  !===Initialization of magnetic field and scalar potential (if present)
722  SUBROUTINE init_maxwell(H_mesh, phi_mesh, time, dt, mu_H_field, mu_phi, &
723  list_mode, hn1, hn, phin1, phin)
724  IMPLICIT NONE
725  TYPE(mesh_type) :: h_mesh, phi_mesh
726  REAL(KIND=8), INTENT(OUT):: time
727  REAL(KIND=8), INTENT(IN) :: dt
728  REAL(KIND=8), DIMENSION(:), INTENT(IN) :: mu_h_field
729  REAL(KIND=8), INTENT(IN) :: mu_phi
730  INTEGER, DIMENSION(:), INTENT(IN) :: list_mode
731  REAL(KIND=8), DIMENSION(:,:,:), INTENT(OUT):: hn, hn1
732  REAL(KIND=8), DIMENSION(:,:,:), INTENT(OUT):: phin, phin1
733  INTEGER :: i, k
734 
735  time = -dt
736  DO k=1,6
737  DO i=1, SIZE(list_mode)
738  hn1(:,k,i) = hexact(h_mesh,k, h_mesh%rr, list_mode(i), mu_h_field, time)
739  IF (inputs%nb_dom_phi>0) THEN
740  IF (k<3) THEN
741  phin1(:,k,i) = phiexact(k, phi_mesh%rr, list_mode(i) , mu_phi, time)
742  ENDIF
743  END IF
744  ENDDO
745  ENDDO
746 
747  time = time + dt
748  DO k=1,6
749  DO i=1, SIZE(list_mode)
750  hn(:,k,i) = hexact(h_mesh,k, h_mesh%rr, list_mode(i), mu_h_field, time)
751  IF (inputs%nb_dom_phi>0) THEN
752  IF (k<3) THEN
753  phin(:,k,i) = phiexact(k, phi_mesh%rr, list_mode(i), mu_phi, time)
754  ENDIF
755  END IF
756  ENDDO
757  ENDDO
758  RETURN
759  END SUBROUTINE init_maxwell
760 
761 !!$ !===Analytical permeability (if needed)
762 !!$ !===This function is not needed unless the flag
763 !!$ !=== ===Use FEM Interpolation for magnetic permeability (true/false)
764 !!$ !===is activated and set to .FALSE. in the data data file. Default is .TRUE.
765 !!$ FUNCTION mu_bar_in_fourier_space(H_mesh,nb,ne,pts,pts_ids) RESULT(vv)
766 !!$ IMPLICIT NONE
767 !!$ TYPE(mesh_type), INTENT(IN) :: H_mesh
768 !!$ REAL(KIND=8), DIMENSION(ne-nb+1) :: vv
769 !!$ INTEGER, INTENT(IN) :: nb, ne
770 !!$ REAL(KIND=8),DIMENSION(2,ne-nb+1),OPTIONAL :: pts
771 !!$ INTEGER, DIMENSION(ne-nb+1), OPTIONAL :: pts_ids
772 !!$
773 !!$ vv = 1.d0
774 !!$ RETURN
775 !!$ END FUNCTION mu_bar_in_fourier_space
776 
777 !!$ !===Analytical mu_in_fourier_space (if needed)
778 !!$ !===This function is not needed unless the flag
779 !!$ !=== ===Use FEM Interpolation for magnetic permeability (true/false)
780 !!$ !===is activated and set to .FALSE. in the data data file. Default is .TRUE.
781 !!$ FUNCTION grad_mu_bar_in_fourier_space(pt,pt_id) RESULT(vv)
782 !!$ IMPLICIT NONE
783 !!$ REAL(KIND=8),DIMENSION(2), INTENT(in):: pt
784 !!$ INTEGER,DIMENSION(1), INTENT(in) :: pt_id
785 !!$ REAL(KIND=8),DIMENSION(2) :: vv
786 !!$
787 !!$ vv=0.d0
788 !!$ RETURN
789 !!$ END FUNCTION grad_mu_bar_in_fourier_space
790 
791 !!$ !===Analytical permeability, mu in real space (if needed)
792 !!$ FUNCTION mu_in_real_space(H_mesh,angles,nb_angles,nb,ne,time) RESULT(vv)
793 !!$ IMPLICIT NONE
794 !!$ TYPE(mesh_type), INTENT(IN) :: H_mesh
795 !!$ REAL(KIND=8), DIMENSION(:), INTENT(IN) :: angles
796 !!$ INTEGER, INTENT(IN) :: nb_angles
797 !!$ INTEGER, INTENT(IN) :: nb, ne
798 !!$ REAL(KIND=8), INTENT(IN) :: time
799 !!$ REAL(KIND=8), DIMENSION(nb_angles,ne-nb+1) :: vv
800 !!$
801 !!$ vv = 1.d0
802 !!$ RETURN
803 !!$ END FUNCTION mu_in_real_space
804 
805  !===Coefficient of Kelvin force: F = coeff(T) * grad(H**2/2)
806  FUNCTION kelvin_force_coeff(temp) RESULT(vv)
807  IMPLICIT NONE
808  REAL(KIND=8) :: temp
809  REAL(KIND=8) :: vv
810 
811  vv = temp**2
812  RETURN
813  END FUNCTION kelvin_force_coeff
814 
815 END MODULE boundary_test_33
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)
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, public kelvin_force_coeff(temp)
real(kind=8) function, public jexact_gauss(TYPE, rr, m, mu_phi, sigma, mu_H, t, mesh_id, opt_B_ext)
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 phiexact(TYPE, rr, m, mu_phi, t)
real(kind=8) function, dimension(size(rr, 2)), public source_in_temperature(TYPE, rr, m, t)
real(kind=8) function, public eexact_gauss(TYPE, rr, m, mu_phi, sigma, mu_H, 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
real(kind=8) function, dimension(size(rr, 2)), public hexact(H_mesh, TYPE, rr, m, mu_H_field, t)
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