21 INTEGER,
DIMENSION(:),
INTENT(IN) :: jjs
22 REAL(KIND=8),
DIMENSION(:),
INTENT(IN) :: us
23 REAL(KIND=8),
DIMENSION(:),
INTENT(INOUT) :: ff
27 DO n = 1,
SIZE(jjs); i = jjs(n)
34 SUBROUTINE moy(mesh,p,RESULT)
41 REAL(KIND=8),
DIMENSION(:) ,
INTENT(IN) :: p
42 REAL(KIND=8) ,
INTENT(OUT) :: result
44 INTEGER :: m, l , i , ni,k
47 TYPE(mesh_type
),
TARGET :: mesh
48 INTEGER,
DIMENSION(:,:),
POINTER :: jj
49 INTEGER,
POINTER :: me
65 DO ni = 1, n_w; i = jj(ni,m)
66 ray = ray + mesh%rr(1,i)*ww(ni,l)
69 result = result + sum(p(jj(:,m)) * ww(:,l))* ray* rj(l,m)
70 vol = vol + ray* rj(l,m)
81 SUBROUTINE qs_00_ssr(mesh, ff, u0) !sans le rayon dans l'integration
90 REAL(KIND=8),
DIMENSION(:),
INTENT(IN) :: ff
91 REAL(KIND=8),
DIMENSION(:),
INTENT(OUT) :: u0
93 INTEGER :: m, l , i , ni
96 TYPE(mesh_type
),
TARGET :: mesh
97 INTEGER,
DIMENSION(:,:),
POINTER :: jj
98 INTEGER,
POINTER :: me
112 DO ni = 1, n_w; i = jj(ni,m)
113 ray = ray + mesh%rr(1,i)*ww(ni,l)
116 fl = sum(ff(jj(:,m)) * ww(:,l)) * rj(l,m)
119 u0(jj(ni,m)) = u0(jj(ni,m)) + ww(ni,l) * fl
127 SUBROUTINE qs_00(mesh, ff, u0) !avec le rayon dans l'integration
135 REAL(KIND=8),
DIMENSION(:),
INTENT(IN) :: ff
136 REAL(KIND=8),
DIMENSION(:),
INTENT(OUT) :: u0
138 INTEGER :: m, l , i , ni
141 TYPE(mesh_type
),
TARGET :: mesh
142 INTEGER,
DIMENSION(:,:),
POINTER :: jj
143 INTEGER,
POINTER :: me
157 DO ni = 1, n_w; i = jj(ni,m)
158 ray = ray + mesh%rr(1,i)*ww(ni,l)
161 fl = sum(ff(jj(:,m)) * ww(:,l)) * rj(l,m)
164 u0(jj(ni,m)) = u0(jj(ni,m)) + ww(ni,l) * fl *ray
181 REAL(KIND=8),
DIMENSION(:),
INTENT(IN) :: ff
182 REAL(KIND=8),
DIMENSION(:),
INTENT(OUT) :: u0
183 REAL(KIND=8),
DIMENSION(:),
INTENT(IN) :: t1
184 REAL(KIND=8),
DIMENSION(:),
INTENT(IN) :: t2
185 REAL(KIND=8),
INTENT(IN) :: dt
186 INTEGER :: m, l , i , ni
190 TYPE(mesh_type
),
TARGET :: mesh
191 INTEGER,
DIMENSION(:,:),
POINTER :: jj
192 INTEGER,
POINTER :: me
206 DO ni = 1, n_w; i = jj(ni,m)
207 ray = ray + mesh%rr(1,i)*ww(ni,l)
210 fl = sum(ff(jj(:,m)) * ww(:,l)) * rj(l,m) &
211 + sum(4 * t1(jj(:,m)) * ww(:,l) - t2(jj(:,m))*ww(:,l))*rj(l,m) /(2* dt)*ray
214 u0(jj(ni,m)) = u0(jj(ni,m)) + ww(ni,l) * fl
231 REAL(KIND=8),
DIMENSION(:),
INTENT(IN) :: ff
232 REAL(KIND=8),
DIMENSION(:),
INTENT(OUT) :: u0
233 REAL(KIND=8),
DIMENSION(:),
INTENT(IN) :: t1
234 REAL(KIND=8),
DIMENSION(:),
INTENT(IN) :: t2
235 REAL(KIND=8),
INTENT(IN) :: dt
236 INTEGER :: m, l , i , ni
240 TYPE(mesh_type
),
TARGET :: mesh
241 INTEGER,
DIMENSION(:,:),
POINTER :: jj
242 INTEGER,
POINTER :: me
256 DO ni = 1, n_w; i = jj(ni,m)
257 ray = ray + mesh%rr(1,i)*ww(ni,l)
260 fl = sum(ff(jj(:,m)) * ww(:,l)) * rj(l,m) &
261 + sum(4 * t1(jj(:,m)) * ww(:,l) - t2(jj(:,m))*ww(:,l))*rj(l,m) /(2* dt)
264 u0(jj(ni,m)) = u0(jj(ni,m)) + ww(ni,l) * fl *ray
281 REAL(KIND=8),
DIMENSION(:),
INTENT(IN) :: ff
282 REAL(KIND=8),
DIMENSION(:),
INTENT(OUT) :: u0
283 REAL(KIND=8),
DIMENSION(:),
INTENT(IN) :: t1
284 REAL(KIND=8),
INTENT(IN) :: dt
285 INTEGER :: m, l , i , ni
289 TYPE(mesh_type
),
TARGET :: mesh
290 INTEGER,
DIMENSION(:,:),
POINTER :: jj
291 INTEGER,
POINTER :: me
305 DO ni = 1, n_w; i = jj(ni,m)
306 ray = ray + mesh%rr(1,i)*ww(ni,l)
309 fl = sum(ff(jj(:,m)) * ww(:,l)) * rj(l,m) &
310 +sum(t1(jj(:,m)) * ww(:,l)) * rj(l,m) / dt
314 u0(jj(ni,m)) = u0(jj(ni,m)) + ww(ni,l) * ray * fl
330 REAL(KIND=8),
DIMENSION(:),
INTENT(IN) :: ff
331 REAL(KIND=8),
DIMENSION(:),
INTENT(OUT) :: u0
332 REAL(KIND=8),
DIMENSION(:),
INTENT(IN) :: t1
333 REAL(KIND=8),
INTENT(IN) :: dt
334 INTEGER :: m, l , i , ni
338 TYPE(mesh_type
),
TARGET :: mesh
339 INTEGER,
DIMENSION(:,:),
POINTER :: jj
340 INTEGER,
POINTER :: me
354 DO ni = 1, n_w; i = jj(ni,m)
355 ray = ray + mesh%rr(1,i)*ww(ni,l)
358 fl = sum(ff(jj(:,m)) * ww(:,l)) * rj(l,m) &
359 +sum(t1(jj(:,m)) * ww(:,l)) * rj(l,m) / dt *ray
363 u0(jj(ni,m)) = u0(jj(ni,m)) + ww(ni,l) * fl
384 REAL(KIND=8),
DIMENSION(:),
INTENT(IN) :: ff
385 REAL(KIND=8),
DIMENSION(:),
INTENT(OUT) :: u0
386 REAL(KIND=8),
DIMENSION(:),
INTENT(IN) :: t1
387 REAL(KIND=8),
DIMENSION(:),
INTENT(IN) :: t2
388 REAL(KIND=8),
INTENT(IN) :: dt
389 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: gg
390 REAL(KIND=8),
INTENT(IN) :: eps
391 INTEGER ,
INTENT(IN) :: mode
392 INTEGER :: m, l , i , ni
393 REAL(KIND=8) :: fs , ft , fexp
396 TYPE(mesh_type
),
TARGET :: mesh
397 INTEGER,
DIMENSION(:,:),
POINTER :: jj
398 INTEGER,
POINTER :: me
412 DO ni = 1, n_w; i = jj(ni,m)
413 ray = ray + mesh%rr(1,i)*ww(ni,l)
416 fs = sum(ff(jj(:,m)) * ww(:,l)) * rj(l,m)
417 ft = sum(4 * t1(jj(:,m)) * ww(:,l) - t2(jj(:,m))*ww(:,l))*rj(l,m) /(2* dt)
418 fexp = eps * mode * sum(gg(3,jj(:,m)) * t1(jj(:,m)) * ww(:,l))*rj(l,m) / ray
422 u0(jj(ni,m)) = u0(jj(ni,m)) + ww(ni,l) * ray *( fs + ft + fexp)
442 REAL(KIND=8),
DIMENSION(:),
INTENT(IN) :: ff
443 REAL(KIND=8),
DIMENSION(:),
INTENT(OUT) :: u0
444 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: t1
445 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: t2
446 REAL(KIND=8),
INTENT(IN) :: dt
447 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: gg
448 REAL(KIND=8),
INTENT(IN) :: eps
449 INTEGER ,
INTENT(IN) :: mode
450 INTEGER :: m, l , i , ni
451 REAL(KIND=8) :: fs , ft , fexp
454 TYPE(mesh_type
),
TARGET :: mesh
455 INTEGER,
DIMENSION(:,:),
POINTER :: jj
456 INTEGER,
POINTER :: me
470 DO ni = 1, n_w; i = jj(ni,m)
471 ray = ray + mesh%rr(1,i)*ww(ni,l)
474 fs = sum(ff(jj(:,m)) * ww(:,l)) * rj(l,m)
476 IF (eps.EQ.(-1.d0))
THEN
477 ft = sum(4 * t1(1,jj(:,m)) * ww(:,l) - t2(1,jj(:,m))*ww(:,l))*rj(l,m) /(2* dt)
478 fexp = eps * mode * sum(gg(3,jj(:,m)) * (2*t1(2,jj(:,m))-t2(2,jj(:,m))) &
479 * ww(:,l))*rj(l,m)/ray
480 ELSEIF (eps.EQ.(1.d0))
THEN
481 ft = sum(4 * t1(2,jj(:,m)) * ww(:,l) - t2(2,jj(:,m))*ww(:,l))*rj(l,m) /(2* dt)
482 fexp = eps * mode * sum(gg(3,jj(:,m)) * (2*t1(1,jj(:,m))-t2(1,jj(:,m))) &
483 * ww(:,l))*rj(l,m)/ray
485 WRITE(*,*)
'probleme ds calcul second membre avec epsilon'
489 u0(jj(ni,m)) = u0(jj(ni,m)) + ww(ni,l) * (fs + ray * (ft + fexp))
508 REAL(KIND=8),
DIMENSION(:),
INTENT(IN) :: ff
509 REAL(KIND=8),
DIMENSION(:),
INTENT(OUT) :: u0
510 REAL(KIND=8),
DIMENSION(:),
INTENT(IN) :: t1
511 REAL(KIND=8),
INTENT(IN) :: dt
512 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN):: gg
513 REAL(KIND=8),
INTENT(IN) :: eps
514 INTEGER ,
INTENT(IN) :: mode
515 INTEGER :: m, l , i , ni
516 REAL(KIND=8) :: fs , ft , fexp
519 TYPE(mesh_type
),
TARGET :: mesh
520 INTEGER,
DIMENSION(:,:),
POINTER :: jj
521 INTEGER,
POINTER :: me
535 DO ni = 1, n_w; i = jj(ni,m)
536 ray = ray + mesh%rr(1,i)*ww(ni,l)
539 fs = sum(ff(jj(:,m)) * ww(:,l)) * rj(l,m)
540 ft = sum(t1(jj(:,m)) * ww(:,l)) * rj(l,m) / dt
541 fexp = eps * mode * sum(gg(3,jj(:,m)) * t1(jj(:,m)) * ww(:,l))*rj(l,m)/ray
544 u0(jj(ni,m)) = u0(jj(ni,m)) + ww(ni,l) * ray * (fs + ft + fexp)
561 REAL(KIND=8),
DIMENSION(:),
INTENT(IN) :: ff
562 REAL(KIND=8),
DIMENSION(:),
INTENT(OUT) :: u0
563 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: t1
564 REAL(KIND=8),
INTENT(IN) :: dt
565 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: gg
566 REAL(KIND=8),
INTENT(IN) :: eps
567 INTEGER ,
INTENT(IN) :: mode
568 INTEGER :: m, l , i , ni
569 REAL(KIND=8) :: fs , ft , fexp
572 TYPE(mesh_type
),
TARGET :: mesh
573 INTEGER,
DIMENSION(:,:),
POINTER :: jj
574 INTEGER,
POINTER :: me
588 DO ni = 1, n_w; i = jj(ni,m)
589 ray = ray + mesh%rr(1,i)*ww(ni,l)
592 IF (eps.EQ.(-1.d0))
THEN
593 fs = sum(ff(jj(:,m)) * ww(:,l)) * rj(l,m)
594 ft = sum(t1(1,jj(:,m)) * ww(:,l)) * rj(l,m) / dt
595 fexp = eps * mode * sum(gg(3,jj(:,m)) * t1(2,jj(:,m)) * ww(:,l))*rj(l,m)/ray
596 ELSEIF (eps.EQ.(1.d0))
THEN
597 fs = sum(ff(jj(:,m)) * ww(:,l)) * rj(l,m)
598 ft = sum(t1(2,jj(:,m)) * ww(:,l)) * rj(l,m) / dt
599 fexp = eps * mode * sum(gg(3,jj(:,m)) * t1(1,jj(:,m))* ww(:,l))*rj(l,m)/ray
601 WRITE(*,*)
'problem calc second membre init avec epsilon'
605 u0(jj(ni,m)) = u0(jj(ni,m)) + ww(ni,l) * (fs + ray *(ft + fexp))
614 SUBROUTINE qs_00_lap(mesh,alpha,mode,ff,V2,V1,dt,u0) !avec le rayon
624 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: ff
625 REAL(KIND=8),
DIMENSION(:,:),
INTENT(OUT) :: u0
626 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: v2
627 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: v1
628 REAL(KIND=8),
INTENT(IN) :: dt , alpha
629 INTEGER ,
INTENT(IN) :: mode
630 INTEGER :: m, l , i , ni , j
631 REAL(KIND=8),
DIMENSION(6) :: fs , ft
636 TYPE(mesh_type
),
TARGET :: mesh
637 INTEGER,
DIMENSION(:,:),
POINTER :: jj
638 INTEGER,
POINTER :: me
652 DO ni = 1, n_w; i = jj(ni,m)
653 ray = ray + mesh%rr(1,i)*ww(ni,l)
658 fs(1) = sum(ff(1,jj(:,m)) * ww(:,l)) * rj(l,m)
659 ft(1) = 1/(2*dt) * sum((4*v2(1,jj(:,m))-v1(1,jj(:,m))) * ww(:,l)) * rj(l,m)
663 fs(2) = sum(ff(2,jj(:,m)) * ww(:,l)) * rj(l,m)
664 ft(2) = 1/(2*dt)* sum((4*v2(2,jj(:,m))-v1(2,jj(:,m))) * ww(:,l)) * rj(l,m)
670 u0(j,jj(ni,m)) = u0(j,jj(ni,m)) + ww(ni,l) * ray * (fs(j) + ft(j))
690 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: ff
691 REAL(KIND=8),
DIMENSION(:,:),
INTENT(OUT) :: u0
692 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: v
693 REAL(KIND=8),
INTENT(IN) :: dt , alpha
694 INTEGER ,
INTENT(IN) :: mode
695 INTEGER :: m, l , i , ni , j
696 REAL(KIND=8),
DIMENSION(6) :: fs , ft
701 TYPE(mesh_type
),
TARGET :: mesh
702 INTEGER,
DIMENSION(:,:),
POINTER :: jj
703 INTEGER,
POINTER :: me
717 DO ni = 1, n_w; i = jj(ni,m)
718 ray = ray + mesh%rr(1,i)*ww(ni,l)
723 fs(1) = sum(ff(2,jj(:,m)) * ww(:,l)) * rj(l,m)
724 ft(1) = 1/(dt) * sum(v(2,jj(:,m)) * ww(:,l)) * rj(l,m)
728 fs(2) = sum(ff(2,jj(:,m)) * ww(:,l)) * rj(l,m)
729 ft(2) = 1/(dt) * sum(v(2,jj(:,m)) * ww(:,l)) * rj(l,m)
735 u0(j,jj(ni,m)) = u0(j,jj(ni,m)) + ww(ni,l) * ray * (fs(j) + ft(j))
749 REAL(KIND=8),
DIMENSION(:,:),
INTENT(OUT) :: rot
750 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: v
751 INTEGER ,
INTENT(IN) :: mode
752 INTEGER :: m, l , i , ni , j
753 REAL(KIND=8),
DIMENSION(6) ::
f
754 TYPE(mesh_type
),
TARGET :: mesh
755 INTEGER,
DIMENSION(:,:),
POINTER :: jj
756 INTEGER,
POINTER :: me
757 REAL(KIND=8) :: ray, rayj
758 REAL(KIND=8),
DIMENSION(6,mesh%gauss%n_w) :: v_loc
759 REAL(KIND=8),
DIMENSION(2,mesh%gauss%n_w) :: dw_loc
760 INTEGER,
DIMENSION(mesh%gauss%n_w) :: j_loc
775 DO ni = 1, n_w; i = jj(ni,m)
776 ray = ray + mesh%rr(1,i)*ww(ni,l)
781 f(1) = ( mode/ray*sum(v_loc(6,:)*ww(:,l))-sum(v_loc(3,:)*dw_loc(2,:)))*rayj
782 f(2) = (-mode/ray*sum(v_loc(5,:)*ww(:,l))-sum(v_loc(4,:)*dw_loc(2,:)))*rayj
783 f(3) = (sum(v_loc(1,:)*dw_loc(2,:))-sum(v_loc(5,:)*dw_loc(1,:)))*rayj
784 f(4) = (sum(v_loc(2,:)*dw_loc(2,:))-sum(v_loc(6,:)*dw_loc(1,:)))*rayj
785 f(5) = (1/ray*sum(v_loc(3,:)*ww(:,l))+sum(v_loc(3,:)*dw_loc(1,:))&
786 -mode/ray*sum(v_loc(2,:)*ww(:,l)))*rayj
787 f(6) = (1/ray*sum(v_loc(4,:)*ww(:,l))+sum(v_loc(4,:)*dw_loc(1,:))&
788 +mode/ray*sum(v_loc(1,:)*ww(:,l)))*rayj
792 rot(j,j_loc(ni)) = rot(j,j_loc(ni)) + ww(ni,l)*
f(j)
813 REAL(KIND=8),
DIMENSION(:,:),
INTENT(OUT) :: u0
814 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: v
815 INTEGER ,
INTENT(IN) :: mode
816 INTEGER :: m, l , i , ni , j
817 REAL(KIND=8),
DIMENSION(2) :: fr , fth , fz
823 TYPE(mesh_type
),
TARGET :: mesh
824 INTEGER,
DIMENSION(:,:),
POINTER :: jj
825 INTEGER,
POINTER :: me
842 DO ni = 1, n_w; i = jj(ni,m)
843 ray = ray + mesh%rr(1,i)*ww(ni,l)
855 fr(1) = -sum(v(1,jj(:,m))*dw(1,:,l,m))* rj(l,m)
856 fth(1) = mode*sum(v(4,jj(:,m))*ww(:,l))* rj(l,m)/ray
857 fz(1) = -sum(v(5,jj(:,m))*dw(2,:,l,m))* rj(l,m)
859 fr(2) = -sum(v(2,jj(:,m))*dw(1,:,l,m))* rj(l,m)
860 fth(2) = - mode*sum(v(3,jj(:,m))*ww(:,l))* rj(l,m)/ray
861 fz(2) = -sum(v(6,jj(:,m))*dw(2,:,l,m))* rj(l,m)
870 u0(j,jj(ni,m)) = u0(j,jj(ni,m)) + ray*ww(ni,l)*(fth(j) + fr(j) + fz(j))
892 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: gg
893 REAL(KIND=8),
DIMENSION(:,:),
INTENT(OUT) :: u0_c
894 INTEGER ,
INTENT(IN) :: mode
895 TYPE(mesh_type
),
TARGET :: uu_mesh, pp_mesh
896 INTEGER,
DIMENSION(:,:),
POINTER :: jj, jj_c
897 INTEGER,
POINTER :: me, nwc
899 REAL(KIND=8),
DIMENSION(pp_mesh%gauss%n_w, uu_mesh%gauss%l_G) :: w_c
900 INTEGER :: m, l, n, k,i,j
902 REAL(KIND=8),
DIMENSION(6,uu_mesh%gauss%n_w) :: ggkn
903 REAL(KIND=8),
DIMENSION(uu_mesh%gauss%k_d,uu_mesh%gauss%n_w) :: dwkn
904 REAL(KIND=8),
DIMENSION(2,uu_mesh%gauss%n_w) :: u0_cn
905 REAL(KIND=8),
DIMENSION(2) :: fr , fth , fz
906 REAL(KIND=8),
DIMENSION(3,2) ::
f
908 INTEGER,
DIMENSION(uu_mesh%gauss%n_w) :: j_loc
909 REAL(KIND=8),
DIMENSION(2) :: smb
912 REAL(KIND=8) :: tps, dummy, tt,tps1
919 nwc => pp_mesh%gauss%n_w
922 w_c(1,l) = ww(1,l) + 0.5*(ww(n_w-1,l) + ww(n_w,l))
923 w_c(2,l) = ww(2,l) + 0.5*(ww(n_w,l) + ww(4,l))
924 w_c(3,l) = ww(3,l) + 0.5*(ww(4,l) + ww(n_w-1,l))
941 DO n = 1, n_w; i = jj(n,m)
942 ray = ray + uu_mesh%rr(1,i)*ww(n,l)
948 f(1,1) = (ray*sum(gg(1,j_loc(:))*dw(1,:,l,m)) &
949 + sum(gg(1,j_loc(:))*ww(:,l)))
950 f(2,1) = mode*sum(gg(4,j_loc(:))*ww(:,l))
951 f(3,1) = sum(gg(5,j_loc(:))*dw(2,:,l,m)) * ray
954 f(1,2) = (ray*sum(gg(2,j_loc(:))*dw(1,:,l,m)) &
955 + sum(gg(2,j_loc(:))*ww(:,l)))
956 f(2,2) = - mode*sum(gg(3,j_loc(:))*ww(:,l))
957 f(3,2) = sum(gg(6,j_loc(:))*dw(2,:,l,m)) * ray
967 u0_c(j,jj_c(n,m)) = u0_c(j,jj_c(n,m)) + w_c(n,l)*(
f(1,j)+
f(2,j)+
f(3,j))
1020 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: gg
1021 REAL(KIND=8),
DIMENSION(:,:),
INTENT(OUT) :: u0_c
1022 INTEGER ,
INTENT(IN) :: mode
1023 TYPE(mesh_type
),
TARGET :: uu_mesh, pp_mesh
1024 INTEGER,
DIMENSION(:,:),
POINTER :: jj, jj_c
1025 INTEGER,
POINTER :: me, nwc
1027 REAL(KIND=8),
DIMENSION(pp_mesh%gauss%n_w, uu_mesh%gauss%l_G) :: w_c
1028 INTEGER :: m, l, n, k,i,j
1030 REAL(KIND=8),
DIMENSION(6,uu_mesh%gauss%n_w) :: ggkn
1031 REAL(KIND=8),
DIMENSION(uu_mesh%gauss%k_d,uu_mesh%gauss%n_w) :: dwkn
1032 REAL(KIND=8),
DIMENSION(2,uu_mesh%gauss%n_w) :: u0_cn
1033 REAL(KIND=8),
DIMENSION(2) :: fr , fth , fz
1034 REAL(KIND=8),
DIMENSION(3,2) ::
f
1036 INTEGER,
DIMENSION(uu_mesh%gauss%n_w) :: j_loc
1037 REAL(KIND=8),
DIMENSION(2) :: smb
1040 REAL(KIND=8) :: tps, dummy, tt,tps1
1047 nwc => pp_mesh%gauss%n_w
1050 w_c(1,l) = ww(1,l) + 0.5*(ww(n_w-1,l) + ww(n_w,l))
1051 w_c(2,l) = ww(2,l) + 0.5*(ww(n_w,l) + ww(4,l))
1052 w_c(3,l) = ww(3,l) + 0.5*(ww(4,l) + ww(n_w-1,l))
1065 DO n = 1, n_w; i = jj(n,m)
1066 ray = ray + uu_mesh%rr(1,i)*ww(n,l)
1072 f(1,1) = (ray*sum(gg(1,j_loc(:))*dw(1,:,l,m)) &
1073 + sum(gg(1,j_loc(:))*ww(:,l)))
1074 f(2,1) = mode*sum(gg(4,j_loc(:))*ww(:,l))
1075 f(3,1) = sum(gg(5,j_loc(:))*dw(2,:,l,m)) * ray
1078 f(1,2) = (ray*sum(gg(2,j_loc(:))*dw(1,:,l,m)) &
1079 + sum(gg(2,j_loc(:))*ww(:,l)))
1080 f(2,2) = - mode*sum(gg(3,j_loc(:))*ww(:,l))
1081 f(3,2) = sum(gg(6,j_loc(:))*dw(2,:,l,m)) * ray
1089 u0_c(j,jj_c(n,m)) = u0_c(j,jj_c(n,m)) + w_c(n,l)*(
f(1,j)+
f(2,j)+
f(3,j))
1109 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: gg
1110 REAL(KIND=8),
DIMENSION(:,:),
INTENT(OUT) :: u0_c
1111 INTEGER ,
INTENT(IN) :: mode
1112 TYPE(mesh_type
),
TARGET :: uu_mesh, pp_mesh
1113 INTEGER,
DIMENSION(:,:),
POINTER :: jj, jj_c
1114 INTEGER,
POINTER :: me, nwc
1116 REAL(KIND=8),
DIMENSION(pp_mesh%gauss%n_w, uu_mesh%gauss%l_G) :: w_c
1117 INTEGER :: m, l, n, k,i,j
1118 REAL(KIND=8),
DIMENSION(3,2) ::
f
1120 INTEGER,
DIMENSION(uu_mesh%gauss%n_w) :: j_loc
1127 nwc => pp_mesh%gauss%n_w
1133 w_c(1,l) = ww(1,l) + 0.5*(ww(n_w-1,l) + ww(n_w,l))
1134 w_c(2,l) = ww(2,l) + 0.5*(ww(n_w,l) + ww(4,l))
1135 w_c(3,l) = ww(3,l) + 0.5*(ww(4,l) + ww(n_w-1,l))
1150 DO n = 1, n_w; i = jj(n,m)
1151 ray = ray + uu_mesh%rr(1,i)*ww(n,l)
1156 f(1,1) = (ray*sum(gg(j_loc,1)*dw(1,:,l,m)) &
1157 + sum(gg(j_loc,1)*ww(:,l)))
1158 f(2,1) = mode*sum(gg(j_loc,4)*ww(:,l))
1159 f(3,1) = sum(gg(j_loc,5)*dw(2,:,l,m)) * ray
1162 f(1,2) = (ray*sum(gg(j_loc,2)*dw(1,:,l,m)) &
1163 + sum(gg(j_loc,2)*ww(:,l)))
1164 f(2,2) =-mode*sum(gg(j_loc,3)*ww(:,l))
1165 f(3,2) = sum(gg(j_loc,6)*dw(2,:,l,m)) * ray
1170 x =
f(1,j)+
f(2,j)+
f(3,j)
1172 u0_c(jj_c(n,m),j) = u0_c(jj_c(n,m),j) + w_c(n,l)*x
1190 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: pp
1191 REAL(KIND=8),
DIMENSION(:,:),
INTENT(OUT) :: u0
1192 INTEGER ,
INTENT(IN) :: mode
1193 INTEGER :: m, l , i , ni , j
1194 REAL(KIND=8),
DIMENSION(6) :: fp
1197 TYPE(mesh_type
),
TARGET :: mesh
1198 INTEGER,
DIMENSION(:,:),
POINTER :: jj
1199 INTEGER,
POINTER :: me
1213 DO ni = 1, n_w; i = jj(ni,m)
1214 ray = ray + mesh%rr(1,i)*ww(ni,l)
1218 fp(1) = sum(pp(1,jj(:,m))*dw(1,:,l,m)) * rj(l,m)
1219 fp(2) = sum(pp(2,jj(:,m))*dw(1,:,l,m)) * rj(l,m)
1222 fp(3) = mode/ray * sum(pp(2,jj(:,m))*ww(:,l)) * rj(l,m)
1223 fp(4) = -mode/ray * sum(pp(1,jj(:,m))*ww(:,l)) * rj(l,m)
1225 fp(5) = sum(pp(1,jj(:,m))*dw(2,:,l,m)) * rj(l,m)
1226 fp(6) = sum(pp(2,jj(:,m))*dw(2,:,l,m)) * rj(l,m)
1230 u0(j,jj(ni,m)) = u0(j,jj(ni,m)) + ww(ni,l) * ray * fp(j)
1246 TYPE(mesh_type
),
TARGET :: mesh
1247 INTEGER ,
INTENT(IN) :: mod_max
1248 REAL(KIND=8),
DIMENSION(2,mesh%np,0:mod_max),
INTENT(IN) :: pp
1249 REAL(KIND=8),
DIMENSION(6,mesh%np,0:mod_max),
INTENT(OUT) :: u0
1250 INTEGER :: m, l , i , ni , j, k
1251 REAL(KIND=8),
DIMENSION(6) :: fp
1255 INTEGER,
DIMENSION(:,:),
POINTER :: jj
1256 INTEGER,
POINTER :: me
1272 DO ni = 1, n_w; i = jj(ni,m)
1273 ray = ray + mesh%rr(1,i)*ww(ni,l)
1279 fp(1) = sum(pp(1,jj(:,m),k)*dw(1,:,l,m)) * rj(l,m)*ray
1280 fp(2) = sum(pp(2,jj(:,m),k)*dw(1,:,l,m)) * rj(l,m)*ray
1282 fp(3) = k* sum(pp(2,jj(:,m),k)*ww(:,l)) * rj(l,m)
1283 fp(4) = -k* sum(pp(1,jj(:,m),k)*ww(:,l)) * rj(l,m)
1285 fp(5) = sum(pp(1,jj(:,m),k)*dw(2,:,l,m)) * rj(l,m)*ray
1286 fp(6) = sum(pp(2,jj(:,m),k)*dw(2,:,l,m)) * rj(l,m)*ray
1290 u0(j,jj(ni,m),k) = u0(j,jj(ni,m),k) + ww(ni,l) * fp(j)
1321 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: ff
1322 REAL(KIND=8),
DIMENSION(:,:),
INTENT(OUT) :: u0
1323 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: v1
1324 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: v2
1325 REAL(KIND=8),
INTENT(IN) :: dt , alpha
1326 INTEGER ,
INTENT(IN) :: mode
1327 INTEGER :: m, l , i , ni , j
1328 REAL(KIND=8),
DIMENSION(6) :: fs , fexp , ft
1335 TYPE(mesh_type
),
TARGET :: mesh
1336 INTEGER,
DIMENSION(:,:),
POINTER :: jj
1337 INTEGER,
POINTER :: me
1354 DO ni = 1, n_w; i = jj(ni,m)
1355 ray = ray + mesh%rr(1,i)*ww(ni,l)
1360 fs(1) = sum(ff(1,jj(:,m)) * ww(:,l)) * rj(l,m)
1361 fexp(1) = -(alpha*2*mode/ray**2)*sum((2*v2(4,jj(:,m))-v1(4,jj(:,m)))* ww(:,l)) * rj(l,m)
1362 ft(1) = 1/(2*dt) * sum((4*v2(1,jj(:,m))-v1(1,jj(:,m))) * ww(:,l)) * rj(l,m)
1366 fs(2) = sum(ff(2,jj(:,m)) * ww(:,l)) * rj(l,m)
1367 fexp(2) = (alpha*2*mode/ray**2)*sum((2*v2(3,jj(:,m))-v1(3,jj(:,m)))* ww(:,l)) * rj(l,m)
1368 ft(2) = 1/(2*dt) * sum((4*v2(2,jj(:,m))-v1(2,jj(:,m))) * ww(:,l)) * rj(l,m)
1372 fs(3) = sum(ff(3,jj(:,m)) * ww(:,l)) * rj(l,m)
1373 fexp(3) = (alpha*2*mode/ray**2)*sum((2*v2(2,jj(:,m))-v1(2,jj(:,m)))* ww(:,l)) * rj(l,m)
1374 ft(3) = 1/(2*dt) * sum((4*v2(3,jj(:,m))-v1(3,jj(:,m))) * ww(:,l)) * rj(l,m)
1378 fs(4) = sum(ff(4,jj(:,m)) * ww(:,l)) * rj(l,m)
1379 fexp(4) = -(alpha*2*mode/ray**2)*sum((2*v2(1,jj(:,m))-v1(1,jj(:,m)))* ww(:,l)) * rj(l,m)
1380 ft(4) = 1/(2*dt) * sum((4*v2(4,jj(:,m))-v1(4,jj(:,m))) * ww(:,l)) * rj(l,m)
1384 fs(5) = sum(ff(5,jj(:,m)) * ww(:,l)) * rj(l,m)
1386 ft(5) = 1/(2*dt) * sum((4*v2(5,jj(:,m))-v1(5,jj(:,m))) * ww(:,l)) * rj(l,m)
1390 fs(6) = sum(ff(6,jj(:,m)) * ww(:,l)) * rj(l,m)
1392 ft(6) = 1/(2*dt) * sum((4*v2(6,jj(:,m))-v1(6,jj(:,m))) * ww(:,l)) * rj(l,m)
1399 u0(j,jj(ni,m)) = u0(j,jj(ni,m)) + ww(ni,l) * ray * (fs(j) + fexp(j) + ft(j))
1424 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: ff
1425 REAL(KIND=8),
DIMENSION(:,:),
INTENT(OUT) :: u0
1426 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: v1
1427 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: v2
1428 REAL(KIND=8),
INTENT(IN) :: dt , alpha
1429 INTEGER ,
INTENT(IN) :: mode
1430 INTEGER :: m, l , i , ni , j
1431 REAL(KIND=8),
DIMENSION(6) :: fs , fexp , ft
1438 TYPE(mesh_type
),
TARGET :: mesh
1439 INTEGER,
DIMENSION(:,:),
POINTER :: jj
1440 INTEGER,
POINTER :: me
1454 DO ni = 1, n_w; i = jj(ni,m)
1455 ray = ray + mesh%rr(1,i)*ww(ni,l)
1460 fs(1) = sum(ff(1,jj(:,m)) * ww(:,l)) * rj(l,m)
1461 fexp(1) = -(alpha*2*mode/ray**2)*sum((2*v2(4,jj(:,m))-v1(4,jj(:,m)))* ww(:,l)) * rj(l,m)
1462 ft(1) = 1/(2*dt) * sum((4*v2(1,jj(:,m))-v1(1,jj(:,m))) * ww(:,l)) * rj(l,m)
1466 fs(2) = sum(ff(2,jj(:,m)) * ww(:,l)) * rj(l,m)
1467 fexp(2) = (alpha*2*mode/ray**2)*sum((2*v2(3,jj(:,m))-v1(3,jj(:,m)))* ww(:,l)) * rj(l,m)
1468 ft(2) = 1/(2*dt) * sum((4*v2(2,jj(:,m))-v1(2,jj(:,m))) * ww(:,l)) * rj(l,m)
1472 fs(3) = sum(ff(3,jj(:,m)) * ww(:,l)) * rj(l,m)
1473 fexp(3) = (alpha*2*mode/ray**2)*sum((2*v2(2,jj(:,m))-v1(2,jj(:,m)))* ww(:,l)) * rj(l,m)
1474 ft(3) = 1/(2*dt) * sum((4*v2(3,jj(:,m))-v1(3,jj(:,m))) * ww(:,l)) * rj(l,m)
1478 fs(4) = sum(ff(4,jj(:,m)) * ww(:,l)) * rj(l,m)
1479 fexp(4) = -(alpha*2*mode/ray**2)*sum((2*v2(1,jj(:,m))-v1(1,jj(:,m)))* ww(:,l)) * rj(l,m)
1480 ft(4) = 1/(2*dt) * sum((4*v2(4,jj(:,m))-v1(4,jj(:,m))) * ww(:,l)) * rj(l,m)
1484 fs(5) = sum(ff(5,jj(:,m)) * ww(:,l)) * rj(l,m)
1486 ft(5) = 1/(2*dt) * sum((4*v2(5,jj(:,m))-v1(5,jj(:,m))) * ww(:,l)) * rj(l,m)
1490 fs(6) = sum(ff(6,jj(:,m)) * ww(:,l)) * rj(l,m)
1492 ft(6) = 1/(2*dt) * sum((4*v2(6,jj(:,m))-v1(6,jj(:,m))) * ww(:,l)) * rj(l,m)
1500 u0(j,jj(ni,m)) = u0(j,jj(ni,m)) + ww(ni,l) * (fs(j) + ray*(fexp(j) + ft(j)))
1524 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: ff
1525 REAL(KIND=8),
DIMENSION(:,:),
INTENT(OUT) :: u0
1526 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: v2
1527 REAL(KIND=8),
INTENT(IN) :: dt , alpha
1528 INTEGER ,
INTENT(IN) :: mode
1529 INTEGER :: m, l , i , ni , j
1530 REAL(KIND=8),
DIMENSION(6) :: fs , fexp , ft
1537 TYPE(mesh_type
),
TARGET :: mesh
1538 INTEGER,
DIMENSION(:,:),
POINTER :: jj
1539 INTEGER,
POINTER :: me
1553 DO ni = 1, n_w; i = jj(ni,m)
1554 ray = ray + mesh%rr(1,i)*ww(ni,l)
1559 fs(1) = sum(ff(1,jj(:,m)) * ww(:,l)) * rj(l,m)
1560 fexp(1) = (-alpha*2*mode/ray**2)*sum(v2(4,jj(:,m))* ww(:,l)) * rj(l,m)
1561 ft(1) = 1/(dt) * sum(v2(1,jj(:,m)) * ww(:,l)) * rj(l,m)
1565 fs(2) = sum(ff(2,jj(:,m)) * ww(:,l)) * rj(l,m)
1566 fexp(2) = (alpha*2*mode/ray**2)*sum(v2(3,jj(:,m))* ww(:,l)) * rj(l,m)
1567 ft(2) = 1/(dt) * sum(v2(2,jj(:,m)) * ww(:,l)) * rj(l,m)
1571 fs(3) = sum(ff(3,jj(:,m)) * ww(:,l)) * rj(l,m)
1572 fexp(3) = (alpha*2*mode/ray**2)*sum(v2(2,jj(:,m))* ww(:,l)) * rj(l,m)
1573 ft(3) = 1/(dt) * sum(v2(3,jj(:,m)) * ww(:,l)) * rj(l,m)
1577 fs(4) = sum(ff(4,jj(:,m)) * ww(:,l)) * rj(l,m)
1578 fexp(4) = (-alpha*2*mode/ray**2)*sum(v2(1,jj(:,m))* ww(:,l)) * rj(l,m)
1579 ft(4) = 1/(dt) * sum(v2(4,jj(:,m)) * ww(:,l)) * rj(l,m)
1583 fs(5) = sum(ff(5,jj(:,m)) * ww(:,l)) * rj(l,m)
1585 ft(5) = 1/(dt) * sum(v2(5,jj(:,m)) * ww(:,l)) * rj(l,m)
1589 fs(6) = sum(ff(6,jj(:,m)) * ww(:,l)) * rj(l,m)
1591 ft(6) = 1/(dt) * sum(v2(6,jj(:,m)) * ww(:,l)) * rj(l,m)
1597 u0(j,jj(ni,m)) = u0(j,jj(ni,m)) + ww(ni,l) * ray * (fs(j) + fexp(j) + ft(j))
1622 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: ff
1623 REAL(KIND=8),
DIMENSION(:,:),
INTENT(OUT) :: u0
1624 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: v
1625 REAL(KIND=8),
INTENT(IN) :: dt , alpha
1626 INTEGER ,
INTENT(IN) :: mode
1627 INTEGER :: m, l , i , ni , j
1628 REAL(KIND=8),
DIMENSION(6) :: fs , fexp , ft
1635 TYPE(mesh_type
),
TARGET :: mesh
1636 INTEGER,
DIMENSION(:,:),
POINTER :: jj
1637 INTEGER,
POINTER :: me
1653 DO ni = 1, n_w; i = jj(ni,m)
1654 ray = ray + mesh%rr(1,i)*ww(ni,l)
1659 fs(1) = sum(ff(1,jj(:,m)) * ww(:,l)) * rj(l,m)
1660 fexp(1) = -(alpha*2*mode/ray**2)*sum(v(4,jj(:,m))* ww(:,l)) * rj(l,m)
1661 ft(1) = 1/(dt) * sum(v(1,jj(:,m)) * ww(:,l)) * rj(l,m)
1665 fs(2) = sum(ff(2,jj(:,m)) * ww(:,l)) * rj(l,m)
1666 fexp(2) = (alpha*2*mode/ray**2)*sum(v(3,jj(:,m))* ww(:,l)) * rj(l,m)
1667 ft(2) = 1/(dt) * sum(v(2,jj(:,m)) * ww(:,l)) * rj(l,m)
1671 fs(3) = sum(ff(3,jj(:,m)) * ww(:,l)) * rj(l,m)
1672 fexp(3) = (alpha*2*mode/ray**2)*sum(v(2,jj(:,m))* ww(:,l)) * rj(l,m)
1673 ft(3) = 1/(dt) * sum(v(3,jj(:,m)) * ww(:,l)) * rj(l,m)
1677 fs(4) = sum(ff(4,jj(:,m)) * ww(:,l)) * rj(l,m)
1678 fexp(4) = -(alpha*2*mode/ray**2)*sum(v(1,jj(:,m))* ww(:,l)) * rj(l,m)
1679 ft(4) = 1/(dt) * sum(v(4,jj(:,m)) * ww(:,l)) * rj(l,m)
1683 fs(5) = sum(ff(5,jj(:,m)) * ww(:,l)) * rj(l,m)
1685 ft(5) = 1/(dt) * sum(v(5,jj(:,m)) * ww(:,l)) * rj(l,m)
1689 fs(6) = sum(ff(6,jj(:,m)) * ww(:,l)) * rj(l,m)
1691 ft(6) = 1/(dt) * sum(v(6,jj(:,m)) * ww(:,l)) * rj(l,m)
1699 u0(j,jj(ni,m)) = u0(j,jj(ni,m)) + ww(ni,l) * (fs(j) + ray*(fexp(j) + ft(j)))
1719 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: ff
1720 REAL(KIND=8),
DIMENSION(:,:),
INTENT(OUT) :: u0
1721 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: v
1722 REAL(KIND=8),
INTENT(IN) :: dt , alpha
1723 INTEGER ,
INTENT(IN) :: mode
1724 INTEGER :: m, l , i , ni , j
1725 REAL(KIND=8),
DIMENSION(6) :: fs , fexp , ft , fv
1726 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: gg
1732 TYPE(mesh_type
),
TARGET :: mesh
1733 INTEGER,
DIMENSION(:,:),
POINTER :: jj
1734 INTEGER,
POINTER :: me
1750 DO ni = 1, n_w; i = jj(ni,m)
1751 ray = ray + mesh%rr(1,i)*ww(ni,l)
1756 fs(1) = sum(ff(1,jj(:,m)) * ww(:,l)) * rj(l,m)
1757 fexp(1) = -(alpha*2*mode/ray**2)*sum(v(4,jj(:,m))* ww(:,l)) * rj(l,m)
1758 ft(1) = 1/(dt) * sum(v(1,jj(:,m)) * ww(:,l)) * rj(l,m)
1759 fv(1) = -mode/ray*sum(gg(3,jj(:,m))*v(2,jj(:,m))* ww(:,l))* rj(l,m)
1763 fs(2) = sum(ff(2,jj(:,m)) * ww(:,l)) * rj(l,m)
1764 fexp(2) = (alpha*2*mode/ray**2)*sum(v(3,jj(:,m))* ww(:,l)) * rj(l,m)
1765 ft(2) = 1/(dt) * sum(v(2,jj(:,m)) * ww(:,l)) * rj(l,m)
1766 fv(2) = mode/ray*sum(gg(3,jj(:,m))*v(1,jj(:,m))* ww(:,l))* rj(l,m)
1770 fs(3) = sum(ff(3,jj(:,m)) * ww(:,l)) * rj(l,m)
1771 fexp(3) = (alpha*2*mode/ray**2)*sum(v(2,jj(:,m))* ww(:,l)) * rj(l,m)
1772 ft(3) = 1/(dt) * sum(v(3,jj(:,m)) * ww(:,l)) * rj(l,m)
1773 fv(3) = -mode/ray*sum(gg(3,jj(:,m))*v(4,jj(:,m))* ww(:,l))* rj(l,m)
1777 fs(4) = sum(ff(4,jj(:,m)) * ww(:,l)) * rj(l,m)
1778 fexp(4) = -(alpha*2*mode/ray**2)*sum(v(1,jj(:,m))* ww(:,l)) * rj(l,m)
1779 ft(4) = 1/(dt) * sum(v(4,jj(:,m)) * ww(:,l)) * rj(l,m)
1780 fv(4) = mode/ray*sum(gg(3,jj(:,m))*v(3,jj(:,m))* ww(:,l))* rj(l,m)
1784 fs(5) = sum(ff(5,jj(:,m)) * ww(:,l)) * rj(l,m)
1786 ft(5) = 1/(dt) * sum(v(5,jj(:,m)) * ww(:,l)) * rj(l,m)
1787 fv(5) = -mode/ray*sum(gg(3,jj(:,m))*v(6,jj(:,m))* ww(:,l))* rj(l,m)
1791 fs(6) = sum(ff(6,jj(:,m)) * ww(:,l)) * rj(l,m)
1793 ft(6) = 1/(dt) * sum(v(6,jj(:,m)) * ww(:,l)) * rj(l,m)
1794 fv(6) = mode/ray*sum(gg(3,jj(:,m))*v(5,jj(:,m))* ww(:,l))* rj(l,m)
1802 u0(j,jj(ni,m)) = u0(j,jj(ni,m)) + ww(ni,l) * (fs(j) + ray*(fexp(j) + ft(j) + fv(j)))
1822 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: ff
1823 REAL(KIND=8),
DIMENSION(:,:),
INTENT(OUT) :: u0
1824 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: v
1825 REAL(KIND=8),
INTENT(IN) :: dt , alpha
1826 INTEGER ,
INTENT(IN) :: mode
1827 INTEGER :: m, l , i , ni , j
1828 REAL(KIND=8),
DIMENSION(6) :: fs , fexp , ft , fv
1829 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: gg
1835 TYPE(mesh_type
),
TARGET :: mesh
1836 INTEGER,
DIMENSION(:,:),
POINTER :: jj
1837 INTEGER,
POINTER :: me
1853 DO ni = 1, n_w; i = jj(ni,m)
1854 ray = ray + mesh%rr(1,i)*ww(ni,l)
1859 fs(1) = sum(ff(1,jj(:,m)) * ww(:,l)) * rj(l,m)
1860 fexp(1) = -(alpha*2*mode/ray**2)*sum(v(4,jj(:,m))* ww(:,l)) * rj(l,m)
1861 ft(1) = 1/(dt) * sum(v(1,jj(:,m)) * ww(:,l)) * rj(l,m)
1862 fv(1) = -mode/ray*sum(gg(3,jj(:,m))*v(2,jj(:,m))* ww(:,l))* rj(l,m)
1866 fs(2) = sum(ff(2,jj(:,m)) * ww(:,l)) * rj(l,m)
1867 fexp(2) = (alpha*2*mode/ray**2)*sum(v(3,jj(:,m))* ww(:,l)) * rj(l,m)
1868 ft(2) = 1/(dt) * sum(v(2,jj(:,m)) * ww(:,l)) * rj(l,m)
1869 fv(2) = mode/ray*sum(gg(3,jj(:,m))*v(1,jj(:,m))* ww(:,l))* rj(l,m)
1873 fs(3) = sum(ff(3,jj(:,m)) * ww(:,l)) * rj(l,m)
1874 fexp(3) = (alpha*2*mode/ray**2)*sum(v(2,jj(:,m))* ww(:,l)) * rj(l,m)
1875 ft(3) = 1/(dt) * sum(v(3,jj(:,m)) * ww(:,l)) * rj(l,m)
1876 fv(3) = -mode/ray*sum(gg(3,jj(:,m))*v(4,jj(:,m))* ww(:,l))* rj(l,m)
1880 fs(4) = sum(ff(4,jj(:,m)) * ww(:,l)) * rj(l,m)
1881 fexp(4) = -(alpha*2*mode/ray**2)*sum(v(1,jj(:,m))* ww(:,l)) * rj(l,m)
1882 ft(4) = 1/(dt) * sum(v(4,jj(:,m)) * ww(:,l)) * rj(l,m)
1883 fv(4) = mode/ray*sum(gg(3,jj(:,m))*v(3,jj(:,m))* ww(:,l))* rj(l,m)
1887 fs(5) = sum(ff(5,jj(:,m)) * ww(:,l)) * rj(l,m)
1889 ft(5) = 1/(dt) * sum(v(5,jj(:,m)) * ww(:,l)) * rj(l,m)
1890 fv(5) = -mode/ray*sum(gg(3,jj(:,m))*v(6,jj(:,m))* ww(:,l))* rj(l,m)
1894 fs(6) = sum(ff(6,jj(:,m)) * ww(:,l)) * rj(l,m)
1896 ft(6) = 1/(dt) * sum(v(6,jj(:,m)) * ww(:,l)) * rj(l,m)
1897 fv(6) = mode/ray*sum(gg(3,jj(:,m))*v(5,jj(:,m))* ww(:,l))* rj(l,m)
1905 u0(j,jj(ni,m)) = u0(j,jj(ni,m)) + ww(ni,l) * ray*(fs(j) + fexp(j) + ft(j) + fv(j))
1923 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: ff
1924 REAL(KIND=8),
DIMENSION(:,:),
INTENT(OUT) :: u0
1925 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: v1
1926 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: v2
1927 REAL(KIND=8),
INTENT(IN) :: dt , alpha
1928 INTEGER ,
INTENT(IN) :: mode
1929 INTEGER :: m, l , i , ni , j
1930 REAL(KIND=8),
DIMENSION(6) :: fs , fexp , ft , fv
1937 TYPE(mesh_type
),
TARGET :: mesh
1938 INTEGER,
DIMENSION(:,:),
POINTER :: jj
1939 INTEGER,
POINTER :: me
1941 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: gg
1958 DO ni = 1, n_w; i = jj(ni,m)
1959 ray = ray + mesh%rr(1,i)*ww(ni,l)
1964 fs(1) = sum(ff(1,jj(:,m)) * ww(:,l)) * rj(l,m)
1965 fexp(1) = -(alpha*2*mode/ray**2)*sum((2*v2(4,jj(:,m))-v1(4,jj(:,m)))* ww(:,l)) * rj(l,m)
1966 ft(1) = 1/(2*dt) * sum((4*v2(1,jj(:,m))-v1(1,jj(:,m))) * ww(:,l)) * rj(l,m)
1967 fv(1) = -mode/ray*sum(gg(3,jj(:,m))*(2*v2(2,jj(:,m)) - v1(2,jj(:,m)))* ww(:,l))*rj(l,m)
1971 fs(2) = sum(ff(2,jj(:,m)) * ww(:,l)) * rj(l,m)
1972 fexp(2) = (alpha*2*mode/ray**2)*sum((2*v2(3,jj(:,m))-v1(3,jj(:,m)))* ww(:,l))*rj(l,m)
1973 ft(2) = 1/(2*dt) * sum((4*v2(2,jj(:,m))-v1(2,jj(:,m))) * ww(:,l)) *rj(l,m)
1974 fv(2) = mode/ray*sum(gg(3,jj(:,m))*(2*v2(1,jj(:,m)) - v1(1,jj(:,m)))* ww(:,l))*rj(l,m)
1978 fs(3) = sum(ff(3,jj(:,m)) * ww(:,l)) * rj(l,m)
1979 fexp(3) = (alpha*2*mode/ray**2)*sum((2*v2(2,jj(:,m))-v1(2,jj(:,m)))* ww(:,l)) * rj(l,m)
1980 ft(3) = 1/(2*dt) * sum((4*v2(3,jj(:,m))-v1(3,jj(:,m))) * ww(:,l)) * rj(l,m)
1981 fv(3) = -mode/ray*sum(gg(3,jj(:,m))*(2*v2(4,jj(:,m)) - v1(4,jj(:,m)))* ww(:,l))*rj(l,m)
1985 fs(4) = sum(ff(4,jj(:,m)) * ww(:,l)) * rj(l,m)
1986 fexp(4) = -(alpha*2*mode/ray**2)*sum((2*v2(1,jj(:,m))-v1(1,jj(:,m)))* ww(:,l)) * rj(l,m)
1987 ft(4) = 1/(2*dt) * sum((4*v2(4,jj(:,m))-v1(4,jj(:,m))) * ww(:,l)) * rj(l,m)
1988 fv(4) = mode/ray*sum(gg(3,jj(:,m))*(2*v2(3,jj(:,m)) - v1(3,jj(:,m)))* ww(:,l))*rj(l,m)
1992 fs(5) = sum(ff(5,jj(:,m)) * ww(:,l)) * rj(l,m)
1994 ft(5) = 1/(2*dt) * sum((4*v2(5,jj(:,m))-v1(5,jj(:,m))) * ww(:,l)) * rj(l,m)
1995 fv(5) = -mode/ray*sum(gg(3,jj(:,m))*(2*v2(6,jj(:,m)) - v1(6,jj(:,m)))* ww(:,l))*rj(l,m)
1999 fs(6) = sum(ff(6,jj(:,m)) * ww(:,l)) * rj(l,m)
2001 ft(6) = 1/(2*dt) * sum((4*v2(6,jj(:,m))-v1(6,jj(:,m))) * ww(:,l)) * rj(l,m)
2002 fv(6) = mode/ray*sum(gg(3,jj(:,m))*(2*v2(5,jj(:,m)) - v1(5,jj(:,m)))* ww(:,l))*rj(l,m)
2009 u0(j,jj(ni,m)) = u0(j,jj(ni,m)) + ww(ni,l) * (fs(j) + ray*(fexp(j) + ft(j) + fv(j)))
2027 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: ff
2028 REAL(KIND=8),
DIMENSION(:,:),
INTENT(OUT) :: u0
2029 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: v1
2030 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: v2
2031 REAL(KIND=8),
INTENT(IN) :: dt , alpha
2032 INTEGER ,
INTENT(IN) :: mode
2033 INTEGER :: m, l , i , ni , j
2034 REAL(KIND=8),
DIMENSION(6) :: fs , fexp , ft , fv
2041 TYPE(mesh_type
),
TARGET :: mesh
2042 INTEGER,
DIMENSION(:,:),
POINTER :: jj
2043 INTEGER,
POINTER :: me
2045 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: gg
2062 DO ni = 1, n_w; i = jj(ni,m)
2063 ray = ray + mesh%rr(1,i)*ww(ni,l)
2068 fs(1) = sum(ff(1,jj(:,m)) * ww(:,l)) * rj(l,m)
2069 fexp(1) = -(alpha*2*mode/ray**2)*sum((2*v2(4,jj(:,m))-v1(4,jj(:,m)))* ww(:,l)) * rj(l,m)
2070 ft(1) = 1/(2*dt) * sum((4*v2(1,jj(:,m))-v1(1,jj(:,m))) * ww(:,l)) * rj(l,m)
2071 fv(1) = -mode/ray*sum(gg(3,jj(:,m))*(2*v2(2,jj(:,m)) - v1(2,jj(:,m)))* ww(:,l))*rj(l,m)
2075 fs(2) = sum(ff(2,jj(:,m)) * ww(:,l)) * rj(l,m)
2076 fexp(2) = (alpha*2*mode/ray**2)*sum((2*v2(3,jj(:,m))-v1(3,jj(:,m)))* ww(:,l))*rj(l,m)
2077 ft(2) = 1/(2*dt) * sum((4*v2(2,jj(:,m))-v1(2,jj(:,m))) * ww(:,l)) *rj(l,m)
2078 fv(2) = mode/ray*sum(gg(3,jj(:,m))*(2*v2(1,jj(:,m)) - v1(1,jj(:,m)))* ww(:,l))*rj(l,m)
2082 fs(3) = sum(ff(3,jj(:,m)) * ww(:,l)) * rj(l,m)
2083 fexp(3) = (alpha*2*mode/ray**2)*sum((2*v2(2,jj(:,m))-v1(2,jj(:,m)))* ww(:,l)) * rj(l,m)
2084 ft(3) = 1/(2*dt) * sum((4*v2(3,jj(:,m))-v1(3,jj(:,m))) * ww(:,l)) * rj(l,m)
2085 fv(3) = -mode/ray*sum(gg(3,jj(:,m))*(2*v2(4,jj(:,m)) - v1(4,jj(:,m)))* ww(:,l))*rj(l,m)
2089 fs(4) = sum(ff(4,jj(:,m)) * ww(:,l)) * rj(l,m)
2090 fexp(4) = -(alpha*2*mode/ray**2)*sum((2*v2(1,jj(:,m))-v1(1,jj(:,m)))* ww(:,l)) * rj(l,m)
2091 ft(4) = 1/(2*dt) * sum((4*v2(4,jj(:,m))-v1(4,jj(:,m))) * ww(:,l)) * rj(l,m)
2092 fv(4) = mode/ray*sum(gg(3,jj(:,m))*(2*v2(3,jj(:,m)) - v1(3,jj(:,m)))* ww(:,l))*rj(l,m)
2096 fs(5) = sum(ff(5,jj(:,m)) * ww(:,l)) * rj(l,m)
2098 ft(5) = 1/(2*dt) * sum((4*v2(5,jj(:,m))-v1(5,jj(:,m))) * ww(:,l)) * rj(l,m)
2099 fv(5) = -mode/ray*sum(gg(3,jj(:,m))*(2*v2(6,jj(:,m)) - v1(6,jj(:,m)))* ww(:,l))*rj(l,m)
2103 fs(6) = sum(ff(6,jj(:,m)) * ww(:,l)) * rj(l,m)
2105 ft(6) = 1/(2*dt) * sum((4*v2(6,jj(:,m))-v1(6,jj(:,m))) * ww(:,l)) * rj(l,m)
2106 fv(6) = mode/ray*sum(gg(3,jj(:,m))*(2*v2(5,jj(:,m)) - v1(5,jj(:,m)))* ww(:,l))*rj(l,m)
2113 u0(j,jj(ni,m)) = u0(j,jj(ni,m)) + ww(ni,l) * ray* (fs(j) + fexp(j) + ft(j) + fv(j))
2136 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: ff
2137 REAL(KIND=8),
DIMENSION(:,:),
INTENT(OUT) :: u0
2138 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: v
2139 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: p
2140 REAL(KIND=8),
INTENT(IN) :: dt , alpha
2141 INTEGER ,
INTENT(IN) :: mode
2142 INTEGER :: m, l , i , ni , j
2143 REAL(KIND=8),
DIMENSION(6) :: fs , fexp , ft , fp
2150 TYPE(mesh_type
),
TARGET :: mesh
2151 INTEGER,
DIMENSION(:,:),
POINTER :: jj
2152 INTEGER,
POINTER :: me
2169 DO ni = 1, n_w; i = jj(ni,m)
2170 ray = ray + mesh%rr(1,i)*ww(ni,l)
2175 fs(1) = sum(ff(1,jj(:,m)) * ww(:,l)) * rj(l,m)
2176 fexp(1) = -(alpha*2*mode/ray**2)*sum(v(4,jj(:,m))* ww(:,l)) * rj(l,m)
2177 ft(1) = 1/(dt) * sum(v(1,jj(:,m)) * ww(:,l)) * rj(l,m)
2179 fp(1) = sum(p(1,jj(:,m))*dw(1,:,l,m)) * rj(l,m)
2182 fs(2) = sum(ff(2,jj(:,m)) * ww(:,l)) * rj(l,m)
2183 fexp(2) = (alpha*2*mode/ray**2)*sum(v(3,jj(:,m))* ww(:,l)) * rj(l,m)
2184 ft(2) = 1/(dt) * sum(v(2,jj(:,m)) * ww(:,l)) * rj(l,m)
2186 fp(2) = sum(p(2,jj(:,m))*dw(1,:,l,m)) * rj(l,m)
2189 fs(3) = sum(ff(3,jj(:,m)) * ww(:,l)) * rj(l,m)
2190 fexp(3) = (alpha*2*mode/ray**2)*sum(v(2,jj(:,m))* ww(:,l)) * rj(l,m)
2191 ft(3) = 1/(dt) * sum(v(3,jj(:,m)) * ww(:,l)) * rj(l,m)
2192 fp(3) = sum(p(2,jj(:,m))*ww(:,l)) * rj(l,m)/ray*mode
2196 fs(4) = sum(ff(4,jj(:,m)) * ww(:,l)) * rj(l,m)
2197 fexp(4) = -(alpha*2*mode/ray**2)*sum(v(1,jj(:,m))* ww(:,l)) * rj(l,m)
2198 ft(4) = 1/(dt) * sum(v(4,jj(:,m)) * ww(:,l)) * rj(l,m)
2199 fp(4) = -sum(p(1,jj(:,m))*ww(:,l)) * rj(l,m)/ray *mode
2203 fs(5) = sum(ff(5,jj(:,m)) * ww(:,l)) * rj(l,m)
2205 ft(5) = 1/(dt) * sum(v(5,jj(:,m)) * ww(:,l)) * rj(l,m)
2207 fp(5) = sum(p(1,jj(:,m))*dw(2,:,l,m)) * rj(l,m)
2210 fs(6) = sum(ff(6,jj(:,m)) * ww(:,l)) * rj(l,m)
2212 ft(6) = 1/(dt) * sum(v(6,jj(:,m)) * ww(:,l)) * rj(l,m)
2214 fp(6) = sum(p(2,jj(:,m))*dw(2,:,l,m)) * rj(l,m)
2221 u0(j,jj(ni,m)) = u0(j,jj(ni,m)) + ww(ni,l) * ray* (fs(j) + fexp(j) + ft(j) + fp(j))
2255 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: ff
2256 REAL(KIND=8),
DIMENSION(:,:),
INTENT(OUT) :: u0
2257 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: v1
2258 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: v2
2259 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: p
2260 REAL(KIND=8),
INTENT(IN) :: dt , alpha
2261 INTEGER ,
INTENT(IN) :: mode
2262 INTEGER :: m, l , i , ni , j
2263 REAL(KIND=8),
DIMENSION(6) :: fs , fexp , ft , fp
2270 TYPE(mesh_type
),
TARGET :: mesh
2271 INTEGER,
DIMENSION(:,:),
POINTER :: jj
2272 INTEGER,
POINTER :: me
2289 DO ni = 1, n_w; i = jj(ni,m)
2290 ray = ray + mesh%rr(1,i)*ww(ni,l)
2295 fs(1) = sum(ff(1,jj(:,m)) * ww(:,l)) * rj(l,m)
2296 fexp(1) = (-alpha*2*mode/ray**2)*sum((2*v2(4,jj(:,m))-v1(4,jj(:,m)))* ww(:,l)) * rj(l,m)
2297 ft(1) = 1/(2*dt) * sum((4*v2(1,jj(:,m))-v1(1,jj(:,m))) * ww(:,l)) * rj(l,m)
2298 fp(1) = -sum(p(1,jj(:,m))*dw(1,:,l,m)) * rj(l,m)
2302 fs(2) = sum(ff(2,jj(:,m)) * ww(:,l)) * rj(l,m)
2303 fexp(2) = (alpha*2*mode/ray**2)*sum((2*v2(3,jj(:,m))-v1(3,jj(:,m)))* ww(:,l)) * rj(l,m)
2304 ft(2) = 1/(2*dt) * sum((4*v2(2,jj(:,m))-v1(2,jj(:,m))) * ww(:,l)) * rj(l,m)
2305 fp(2) = -sum(p(2,jj(:,m))*dw(1,:,l,m)) * rj(l,m)
2309 fs(3) = sum(ff(3,jj(:,m)) * ww(:,l)) * rj(l,m)
2310 fexp(3) = (alpha*2*mode/ray**2)*sum((2*v2(2,jj(:,m))-v1(2,jj(:,m)))* ww(:,l)) * rj(l,m)
2311 ft(3) = 1/(2*dt) * sum((4*v2(3,jj(:,m))-v1(3,jj(:,m))) * ww(:,l)) * rj(l,m)
2312 fp(3) = -sum(p(2,jj(:,m))*ww(:,l)) * rj(l,m)/ray*mode
2316 fs(4) = sum(ff(4,jj(:,m)) * ww(:,l)) * rj(l,m)
2317 fexp(4) = (-alpha*2*mode/ray**2)*sum((2*v2(1,jj(:,m))-v1(1,jj(:,m)))* ww(:,l)) * rj(l,m)
2318 ft(4) = 1/(2*dt) * sum((4*v2(4,jj(:,m))-v1(4,jj(:,m))) * ww(:,l)) * rj(l,m)
2319 fp(4) = sum(p(1,jj(:,m))*ww(:,l)) * rj(l,m)/ray *mode
2323 fs(5) = sum(ff(5,jj(:,m)) * ww(:,l)) * rj(l,m)
2325 ft(5) = 1/(2*dt) * sum((4*v2(5,jj(:,m))-v1(5,jj(:,m))) * ww(:,l)) * rj(l,m)
2326 fp(5) = -sum(p(1,jj(:,m))*dw(2,:,l,m)) * rj(l,m)
2330 fs(6) = sum(ff(6,jj(:,m)) * ww(:,l)) * rj(l,m)
2332 ft(6) = 1/(2*dt) * sum((4*v2(6,jj(:,m))-v1(6,jj(:,m))) * ww(:,l)) * rj(l,m)
2333 fp(6) = -sum(p(2,jj(:,m))*dw(2,:,l,m)) * rj(l,m)
2339 u0(j,jj(ni,m)) = u0(j,jj(ni,m)) + ww(ni,l) * ray * (fs(j) + fexp(j) + ft(j) + fp(j))
2356 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: ff
2357 REAL(KIND=8),
DIMENSION(:,:),
INTENT(OUT) :: u0
2358 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: v1
2359 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: v2
2360 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: p
2361 REAL(KIND=8),
INTENT(IN) :: dt , alpha
2362 INTEGER ,
INTENT(IN) :: mode
2363 INTEGER :: m, l , i , ni , j
2364 REAL(KIND=8),
DIMENSION(6) :: fs , ft , fp
2371 TYPE(mesh_type
),
TARGET :: mesh
2372 INTEGER,
DIMENSION(:,:),
POINTER :: jj
2373 INTEGER,
POINTER :: me
2390 DO ni = 1, n_w; i = jj(ni,m)
2391 ray = ray + mesh%rr(1,i)*ww(ni,l)
2396 fs(1) = sum(ff(1,jj(:,m)) * ww(:,l)) * rj(l,m)
2397 ft(1) = 1/(2*dt) * sum((4*v2(1,jj(:,m))-v1(1,jj(:,m))) * ww(:,l)) * rj(l,m)
2398 fp(1) = -sum(p(1,jj(:,m))*dw(1,:,l,m)) * rj(l,m)
2402 fs(2) = sum(ff(2,jj(:,m)) * ww(:,l)) * rj(l,m)
2403 ft(2) = 1/(2*dt) * sum((4*v2(2,jj(:,m))-v1(2,jj(:,m))) * ww(:,l)) * rj(l,m)
2404 fp(2) = -sum(p(2,jj(:,m))*dw(1,:,l,m)) * rj(l,m)
2408 fs(3) = sum(ff(3,jj(:,m)) * ww(:,l)) * rj(l,m)
2409 ft(3) = 1/(2*dt) * sum((4*v2(3,jj(:,m))-v1(3,jj(:,m))) * ww(:,l)) * rj(l,m)
2410 fp(3) = -sum(p(2,jj(:,m))*ww(:,l)) * rj(l,m)/ray*mode
2414 fs(4) = sum(ff(4,jj(:,m)) * ww(:,l)) * rj(l,m)
2415 ft(4) = 1/(2*dt) * sum((4*v2(4,jj(:,m))-v1(4,jj(:,m))) * ww(:,l)) * rj(l,m)
2416 fp(4) = sum(p(1,jj(:,m))*ww(:,l)) * rj(l,m)/ray *mode
2420 fs(5) = sum(ff(5,jj(:,m)) * ww(:,l)) * rj(l,m)
2421 ft(5) = 1/(2*dt) * sum((4*v2(5,jj(:,m))-v1(5,jj(:,m))) * ww(:,l)) * rj(l,m)
2422 fp(5) = -sum(p(1,jj(:,m))*dw(2,:,l,m)) * rj(l,m)
2426 fs(6) = sum(ff(6,jj(:,m)) * ww(:,l)) * rj(l,m)
2427 ft(6) = 1/(2*dt) * sum((4*v2(6,jj(:,m))-v1(6,jj(:,m))) * ww(:,l)) * rj(l,m)
2428 fp(6) = -sum(p(2,jj(:,m))*dw(2,:,l,m)) * rj(l,m)
2434 u0(j,jj(ni,m)) = u0(j,jj(ni,m)) + ww(ni,l) * (fs(j) + ray*(ft(j) + fp(j)))
2453 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: ff
2454 REAL(KIND=8),
DIMENSION(:,:),
INTENT(OUT) :: u0
2455 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: v1
2456 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: v2
2457 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: p
2458 REAL(KIND=8),
INTENT(IN) :: dt , alpha
2459 INTEGER ,
INTENT(IN) :: mode
2460 INTEGER :: m, l , i , ni , j
2461 REAL(KIND=8),
DIMENSION(6) :: fs , ft , fp
2469 TYPE(mesh_type
),
TARGET :: mesh
2470 INTEGER,
DIMENSION(:,:),
POINTER :: jj
2471 INTEGER,
POINTER :: me
2488 DO ni = 1, n_w; i = jj(ni,m)
2489 ray = ray + mesh%rr(1,i)*ww(ni,l)
2494 fs(1) = sum(ff(1,jj(:,m)) * ww(:,l)) * rj(l,m)
2495 ft(1) = 1/(2*dt) * sum((4*v2(1,jj(:,m))-v1(1,jj(:,m))) * ww(:,l)) * rj(l,m)
2496 fp(1) = -sum(p(1,jj(:,m))*dw(1,:,l,m)) * rj(l,m)
2500 fs(2) = sum(ff(2,jj(:,m)) * ww(:,l)) * rj(l,m)
2501 ft(2) = 1/(2*dt) * sum((4*v2(2,jj(:,m))-v1(2,jj(:,m))) * ww(:,l)) * rj(l,m)
2502 fp(2) = -sum(p(2,jj(:,m))*dw(1,:,l,m)) * rj(l,m)
2507 fs(3) = sum(ff(3,jj(:,m)) * ww(:,l)) * rj(l,m)
2508 ft(3) = 1/(2*dt) * sum((4*v2(3,jj(:,m))-v1(3,jj(:,m))) * ww(:,l)) * rj(l,m)
2509 fp(3) = -sum(p(2,jj(:,m))*ww(:,l)) * rj(l,m)/ray*mode
2514 fs(4) = sum(ff(4,jj(:,m)) * ww(:,l)) * rj(l,m)
2515 ft(4) = 1/(2*dt) * sum((4*v2(4,jj(:,m))-v1(4,jj(:,m))) * ww(:,l)) * rj(l,m)
2516 fp(4) = sum(p(1,jj(:,m))*ww(:,l)) * rj(l,m)/ray *mode
2520 fs(5) = sum(ff(5,jj(:,m)) * ww(:,l)) * rj(l,m)
2521 ft(5) = 1/(2*dt) * sum((4*v2(5,jj(:,m))-v1(5,jj(:,m))) * ww(:,l)) * rj(l,m)
2522 fp(5) = -sum(p(1,jj(:,m))*dw(2,:,l,m)) * rj(l,m)
2526 fs(6) = sum(ff(6,jj(:,m)) * ww(:,l)) * rj(l,m)
2527 ft(6) = 1/(2*dt) * sum((4*v2(6,jj(:,m))-v1(6,jj(:,m))) * ww(:,l)) * rj(l,m)
2528 fp(6) = -sum(p(2,jj(:,m))*dw(2,:,l,m)) * rj(l,m)
2534 u0(j,jj(ni,m)) = u0(j,jj(ni,m)) + ww(ni,l) * ray * (fs(j) + ft(j) + fp(j))
2545 SUBROUTINE qs_00_ns_inline(mesh,mode,m_max,ff,V1,V2,P,dt,u0,meth,tps_sft)
2557 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: ff
2558 REAL(KIND=8),
DIMENSION(:,:),
INTENT(OUT) :: u0
2559 TYPE(mesh_type
),
TARGET :: mesh
2560 INTEGER ,
INTENT(IN) :: m_max
2561 REAL(KIND=8),
DIMENSION(6,mesh%np,0:m_max),
INTENT(IN) :: v1
2562 REAL(KIND=8),
DIMENSION(6,mesh%np,0:m_max),
INTENT(IN) :: v2
2563 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: p
2564 REAL(KIND=8),
INTENT(IN) :: dt
2565 INTEGER ,
INTENT(IN) :: mode
2566 CHARACTER(len=200),
INTENT(IN) :: meth
2567 REAL(KIND=8),
OPTIONAL,
INTENT(INOUT):: tps_sft
2569 INTEGER :: m, l , i , ni , j
2570 REAL(KIND=8),
DIMENSION(6) :: fs , ft , fp, fnl, smb
2571 REAL(KIND=8),
ALLOCATABLE,
DIMENSION(:,:,:,:),
SAVE :: w, rotv, prodi
2573 INTEGER,
ALLOCATABLE,
DIMENSION(:),
SAVE :: j_loc
2574 REAL(KIND=8),
ALLOCATABLE,
DIMENSION(:,:),
SAVE :: dw_loc
2575 REAL(KIND=8),
DIMENSION(6,mesh%np) :: v1m, v2m, vt
2576 REAL(KIND=8),
DIMENSION(6,mesh%np,0:m_max) :: vs
2577 LOGICAL,
SAVE :: once = .true.
2578 LOGICAL,
SAVE :: once_sft = .true.
2581 REAL(KIND=8) :: tps, dummy, tt,tps1
2594 INTEGER,
DIMENSION(:,:),
POINTER :: jj
2595 INTEGER,
POINTER :: me
2599 REAL(KIND=8),
ALLOCATABLE,
DIMENSION(:,:),
SAVE :: aps, asp
2602 REAL(KIND=8),
DIMENSION(0:2*m_max) :: vr_p, vth_p, vz_p, rotr_p, rotth_p, rotz_p
2603 REAL(KIND=8),
DIMENSION(3,0:2*m_max) :: prod
2621 ALLOCATE(asp(0:2*m_max, 0:2*m_max))
2622 ALLOCATE(aps(0:2*m_max, 0:2*m_max))
2628 IF(ki < m_max+1)
THEN
2629 asp(li, ki) = cos(li*2.d0*pi/(2*m_max+1)*ki)
2631 asp(li, ki) = sin(li*2.d0*pi/(2*m_max+1)*(ki-m_max))
2640 IF(ki < m_max+1)
THEN
2644 aps(ki, li) = cos(li*2.d0*pi/(2*m_max+1)*ki)
2647 aps(ki, li) = sin(li*2.d0*pi/(2*m_max+1)*(ki-m_max))
2649 aps(ki, li) = aps(ki, li) * 2.d0/(float(2*m_max)+1.d0)
2652 WRITE(*,*)
'calcul matrices sft done'
2654 ALLOCATE(w(3,0:2*m_max,l_g,me))
2655 ALLOCATE(rotv(3,0:2*m_max,l_g,me))
2656 ALLOCATE(prodi(3,0:2*m_max,l_g,me))
2657 ALLOCATE(j_loc(mesh%gauss%n_w))
2658 ALLOCATE(dw_loc(mesh%gauss%k_d,mesh%gauss%n_w))
2679 dw_loc = dw(:,:,l,m)
2682 ray = sum(mesh%rr(1,j_loc)*ww(:,l))
2685 fs(j) = sum(ff(j,j_loc) * ww(:,l))
2686 ft(j) = 1/(2*dt)*sum((4*v2m(j,j_loc)-v1m(j,j_loc)) * ww(:,l))
2689 fp(1) = -sum(p(1,j_loc)*dw_loc(1,:))
2690 fp(2) = -sum(p(2,j_loc)*dw_loc(1,:))
2691 fp(3) = -sum(p(2,j_loc)*ww(:,l))/ray*mode
2692 fp(4) = sum(p(1,j_loc)*ww(:,l))/ray *mode
2693 fp(5) = -sum(p(1,j_loc)*dw_loc(2,:))
2694 fp(6) = -sum(p(2,j_loc)*dw_loc(2,:))
2705 rotv(1,j,l,m) = j/ray*sum(vs(6,j_loc,j)*ww(:,l)) &
2706 -sum(vs(3,j_loc,j)*dw_loc(2,:))
2707 rotv(2,j,l,m) = sum(vs(1,j_loc,j)*dw_loc(2,:)) &
2708 -sum(vs(5,j_loc,j)*dw_loc(1,:))
2709 rotv(3,j,l,m) = 1/ray*sum(vs(3,j_loc,j)*ww(:,l))+sum(vs(3,j_loc,j) &
2710 *dw_loc(1,:))-j/ray*sum(vs(2,j_loc,j)*ww(:,l))
2713 rotv(1,j+m_max,l,m) = -j/ray*sum(vs(5,j_loc,j)*ww(:,l)) &
2714 -sum(vs(4,j_loc,j)*dw_loc(2,:))
2715 rotv(2,j+m_max,l,m) = sum(vs(2,j_loc,j)*dw_loc(2,:)) &
2716 -sum(vs(6,j_loc,j)*dw_loc(1,:))
2717 rotv(3,j+m_max,l,m) = 1/ray*sum(vs(4,j_loc,j)*ww(:,l))+sum(vs(4,j_loc,j) &
2718 *dw_loc(1,:))+j/ray*sum(vs(1,j_loc,j)*ww(:,l))
2721 w(1,j,l,m) = sum(vs(1,j_loc,j)*ww(:,l))
2722 w(2,j,l,m) = sum(vs(3,j_loc,j)*ww(:,l))
2723 w(3,j,l,m) = sum(vs(5,j_loc,j)*ww(:,l))
2725 w(1,j+m_max,l,m) = sum(vs(2,j_loc,j)*ww(:,l))
2726 w(2,j+m_max,l,m) = sum(vs(4,j_loc,j)*ww(:,l))
2727 w(3,j+m_max,l,m) = sum(vs(6,j_loc,j)*ww(:,l))
2735 IF (meth ==
'sum')
THEN
2737 vr_p(li) = sum(asp(li,:)*w(1,:,l,m))
2738 vth_p(li) = sum(asp(li,:)*w(2,:,l,m))
2739 vz_p(li) = sum(asp(li,:)*w(3,:,l,m))
2740 rotr_p(li) = sum(asp(li,:)*rotv(1,:,l,m))
2741 rotth_p(li) = sum(asp(li,:)*rotv(2,:,l,m))
2742 rotz_p(li) = sum(asp(li,:)*rotv(3,:,l,m))
2744 ELSEIF (meth ==
'matmul')
THEN
2745 vr_p = matmul(asp,w(1,:,l,m) )
2746 vth_p = matmul(asp,w(2,:,l,m))
2747 vz_p = matmul(asp,w(3,:,l,m))
2749 rotr_p = matmul(asp,rotv(1,:,l,m))
2750 rotth_p = matmul(asp,rotv(2,:,l,m))
2751 rotz_p = matmul(asp,rotv(3,:,l,m))
2754 prod(1,:) = rotz_p(:) * vth_p(:) - rotth_p(:) * vz_p(:)
2755 prod(2,:) = rotr_p(:) * vz_p(:) - rotz_p(:) * vr_p(:)
2756 prod(3,:) = rotth_p(:) * vr_p(:) - rotr_p(:) * vth_p(:)
2758 IF (meth ==
'sum')
THEN
2760 prodi(1,li,l,m) = sum(aps(li,:)*prod(1,:))
2761 prodi(2,li,l,m) = sum(aps(li,:)*prod(2,:))
2762 prodi(3,li,l,m) = sum(aps(li,:)*prod(3,:))
2764 ELSEIF (meth ==
'matmul')
THEN
2765 prodi(1,:,l,m) = matmul(aps, prod(1,:))
2766 prodi(2,:,l,m) = matmul(aps, prod(2,:))
2767 prodi(3,:,l,m) = matmul(aps, prod(3,:))
2770 IF (present(tps_sft))
THEN
2771 tps_sft = tps_sft +
user_time(dummy) - tt
2778 fnl(1) = prodi(1,mode,l,m)
2779 fnl(3) = prodi(2,mode,l,m)
2780 fnl(5) = prodi(3,mode,l,m)
2782 fnl(2) = prodi(1,mode+m_max,l,m)
2783 fnl(4) = prodi(2,mode+m_max,l,m)
2784 fnl(6) = prodi(3,mode+m_max,l,m)
2787 smb = (fnl+ft+fp+fs)*ray*rj(l,m)
2793 u0(j,jj(ni,m)) = u0(j,jj(ni,m)) + ww(ni,l)*smb(j)
2814 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: ff
2815 REAL(KIND=8),
DIMENSION(:,:),
INTENT(OUT) :: u0
2816 TYPE(mesh_type
),
TARGET :: mesh
2817 REAL(KIND=8),
DIMENSION(6,mesh%np),
INTENT(IN) :: v1m, v2m
2818 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: p
2819 REAL(KIND=8),
INTENT(IN) :: dt
2820 INTEGER ,
INTENT(IN) :: mode
2822 INTEGER :: m, l , i , ni , j
2823 REAL(KIND=8),
DIMENSION(6) :: fs , ft , fp, fnl, smb
2825 INTEGER,
DIMENSION(mesh%gauss%n_w) :: j_loc
2826 REAL(KIND=8),
DIMENSION(mesh%gauss%k_d,mesh%gauss%n_w) :: dw_loc
2827 LOGICAL,
SAVE :: once = .true.
2828 LOGICAL,
SAVE :: once_sft = .true.
2834 INTEGER,
DIMENSION(:,:),
POINTER :: jj
2835 INTEGER,
POINTER :: me
2847 dw_loc = dw(:,:,l,m)
2850 ray = sum(mesh%rr(1,j_loc)*ww(:,l))
2855 fs(1) = sum(ff(1,j_loc) * ww(:,l))
2856 ft(1) = 1/(2*dt) * sum((4*v2m(1,j_loc)-v1m(1,j_loc)) * ww(:,l))
2857 fp(1) = -sum(p(1,j_loc)*dw_loc(1,:))
2861 fs(2) = sum(ff(2,j_loc) * ww(:,l))
2862 ft(2) = 1/(2*dt) * sum((4*v2m(2,j_loc)-v1m(2,j_loc)) * ww(:,l))
2863 fp(2) = -sum(p(2,j_loc)*dw_loc(1,:))
2867 fs(3) = sum(ff(3,j_loc) * ww(:,l))
2868 ft(3) = 1/(2*dt) * sum((4*v2m(3,j_loc)-v1m(3,j_loc)) * ww(:,l))
2869 fp(3) = -sum(p(2,j_loc)*ww(:,l))/ray*mode
2873 fs(4) = sum(ff(4,j_loc) * ww(:,l))
2874 ft(4) = 1/(2*dt) * sum((4*v2m(4,j_loc)-v1m(4,j_loc)) * ww(:,l))
2875 fp(4) = sum(p(1,j_loc)*ww(:,l))/ray *mode
2879 fs(5) = sum(ff(5,j_loc) * ww(:,l))
2880 ft(5) = 1/(2*dt) * sum((4*v2m(5,j_loc)-v1m(5,j_loc)) * ww(:,l))
2881 fp(5) = -sum(p(1,j_loc)*dw_loc(2,:))
2885 fs(6) = sum(ff(6,j_loc) * ww(:,l))
2886 ft(6) = 1/(2*dt) * sum((4*v2m(6,j_loc)-v1m(6,j_loc)) * ww(:,l))
2887 fp(6) = -sum(p(2,j_loc)*dw_loc(2,:))
2892 smb = (ft+fp+fs)*ray*rj(l,m)
2898 u0(j,jj(ni,m)) = u0(j,jj(ni,m)) + ww(ni,l)*smb(j)
2916 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: ff, rotv_v
2917 REAL(KIND=8),
DIMENSION(:,:),
INTENT(OUT):: u0
2918 TYPE(mesh_type
),
TARGET :: mesh
2919 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: v1m
2920 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: p
2921 REAL(KIND=8),
INTENT(IN) :: dt
2922 INTEGER,
INTENT(IN) :: mode
2923 INTEGER :: m, l , i , ni , j, index
2924 REAL(KIND=8),
DIMENSION(6) :: fs , ft , fp, smb
2925 INTEGER,
DIMENSION(mesh%gauss%n_w) :: j_loc
2926 REAL(KIND=8),
DIMENSION(mesh%gauss%k_d,mesh%gauss%n_w) :: dw_loc
2932 INTEGER,
DIMENSION(:,:),
POINTER :: jj
2933 INTEGER,
POINTER :: me
2948 dw_loc = dw(:,:,l,m)
2951 ray = sum(mesh%rr(1,j_loc)*ww(:,l))
2956 fs(1) = sum(ff(j_loc,1) * ww(:,l))
2957 ft(1) = sum(v1m(j_loc,1) * ww(:,l))
2958 fp(1) = -sum(p(j_loc,1)*dw_loc(1,:))
2962 fs(2) = sum(ff(j_loc,2) * ww(:,l))
2963 ft(2) = sum(v1m(j_loc,2) * ww(:,l))
2964 fp(2) = -sum(p(j_loc,2)*dw_loc(1,:))
2968 fs(3) = sum(ff(j_loc,3) * ww(:,l))
2969 ft(3) = sum(v1m(j_loc,3) * ww(:,l))
2970 fp(3) = -sum(p(j_loc,2)*ww(:,l))/ray*mode
2974 fs(4) = sum(ff(j_loc,4) * ww(:,l))
2975 ft(4) = sum(v1m(j_loc,4) * ww(:,l))
2976 fp(4) = sum(p(j_loc,1)*ww(:,l))/ray *mode
2980 fs(5) = sum(ff(j_loc,5) * ww(:,l))
2981 ft(5) = sum(v1m(j_loc,5) * ww(:,l))
2982 fp(5) = -sum(p(j_loc,1)*dw_loc(2,:))
2986 fs(6) = sum(ff(j_loc,6) * ww(:,l))
2987 ft(6) = sum(v1m(j_loc,6) * ww(:,l))
2988 fp(6) = -sum(p(j_loc,2)*dw_loc(2,:))
2993 smb = (ft+fp+fs-rotv_v(:,index))*ray*rj(l,m)
2998 u0(jj(ni,m),j) = u0(jj(ni,m),j) + ww(ni,l)*smb(j)
3016 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: ff, rotv_v
3017 REAL(KIND=8),
DIMENSION(:,:),
INTENT(OUT):: u0
3018 TYPE(mesh_type
),
TARGET :: mesh
3019 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: v1m
3020 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: p
3021 REAL(KIND=8),
INTENT(IN) :: dt
3022 INTEGER,
INTENT(IN) :: mode
3023 INTEGER :: m, l , i , ni , j, index
3024 REAL(KIND=8),
DIMENSION(6) :: fs , ft , fp, smb
3025 INTEGER,
DIMENSION(mesh%gauss%n_w) :: j_loc
3026 REAL(KIND=8),
DIMENSION(mesh%gauss%k_d,mesh%gauss%n_w) :: dw_loc
3032 INTEGER,
DIMENSION(:,:),
POINTER :: jj
3033 INTEGER,
POINTER :: me
3048 dw_loc = dw(:,:,l,m)
3051 ray = sum(mesh%rr(1,j_loc)*ww(:,l))
3056 fs(1) = sum(ff(j_loc,1) * ww(:,l))
3057 ft(1) = sum(v1m(j_loc,1) * ww(:,l))
3058 fp(1) = -sum(p(j_loc,1)*dw_loc(1,:))
3062 fs(2) = sum(ff(j_loc,2) * ww(:,l))
3063 ft(2) = sum(v1m(j_loc,2) * ww(:,l))
3064 fp(2) = -sum(p(j_loc,2)*dw_loc(1,:))
3068 fs(3) = sum(ff(j_loc,3) * ww(:,l))
3069 ft(3) = sum(v1m(j_loc,3) * ww(:,l))
3070 fp(3) = -sum(p(j_loc,2)*ww(:,l))/ray*mode
3074 fs(4) = sum(ff(j_loc,4) * ww(:,l))
3075 ft(4) = sum(v1m(j_loc,4) * ww(:,l))
3076 fp(4) = sum(p(j_loc,1)*ww(:,l))/ray *mode
3080 fs(5) = sum(ff(j_loc,5) * ww(:,l))
3081 ft(5) = sum(v1m(j_loc,5) * ww(:,l))
3082 fp(5) = -sum(p(j_loc,1)*dw_loc(2,:))
3086 fs(6) = sum(ff(j_loc,6) * ww(:,l))
3087 ft(6) = sum(v1m(j_loc,6) * ww(:,l))
3088 fp(6) = -sum(p(j_loc,2)*dw_loc(2,:))
3093 smb = (ft+fp+fs-rotv_v(index,:))*ray*rj(l,m)
3098 u0(jj(ni,m),j) = u0(jj(ni,m),j) + ww(ni,l)*smb(j)
3108 SUBROUTINE qs_ns_stab_new(mesh,mode,ff,vel_tot,V1m,vit,P,dudt,phalf,nlhalf,dt,u0,rotv_v)
3118 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: ff, rotv_v, dudt, nlhalf
3119 REAL(KIND=8),
DIMENSION(:),
INTENT(IN) :: vel_tot
3120 REAL(KIND=8),
DIMENSION(:,:),
INTENT(OUT):: u0
3121 TYPE(mesh_type
),
TARGET :: mesh
3122 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: v1m, vit
3123 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: p, phalf
3124 REAL(KIND=8),
INTENT(IN) :: dt
3125 INTEGER,
INTENT(IN) :: mode
3126 INTEGER :: m, l , i , ni , j, index, index2, type, k
3127 REAL(KIND=8),
DIMENSION(6) :: fs , ft , fp, smb, fv, mult
3128 REAL(KIND=8),
DIMENSION(2,6,mesh%gauss%l_G) :: grad
3129 REAL(KIND=8),
DIMENSION(6,mesh%gauss%l_G) :: vitloc
3130 INTEGER,
DIMENSION(mesh%gauss%n_w) :: j_loc
3131 REAL(KIND=8),
DIMENSION(mesh%gauss%k_d,mesh%gauss%n_w) :: dw_loc
3132 REAL(KIND=8),
DIMENSION(mesh%me) :: visc_plot
3137 INTEGER,
DIMENSION(:,:),
POINTER :: jj
3138 INTEGER,
POINTER :: me
3139 REAL(KIND=8),
DIMENSION(6) :: visc1, visc2
3140 REAL(KIND=8) :: ray, div1, div2, hloc, cfl, vloc, normal_vit
3141 REAL(KIND=8),
SAVE :: coeff_ed_st, r_eff, visc_eff, surf, nu_loc, coeff_visc_ordre_un
3142 LOGICAL,
SAVE :: once = .true.
3143 REAL(KIND=8),
ALLOCATABLE,
DIMENSION(:),
SAVE :: h
3145 REAL(KIND=8),
DIMENSION(mesh%gauss%n_w,mesh%gauss%l_Gs, 2) :: dwni_loc
3146 INTEGER,
DIMENSION(mesh%gauss%n_w,2) :: jji_loc
3147 REAL(KIND=8),
DIMENSION(mesh%gauss%n_w,2) :: u0loci, uloci
3148 REAL(KIND=8),
DIMENSION(6,mesh%me) :: viscosity
3149 REAL(KIND=8),
DIMENSION(6) :: norm_vit
3150 REAL(KIND=8) :: dul, h2, type_fe, ed_st
3151 INTEGER :: ms, ls, cotei
3163 ray = sum(mesh%rr(1,mesh%jj(:,m))*mesh%gauss%ww(:,l))
3164 surf = surf + rj(l,m)*ray
3168 IF (mesh%gauss%n_w==3)
THEN
3173 coeff_ed_st = inputs%coeff3*0.02d0/type_fe
3174 coeff_visc_ordre_un = inputs%coeff4*0.25d0
3175 ALLOCATE(h(mesh%mi))
3177 h(ms) = sum(mesh%gauss%rji(:,ms))/type_fe
3184 normal_vit = maxval(vel_tot)
3186 norm_vit(type) = sum(abs(vit(:,type)))/mesh%np + 1.d-14
3194 IF (inputs%LES)
THEN
3197 hloc = sqrt(sum(rj(:,m)))
3198 vloc = maxval(vel_tot(j_loc))
3199 cfl = max(vloc*dt/hloc,cfl)
3204 dw_loc = dw(:,:,l,m)
3207 ray = sum(mesh%rr(1,j_loc)*ww(:,l))
3210 ft(1) = sum(dudt(j_loc,1) * ww(:,l))
3211 fp(1) = sum(phalf(j_loc,1)*dw_loc(1,:))
3213 ft(2) = sum(dudt(j_loc,2) * ww(:,l))
3214 fp(2) = sum(phalf(j_loc,2)*dw_loc(1,:))
3216 ft(3) = sum(dudt(j_loc,3) * ww(:,l))
3217 fp(3) = sum(phalf(j_loc,2)*ww(:,l))*mode/ray
3219 ft(4) = sum(dudt(j_loc,4) * ww(:,l))
3220 fp(4) = -sum(phalf(j_loc,1)*ww(:,l))*mode/ray
3222 ft(5) = sum(dudt(j_loc,5) * ww(:,l))
3223 fp(5) = sum(phalf(j_loc,1)*dw_loc(2,:))
3225 ft(6) = sum(dudt(j_loc,6) * ww(:,l))
3226 fp(6) = sum(phalf(j_loc,2)*dw_loc(2,:))
3229 visc1= max(visc1,abs(ft+fp+nlhalf(index2,:)))
3234 grad(k,type,l) = sum(vit(j_loc,type)*dw_loc(k,:))
3236 vitloc(type,l) = sum(vit(j_loc,type)*ww(:,l))
3240 div1 = abs(grad(1,1,l) + vitloc(1,l)/ray + grad(2,5,l) + vitloc(4,l)*mode/ray)
3241 div2 = abs(grad(1,2,l) + vitloc(2,l)/ray + grad(2,6,l) - vitloc(3,l)*mode/ray)
3242 visc2(1) = max(visc2(1),div1)
3243 visc2(2) = max(visc2(1),div2)
3245 visc2(4) = visc2(1); visc2(5) = visc2(1)
3246 visc2(3) = visc2(2); visc2(6) = visc2(2)
3252 visc1(type) = max(visc1(type),2*visc2(type)*normal_vit)/norm_vit(type)
3253 visc_eff = hloc*min(coeff_visc_ordre_un*normal_vit,inputs%coeff2*hloc*visc1(type))
3254 nu_loc = nu_loc + visc_eff
3256 viscosity(type,m) = inputs%coeff1*hloc - visc_eff
3259 r_eff = r_eff + (nu_loc + dt**2*hloc)*hloc**2*ray
3260 visc_plot(m) = (nu_loc/6)/(coeff_visc_ordre_un*hloc*normal_vit)
3266 hloc = sqrt(sum(rj(:,m)))
3267 vloc = maxval(vel_tot(j_loc))
3268 cfl = max(vloc*dt/hloc,cfl)
3276 mult(:)= viscosity(:,m)
3280 dw_loc = dw(:,:,l,m)
3283 grad(k,type,l) = sum(vit(j_loc,type)*dw_loc(k,:))
3285 vitloc(type,l) = sum(vit(j_loc,type)*ww(:,l))
3289 ray = sum(mesh%rr(1,j_loc)*ww(:,l))
3293 fs(1) = sum(ff(j_loc,1) * ww(:,l))
3294 ft(1) = sum(v1m(j_loc,1) * ww(:,l))
3295 fp(1) = -sum(p(j_loc,1)*dw_loc(1,:))
3296 fv(1) = ((mode*vitloc(1,l)+vitloc(4,l))*mode +mode*vitloc(4,l)+vitloc(1,l))/ray**2
3298 fs(2) = sum(ff(j_loc,2) * ww(:,l))
3299 ft(2) = sum(v1m(j_loc,2) * ww(:,l))
3300 fp(2) = -sum(p(j_loc,2)*dw_loc(1,:))
3301 fv(2) = ((mode*vitloc(2,l)-vitloc(3,l))*mode -mode*vitloc(3,l)+vitloc(2,l))/ray**2
3303 fs(3) = sum(ff(j_loc,3) * ww(:,l))
3304 ft(3) = sum(v1m(j_loc,3) * ww(:,l))
3305 fp(3) = -sum(p(j_loc,2)*ww(:,l))/ray*mode
3306 fv(3) = (-mode*vitloc(2,l)+vitloc(3,l) +(mode*vitloc(3,l)-vitloc(2,l))*mode)/ray**2
3308 fs(4) = sum(ff(j_loc,4) * ww(:,l))
3309 ft(4) = sum(v1m(j_loc,4) * ww(:,l))
3310 fp(4) = sum(p(j_loc,1)*ww(:,l))/ray *mode
3311 fv(4) = (mode*vitloc(1,l)+vitloc(4,l) +(mode*vitloc(4,l)+vitloc(1,l))*mode)/ray**2
3313 fs(5) = sum(ff(j_loc,5) * ww(:,l))
3314 ft(5) = sum(v1m(j_loc,5) * ww(:,l))
3315 fp(5) = -sum(p(j_loc,1)*dw_loc(2,:))
3316 fv(5) = vitloc(5,l)*(mode/ray)**2
3318 fs(6) = sum(ff(j_loc,6) * ww(:,l))
3319 ft(6) = sum(v1m(j_loc,6) * ww(:,l))
3320 fp(6) = -sum(p(j_loc,2)*dw_loc(2,:))
3321 fv(6) = vitloc(6,l)*(mode/ray)**2
3327 smb = (ft+fp+fs+fv-rotv_v(index,:))*ray*rj(l,m)
3329 grad(:,type,l) = mult(type)*grad(:,type,l)*ray*rj(l,m)
3334 u0(j_loc(ni),j) = u0(j_loc(ni),j) + ww(ni,l)*smb(j) + sum(dw_loc(:,ni)*grad(:,j,l))
3342 WRITE(*,*)
' CFL = ', cfl,
' inputs%LES ', inputs%LES
3343 IF (inputs%LES)
THEN
3344 WRITE(*,*)
' R_eff for mode', mode, normal_vit*6*surf/r_eff
3350 IF (inputs%LES)
THEN
3351 ed_st = normal_vit*coeff_ed_st
3353 dwni_loc = mesh%gauss%dwni(:,:,:,ms)
3354 jji_loc = mesh%jji(:,:,ms)
3355 h2 = -ed_st*h(ms)**2
3357 j_loc(1:n_ws) = mesh%jjsi(1:n_ws,ms)
3360 uloci(:,cotei) = vit(jji_loc(:,cotei),type)
3364 DO ls = 1, mesh%gauss%l_Gs
3366 ray = mesh%rr(1,j_loc(3))
3369 dul = sum(dwni_loc(:,ls,:)*uloci)*mesh%gauss%rji(ls,ms)*h2*ray
3371 DO ni = 1, mesh%gauss%n_w
3372 u0loci(ni, cotei) = u0loci(ni, cotei) + dwni_loc(ni,ls,cotei)*dul
3378 DO ni = 1, mesh%gauss%n_w
3379 u0(jji_loc(ni,cotei),type) = u0(jji_loc(ni,cotei),type) + u0loci(ni,cotei)
3389 SUBROUTINE qs_ns_stab_2010(mesh, pp_mesh, mode,ff,vel_tot,V1m,vit,P,dudt,phalf,nlhalf,dt,u0,rotv_v)
3398 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: ff, rotv_v, dudt, nlhalf
3399 REAL(KIND=8),
DIMENSION(:),
INTENT(IN) :: vel_tot
3400 REAL(KIND=8),
DIMENSION(:,:),
INTENT(OUT):: u0
3401 TYPE(mesh_type
),
TARGET :: mesh, pp_mesh
3402 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: v1m, vit
3403 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: p, phalf
3404 REAL(KIND=8),
INTENT(IN) :: dt
3405 INTEGER,
INTENT(IN) :: mode
3406 INTEGER :: m, l , i , ni , j, index, index2, type, k
3407 REAL(KIND=8),
DIMENSION(6) :: fs , ft , fp, smb, fv, mult
3408 REAL(KIND=8),
DIMENSION(2,6,mesh%gauss%l_G) :: grad
3409 REAL(KIND=8),
DIMENSION(6,mesh%gauss%l_G) :: vitloc
3410 INTEGER,
DIMENSION(mesh%gauss%n_w) :: j_loc
3411 REAL(KIND=8),
DIMENSION(mesh%gauss%k_d,mesh%gauss%n_w) :: dw_loc
3412 REAL(KIND=8),
DIMENSION(mesh%me) :: visc_plot
3417 INTEGER,
DIMENSION(:,:),
POINTER :: jj
3418 INTEGER,
POINTER :: me
3419 REAL(KIND=8),
DIMENSION(6) :: visc1, visc2
3420 REAL(KIND=8) :: ray, div1, div2, hloc, cfl, vloc, normal_vit
3421 REAL(KIND=8),
SAVE :: coeff_ed_st, r_eff, visc_eff, surf, nu_loc
3422 LOGICAL,
SAVE :: once = .true.
3423 REAL(KIND=8),
ALLOCATABLE,
DIMENSION(:),
SAVE :: h
3425 REAL(KIND=8),
DIMENSION(mesh%gauss%n_w,mesh%gauss%l_Gs, 2) :: dwni_loc
3426 INTEGER,
DIMENSION(mesh%gauss%n_w,2) :: jji_loc
3427 REAL(KIND=8),
DIMENSION(mesh%gauss%n_w,2) :: u0loci, uloci
3428 REAL(KIND=8),
DIMENSION(6,mesh%me) :: viscosity
3429 REAL(KIND=8),
DIMENSION(6,pp_mesh%np) :: viscosity_clement
3430 REAL(KIND=8),
DIMENSION(mesh%np,6) :: visc_p2
3431 REAL(KIND=8),
DIMENSION(6) :: norm_vit
3432 REAL(KIND=8) :: dul, h2, type_fe, ed_st
3433 INTEGER :: ms, ls, cotei
3443 IF (mesh%gauss%n_w==3)
THEN
3448 coeff_ed_st = 0.05d0/type_fe
3450 ALLOCATE(h(mesh%mi))
3452 h(ms) = sum(mesh%gauss%rji(:,ms))/type_fe
3459 normal_vit = maxval(vel_tot)
3461 norm_vit(type) = sum(abs(vit(:,type)))/mesh%np + 1.d-14
3471 hloc = sqrt(sum(rj(:,m)))
3473 vloc = maxval(vel_tot(j_loc))
3474 cfl = max(vloc*dt/hloc,cfl)
3479 dw_loc = dw(:,:,l,m)
3482 ray = sum(mesh%rr(1,j_loc)*ww(:,l))
3485 ft(1) = sum(dudt(j_loc,1) * ww(:,l))
3486 fp(1) = sum(phalf(j_loc,1)*dw_loc(1,:))
3488 ft(2) = sum(dudt(j_loc,2) * ww(:,l))
3489 fp(2) = sum(phalf(j_loc,2)*dw_loc(1,:))
3491 ft(3) = sum(dudt(j_loc,3) * ww(:,l))
3492 fp(3) = sum(phalf(j_loc,2)*ww(:,l))*mode/ray
3494 ft(4) = sum(dudt(j_loc,4) * ww(:,l))
3495 fp(4) = -sum(phalf(j_loc,1)*ww(:,l))*mode/ray
3497 ft(5) = sum(dudt(j_loc,5) * ww(:,l))
3498 fp(5) = sum(phalf(j_loc,1)*dw_loc(2,:))
3500 ft(6) = sum(dudt(j_loc,6) * ww(:,l))
3501 fp(6) = sum(phalf(j_loc,2)*dw_loc(2,:))
3504 visc1= max(visc1,abs(ft+fp+nlhalf(index2,:)))
3509 grad(k,type,l) = sum(vit(j_loc,type)*dw_loc(k,:))
3511 vitloc(type,l) = sum(vit(j_loc,type)*ww(:,l))
3515 div1 = abs(grad(1,1,l) + vitloc(1,l)/ray + grad(2,5,l) + vitloc(4,l)*mode/ray)
3516 div2 = abs(grad(1,2,l) + vitloc(2,l)/ray + grad(2,6,l) - vitloc(3,l)*mode/ray)
3517 visc2(1) = div1; visc2(4) = div1; visc2(5) = div1
3518 visc2(2) = div2; visc2(3) = div2; visc2(6) = div2
3524 visc1(type) = max(visc1(type),visc2(type)*normal_vit)/norm_vit(type)
3526 visc_eff = hloc*min(0.25d0*normal_vit,inputs%coeff2*hloc*visc1(type))
3527 nu_loc = nu_loc + visc_eff
3529 viscosity(type,m) = inputs%coeff1*hloc - visc_eff
3536 r_eff = r_eff + (nu_loc + dt**2*hloc)*hloc**2
3537 visc_plot(m) = (nu_loc/6)/(0.25*hloc*normal_vit)
3542 CALL
clement_c(pp_mesh,viscosity(type,:),viscosity_clement(type,:))
3543 CALL
inject_clement(pp_mesh%jj, mesh%jj, viscosity_clement(type,:), visc_p2(:,type))
3557 dw_loc = dw(:,:,l,m)
3559 mult(type)= sum(visc_p2(j_loc(:),type)*ww(:,l))
3561 grad(k,type,l) = sum(vit(j_loc,type)*dw_loc(k,:))
3563 vitloc(type,l) = sum(vit(j_loc,type)*ww(:,l))
3567 ray = sum(mesh%rr(1,j_loc)*ww(:,l))
3571 fs(1) = sum(ff(j_loc,1) * ww(:,l))
3572 ft(1) = sum(v1m(j_loc,1) * ww(:,l))
3573 fp(1) = -sum(p(j_loc,1)*dw_loc(1,:))
3574 fv(1) = ((mode*vitloc(1,l)+vitloc(4,l))*mode +mode*vitloc(4,l)+vitloc(1,l))/ray**2
3576 fs(2) = sum(ff(j_loc,2) * ww(:,l))
3577 ft(2) = sum(v1m(j_loc,2) * ww(:,l))
3578 fp(2) = -sum(p(j_loc,2)*dw_loc(1,:))
3579 fv(2) = ((mode*vitloc(2,l)-vitloc(3,l))*mode -mode*vitloc(3,l)+vitloc(2,l))/ray**2
3581 fs(3) = sum(ff(j_loc,3) * ww(:,l))
3582 ft(3) = sum(v1m(j_loc,3) * ww(:,l))
3583 fp(3) = -sum(p(j_loc,2)*ww(:,l))/ray*mode
3584 fv(3) = (-mode*vitloc(2,l)+vitloc(3,l) +(mode*vitloc(3,l)-vitloc(2,l))*mode)/ray**2
3586 fs(4) = sum(ff(j_loc,4) * ww(:,l))
3587 ft(4) = sum(v1m(j_loc,4) * ww(:,l))
3588 fp(4) = sum(p(j_loc,1)*ww(:,l))/ray *mode
3589 fv(4) = (mode*vitloc(1,l)+vitloc(4,l) +(mode*vitloc(4,l)+vitloc(1,l))*mode)/ray**2
3591 fs(5) = sum(ff(j_loc,5) * ww(:,l))
3592 ft(5) = sum(v1m(j_loc,5) * ww(:,l))
3593 fp(5) = -sum(p(j_loc,1)*dw_loc(2,:))
3594 fv(5) = vitloc(5,l)*(mode/ray)**2
3596 fs(6) = sum(ff(j_loc,6) * ww(:,l))
3597 ft(6) = sum(v1m(j_loc,6) * ww(:,l))
3598 fp(6) = -sum(p(j_loc,2)*dw_loc(2,:))
3599 fv(6) = vitloc(6,l)*(mode/ray)**2
3606 smb = (ft+fp+fs+fv-rotv_v(index,:))*ray*rj(l,m)
3608 grad(:,type,l) = mult(type)*grad(:,type,l)*ray*rj(l,m)
3613 u0(j_loc(ni),j) = u0(j_loc(ni),j) + ww(ni,l)*smb(j) + sum(dw_loc(:,ni)*grad(:,j,l))
3620 WRITE(*,*)
' CFL = ', cfl,
' R_eff for mode', mode, 6*surf/r_eff
3625 ed_st = normal_vit*coeff_ed_st
3627 dwni_loc = mesh%gauss%dwni(:,:,:,ms)
3628 jji_loc = mesh%jji(:,:,ms)
3629 h2 = -ed_st*h(ms)**2
3631 j_loc(1:n_ws) = mesh%jjsi(1:n_ws,ms)
3634 uloci(:,cotei) = vit(jji_loc(:,cotei),type)
3638 DO ls = 1, mesh%gauss%l_Gs
3641 ray = mesh%rr(1,j_loc(3))
3644 dul = sum(dwni_loc(:,ls,:)*uloci)*mesh%gauss%rji(ls,ms)*h2*ray
3646 DO ni = 1, mesh%gauss%n_w
3647 u0loci(ni, cotei) = u0loci(ni, cotei) + dwni_loc(ni,ls,cotei)*dul
3653 DO ni = 1, mesh%gauss%n_w
3654 u0(jji_loc(ni,cotei),type) = u0(jji_loc(ni,cotei),type) + u0loci(ni,cotei)
3672 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: ff, rotv_v
3673 REAL(KIND=8),
DIMENSION(:),
INTENT(IN) :: vel_tot
3674 REAL(KIND=8),
DIMENSION(:,:),
INTENT(OUT):: u0
3675 TYPE(mesh_type
),
TARGET :: mesh
3676 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: v1m, vit
3677 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: p
3678 REAL(KIND=8),
INTENT(IN) :: dt
3679 INTEGER,
INTENT(IN) :: mode
3680 INTEGER :: m, l , i , ni , j, index, type, k
3681 REAL(KIND=8),
DIMENSION(6) :: fs , ft , fp, smb, fv, mult
3682 REAL(KIND=8),
DIMENSION(2,6,mesh%gauss%l_G) :: grad
3683 REAL(KIND=8),
DIMENSION(6,mesh%gauss%l_G) :: vitloc
3684 INTEGER,
DIMENSION(mesh%gauss%n_w) :: j_loc
3685 REAL(KIND=8),
DIMENSION(mesh%gauss%k_d,mesh%gauss%n_w) :: dw_loc
3691 INTEGER,
DIMENSION(:,:),
POINTER :: jj
3692 INTEGER,
POINTER :: me
3693 REAL(KIND=8) :: ray, div1, div2, vit1=0.005d0, vit2=.075d0, visc1, visc2, hloc, cfl, vloc
3694 REAL(KIND=8),
SAVE :: r_eff, visc_eff
3695 LOGICAL,
SAVE :: once = .true.
3711 hloc = sqrt(sum(rj(:,m)))
3712 vloc = min(maxval(vel_tot(j_loc)),1.d0)
3714 cfl = max(vloc*dt/hloc,cfl)
3718 dw_loc = dw(:,:,l,m)
3721 ray = sum(mesh%rr(1,j_loc)*ww(:,l))
3726 grad(k,type,l) = sum(vit(j_loc,type)*dw_loc(k,:))
3728 vitloc(type,l) = sum(vit(j_loc,type)*ww(:,l))
3732 div1 = grad(1,1,l) + vitloc(1,l)/ray + grad(2,5,l) + vitloc(4,l)*mode/ray
3733 div2 = grad(1,2,l) + vitloc(2,l)/ray + grad(2,6,l) - vitloc(3,l)*mode/ray
3734 visc1 = max(visc1,abs(div1))
3735 visc2 = max(visc2,abs(div2))
3737 visc1 = max(visc1,visc2)
3740 visc_eff = hloc*min(inputs%coeff1,inputs%coeff2*visc1)
3747 r_eff = r_eff + visc_eff*sum(rj(:,m))
3752 dw_loc = dw(:,:,l,m)
3755 ray = sum(mesh%rr(1,j_loc)*ww(:,l))
3760 fs(1) = sum(ff(j_loc,1) * ww(:,l))
3761 ft(1) = sum(v1m(j_loc,1) * ww(:,l))
3762 fp(1) = -sum(p(j_loc,1)*dw_loc(1,:))
3763 fv(1) = ((mode*vitloc(1,l)+vitloc(4,l))*mode +mode*vitloc(4,l)+vitloc(1,l))/ray**2
3766 fs(2) = sum(ff(j_loc,2) * ww(:,l))
3767 ft(2) = sum(v1m(j_loc,2) * ww(:,l))
3768 fp(2) = -sum(p(j_loc,2)*dw_loc(1,:))
3769 fv(2) = ((mode*vitloc(2,l)-vitloc(3,l))*mode -mode*vitloc(3,l)+vitloc(2,l))/ray**2
3772 fs(3) = sum(ff(j_loc,3) * ww(:,l))
3773 ft(3) = sum(v1m(j_loc,3) * ww(:,l))
3774 fp(3) = -sum(p(j_loc,2)*ww(:,l))/ray*mode
3775 fv(3) = (-mode*vitloc(2,l)+vitloc(3,l) +(mode*vitloc(3,l)-vitloc(2,l))*mode)/ray**2
3778 fs(4) = sum(ff(j_loc,4) * ww(:,l))
3779 ft(4) = sum(v1m(j_loc,4) * ww(:,l))
3780 fp(4) = sum(p(j_loc,1)*ww(:,l))/ray *mode
3781 fv(4) = (mode*vitloc(1,l)+vitloc(4,l) +(mode*vitloc(4,l)+vitloc(1,l))*mode)/ray**2
3784 fs(5) = sum(ff(j_loc,5) * ww(:,l))
3785 ft(5) = sum(v1m(j_loc,5) * ww(:,l))
3786 fp(5) = -sum(p(j_loc,1)*dw_loc(2,:))
3787 fv(5) = vitloc(5,l)*(mode/ray)**2
3790 fs(6) = sum(ff(j_loc,6) * ww(:,l))
3791 ft(6) = sum(v1m(j_loc,6) * ww(:,l))
3792 fp(6) = -sum(p(j_loc,2)*dw_loc(2,:))
3793 fv(6) = vitloc(6,l)*(mode/ray)**2
3797 mult(1) = visc1; mult(4) = visc1; mult(5) = visc1
3798 mult(2) = visc2; mult(3) = visc2; mult(6) = visc2
3801 smb = (ft+fp+fs+fv-rotv_v(index,:))*ray*rj(l,m)
3803 grad(:,type,l) = mult(type)*grad(:,type,l)*ray*rj(l,m)
3808 u0(jj(ni,m),j) = u0(jj(ni,m),j) + ww(ni,l)*smb(j) + sum(dw_loc(:,ni)*grad(:,j,l))
3815 WRITE(*,*)
' CFL = ', cfl,
' R_eff ', 1/r_eff
3832 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: ff
3833 REAL(KIND=8),
DIMENSION(:,:),
INTENT(OUT) :: u0
3834 TYPE(mesh_type
),
TARGET :: mesh
3835 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: v1m
3836 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: p
3837 REAL(KIND=8),
INTENT(IN) :: dt
3838 INTEGER ,
INTENT(IN) :: mode
3840 INTEGER :: m, l , i , ni , j
3841 REAL(KIND=8),
DIMENSION(6) :: fs , ft , fp, fnl, smb
3843 INTEGER,
DIMENSION(mesh%gauss%n_w) :: j_loc
3844 REAL(KIND=8),
DIMENSION(mesh%gauss%k_d,mesh%gauss%n_w) :: dw_loc
3845 LOGICAL,
SAVE :: once = .true.
3846 LOGICAL,
SAVE :: once_sft = .true.
3849 REAL(KIND=8) :: tps, dummy, tt
3855 INTEGER,
DIMENSION(:,:),
POINTER :: jj
3856 INTEGER,
POINTER :: me
3875 dw_loc = dw(:,:,l,m)
3878 ray = sum(mesh%rr(1,j_loc)*ww(:,l))
3883 fs(1) = sum(ff(j_loc,1) * ww(:,l))
3884 ft(1) = 1/(2*dt) * sum(v1m(j_loc,1) * ww(:,l))
3885 fp(1) = -sum(p(j_loc,1)*dw_loc(1,:))
3889 fs(2) = sum(ff(j_loc,2) * ww(:,l))
3890 ft(2) = 1/(2*dt) * sum(v1m(j_loc,2) * ww(:,l))
3891 fp(2) = -sum(p(j_loc,2)*dw_loc(1,:))
3895 fs(3) = sum(ff(j_loc,3) * ww(:,l))
3896 ft(3) = 1/(2*dt) * sum(v1m(j_loc,3) * ww(:,l))
3897 fp(3) = -sum(p(j_loc,2)*ww(:,l))/ray*mode
3901 fs(4) = sum(ff(j_loc,4) * ww(:,l))
3902 ft(4) = 1/(2*dt) * sum(v1m(j_loc,4) * ww(:,l))
3903 fp(4) = sum(p(j_loc,1)*ww(:,l))/ray *mode
3907 fs(5) = sum(ff(j_loc,5) * ww(:,l))
3908 ft(5) = 1/(2*dt) * sum(v1m(j_loc,5) * ww(:,l))
3909 fp(5) = -sum(p(j_loc,1)*dw_loc(2,:))
3913 fs(6) = sum(ff(j_loc,6) * ww(:,l))
3914 ft(6) = 1/(2*dt) * sum(v1m(j_loc,6) * ww(:,l))
3915 fp(6) = -sum(p(j_loc,2)*dw_loc(2,:))
3920 smb = (ft+fp+fs)*ray*rj(l,m)
3926 u0(jj(ni,m),j) = u0(jj(ni,m),j) + ww(ni,l)*smb(j)
3945 TYPE(mesh_type
),
TARGET :: mesh
3946 INTEGER,
INTENT(IN) :: mode
3947 REAL(KIND=8),
DIMENSION(mesh%np,6),
INTENT(IN) :: v
3948 REAL(KIND=8),
DIMENSION(:,:),
INTENT(OUT) :: u0
3950 INTEGER :: m, l , i , ni , j
3951 REAL(KIND=8),
DIMENSION(6) :: frot
3954 REAL(KIND=8) :: tps, dummy, tt
3955 INTEGER,
DIMENSION(:,:),
POINTER :: jj
3956 INTEGER,
POINTER :: me
3972 DO ni = 1, n_w; i = jj(ni,m)
3973 ray = ray + mesh%rr(1,i)*ww(ni,l)
3995 frot(1) = mode/ray*sum(v(jj(:,m),6)*ww(:,l)) &
3996 -sum(v(jj(:,m),3)*dw(2,:,l,m))
3997 frot(3) = sum(v(jj(:,m),1)*dw(2,:,l,m)) &
3998 -sum(v(jj(:,m),5)*dw(1,:,l,m))
3999 frot(5) = 1/ray*sum(v(jj(:,m),3)*ww(:,l))+sum(v(jj(:,m),3) &
4000 *dw(1,:,l,m))-mode/ray*sum(v(jj(:,m),2)*ww(:,l))
4003 frot(2) = -mode/ray*sum(v(jj(:,m),5)*ww(:,l)) &
4004 -sum(v(jj(:,m),4)*dw(2,:,l,m))
4005 frot(4) = sum(v(jj(:,m),2)*dw(2,:,l,m)) &
4006 -sum(v(jj(:,m),6)*dw(1,:,l,m))
4007 frot(6) = 1/ray*sum(v(jj(:,m),4)*ww(:,l))+sum(v(jj(:,m),4) &
4008 *dw(1,:,l,m))+mode/ray*sum(v(jj(:,m),1)*ww(:,l))
4010 frot = frot * rj(l,m)
4016 u0(jj(ni,m),j) = u0(jj(ni,m),j) + ww(ni,l) * ray * frot(j)
4030 REAL(KIND=8),
DIMENSION(:),
INTENT(OUT) :: visc
4031 TYPE(mesh_type
) :: mesh
4032 REAL(KIND=8),
DIMENSION(mesh%me) :: vint
4033 REAL(KIND=8),
ALLOCATABLE,
DIMENSION(:),
SAVE :: surf, sloc
4034 LOGICAL,
SAVE :: once=.true.
4035 INTEGER :: m, n, mp, it
4039 ALLOCATE(surf(mesh%me),sloc(mesh%me))
4041 surf(m) = sum(mesh%gauss%rj(:,m))
4046 IF (mesh%neigh(n,m)==0) cycle
4047 mp = mesh%neigh(n,m)
4048 sloc(m) = sloc(m) + surf(m)+surf(mp)
4056 IF (mesh%neigh(n,m)==0) cycle
4057 mp = mesh%neigh(n,m)
4058 vint(m) = vint(m)+ visc(m)*surf(m) + visc(mp)*surf(mp)
4060 vint(m) = vint(m)/sloc(m)
4071 REAL(KIND=8),
DIMENSION(:) :: phi, phi_s
4073 INTEGER :: m, l, k, i
4074 REAL(KIND=8),
DIMENSION(mesh%np) :: a0, b1, b2
4075 REAL(KIND=8),
DIMENSION(:),
ALLOCATABLE,
SAVE :: s, a11, a22, a12, oneoverdet
4076 REAL(KIND=8),
DIMENSION(:,:),
ALLOCATABLE,
SAVE :: rg
4077 INTEGER,
DIMENSION(mesh%gauss%n_w) :: j_loc
4078 REAL(KIND=8),
DIMENSION(mesh%gauss%n_w) :: phi_loc
4079 REAL(KIND=8),
DIMENSION(mesh%gauss%k_d) :: r_loc
4080 REAL(KIND=8) :: x1l, x2l, phil, a1, a2, det
4082 LOGICAL,
SAVE :: once=.true.
4084 IF (mesh%gauss%n_w/=3)
THEN
4085 WRITE(*,*)
' BUG CLEMENT: Only P1 programmed yet'
4090 IF (
SIZE(phi)==mesh%me) const=.true.
4094 ALLOCATE(rg(mesh%gauss%k_d,mesh%np), s(mesh%np))
4098 j_loc = mesh%jj(:,m)
4099 s(j_loc(:)) = s(j_loc(:)) + sum(mesh%gauss%rj(:,m))
4100 DO k = 1, mesh%gauss%k_d
4102 DO l = 1, mesh%gauss%l_G
4103 r_loc(k) = r_loc(k) + sum(mesh%rr(k,j_loc(:))*mesh%gauss%ww(:,l))*mesh%gauss%rj(l,m)
4105 rg(k,j_loc(:)) = rg(k,j_loc(:)) + r_loc(k)
4109 rg(:,i) = rg(:,i)/s(i)
4111 ALLOCATE(a11(mesh%np), a22(mesh%np), a12(mesh%np), oneoverdet(mesh%np))
4116 j_loc = mesh%jj(:,m)
4117 DO l = 1, mesh%gauss%l_G
4118 x1l = sum(mesh%rr(1,j_loc(:))*mesh%gauss%ww(:,l))
4119 x2l = sum(mesh%rr(2,j_loc(:))*mesh%gauss%ww(:,l))
4120 a11(j_loc(:)) = a11(j_loc(:)) + (x1l-rg(1,j_loc(:)))**2*mesh%gauss%rj(l,m)
4121 a22(j_loc(:)) = a22(j_loc(:)) + (x2l-rg(2,j_loc(:)))**2*mesh%gauss%rj(l,m)
4122 a12(j_loc(:)) = a12(j_loc(:)) + (x1l-rg(1,j_loc(:)))*(x2l-rg(2,j_loc(:)))*mesh%gauss%rj(l,m)
4125 oneoverdet = 1.d0/(a11*a22-a12*a12)
4133 j_loc = mesh%jj(:,m)
4139 DO l = 1, mesh%gauss%l_G
4140 x1l = sum(mesh%rr(1,j_loc(:))*mesh%gauss%ww(:,l))
4141 x2l = sum(mesh%rr(2,j_loc(:))*mesh%gauss%ww(:,l))
4142 phil = sum(phi_loc*mesh%gauss%ww(:,l))
4143 b1(j_loc(:)) = b1(j_loc(:)) + (x1l-rg(1,j_loc(:)))*phil*mesh%gauss%rj(l,m)
4144 b2(j_loc(:)) = b2(j_loc(:)) + (x2l-rg(2,j_loc(:)))*phil*mesh%gauss%rj(l,m)
4145 a0(j_loc) = a0(j_loc) + phil*mesh%gauss%rj(l,m)
4152 a1 =(b1(i)*a22(i)-b2(i)*a12(i))*oneoverdet(i)
4153 a2 =(b2(i)*a11(i)-b1(i)*a12(i))*oneoverdet(i)
4154 phi_s(i) = a0(i)+a1*(mesh%rr(1,i)-rg(1,i))+a2*(mesh%rr(2,i)-rg(2,i))
4164 INTEGER,
DIMENSION(:,:),
INTENT(IN) :: jj_c, jj_f
4165 REAL(KIND=8),
DIMENSION(:),
INTENT(IN) :: pp_c
4166 REAL(KIND=8),
DIMENSION(:),
INTENT(OUT) :: pp_f
4167 REAL(KIND=8) :: half = 0.5
4170 DO m = 1,
SIZE(jj_f,2)
4171 pp_f(jj_f(1:3,m)) = pp_c(jj_c(:,m))
4172 pp_f(jj_f(4,m)) = (pp_c(jj_c(2,m)) + pp_c(jj_c(3,m)))*half
4173 pp_f(jj_f(5,m)) = (pp_c(jj_c(3,m)) + pp_c(jj_c(1,m)))*half
4174 pp_f(jj_f(6,m)) = (pp_c(jj_c(1,m)) + pp_c(jj_c(2,m)))*half
subroutine qs_ns_stab_2008(mesh, mode, ff, vel_tot, V1m, vit, P, dt, u0, rotv_v)
subroutine qs_00_stokes_3d(mesh, alpha, mode, ff, V1, V2, P, dt, u0)
subroutine qs_ns_stab_new(mesh, mode, ff, vel_tot, V1m, vit, P, dudt, phalf, nlhalf, dt, u0, rotv_v)
subroutine qs_00_adv_diff_vect_3d(mesh, alpha, mode, gg, ff, V1, V2, dt, u0)
subroutine qs_00_inst_init_3d(mesh, ff, T1, dt, u0)
subroutine qs_01_rot(mesh, mode, V, u0)
subroutine qs_01_div_hybrid(uu_mesh, pp_mesh, mode, gg, u0_c)
subroutine qs_00_inst_vect_3d_init_ssr(mesh, alpha, mode, ff, V, dt, u0)
subroutine qs_00_stokes_3d_new_ssr(mesh, alpha, mode, ff, V1, V2, P, dt, u0)
subroutine qs_00_inst_vect_3d_ssr(mesh, alpha, mode, ff, V1, V2, dt, u0)
subroutine qs_01_grad(mesh, mode, pp, u0)
subroutine qs_00_inst_vect_3d_init(mesh, alpha, mode, ff, V2, dt, u0)
subroutine qs_ns_stab_2010(mesh, pp_mesh, mode, ff, vel_tot, V1m, vit, P, dudt, phalf, nlhalf, dt, u0, rotv_v)
subroutine qs_00_adv_diff_vect_3d_init(mesh, alpha, mode, gg, ff, V, dt, u0)
subroutine qs_00_inst_init_3d_ssr(mesh, ff, T1, dt, u0)
subroutine qs_navier_stokes_2006(mesh, mode, ff, V1m, P, dt, u0, rotv_v)
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 f with f(r,\theta, z)\f $the cylindrical coordinates
subroutine qs_00_ssr(mesh, ff, u0)
subroutine moy(mesh, p, RESULT)
subroutine qs_00_adv_diff_init_3d(eps, mode, mesh, ff, T1, dt, gg, u0)
subroutine qs_01_grad_gl(mesh, mod_max, pp, u0)
subroutine qs_00_adv_diff_vect_3d_init_ssr(mesh, alpha, mode, gg, ff, V, dt, u0)
subroutine inject_clement(jj_c, jj_f, pp_c, pp_f)
subroutine qs_00_adv_diff_3d_ssr(eps, mode, mesh, ff, T1, T2, gg, dt, u0)
subroutine qs_00_lap_init(mesh, alpha, mode, ff, V, dt, u0)
subroutine qs_stokes_2006(mesh, mode, ff, V1m, P, dt, u0)
subroutine clement_c(mesh, phi, phi_s)
subroutine qs_ns_2006(mesh, mode, ff, V1m, P, dt, u0, rotv_v)
subroutine qs_00_adv_diff_3d(eps, mode, mesh, ff, T1, T2, gg, dt, u0)
subroutine qs_00_stokes_3d_init(mesh, alpha, mode, ff, P, V, dt, u0)
subroutine qs_00_stokes_3d_new(mesh, alpha, mode, ff, V1, V2, P, dt, u0)
subroutine qs_00(mesh, ff, u0)
subroutine qs_stokes(mesh, mode, ff, V1m, V2m, P, dt, u0)
subroutine qs_01_div_hybrid_2006(uu_mesh, pp_mesh, mode, gg, u0_c)
subroutine qs_00_div(mesh, mode, V, u0)
subroutine qs_00_lap(mesh, alpha, mode, ff, V2, V1, dt, u0)
real(kind=8) function user_time()
subroutine qs_00_rot(mesh, mode, V, Rot)
subroutine qs_00_adv_diff_vect_3d_ssr(mesh, alpha, mode, gg, ff, V1, V2, dt, u0)
subroutine qs_00_ns_inline(mesh, mode, m_max, ff, V1, V2, P, dt, u0, meth, tps_sft)
subroutine qs_00_inst_3d(mesh, ff, T1, T2, dt, u0)
subroutine average(mesh, visc)
subroutine qs_00_inst_vect_3d(mesh, alpha, mode, ff, V1, V2, dt, u0)
subroutine dirichlet(jjs, us, ff)
subroutine qs_01_div_hybrid_old(uu_mesh, pp_mesh, mode, gg, u0_c)
subroutine qs_00_adv_diff_init_3d_ssr(eps, mode, mesh, ff, T1, dt, gg, u0)
subroutine qs_00_inst_3d_ssr(mesh, ff, T1, T2, dt, u0)