33 REAL(KIND=8),
PRIVATE :: test11_Ri=7/13d0, test11_Ro=20/13d0, test11_Rossby=1.d-4
34 REAL(KIND=8) :: pi=ACOS(-1.d0)
43 un_m1, un, pn_m1, pn, phin_m1, phin)
46 REAL(KIND=8),
INTENT(OUT):: time
47 REAL(KIND=8),
INTENT(IN) :: dt
48 INTEGER,
DIMENSION(:),
INTENT(IN) :: list_mode
49 REAL(KIND=8),
DIMENSION(:,:,:),
INTENT(OUT):: un_m1, un
50 REAL(KIND=8),
DIMENSION(:,:,:),
INTENT(OUT):: pn_m1, pn, phin_m1, phin
52 REAL(KIND=8),
DIMENSION(mesh_c%np) :: pn_m2
55 DO i= 1,
SIZE(list_mode)
59 un_m1(:,j,i) =
vv_exact(j,mesh_f%rr,mode,time-dt)
60 un(:,j,i) =
vv_exact(j,mesh_f%rr,mode,time)
64 pn_m2(:) =
pp_exact(j,mesh_c%rr,mode,time-2*dt)
65 pn_m1(:,j,i) =
pp_exact(j,mesh_c%rr,mode,time-dt)
66 pn(:,j,i) =
pp_exact(j,mesh_c%rr,mode,time)
67 phin_m1(:,j,i) = pn_m1(:,j,i) - pn_m2(:)
68 phin(:,j,i) = pn(:,j,i) - pn_m1(:,j,i)
77 REAL(KIND=8),
INTENT(OUT):: time
78 REAL(KIND=8),
INTENT(IN) :: dt
79 INTEGER,
DIMENSION(:),
INTENT(IN) :: list_mode
80 REAL(KIND=8),
DIMENSION(:,:,:),
INTENT(OUT):: tempn_m1, tempn
84 DO i= 1,
SIZE(list_mode)
121 INTEGER ,
INTENT(IN) :: type
122 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: rr
123 INTEGER ,
INTENT(IN) :: mode, i
124 REAL(KIND=8),
INTENT(IN) :: time
125 REAL(KIND=8),
INTENT(IN) :: re
126 CHARACTER(LEN=2),
INTENT(IN) :: ty
127 REAL(KIND=8),
DIMENSION(:,:,:),
OPTIONAL,
INTENT(IN) :: opt_density
128 REAL(KIND=8),
DIMENSION(:,:,:),
OPTIONAL,
INTENT(IN) :: opt_tempn
129 REAL(KIND=8),
DIMENSION(SIZE(rr,2)) :: vv
130 REAL(KIND=8),
DIMENSION(SIZE(rr,2)) :: r,
z
132 CHARACTER(LEN=2) :: np
134 IF (present(opt_density)) CALL
error_petsc(
'density should not be present for test 11')
140 vv = inputs%gravity_coefficient*opt_tempn(:,1,i)*r
141 ELSE IF (type==2)
THEN
142 vv = inputs%gravity_coefficient*opt_tempn(:,2,i)*r
143 ELSE IF (type==5)
THEN
144 vv = inputs%gravity_coefficient*opt_tempn(:,1,i)*
z
145 ELSE IF (type==6)
THEN
146 vv = inputs%gravity_coefficient*opt_tempn(:,2,i)*
z
153 n=mode; r=time; r=re; np=ty
160 INTEGER ,
INTENT(IN) :: type
161 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: rr
162 INTEGER ,
INTENT(IN) :: m
163 REAL(KIND=8),
INTENT(IN) ::
t
164 REAL(KIND=8),
DIMENSION(SIZE(rr,2)) :: vv
172 n=type; r=rr(1,1); n=m; r=
t
193 INTEGER ,
INTENT(IN) :: type
194 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: rr
195 INTEGER,
INTENT(IN) :: m
196 REAL(KIND=8),
INTENT(IN) ::
t
197 REAL(KIND=8),
DIMENSION(SIZE(rr,2)) :: vv
205 n=type; r=rr(1,1); n=m; r=
t
228 INTEGER ,
INTENT(IN) :: type
229 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: rr
230 INTEGER ,
INTENT(IN) :: m
231 REAL(KIND=8),
INTENT(IN) ::
t
232 REAL(KIND=8),
DIMENSION(SIZE(rr,2)) :: vv
240 n=type; r=rr(1,1); n=m; r=
t
247 INTEGER ,
INTENT(IN) :: type
248 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: rr
249 INTEGER ,
INTENT(IN) :: m
250 REAL(KIND=8),
INTENT(IN) ::
t
251 REAL(KIND=8),
DIMENSION(SIZE(rr,2)) :: vv
252 REAL(KIND=8),
DIMENSION(SIZE(rr,2)) :: r,
z, rho, a
254 REAL(KIND=8) ::
theta, x
259 IF (m==0 .OR. m==4)
THEN
261 rho(n)=sqrt(rr(1,n)**2+rr(2,n)**2)
262 theta=atan2(rr(1,n),rr(2,n))
263 x=2*rho(n) - test11_ri - test11_ro
264 a(n)=(21/sqrt(17920*pi))*(1-3*x**2+3*x**4-x**6)*sin(
theta)**4
266 IF (m==0 .AND. type==1)
THEN
267 vv=( test11_ri*test11_ro/rho)- test11_ri
268 ELSE IF (m==4 .AND. type==1)
THEN
325 INTEGER ,
INTENT(IN) :: type, n_start
326 INTEGER,
INTENT(IN) :: mode
327 REAL(KIND=8),
INTENT(IN) ::
t
328 REAL(KIND=8),
DIMENSION(H_Mesh%np) :: vv
336 n=h_mesh%np; r=
t; n=type; n=mode; n=n_start
369 FUNCTION hexact(H_mesh,TYPE, rr, m, mu_H_field, t) RESULT(vv)
372 INTEGER ,
INTENT(IN) :: type
373 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: rr
374 INTEGER ,
INTENT(IN) :: m
375 REAL(KIND=8),
INTENT(IN) ::
t
376 REAL(KIND=8),
DIMENSION(:),
INTENT(IN) :: mu_h_field
377 REAL(KIND=8),
DIMENSION(SIZE(rr,2)) :: vv
385 n=h_mesh%np; n=type; r=rr(1,1); n=m; r=
t; r=mu_h_field(1)
390 FUNCTION phiexact(TYPE, rr, m, mu_phi,t) RESULT(vv)
392 INTEGER ,
INTENT(IN) :: type
393 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: rr
394 INTEGER ,
INTENT(IN) :: m
395 REAL(KIND=8),
INTENT(IN) :: mu_phi,
t
396 REAL(KIND=8),
DIMENSION(SIZE(rr,2)) :: vv
401 CALL
error_petsc(
'Phiexact: should not be called for this test')
405 n=type; r=rr(1,1); n=m; r=mu_phi; r=
t
410 FUNCTION jexact_gauss(TYPE, rr, m, mu_phi, sigma, mu_H, t, mesh_id, opt_B_ext) RESULT(vv)
412 INTEGER ,
INTENT(IN) :: type
413 REAL(KIND=8),
DIMENSION(:),
INTENT(IN) :: rr
414 INTEGER ,
INTENT(IN) :: m
415 REAL(KIND=8),
INTENT(IN) :: mu_phi, sigma, mu_h,
t
416 INTEGER ,
INTENT(IN) :: mesh_id
417 REAL(KIND=8),
DIMENSION(6),
OPTIONAL,
INTENT(IN) :: opt_b_ext
426 r=rr(1); r=mu_phi; r=sigma; r=mu_h; r=
t; n=type; n=m; n=mesh_id
427 IF (present(opt_b_ext)) r=opt_b_ext(1)
434 INTEGER,
INTENT(IN) :: type
435 REAL(KIND=8),
DIMENSION(:),
INTENT(IN) :: rr
436 INTEGER,
INTENT(IN) :: m
437 REAL(KIND=8),
INTENT(IN) :: mu_phi, sigma, mu_h,
t
446 r=rr(1); r=mu_phi; r=sigma; r=mu_h; r=
t; n=type; n=m
451 SUBROUTINE init_maxwell(H_mesh, phi_mesh, time, dt, mu_H_field, mu_phi, &
452 list_mode, hn1, hn, phin1, phin)
455 REAL(KIND=8),
INTENT(OUT):: time
456 REAL(KIND=8),
INTENT(IN) :: dt
457 REAL(KIND=8),
DIMENSION(:),
INTENT(IN) :: mu_h_field
458 REAL(KIND=8),
INTENT(IN) :: mu_phi
459 INTEGER,
DIMENSION(:),
INTENT(IN) :: list_mode
460 REAL(KIND=8),
DIMENSION(:,:,:),
INTENT(OUT):: hn, hn1
461 REAL(KIND=8),
DIMENSION(:,:,:),
INTENT(OUT):: phin, phin1
463 REAL(KIND=8) :: rho, normalization
464 REAL(KIND=8),
DIMENSION(SIZE(H_mesh%rr,2)) :: brho, btheta, bphi,
theta
467 DO i = 1,
SIZE(list_mode)
468 IF (list_mode(i) /= 0) cycle
469 normalization = 2*sqrt(test11_rossby)
470 DO n=1,
SIZE(h_mesh%rr,2)
471 rho = sqrt(h_mesh%rr(1,n)**2+h_mesh%rr(2,n)**2)
472 theta(n) = atan2(h_mesh%rr(1,n),h_mesh%rr(2,n))
473 brho(n) = cos(
theta(n))*5.d0/(sqrt(2d0)*8d0)*(-48d0*test11_ri*test11_ro &
474 + (4*test11_ro+test11_ri*(4+3*test11_ro))*6*rho &
475 -4*(4+3*(test11_ri+test11_ro))*rho**2+9*rho**3)/rho
476 btheta(n)= sin(
theta(n))*(-15d0/(sqrt(2d0)*4d0))*((rho-test11_ri)*(rho-test11_ro)*(3*rho-4))/rho
477 bphi(n) = sin(2*
theta(n))*(15d0/(sqrt(2d0)*8d0))*sin(pi*(rho-test11_ri))
479 hn(:,1,i) = normalization*(brho*sin(
theta) + btheta*cos(
theta))
480 hn(:,3,i) = normalization*bphi
481 hn(:,5,i) = normalization*(brho*cos(
theta) - btheta*sin(
theta))
490 rho=dt; rho=mu_h_field(1); rho=mu_phi; n=phi_mesh%np;
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)
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 theta
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 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)
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