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)
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)) :: lap, temp,gpress
129 REAL(KIND=8),
DIMENSION(SIZE(rr,2)) :: r,
z
130 REAL(KIND=8),
DIMENSION(SIZE(rr,2),6) :: rot, v, vt
133 CHARACTER(LEN=2) :: np
135 IF (present(opt_density)) CALL
error_petsc(
'density should not be present for test 2')
136 IF (present(opt_tempn)) CALL
error_petsc(
'temperature should not be present for test 2')
144 IF (m==0 .OR. m==2)
THEN
146 rot(:,1) = r*(4.d0*cos(k*
z)+1)
149 rot(:,4) = r*(k**2*r**2*cos(k*
z)-8.d0*cos(k*
z)-2.d0)
150 rot(:,5) = r*(-8.d0-k*r*sin(k*
z))
157 vt(:,2) = -r**2*(1.d0-k*r*sin(k*
z))
159 vt(:,6) = r**2*(4.d0*cos(k*
z)+1.d0)
163 IF (mod(type,2)==0)
THEN
165 ELSEIF (
TYPE == 1) then
166 vv = -0.5d0*(v(:,3)*rot(:,5)-v(:,5)*rot(:,3) &
167 + v(:,4)*rot(:,6)-v(:,6)*rot(:,4))
168 ELSEIF (
TYPE == 3) then
169 vv = -0.5d0*(v(:,5)*rot(:,1)-v(:,1)*rot(:,5) &
170 + v(:,6)*rot(:,2)-v(:,2)*rot(:,6))
171 ELSEIF (
TYPE == 5) then
172 vv = -0.5d0*(v(:,1)*rot(:,3)-v(:,3)*rot(:,1) &
173 + v(:,2)*rot(:,4)-v(:,4)*rot(:,2))
180 vv = -0.5d0*(v(:,3)*rot(:,5)-v(:,5)*rot(:,3) &
181 - (v(:,4)*rot(:,6)-v(:,6)*rot(:,4)))
182 ELSE IF (
TYPE == 2) then
183 vv = -0.5d0*(v(:,3)*rot(:,6)-v(:,5)*rot(:,4) &
184 + (v(:,4)*rot(:,5)-v(:,6)*rot(:,3)))
185 ELSEIF (
TYPE == 3) then
186 vv = -0.5d0*(v(:,5)*rot(:,1)-v(:,1)*rot(:,5) &
187 - (v(:,6)*rot(:,2)-v(:,2)*rot(:,6)))
188 ELSEIF (
TYPE == 4) then
189 vv = -0.5d0*(v(:,5)*rot(:,2)-v(:,1)*rot(:,6) &
190 + (v(:,6)*rot(:,1)-v(:,2)*rot(:,5)))
191 ELSEIF (
TYPE == 5) then
192 vv = -0.5d0*(v(:,1)*rot(:,3)-v(:,3)*rot(:,1) &
193 - (v(:,2)*rot(:,4)-v(:,4)*rot(:,2)))
194 ELSEIF (
TYPE == 6) then
195 vv = -0.5d0*(v(:,1)*rot(:,4)-v(:,3)*rot(:,2) &
196 + (v(:,2)*rot(:,3)-v(:,4)*rot(:,1)))
211 vt(:,2) = -r**2*(1.d0-k*r*sin(k*
z))
213 vt(:,6) = r**2*(4.d0*cos(k*
z)+1.d0)
218 gpress(:) = 2*r*cos(k*
z)
219 ELSEIF (
TYPE == 2) then
220 lap(:) = (-8+7*k*r*sin(k*
z)-k**3*r**3*sin(k*
z))
223 ELSEIF (
TYPE == 3) then
224 lap(:) = (-8+2*k*r*sin(k*
z))
227 ELSEIF (
TYPE == 4) then
230 gpress(:) = -r*cos(k*
z)
231 ELSEIF (
TYPE == 5) then
234 gpress(:) = -r**2*k*sin(k*
z)
235 ELSEIF (
TYPE == 6) then
236 lap(:) = (3+12*cos(k*
z)-4*k**2*r**2*cos(k*
z))
241 vv(:) = -sin(
t)*temp + cos(
t)*(-lap/re) + gpress*cos(
t)
281 INTEGER ,
INTENT(IN) :: type
282 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: rr
283 INTEGER,
INTENT(IN) :: m
284 REAL(KIND=8),
INTENT(IN) ::
t
285 REAL(KIND=8),
DIMENSION(SIZE(rr,2)) :: vv
287 REAL(KIND=8),
DIMENSION(SIZE(rr,2)) :: r,
z
300 ELSEIF (
TYPE == 2 .AND. m /= 0)
THEN
301 vv(:) = -r**2*(1-k*r*sin(k*
z))
302 ELSEIF (
TYPE == 3) then
304 ELSEIF (
TYPE == 4 .AND. m /= 0)
THEN
306 ELSEIF (
TYPE == 5) then
308 ELSEIF (
TYPE == 6 .AND. m /= 0)
THEN
309 vv(:) = r**2*(4*cos(k*
z)+1)
312 vv(:) = vv(:) * cos(
t)
335 INTEGER ,
INTENT(IN) :: type
336 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: rr
337 INTEGER ,
INTENT(IN) :: m
338 REAL(KIND=8),
INTENT(IN) ::
t
339 REAL(KIND=8),
DIMENSION(SIZE(rr,2)) :: vv
349 vv(:)= rr(1,:)**2*cos(k*rr(2,:))*cos(
t)
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(size(rr, 2)), public vv_exact(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)
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