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