9 FUNCTION dot_product_sf(communicator, mesh, list_mode, v, w) RESULT(norm)
13 #include "petsc/finclude/petsc.h" 17 INTEGER,
DIMENSION(:),
INTENT(IN) :: list_mode
18 REAL(KIND=8),
DIMENSION(:,:,:),
INTENT(IN) :: v, w
19 REAL(KIND=8) :: norm_loc, norm_tot, norm
22 mpi_comm,
DIMENSION(2) :: communicator
32 CALL mpi_allreduce(norm_loc,norm_tot,1,mpi_double_precision, mpi_sum, communicator(2), code)
33 CALL mpi_allreduce(norm_tot,norm,1,mpi_double_precision, mpi_sum, communicator(1), code)
39 FUNCTION norm_sf(communicator, norm_type, mesh, list_mode, v) RESULT(norm)
43 #include "petsc/finclude/petsc.h" 47 INTEGER,
DIMENSION(:),
INTENT(IN) :: list_mode
48 REAL(KIND=8),
DIMENSION(:,:,:),
INTENT(IN) :: v
49 CHARACTER(*),
INTENT(IN) :: norm_type
50 REAL(KIND=8) :: norm_loc, norm_tot, norm
51 INTEGER :: deb, fin, code
53 mpi_comm,
DIMENSION(2) :: communicator
63 IF (norm_type(deb:fin)==
'L2')
THEN 65 ELSE IF (norm_type(deb:fin)==
'H1')
THEN 68 ELSE IF (norm_type(deb:fin)==
'sH1')
THEN 70 ELSE IF (norm_type(deb:fin)==
'div')
THEN 72 ELSE IF (norm_type(deb:fin)==
'curl')
THEN 75 WRITE(*,*)
' BUG in norm, norm_type not programmed yet' , norm_type(deb:fin)
79 CALL mpi_allreduce(norm_loc**2,norm_tot,1,mpi_double_precision, mpi_sum, communicator(2), code)
80 CALL mpi_allreduce(norm_tot,norm,1,mpi_double_precision, mpi_sum, communicator(1), code)
82 norm = sqrt(norm*2*pi)
87 FUNCTION norm_s(communicator, norm_type, mesh, list_mode, v) RESULT(norm)
91 #include "petsc/finclude/petsc.h" 95 INTEGER,
DIMENSION(:),
INTENT(IN) :: list_mode
96 REAL(KIND=8),
DIMENSION(:,:,:),
INTENT(IN) :: v
97 CHARACTER(*),
INTENT(IN) :: norm_type
98 REAL(KIND=8) :: norm_loc, norm_tot, norm
99 INTEGER :: deb, fin, code
101 mpi_comm,
DIMENSION(2) :: communicator
111 IF (norm_type(deb:fin)==
'L2')
THEN 113 ELSE IF (norm_type(deb:fin)==
'H1')
THEN 116 ELSE IF (norm_type(deb:fin)==
'sH1')
THEN 118 ELSE IF (norm_type(deb:fin)==
'div')
THEN 120 ELSE IF (norm_type(deb:fin)==
'curl')
THEN 123 WRITE(*,*)
' BUG in norm, norm_type not programmed yet' , norm_type(deb:fin)
127 CALL mpi_allreduce(norm_loc**2,norm,1,mpi_double_precision, mpi_sum, communicator(1), code)
129 norm = sqrt(norm*2*pi)
138 #include "petsc/finclude/petsc.h" 142 INTEGER,
DIMENSION(:),
INTENT(IN) :: list_mode
143 REAL(KIND=8),
DIMENSION(:,:,:),
INTENT(IN) :: v
144 REAL(KIND=8) :: norm_loc, norm_tot, norm
146 INTEGER,
DIMENSION(1) :: zero_mode
148 mpi_comm,
DIMENSION(2) :: communicator
155 IF (minval(list_mode)==0)
THEN 156 zero_mode = minloc(list_mode)
162 CALL mpi_allreduce(norm_loc,norm,1,mpi_double_precision, mpi_sum, communicator(1), code)
168 #include "petsc/finclude/petsc.h" 171 REAL(KIND=8),
INTENT(IN) :: norm_loc
172 REAL(KIND=8),
INTENT(OUT) :: norm_tot
175 mpi_comm :: communicator
178 CALL mpi_allreduce(norm_loc,norm_tot,1,mpi_double_precision, mpi_sum, communicator, code)
180 norm_tot = norm_tot*(2*pi)
187 #include "petsc/finclude/petsc.h" 191 INTEGER,
DIMENSION(:),
INTENT(IN) :: list_mode
192 REAL(KIND=8),
DIMENSION(:,:,:),
INTENT(IN) :: v
193 REAL(KIND=8) :: norm_loc, norm
194 mpi_comm :: communicator
208 #include "petsc/finclude/petsc.h" 212 INTEGER,
DIMENSION(:),
INTENT(IN) :: list_mode
213 REAL(KIND=8),
DIMENSION(:,:,:),
INTENT(IN) :: v
214 REAL(KIND=8) :: norm_loc, norm
215 mpi_comm :: communicator
222 FUNCTION norme_div_par(communicator, H_mesh, list_mode, v) RESULT(norm)
225 #include "petsc/finclude/petsc.h" 229 INTEGER,
DIMENSION(:),
INTENT(IN) :: list_mode
230 REAL(KIND=8),
DIMENSION(:,:,:),
INTENT(IN) :: v
231 REAL(KIND=8) :: norm_loc, norm
232 mpi_comm :: communicator
234 norm_loc =
norme_div(h_mesh, list_mode, v)
239 FUNCTION norme_curl_par(communicator, H_mesh, list_mode, v) RESULT(norm)
242 #include "petsc/finclude/petsc.h" 246 INTEGER,
DIMENSION(:),
INTENT(IN) :: list_mode
247 REAL(KIND=8),
DIMENSION(:,:,:),
INTENT(IN) :: v
248 REAL(KIND=8) :: norm_loc, norm
249 mpi_comm :: communicator
258 #include "petsc/finclude/petsc.h" 262 INTEGER,
DIMENSION(:),
INTENT(IN) :: list_mode
263 REAL(KIND=8),
DIMENSION(:,:,:),
INTENT(IN) :: field
264 REAL(KIND=8),
DIMENSION(3),
INTENT(OUT) :: moments_out
266 REAL(KIND=8),
DIMENSION(3) :: moments_loc
267 REAL(KIND=8) :: m_x, m_y, ray, zed, urc, urs, utc, uts, uzc, uzs
268 INTEGER :: m, l, code, i
272 DO i = 1,
SIZE(list_mode)
273 IF (list_mode(i)==0)
THEN 275 DO l = 1, mesh%gauss%l_G
276 ray = sum(mesh%rr(1,mesh%jj(:,m))*mesh%gauss%ww(:,l))
277 utc = sum(field(mesh%jj(:,m),3,i)*mesh%gauss%ww(:,l))
279 moments_loc(3) = moments_loc(3) + ray**2*utc*mesh%gauss%rj(l,m)
282 ELSE IF (list_mode(i)==1)
THEN 284 DO l = 1, mesh%gauss%l_G
285 ray = sum(mesh%rr(1,mesh%jj(:,m))*mesh%gauss%ww(:,l))
286 zed = sum(mesh%rr(2,mesh%jj(:,m))*mesh%gauss%ww(:,l))
288 urc = sum(field(mesh%jj(:,m),1,i)*mesh%gauss%ww(:,l))
289 urs = sum(field(mesh%jj(:,m),2,i)*mesh%gauss%ww(:,l))
290 utc = sum(field(mesh%jj(:,m),3,i)*mesh%gauss%ww(:,l))
291 uts = sum(field(mesh%jj(:,m),4,i)*mesh%gauss%ww(:,l))
292 uzc = sum(field(mesh%jj(:,m),5,i)*mesh%gauss%ww(:,l))
293 uzs = sum(field(mesh%jj(:,m),6,i)*mesh%gauss%ww(:,l))
295 m_x = ray*(-zed*(urs+utc)+ray*uzs)
296 m_y = ray*(zed*(urc-uts)-ray*uzc)
298 moments_loc(1) = moments_loc(1) + m_x*mesh%gauss%rj(l,m)
299 moments_loc(2) = moments_loc(2) + m_y*mesh%gauss%rj(l,m)
305 CALL mpi_allreduce(moments_loc, moments_out, 3, mpi_double_precision, mpi_sum, petsc_comm_world, code)
306 moments_out(1) = acos(-1.d0)*moments_out(1)
307 moments_out(2) = acos(-1.d0)*moments_out(2)
308 moments_out(3) = 2*acos(-1.d0)*moments_out(3)
312 SUBROUTINE moments(communicator, H_mesh, list_mode, v, dipole_out, quadripole_out)
314 #include "petsc/finclude/petsc.h" 318 INTEGER,
DIMENSION(:),
INTENT(IN) :: list_mode
319 REAL(KIND=8),
DIMENSION(:,:,:),
INTENT(IN) :: v
320 REAL(KIND=8),
DIMENSION(3),
INTENT(OUT):: dipole_out
321 REAL(KIND=8),
DIMENSION(3,3),
INTENT(OUT):: quadripole_out
322 REAL(KIND=8),
DIMENSION(3) :: dipole
323 REAL(KIND=8),
DIMENSION(3,3) :: quadripole
324 REAL(KIND=8) :: jr, ray, zed, pi
325 INTEGER :: m_max_c, mode, k, m, l, ni, i
327 REAL(KIND=8),
DIMENSION(6) :: c
329 mpi_comm :: communicator
332 m_max_c =
SIZE(list_mode)
336 IF (
SIZE(v,1)/=h_mesh%np .OR.
SIZE(v,2)/=6 .OR.
SIZE(v,3)/=m_max_c )
THEN 337 WRITE(*,*)
' BUG in MOMENTS',
SIZE(v,1), h_mesh%np
341 DO l = 1, h_mesh%gauss%l_G
346 DO ni = 1, h_mesh%gauss%n_w; i = h_mesh%jj(ni,m)
347 ray = ray + h_mesh%rr(1,i)*h_mesh%gauss%ww(ni,l)
348 zed = zed + h_mesh%rr(2,i)*h_mesh%gauss%ww(ni,l)
350 jr = ray * h_mesh%gauss%rj(l,m)
354 IF (mode /=0 .AND. mode /=1 .AND. mode /=2) cycle
358 DO ni = 1,h_mesh%gauss%n_w; i = h_mesh%jj(ni,m)
360 c(1) = c(1) + ( mode/ray*v(i,6,k)*h_mesh%gauss%ww(ni,l) &
361 - v(i,3,k)*h_mesh%gauss%dw(2,ni,l,m))
362 c(2) = c(2) + (-mode/ray*v(i,5,k)*h_mesh%gauss%ww(ni,l) &
363 - v(i,4,k)*h_mesh%gauss%dw(2,ni,l,m))
365 c(3) = c(3) + (v(i,1,k)*h_mesh%gauss%dw(2,ni,l,m) &
366 - v(i,5,k)*h_mesh%gauss%dw(1,ni,l,m))
367 c(4) = c(4) + (v(i,2,k)*h_mesh%gauss%dw(2,ni,l,m) &
368 - v(i,6,k)*h_mesh%gauss%dw(1,ni,l,m))
370 c(5) = c(5) + (v(i,3,k)*h_mesh%gauss%dw(1,ni,l,m) &
371 + v(i,3,k)*h_mesh%gauss%ww(ni,l)/ray &
372 - mode/ray*v(i,2,k)*h_mesh%gauss%ww(ni,l))
373 c(6) = c(6) + (v(i,4,k)*h_mesh%gauss%dw(1,ni,l,m) &
374 + v(i,4,k)*h_mesh%gauss%ww(ni,l)/ray &
375 + mode/ray*v(i,1,k)*h_mesh%gauss%ww(ni,l))
380 dipole(3) = dipole(3) + 2*pi*ray*c(3)*jr
381 quadripole(1,1) = quadripole(1,1) + pi*(-zed*ray*c(3))*jr
382 quadripole(1,2) = quadripole(1,2) + pi*(-zed*ray*c(1)+ray*ray*c(5))*jr
383 quadripole(2,1) = quadripole(2,1) + pi*( zed*ray*c(1)-ray*ray*c(5))*jr
384 quadripole(2,2) = quadripole(2,2) + pi*(-zed*ray*c(3))*jr
385 quadripole(3,3) = quadripole(3,3)+2*pi*( zed*ray*c(3))*jr
386 ELSE IF (mode == 1)
THEN 387 dipole(1) = dipole(1) + pi*(-zed*c(3) -zed*c(2) +ray*c(6))*jr
388 dipole(2) = dipole(2) + pi*(-zed*c(4) +zed*c(1) -ray*c(5))*jr
389 quadripole(1,3) = quadripole(1,3) + pi*(-zed*c(3)-zed*c(2)+ray*c(6))*zed*jr
390 quadripole(2,3) = quadripole(2,3) + pi*(-zed*c(4)+zed*c(1)-ray*c(5))*zed*jr
391 quadripole(3,1) = quadripole(3,1) + pi*(ray*ray*c(3))*jr
392 quadripole(3,2) = quadripole(3,2) + pi*(ray*ray*c(4))*jr
393 ELSE IF (mode == 2)
THEN 394 quadripole(1,1) = quadripole(1,1) + pi*(-zed*c(3)-zed*c(2)+ray*c(6))*ray*jr/2
395 quadripole(1,2) = quadripole(1,2) + pi*(-zed*c(4)+zed*c(1)-ray*c(5))*ray*jr/2
396 quadripole(2,1) = quadripole(2,1) + pi*(-zed*c(4)+zed*c(1)-ray*c(5))*ray*jr/2
397 quadripole(2,2) = quadripole(2,2) + pi*( zed*c(3)+zed*c(2)-ray*c(6))*ray*jr/2
405 CALL mpi_allreduce(dipole, dipole_out, 3,mpi_double_precision, &
406 mpi_sum, communicator, code)
407 CALL mpi_allreduce(quadripole,quadripole_out,9,mpi_double_precision, &
408 mpi_sum, communicator, code)
integer function, public last_of_string(string)
subroutine moments(communicator, H_mesh, list_mode, v, dipole_out, quadripole_out)
real(kind=8) function dot_product_champ(mesh, list_mode, v, w)
real(kind=8) function norme_l2_champ(mesh, list_mode, v)
real(kind=8) function norme_l1_one_mode(mesh, mode_idx, v)
real(kind=8) function norme_h1_champ(mesh, list_mode, v)
real(kind=8) function norme_div(H_mesh, list_mode, v)
real(kind=8) function norme_curl(H_mesh, list_mode, v)
subroutine integration_mode(communicator, norm_loc, norm_tot)
real(kind=8) function norm_s(communicator, norm_type, mesh, list_mode, v)
integer function, public start_of_string(string)
subroutine angular_momentum(mesh, list_mode, field, moments_out)
real(kind=8) function norme_h1_champ_par(communicator, mesh, list_mode, v)
real(kind=8) function norme_div_par(communicator, H_mesh, list_mode, v)
real(kind=8) function norme_l2_champ_par(communicator, mesh, list_mode, v)
real(kind=8) function norme_curl_par(communicator, H_mesh, list_mode, v)
real(kind=8) function norm_sf(communicator, norm_type, mesh, list_mode, v)
real(kind=8) function dot_product_sf(communicator, mesh, list_mode, v, w)
real(kind=8) function norm_s_l1_zero_mode(communicator, mesh, list_mode, v)