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