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)
381 #include "petsc/finclude/petsc.h" 384 INTEGER ,
INTENT(IN) :: n_bord
385 TYPE(
dyn_int_line),
DIMENSION(:),
INTENT(IN) :: list, perlist
387 INTEGER,
PARAMETER :: nmaxcols = 300
389 INTEGER,
DIMENSION(nmaxcols) :: cols
390 REAL(KIND=8),
DIMENSION(nmaxcols) :: vals
391 INTEGER,
DIMENSION(:),
ALLOCATABLE :: n_cols_i
392 INTEGER,
DIMENSION(1) :: idxn
393 INTEGER,
DIMENSION(:,:),
ALLOCATABLE :: jdxn
394 REAL(KIND=8),
DIMENSION(:,:),
ALLOCATABLE :: vals_pi
395 INTEGER :: n, l, i, pi, n_D, k
397 petscerrorcode :: ierr
399 WRITE(*,*)
'Entering periodic_matrix_petsc' 401 CALL matsetoption (matrix, mat_row_oriented, petsc_false, ierr)
402 CALL matsetoption (matrix, mat_keep_nonzero_pattern, petsc_true, ierr)
404 DO k = 1,
SIZE(la%loc_to_glob,1)
406 n_d =
SIZE(list(n)%DIL)
408 ALLOCATE(jdxn(n_d,nmaxcols), vals_pi(n_d,nmaxcols), n_cols_i(n_d))
413 DO l = 1,
SIZE(list(n)%DIL)
414 idxn(1) = la%loc_to_glob(k,list(n)%DIL(l))
415 CALL matgetrow(matrix, idxn(1)-1,ncols,cols,vals,ierr)
417 jdxn(l,1:ncols) = cols(1:ncols)
418 vals_pi(l,1:ncols) = vals(1:ncols)
419 CALL matrestorerow(matrix, idxn(1)-1, ncols,cols, vals, ierr)
423 idxn(1) = la%loc_to_glob(k,perlist(n)%DIL(l))-1
424 CALL matsetvalues(matrix, 1,idxn, n_cols_i(l), jdxn(l,1:n_cols_i(l)), &
425 vals_pi(l:l,1:n_cols_i(l)), add_values, ierr)
427 DEALLOCATE(jdxn, vals_pi, n_cols_i)
431 CALL matassemblybegin(matrix,mat_final_assembly,ierr)
432 CALL matassemblyend(matrix,mat_final_assembly,ierr)
435 n_d =
SIZE(list(n)%DIL)
438 CALL matzerorows(matrix, n_d, la%loc_to_glob(k,list(n)%DIL(:))-1, 1.d0, &
439 petsc_null_vec, petsc_null_vec, ierr)
445 DO l = 1,
SIZE(list(n)%DIL)
446 i = la%loc_to_glob(k,list(n)%DIL(l))
447 pi = la%loc_to_glob(k,perlist(n)%DIL(l))
448 CALL matsetvalue(matrix, i-1, pi-1, -1.d0, insert_values, ierr)
451 CALL matassemblybegin(matrix,mat_final_assembly,ierr)
452 CALL matassemblyend(matrix,mat_final_assembly,ierr)
461 #include "petsc/finclude/petsc.h" 464 INTEGER ,
INTENT(IN) :: n_bord
465 TYPE(
dyn_int_line),
DIMENSION(:),
INTENT(IN) :: list, perlist
467 INTEGER,
DIMENSION(:),
ALLOCATABLE :: idxn, jdxn
468 REAL(KIND=8),
DIMENSION(:),
ALLOCATABLE :: vals, bs
471 petscerrorcode :: ierr
474 DO k = 1,
SIZE(la%loc_to_glob,1)
476 n_d =
SIZE(list(n)%DIL)
477 ALLOCATE(idxn(n_d), vals(n_d), jdxn(n_d), bs(n_d))
478 idxn = la%loc_to_glob(k,list(n)%DIL(:))-1
479 jdxn = la%loc_to_glob(k,perlist(n)%DIL(:))-1
480 CALL vecgetvalues(v_rhs, n_d, idxn, vals, ierr)
481 CALL vecassemblybegin(v_rhs,ierr)
482 CALL vecassemblyend(v_rhs,ierr)
485 CALL vecsetvalues(v_rhs, n_d, jdxn, vals, add_values, ierr)
486 CALL vecassemblybegin(v_rhs,ierr)
487 CALL vecassemblyend(v_rhs,ierr)
488 CALL vecsetvalues(v_rhs, n_d, idxn, bs, insert_values, ierr)
489 CALL vecassemblybegin(v_rhs,ierr)
490 CALL vecassemblyend(v_rhs,ierr)
492 IF (
ALLOCATED(idxn))
DEALLOCATE(idxn, jdxn, vals, bs)
subroutine, public prep_periodic_bloc(my_periodic, mesh, periodic, nb_bloc)
subroutine, public prep_periodic_h_p_phi_bc(my_periodic, H_mesh, pmag_mesh, phi_mesh, H_p_phi_per)
subroutine list_periodic(np, jjs, sides, rr, side1, side2, e, list_out, perlist_out)
subroutine, public periodic_matrix_petsc(n_bord, list, perlist, matrix, LA)
subroutine, public prep_periodic(my_periodic, mesh, periodic)
subroutine, public periodic_rhs_petsc(n_bord, list, perlist, v_rhs, LA)