27 INTEGER,
DIMENSION(:),
POINTER :: list_loc, perlist_loc, list_dom, perlist_dom
28 INTEGER :: n, side1, side2, n_b, nx, i
29 REAL(KIND=8),
DIMENSION(:),
POINTER :: e
31 WRITE (*,*)
'Loading periodic-data file ...'
33 IF (mesh%np == 0)
THEN
34 WRITE(*,*)
'no mesh on this proc'
38 ALLOCATE(e(
SIZE(my_periodic%vect_e,1)))
40 periodic%n_bord = my_periodic%nb_periodic_pairs
42 WRITE(*,*)
'PREP_MESH_PERIODIC: trop de bords periodiques'
48 side1 = my_periodic%list_periodic(1,n)
49 side2 = my_periodic%list_periodic(2,n)
50 e = my_periodic%vect_e(:,n)
52 CALL
list_periodic(mesh%np, mesh%jjs, mesh%sides, mesh%rr, side1, side2, e, &
53 list_loc, perlist_loc)
56 ALLOCATE(list_dom(n_b), perlist_dom(n_b))
59 IF (max(list_loc(i),perlist_loc(i)) .LE. mesh%dom_np)
THEN
61 list_dom(nx) = list_loc(i)
62 perlist_dom(nx) = perlist_loc(i)
63 ELSE IF (min(list_loc(i),perlist_loc(i)) .LE. mesh%dom_np)
THEN
64 WRITE(*,*)
'BUG in prep_periodic'
68 IF (n_b /= nx)
WRITE(*,*)
'WARNING, I have removed', n_b-nx,
' periodic pairs in prep_periodic'
72 periodic%list(n)%DIL = list_dom(1:n_b)
73 periodic%perlist(n)%DIL = perlist_dom(1:n_b)
75 DEALLOCATE(list_loc,perlist_loc, list_dom, perlist_dom)
80 WRITE (*,*)
'Treatment of periodic-data done'
92 INTEGER,
INTENT(IN) :: nb_bloc
93 INTEGER,
DIMENSION(:),
POINTER :: list_loc, perlist_loc, list_dom, perlist_dom
94 INTEGER :: n, side1, side2, nsize, n_b
95 INTEGER :: k, k_deb, k_fin, nx, i
96 REAL(KIND=8),
DIMENSION(2) :: e
98 WRITE (*,*)
'Loading periodic-data file ...'
100 IF (mesh%np == 0)
THEN
101 WRITE(*,*)
'no mesh on this proc'
105 periodic%n_bord = my_periodic%nb_periodic_pairs
107 WRITE(*,*)
'PREP_MESH_PERIODIC: trop de bords periodiques'
113 side1 = my_periodic%list_periodic(1,n)
114 side2 = my_periodic%list_periodic(2,n)
115 e = my_periodic%vect_e(:,n)
117 CALL
list_periodic(mesh%np, mesh%jjs, mesh%sides, mesh%rr, side1, side2, e, &
118 list_loc, perlist_loc)
121 n_b =
SIZE(perlist_loc)
122 ALLOCATE(list_dom(n_b), perlist_dom(n_b))
125 IF (max(list_loc(i),perlist_loc(i)) .LE. mesh%dom_np)
THEN
127 list_dom(nx) = list_loc(i)
128 perlist_dom(nx) = perlist_loc(i)
129 ELSE IF (min(list_loc(i),perlist_loc(i)) .LE. mesh%dom_np)
THEN
130 WRITE(*,*)
'BUG in prep_periodic_bloc'
134 IF (n_b /= nx)
WRITE(*,*)
'WARNING, I have removed', n_b-nx,
' periodic pairs in prep_periodic_bloc'
144 periodic%list(n)%DIL(k_deb:k_fin) = list_dom(1:n_b) + (k-1)*mesh%dom_np
145 periodic%perlist(n)%DIL(k_deb:k_fin) = perlist_dom(1:n_b) + (k-1)*mesh%dom_np
148 DEALLOCATE(list_loc,perlist_loc, list_dom, perlist_dom)
152 WRITE (*,*)
'Treatment of periodic-data done'
163 TYPE(mesh_type) :: h_mesh, pmag_mesh, phi_mesh
165 INTEGER,
DIMENSION(:),
POINTER :: b_list_loc, b_perlist_loc, &
166 e_list_loc, e_perlist_loc, p_list_loc, p_perlist_loc
167 INTEGER,
DIMENSION(:),
POINTER :: b_list_dom, b_perlist_dom, &
168 e_list_dom, e_perlist_dom, p_list_dom, p_perlist_dom
169 INTEGER :: n, side1, side2, nsize, n_b, n_e, n_p, nx, i
170 REAL(KIND=8),
DIMENSION(2) :: e
172 WRITE (*,*)
'Loading periodic-data file ...'
174 h_p_phi_per%n_bord = my_periodic%nb_periodic_pairs
175 IF (h_p_phi_per%n_bord .GT. 20)
THEN
176 WRITE(*,*)
'PREP_MESH_PERIODIC: trop de bords periodiques'
180 DO n= 1, h_p_phi_per%n_bord
182 side1 = my_periodic%list_periodic(1,n)
183 side2 = my_periodic%list_periodic(2,n)
184 e = my_periodic%vect_e(:,n)
187 IF (h_mesh%np == 0)
THEN
190 CALL
list_periodic(h_mesh%np, h_mesh%jjs, h_mesh%sides, h_mesh%rr, side1, side2, e, &
191 b_list_loc, b_perlist_loc)
192 n_b =
SIZE(b_list_loc)
193 ALLOCATE(b_list_dom(n_b), b_perlist_dom(n_b))
196 IF (max(b_list_loc(i),b_perlist_loc(i)) .LE. h_mesh%dom_np)
THEN
198 b_list_dom(nx) = b_list_loc(i)
199 b_perlist_dom(nx) = b_perlist_loc(i)
200 ELSE IF (min(b_list_loc(i),b_perlist_loc(i)) .LE. h_mesh%dom_np)
THEN
201 WRITE(*,*)
'BUG in prep_periodic_H_p_phi_bc (H)'
205 IF (n_b /= nx)
WRITE(*,*)
'WARNING, I have removed', n_b-nx,
' periodic pairs on H'
209 IF (pmag_mesh%np == 0)
THEN
212 CALL
list_periodic(pmag_mesh%np, pmag_mesh%jjs, pmag_mesh%sides, pmag_mesh%rr, side1, side2, e, &
213 p_list_loc, p_perlist_loc)
214 n_p =
SIZE(p_list_loc)
215 ALLOCATE(p_list_dom(n_p), p_perlist_dom(n_p))
218 IF (max(p_list_loc(i),p_perlist_loc(i)) .LE. pmag_mesh%dom_np)
THEN
220 p_list_dom(nx) = p_list_loc(i)
221 p_perlist_dom(nx) = p_perlist_loc(i)
222 ELSE IF (min(p_list_loc(i),p_perlist_loc(i)) .LE. pmag_mesh%dom_np)
THEN
223 WRITE(*,*)
'BUG in prep_periodic_H_p_phi_bc (pmag) '
227 IF (n_p /= nx)
WRITE(*,*)
'WARNING, I have removed', n_p-nx,
' periodic pairs on pmag'
231 IF (phi_mesh%np==0)
THEN
234 CALL
list_periodic(phi_mesh%np, phi_mesh%jjs, phi_mesh%sides, phi_mesh%rr, side1, side2, e, &
235 e_list_loc, e_perlist_loc)
236 n_e =
SIZE(e_list_loc)
237 ALLOCATE(e_list_dom(n_e), e_perlist_dom(n_e))
240 IF (max(e_list_loc(i),e_perlist_loc(i)) .LE. phi_mesh%dom_np)
THEN
242 e_list_dom(nx) = e_list_loc(i)
243 e_perlist_dom(nx) = e_perlist_loc(i)
244 ELSE IF (min(e_list_loc(i),e_perlist_loc(i)) .LE. phi_mesh%dom_np)
THEN
245 WRITE(*,*)
'BUG in prep_periodic_H_p_phi_bc (phi) '
249 IF (n_e /= nx)
WRITE(*,*)
'WARNING, I have removed', n_e-nx,
' periodic pairs on phi'
256 nsize = 3*n_b + n_p + n_e
258 ALLOCATE(h_p_phi_per%list(n)%DIL(nsize), h_p_phi_per%perlist(n)%DIL(nsize))
262 h_p_phi_per%list(n)%DIL(1:n_b) = b_list_dom(1:n_b)
263 h_p_phi_per%perlist(n)%DIL(1:n_b) = b_perlist_dom(1:n_b)
265 h_p_phi_per%list(n)%DIL(n_b+1:2*n_b) = b_list_dom(1:n_b) + h_mesh%dom_np
266 h_p_phi_per%perlist(n)%DIL(n_b+1:2*n_b) = b_perlist_dom(1:n_b) + h_mesh%dom_np
268 h_p_phi_per%list(n)%DIL(2*n_b+1:3*n_b) = b_list_dom(1:n_b) + 2*h_mesh%dom_np
269 h_p_phi_per%perlist(n)%DIL(2*n_b+1:3*n_b) = b_perlist_dom(1:n_b) + 2*h_mesh%dom_np
274 h_p_phi_per%list(n)%DIL(3*n_b+1:3*n_b+n_p) = p_list_dom(1:n_p) + 3*h_mesh%dom_np
275 h_p_phi_per%perlist(n)%DIL(3*n_b+1:3*n_b+n_p) = p_perlist_dom(1:n_p) + 3*h_mesh%dom_np
279 h_p_phi_per%list(n)%DIL(3*n_b+n_p+1:) = e_list_dom(1:n_e) + 3*h_mesh%dom_np + pmag_mesh%dom_np
280 h_p_phi_per%perlist(n)%DIL(3*n_b+n_p+1:) = e_perlist_dom(1:n_e) + 3*h_mesh%dom_np + pmag_mesh%dom_np
283 IF (
ASSOCIATED(b_list_loc)) nullify(b_list_loc, b_perlist_loc, b_list_dom, b_perlist_dom)
284 IF (
ASSOCIATED(p_list_loc)) nullify(p_list_loc, p_perlist_loc, p_list_dom, p_perlist_dom)
285 IF (
ASSOCIATED(e_list_loc)) nullify(e_list_loc, e_perlist_loc, e_list_dom, e_perlist_dom)
289 WRITE (*,*)
'Treatment of periodic-data done'
293 SUBROUTINE list_periodic(np, jjs, sides, rr, side1, side2, e, list_out, perlist_out)
296 INTEGER,
INTENT(IN) :: np
297 INTEGER,
DIMENSION(:,:),
INTENT(IN) :: jjs
298 INTEGER,
DIMENSION(:),
INTENT(IN) :: sides
299 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: rr
300 INTEGER,
INTENT(IN) :: side1, side2
301 REAL(KIND=8),
DIMENSION(:),
INTENT(IN) :: e
302 INTEGER,
DIMENSION(:),
POINTER :: list_out, perlist_out
303 INTEGER,
DIMENSION(:),
ALLOCATABLE :: list, perlist
304 LOGICAL,
DIMENSION(np) :: virgin
305 REAL(KIND=8),
DIMENSION(SIZE(rr,1)) :: ri
306 INTEGER :: ms, ns, i, j, long, inter
307 REAL(KIND=8) :: r, epsilon = 1.d-9
310 IF (
ALLOCATED(list))
DEALLOCATE(list)
311 IF (
ALLOCATED(perlist))
DEALLOCATE(perlist)
313 ALLOCATE (list(np), perlist(np))
317 DO ms = 1,
SIZE(sides)
319 IF (sides(ms) .EQ. side1)
THEN
320 DO ns = 1,
SIZE(jjs,1)
321 IF (virgin(jjs(ns,ms)))
THEN
324 virgin(jjs(ns,ms)) = .false.
327 ELSE IF (sides(ms) .EQ. side2)
THEN
328 DO ns = 1,
SIZE(jjs,1)
329 IF (virgin(jjs(ns,ms)))
THEN
331 perlist(j) = jjs(ns,ms)
332 virgin(jjs(ns,ms)) = .false.
341 WRITE(*,*)
' FEM_PERIODIC: side1 and side2 have', &
342 ' different numbers of points'
348 ri = rr(:,list(i))+e(:)
352 r = sum(abs(ri - rr(:,perlist(j))))
354 IF (r .LE. epsilon )
THEN
356 perlist(i) = perlist(j)
363 WRITE(*,*)
' BUG dans data_periodic ou le maillage:', &
364 ' side1 + e /= side2'
365 WRITE(*,*)
' i = ', i
370 ALLOCATE (list_out(long))
371 list_out(1:long) = list(1:long)
372 ALLOCATE (perlist_out(long))
373 perlist_out(1:long) = perlist(1:long)
382 INTEGER ,
INTENT(IN) :: n_bord
383 TYPE(dyn_int_line),
DIMENSION(:),
INTENT(IN) :: list, perlist
385 INTEGER,
PARAMETER :: nmaxcols = 300
387 INTEGER,
DIMENSION(nmaxcols) :: cols
388 REAL(KIND=8),
DIMENSION(nmaxcols) :: vals
389 INTEGER,
DIMENSION(:),
ALLOCATABLE :: n_cols_i
390 INTEGER,
DIMENSION(1) :: idxn
391 INTEGER,
DIMENSION(:,:),
ALLOCATABLE :: jdxn
392 REAL(KIND=8),
DIMENSION(:,:),
ALLOCATABLE :: vals_pi
393 INTEGER :: n, l, i, pi, n_d, k
394 #include "petsc/finclude/petsc.h"
396 petscerrorcode :: ierr
398 WRITE(*,*)
'Entering periodic_matrix_petsc'
400 CALL matsetoption(matrix, mat_row_oriented, petsc_false, ierr)
401 CALL matsetoption(matrix, mat_keep_nonzero_pattern, petsc_true, ierr)
403 DO k = 1,
SIZE(la%loc_to_glob,1)
405 n_d =
SIZE(list(n)%DIL)
407 ALLOCATE(jdxn(n_d,nmaxcols), vals_pi(n_d,nmaxcols), n_cols_i(n_d))
412 DO l = 1,
SIZE(list(n)%DIL)
413 idxn(1) = la%loc_to_glob(k,list(n)%DIL(l))
414 CALL matgetrow(matrix, idxn(1)-1,ncols,cols,vals,ierr)
416 jdxn(l,1:ncols) = cols(1:ncols)
417 vals_pi(l,1:ncols) = vals(1:ncols)
418 CALL matrestorerow(matrix, idxn(1)-1, ncols,cols, vals, ierr)
422 idxn(1) = la%loc_to_glob(k,perlist(n)%DIL(l))-1
423 CALL matsetvalues(matrix, 1,idxn, n_cols_i(l), jdxn(l,1:n_cols_i(l)), &
424 vals_pi(l:l,1:n_cols_i(l)), add_values, ierr)
426 DEALLOCATE(jdxn, vals_pi, n_cols_i)
430 CALL matassemblybegin(matrix,mat_final_assembly,ierr)
431 CALL matassemblyend(matrix,mat_final_assembly,ierr)
434 n_d =
SIZE(list(n)%DIL)
435 CALL matzerorows(matrix, n_d, la%loc_to_glob(k,list(n)%DIL(:))-1, 1.d0, &
436 petsc_null_object, petsc_null_object, ierr)
442 DO l = 1,
SIZE(list(n)%DIL)
443 i = la%loc_to_glob(k,list(n)%DIL(l))
444 pi = la%loc_to_glob(k,perlist(n)%DIL(l))
445 CALL matsetvalue(matrix, i-1, pi-1, -1.d0, insert_values, ierr)
448 CALL matassemblybegin(matrix,mat_final_assembly,ierr)
449 CALL matassemblyend(matrix,mat_final_assembly,ierr)
459 INTEGER ,
INTENT(IN) :: n_bord
460 TYPE(dyn_int_line),
DIMENSION(:),
INTENT(IN) :: list, perlist
462 INTEGER,
DIMENSION(:),
ALLOCATABLE :: idxn, jdxn
463 REAL(KIND=8),
DIMENSION(:),
ALLOCATABLE :: vals, bs
465 #include "petsc/finclude/petsc.h"
467 petscerrorcode :: ierr
470 DO k = 1,
SIZE(la%loc_to_glob,1)
472 n_d =
SIZE(list(n)%DIL)
473 ALLOCATE(idxn(n_d), vals(n_d), jdxn(n_d), bs(n_d))
474 idxn = la%loc_to_glob(k,list(n)%DIL(:))-1
475 jdxn = la%loc_to_glob(k,perlist(n)%DIL(:))-1
476 CALL vecgetvalues(v_rhs, n_d, idxn, vals, ierr)
477 CALL vecassemblybegin(v_rhs,ierr)
478 CALL vecassemblyend(v_rhs,ierr)
481 CALL vecsetvalues(v_rhs, n_d, jdxn, vals, add_values, ierr)
482 CALL vecassemblybegin(v_rhs,ierr)
483 CALL vecassemblyend(v_rhs,ierr)
484 CALL vecsetvalues(v_rhs, n_d, idxn, bs, insert_values, ierr)
485 CALL vecassemblybegin(v_rhs,ierr)
486 CALL vecassemblyend(v_rhs,ierr)
488 IF (
ALLOCATED(idxn))
DEALLOCATE(idxn, jdxn, vals, bs)
subroutine list_periodic(np, jjs, sides, rr, side1, side2, e, list_out, perlist_out)
subroutine, public prep_periodic_bloc(my_periodic, mesh, periodic, nb_bloc)
subroutine, public periodic_rhs_petsc(n_bord, list, perlist, v_rhs, LA)
subroutine, public prep_periodic_h_p_phi_bc(my_periodic, H_mesh, pmag_mesh, phi_mesh, H_p_phi_per)
subroutine, public periodic_matrix_petsc(n_bord, list, perlist, matrix, LA)
subroutine, public prep_periodic(my_periodic, mesh, periodic)