9 INTEGER,
DIMENSION(:),
ALLOCATABLE,
PUBLIC :: vv_mz_LA, H_mz_LA
10 LOGICAL,
PRIVATE :: need_sym=.FALSE.
13 #include "petsc/finclude/petsc.h"
23 INTEGER,
DIMENSION(:) :: ltg_la
24 INTEGER,
DIMENSION(mesh_loc%np) :: i_d
25 INTEGER :: i, j, m, dom, n_w
26 REAL(KIND=8) :: xx, zz, epsilon
31 n_w =
SIZE(mesh_loc%jj,1)
34 i_d(mesh_loc%jj(:,m)) = mesh_loc%i_d(m)
39 zz = -mesh_loc%rr(2,i)
41 DO m = 1, mesh_glob%me
42 IF (mesh_glob%i_d(m) /= dom) cycle
44 IF (abs(xx-mesh_glob%rr(1,mesh_glob%jj(j,m)))+abs(zz-mesh_glob%rr(2,mesh_glob%jj(j,m))) .LT. epsilon)
THEN
45 ltg_la(i) = mesh_glob%jj(j,m)
52 IF (minval(ltg_la)<0 .AND. (mesh_loc%me/=0))
THEN
53 WRITE(*,*)
'bug in symmetric_points'
58 SUBROUTINE symm_champ(communicator, vv_in, mesh, vv_out, if_u_h)
65 REAL(KIND=8),
DIMENSION(:,:,:),
INTENT(IN) :: vv_in
66 REAL(KIND=8),
DIMENSION(:,:,:),
INTENT(OUT) :: vv_out
67 CHARACTER(len=1),
INTENT(IN) :: if_u_h
68 INTEGER,
POINTER,
DIMENSION(:) :: ifrom
71 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: ix
73 LOGICAL,
SAVE :: once_u=.true.
74 LOGICAL,
SAVE :: once_h=.true.
75 #include "petsc/finclude/petscvec.h90"
76 vec,
SAVE :: vv_mz, vv_mz_ghost, h_mz, h_mz_ghost
77 petscerrorcode :: ierr
78 mpi_comm :: communicator
79 petscscalar,
POINTER :: x_loc(:)
82 IF (mesh%me ==0)
RETURN
85 IF (.NOT. need_sym)
RETURN
87 IF ((once_u) .AND. (if_u_h==
'u'))
THEN
88 ALLOCATE(la_u%loc_to_glob(1,mesh%np))
89 la_u%loc_to_glob(1,:) = mesh%loc_to_glob
93 CALL veccreateghost(communicator, n, &
94 petsc_determine,
SIZE(ifrom), ifrom, vv_mz, ierr)
95 CALL vecghostgetlocalform(vv_mz, vv_mz_ghost, ierr)
99 IF ((once_h) .AND. (if_u_h==
'h'))
THEN
100 ALLOCATE(la_hh%loc_to_glob(1,mesh%np))
101 la_hh%loc_to_glob(1,:) = mesh%loc_to_glob
105 CALL veccreateghost(communicator, n, &
106 petsc_determine,
SIZE(ifrom), ifrom, h_mz, ierr)
107 CALL vecghostgetlocalform(h_mz, h_mz_ghost, ierr)
112 ALLOCATE(ix(mesh%np))
113 IF (if_u_h ==
'u')
THEN
115 DO i = 1,
SIZE(vv_in,2)
116 DO j = 1,
SIZE(vv_in,3)
117 CALL veczeroentries(vv_mz, ierr)
118 CALL vecsetvalues(vv_mz, mesh%np, ix, vv_in(:,i,j), insert_values, ierr)
119 CALL vecassemblybegin(vv_mz, ierr)
120 CALL vecassemblyend(vv_mz, ierr)
121 CALL vecghostupdatebegin(vv_mz,insert_values, scatter_forward,ierr)
122 CALL vecghostupdateend(vv_mz,insert_values, scatter_forward, ierr)
123 CALL vecgetarrayf90(vv_mz_ghost, x_loc, ierr)
124 vv_out(:,i,j) = x_loc(:)
125 CALL vecrestorearrayf90(vv_mz_ghost, x_loc, ierr)
129 ELSE IF (if_u_h==
'h')
THEN
131 DO i = 1,
SIZE(vv_in,2)
132 DO j = 1,
SIZE(vv_in,3)
133 CALL veczeroentries(h_mz, ierr)
134 CALL vecsetvalues(h_mz, mesh%np, ix, vv_in(:,i,j), insert_values, ierr)
135 CALL vecassemblybegin(h_mz, ierr)
136 CALL vecassemblyend(h_mz, ierr)
137 CALL vecghostupdatebegin(h_mz,insert_values, scatter_forward,ierr)
138 CALL vecghostupdateend(h_mz,insert_values, scatter_forward, ierr)
139 CALL vecgetarrayf90(h_mz_ghost, x_loc, ierr)
140 vv_out(:,i,j) = x_loc(:)
141 CALL vecrestorearrayf90(h_mz_ghost, x_loc, ierr)
158 INTEGER,
DIMENSION(:),
INTENT(IN) :: list_mode
159 REAL(KIND=8),
DIMENSION(:,:,:) ,
INTENT(IN) :: v
160 REAL(KIND=8),
DIMENSION(mesh%np, SIZE(v,2),SIZE(list_mode)) :: vv
161 REAL(KIND=8),
DIMENSION(SIZE(list_mode)) ,
INTENT(OUT):: e_mode, e_mode_sym, e_mode_anti
162 CHARACTER(len=1),
INTENT(IN) :: if_u_h
163 REAL(KIND=8),
DIMENSION(3) :: type_sym
165 REAL(KIND=8),
DIMENSION(mesh%np,size(v,2),size(list_mode)) :: champ_anti, champ_sym
167 INTEGER,
DIMENSION(:,:),
POINTER :: jj
168 INTEGER,
POINTER :: me
171 #include "petsc/finclude/petsc.h"
172 mpi_comm :: communicator
178 m_max_c =
size(list_mode)
183 CALL
symm_champ(communicator, v, mesh, vv, if_u_h)
186 DO i=1,
size(list_mode)
187 IF (mod(list_mode(i),2) == 0)
THEN
196 champ_anti(:,1:2,i) = 0.5d0*(v(:,1:2,i) - type_sym(1)*vv(:,1:2,i))
197 champ_anti(:,3:4,i) = 0.5d0*(v(:,3:4,i) - type_sym(2)*vv(:,3:4,i))
198 champ_anti(:,5:6,i) = 0.5d0*(v(:,5:6,i) - type_sym(3)*vv(:,5:6,i))
199 champ_sym(:,1:2,i) = 0.5d0*(v(:,1:2,i) + type_sym(1)*vv(:,1:2,i))
200 champ_sym(:,3:4,i) = 0.5d0*(v(:,3:4,i) + type_sym(2)*vv(:,3:4,i))
201 champ_sym(:,5:6,i) = 0.5d0*(v(:,5:6,i) + type_sym(3)*vv(:,5:6,i))
203 write(*,*)
'===Compute anti-symmetric field'
207 e_mode(i) = 0.5*(
norme_l2_champ_par(communicator, mesh, list_mode(i:i), v(:,:,i:i)))**2
208 e_mode_anti(i) = 0.5*(
norme_l2_champ_par(communicator, mesh, list_mode(i:i), champ_anti(:,:,i:i)))**2
209 e_mode_sym(i) = 0.5*(
norme_l2_champ_par(communicator, mesh, list_mode(i:i), champ_sym(:,:,i:i)))**2
216 SUBROUTINE val_ener_sym_glob(communicator, mesh, list_mode, v, e_mode, e_mode_sym, e_mode_anti, type_sym, if_u_h)
226 INTEGER,
DIMENSION(:),
INTENT(IN) :: list_mode
227 REAL(KIND=8),
DIMENSION(:,:,:) ,
INTENT(IN) :: v
228 REAL(KIND=8),
DIMENSION(mesh%np, SIZE(v,2),SIZE(list_mode)) :: vv
229 REAL(KIND=8),
DIMENSION(SIZE(list_mode)) ,
INTENT(OUT):: e_mode, e_mode_sym, e_mode_anti
230 REAL(KIND=8),
DIMENSION(3),
INTENT(IN) :: type_sym
231 CHARACTER(len=1),
INTENT(IN) :: if_u_h
233 REAL(KIND=8),
DIMENSION(mesh%np,size(v,2),size(list_mode)) :: champ_anti, champ_sym
235 INTEGER,
DIMENSION(:,:),
POINTER :: jj
236 INTEGER,
POINTER :: me
239 #include "petsc/finclude/petsc.h"
240 mpi_comm :: communicator
246 m_max_c =
size(list_mode)
251 CALL
symm_champ(communicator, v, mesh, vv, if_u_h)
254 champ_anti(:,1:2,:) = 0.5d0*(v(:,1:2,:) - type_sym(1)*vv(:,1:2,:))
255 champ_anti(:,3:4,:) = 0.5d0*(v(:,3:4,:) - type_sym(2)*vv(:,3:4,:))
256 champ_anti(:,5:6,:) = 0.5d0*(v(:,5:6,:) - type_sym(3)*vv(:,5:6,:))
257 champ_sym(:,1:2,:) = 0.5d0*(v(:,1:2,:) + type_sym(1)*vv(:,1:2,:))
258 champ_sym(:,3:4,:) = 0.5d0*(v(:,3:4,:) + type_sym(2)*vv(:,3:4,:))
259 champ_sym(:,5:6,:) = 0.5d0*(v(:,5:6,:) + type_sym(3)*vv(:,5:6,:))
260 write(*,*)
'===Compute anti-symmetric field'
264 e_mode(i) = 0.5*(
norme_l2_champ_par(communicator, mesh, list_mode(i:i), v(:,:,i:i)))**2
265 e_mode_anti(i) = 0.5*(
norme_l2_champ_par(communicator, mesh, list_mode(i:i), champ_anti(:,:,i:i)))**2
266 e_mode_sym(i) = 0.5*(
norme_l2_champ_par(communicator, mesh, list_mode(i:i), champ_sym(:,:,i:i)))**2
272 SUBROUTINE val_ener_sym_rpi(communicator, mesh, list_mode, v, e_mode, e_mode_sym, e_mode_anti, type_sym, if_u_h)
287 INTEGER,
DIMENSION(:),
INTENT(IN) :: list_mode
288 REAL(KIND=8),
DIMENSION(:,:,:) ,
INTENT(IN) :: v
289 REAL(KIND=8),
DIMENSION(mesh%np, SIZE(v,2),SIZE(list_mode)) :: vv
290 REAL(KIND=8),
DIMENSION(SIZE(list_mode)) ,
INTENT(OUT):: e_mode, e_mode_sym, e_mode_anti
291 REAL(KIND=8),
DIMENSION(6),
INTENT(IN) :: type_sym
292 CHARACTER(len=1),
INTENT(IN) :: if_u_h
294 REAL(KIND=8),
DIMENSION(mesh%np,size(v,2),size(list_mode)) :: champ_anti, champ_sym
296 INTEGER,
DIMENSION(:,:),
POINTER :: jj
297 INTEGER,
POINTER :: me
300 #include "petsc/finclude/petsc.h"
301 mpi_comm :: communicator
307 m_max_c =
size(list_mode)
312 CALL
symm_champ(communicator, v, mesh, vv, if_u_h)
316 champ_anti(:,k,:) = 0.5d0*(v(:,k,:) - type_sym(k)*vv(:,k,:))
317 champ_sym(:,k,:) = 0.5d0*(v(:,k,:) + type_sym(k)*vv(:,k,:))
319 write(*,*)
'===Compute anti-symmetric field'
323 e_mode(i) = 0.5*(
norme_l2_champ_par(communicator, mesh, list_mode(i:i), v(:,:,i:i)))**2
324 e_mode_anti(i) = 0.5*(
norme_l2_champ_par(communicator, mesh, list_mode(i:i), champ_anti(:,:,i:i)))**2
325 e_mode_sym(i) = 0.5*(
norme_l2_champ_par(communicator, mesh, list_mode(i:i), champ_sym(:,:,i:i)))**2
339 INTEGER,
DIMENSION(:),
INTENT(IN) :: list_mode
340 REAL(KIND=8),
DIMENSION(:,:,:) ,
INTENT(IN) :: v
341 REAL(KIND=8),
DIMENSION(mesh%np, SIZE(v,2),SIZE(list_mode)) :: vv_n, vv_s
342 REAL(KIND=8),
DIMENSION(SIZE(list_mode)) ,
INTENT(OUT):: e_north, e_south
343 REAL(KIND=8),
DIMENSION(SIZE(list_mode)),
OPTIONAL :: e_tot
345 REAL(KIND=8) :: ray, pi
346 REAL(KIND=8),
DIMENSION(2) :: e_ns
347 REAL(KIND=8) :: epsilon = 1.d-10
348 INTEGER :: i, k, j, l, m
351 #include "petsc/finclude/petsc.h"
352 mpi_comm :: communicator
354 m_max_c =
size(list_mode)
361 DO i = 1,
SIZE(list_mode)
364 IF (minval(mesh%rr(2,mesh%jj(:,m))) > -epsilon)
THEN
366 ELSE IF (maxval(mesh%rr(2,mesh%jj(:,m))) < epsilon)
THEN
369 WRITE(*,*)
'Attention, element a cheval entre nord et sud'
372 DO l = 1, mesh%gauss%l_G
373 ray = sum(mesh%rr(1,mesh%jj(:,m))*mesh%gauss%ww(:,l))
375 e_ns(j) = e_ns(j) + ray*sum(v(mesh%jj(:,m),k,i)* mesh%gauss%ww(:,l))**2*mesh%gauss%rj(l,m)
380 IF (list_mode(i) /= 0)
THEN
384 CALL mpi_allreduce(e_ns(1),e_north(i),1,mpi_double_precision, mpi_sum, communicator, code)
385 CALL mpi_allreduce(e_ns(2),e_south(i),1,mpi_double_precision, mpi_sum, communicator, code)
386 IF (present(e_tot))
THEN
387 e_tot(i) = 0.5d0*(
norme_l2_champ_par(communicator, mesh, list_mode(i:i), v(:,:,i:i)))**2
406 INTEGER,
DIMENSION(:),
INTENT(IN) :: list_mode
407 REAL(KIND=8),
INTENT(IN) :: eps_sa
408 REAL(KIND=8),
DIMENSION(:,:,:) ,
INTENT(IN) :: v
409 CHARACTER(len=1),
INTENT(IN) :: if_u_h
410 REAL(KIND=8),
DIMENSION(mesh%np, SIZE(v,2),SIZE(list_mode)) :: vv
411 REAL(KIND=8),
DIMENSION(3) :: type_sym
413 REAL(KIND=8),
DIMENSION(mesh%np,size(v,2),size(list_mode)),
INTENT(OUT) :: v_out
415 INTEGER,
DIMENSION(:,:),
POINTER :: jj
416 INTEGER,
POINTER :: me
419 #include "petsc/finclude/petsc.h"
420 mpi_comm :: communicator
426 m_max_c =
size(list_mode)
429 CALL
symm_champ(communicator, v, mesh, vv, if_u_h)
432 DO i=1,
size(list_mode)
433 IF (mod(list_mode(i),2) == 0)
THEN
442 v_out(:,1:2,i) = 0.5d0*(v(:,1:2,i) - eps_sa*type_sym(1)*vv(:,1:2,i))
443 v_out(:,3:4,i) = 0.5d0*(v(:,3:4,i) - eps_sa*type_sym(2)*vv(:,3:4,i))
444 v_out(:,5:6,i) = 0.5d0*(v(:,5:6,i) - eps_sa*type_sym(3)*vv(:,5:6,i))
446 write(*,*)
'===Compute anti-symmetric field'
subroutine, public val_ener_sym_glob(communicator, mesh, list_mode, v, e_mode, e_mode_sym, e_mode_anti, type_sym, if_u_h)
subroutine, public create_my_ghost(mesh, LA, ifrom)
subroutine, public val_ener_north_south(communicator, mesh, list_mode, v, e_north, e_south, e_tot)
subroutine, public symm_champ(communicator, vv_in, mesh, vv_out, if_u_h)
subroutine, public champ_total_anti_sym(communicator, mesh, list_mode, eps_sa, v, v_out, if_u_h)
subroutine, public val_ener_sym_centrale(communicator, mesh, list_mode, v, e_mode, e_mode_sym, e_mode_anti, if_u_h)
subroutine, public symmetric_points(mesh_loc, mesh_glob, ltg_LA)
real(kind=8) function norme_l2_champ_par(communicator, mesh, list_mode, v)
subroutine, public val_ener_sym_rpi(communicator, mesh, list_mode, v, e_mode, e_mode_sym, e_mode_anti, type_sym, if_u_h)