41 un_m1, un, pn_m1, pn, phin_m1, phin)
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
50 REAL(KIND=8),
DIMENSION(mesh_c%np) :: pn_m2
53 DO i= 1,
SIZE(list_mode)
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)
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)
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
82 DO i= 1,
SIZE(list_mode)
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
130 IF (present(opt_density)) CALL
error_petsc(
'density should not be present for test 31')
132 gamma = inputs%gravity_coefficient
138 ELSE IF (type==6)
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)
161 ELSE IF (type==2)
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)
169 ELSE IF (type==3)
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)
182 ELSE IF (type==4)
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
189 ELSE IF (type==5)
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
205 ELSE IF (type==6)
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
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)
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
238 REAL(KIND=8) :: r0 = 0.5d0
239 REAL(KIND=8) :: pi = 3.1415926535897932d0
246 kappa(i) = inputs%temperature_diffusivity(1)
248 kappa(i) = inputs%temperature_diffusivity(2)
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)
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)
270 vv(i) = vv(i) + (-3*pi*r(i)*(r(i) - r0)**4*cos(
t)**2*sin(4*pi*
z(i)))/2.d0
276 vv(i) = vv(i) - 2*pi*r(i)*(r(i) - r0)**4*cos(
t)**2*sin(4*pi*
z(i))
282 vv(i) = vv(i) - 0.5*pi*(r(i)*(r(i) - r0)**4*cos(
t)**2*sin(4*pi*
z(i)))
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
320 IF ((m==0).OR.(m==1))
THEN
321 vv = -2*pi*(r-r0)**2*cos(
t)*cos(2*pi*
z)
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)
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)
337 ELSE IF (type==6)
THEN
339 vv = (r-r0)*cos(
t)*sin(2*pi*
z)/r * (r-r0)
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
378 IF ((type==1).AND.(m==1))
THEN
379 vv = r**3*sin(2*pi*
z)*cos(
t)
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
400 IF ((type==1).AND.((m==0).OR.(m==1)))
THEN
401 vv = r**2*(r-r0)**2*sin(2*pi*
z)*cos(
t)
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
459 n=h_mesh%np; r=
t; n=type; n=mode; n=n_start
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
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)
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