-
Notifications
You must be signed in to change notification settings - Fork 1
/
wwm_petsc_parallel.F90
579 lines (513 loc) · 22.5 KB
/
wwm_petsc_parallel.F90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
#include "wwm_functions.h"
#ifdef PETSC
# ifdef MPI_PARALL_GRID
MODULE PETSC_PARALLEL
implicit none
#include "finclude/petscsysdef.h"
#include "finclude/petscaodef.h"
#include "finclude/petscisdef.h"
#include "finclude/petscvecdef.h"
#include "finclude/petscmatdef.h"
#include "finclude/petsckspdef.h"
!> IA in CSR format in petsc local order
PetscInt, allocatable :: IA_petsc(:)
!> colum index in CSR format in petsc local order
PetscInt, allocatable :: JA_petsc(:)
!> Matrix values in CSR format
PetscScalar, allocatable :: ASPAR_petsc(:)
!> offdiagonal IA CSR format in petsc local order
PetscInt, allocatable :: oIA_petsc(:)
!> offdiagonal colum index in CSR fromat in petsc GLOABL order?? petsc doku said this should in local order
PetscInt, allocatable :: oJA_petsc(:)
!> offdiagonal submatrix values in CSR format
PetscScalar, allocatable :: oASPAR_petsc(:)
!>
PetscInt, allocatable :: CSR_App2PetscLUT(:), o_CSR_App2PetscLUT(:)
contains
!> - Initialize Petsc
!> - create the mappings
!> - create the matrix and vectors
!> - create the solver and preconditioner
SUBROUTINE PETSC_INIT_PARALLEL
USE DATAPOOL, only: MNP, CCON, NNZ, DBG
use datapool, only: comm, np_global
! np_global - # nodes gloabl
! np - # nodes local non augmented
! npg - # ghost
! npa - # nodes aufmented
! nea - # elemnts augmented
! int::iplg(ip) global element index. local to global LUT
! llsit_type::ipgl(ipgb) ipgb is a global node. global to local LUT
! int::nnp(ip) total # of surrounding nodes for node ip
! int::inp(ip, 1:nnp(ip)) list of surrounding nodes
use petscpool
use petscsys
! use petscao
implicit none
integer :: ierr
call MPI_Comm_rank(comm, rank, ierr)
call MPI_Comm_size(comm, nProcs, ierr)
#ifdef PETSC_DEBUG
call PetscPrintf(PETSC_COMM_WORLD, "PETSC_INIT_PARALLEL\n", petscErr);CHKERRQ(petscErr);
#endif
call createMappings()
call createMatrix()
! create X vector
call VecCreateGhost(PETSC_COMM_WORLD, &
& nNodesWithoutInterfaceGhosts, np_global, nghost, onlyGhosts, &
& myX, petscErr);CHKERRQ(petscErr)
! create B vector
call VecCreateGhost(PETSC_COMM_WORLD, &
& nNodesWithoutInterfaceGhosts, np_global, nghost, onlyGhosts, &
& myB, petscErr);CHKERRQ(petscErr)
call createSolver()
#ifdef PETSC_DEBUG
if(rank == 0) call printSolverTolerance(solver)
if(rank == 0) call printKSPType(solver)
#endif
END SUBROUTINE
!> create PETSC matrix which uses fortran arrays
SUBROUTINE createMatrix()
use datapool, only: np_global
use petscpool
use petscsys
use petscmat
implicit none
call createCSR_petsc()
call MatCreateMPIAIJWithSplitArrays(PETSC_COMM_WORLD, &
& nNodesWithoutInterfaceGhosts, nNodesWithoutInterfaceGhosts, &
& np_global, np_global, &
& IA_petsc, JA_petsc, ASPAR_petsc, &
& oIA_petsc, oJA_petsc, oASPAR_petsc, &
& matrix, petscErr);CHKERRQ(petscErr);
! indicates that any add or insertion that would generate a new entry in the nonzero structure instead produces an error.
call MatSetOption(matrix, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE, petscErr);CHKERRQ(petscErr)
! indicates that any add or insertion that would generate a new entry that has not been preallocated will instead produce an error
call MatSetOption(matrix, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE, petscErr);CHKERRQ(petscErr)
! indicates entries destined for other processors should be dropped, rather than stashed
call MatSetOption(matrix, MAT_IGNORE_OFF_PROC_ENTRIES, PETSC_TRUE, petscErr);CHKERRQ(petscErr)
! you know each process will only zero its own rows.
! This avoids all reductions in the zero row routines and thus improves performance for very large process counts.
call MatSetOption(matrix, MAT_NO_OFF_PROC_ZERO_ROWS, PETSC_TRUE, petscErr);CHKERRQ(petscErr)
end SUBROUTINE
! 1. create IA JA ASPAR petsc arrays
SUBROUTINE createCSR_petsc()
use datapool, only: NNZ, MNE, INE, MNP, DBG, iplg, JA
use petscpool
use algorithm, only: bubbleSort, genericData
implicit none
! max number of adj nodes per node
integer :: maxNumConnNode = 0
integer :: istat
! running variable node number
integer :: IP = 0
! node number in petsc order
integer :: IP_petsc = 0
! running variable
integer :: i = 0, j = 0, o_j = 0
! number of nonzero without interface and ghosts
integer :: nnz_new = 0
! number of nonzeros in the offdiagonal submatrix without interface and ghosts
integer :: o_nnz_new = 0
type(genericData), allocatable :: toSort(:)
integer :: nToSort = 0
type(genericData), allocatable :: o_toSort(:)
integer :: o_nToSort = 0
! calc max number of adj nodes per node
maxNumConnNode = 0
do IP = 1, MNP
if(IA_P(IP+1) - IA_P(IP)-1 > maxNumConnNode) then
maxNumConnNode = IA_P(IP+1) - IA_P(IP)-1
end if
end do
! calc NNZ and offdiagonal NNZ
! iterate over all petsc rows and the nodes in this row
! if one is a ghost or interface nodes, increase offdiagonal NNZ
nnz_new = 0
o_nnz_new = 0
do IP_petsc = 1, nNodesWithoutInterfaceGhosts
IP = PLO2ALO(IP_petsc-1)+1
do i = IA_P(IP)+1, IA_P(IP+1)
if(ALOold2ALO(JA(i)) .eq. -999) then
o_nnz_new = o_nnz_new + 1
else
nnz_new = nnz_new + 1
endif
end do
end do
! write(DBG%FHNDL,*) rank, "nnz_new", nnz_new, " old", NNZ
! write(DBG%FHNDL,*) rank, "o_nnz_new", o_nnz_new
! we have now for every node their connected nodes
! iterate over connNode array to create IA and JA
allocate(IA_petsc(nNodesWithoutInterfaceGhosts+1), &
& JA_petsc(nnz_new), &
& ASPAR_petsc(nnz_new), &
& oIA_petsc(nNodesWithoutInterfaceGhosts+1), &
& oJA_petsc(o_nnz_new), &
& oASPAR_petsc(o_nnz_new), &
! +1 because we have to store the diagonal node number too
& toSort(maxNumConnNode+1), &
& o_toSort(maxNumConnNode+1), &
& stat=istat)
if(istat /= 0) then
write(DBG%FHNDL,*) __FILE__, " Line", __LINE__
stop 'wwm_petsc_parallel l.171'
endif
allocate(CSR_App2PetscLUT(NNZ), o_CSR_App2PetscLUT(NNZ), stat=istat)
if(istat /= 0) then
write(DBG%FHNDL,*) __FILE__, " Line", __LINE__
stop 'wwm_petsc_parallel l.178'
endif
CSR_App2PetscLUT = -999
o_CSR_App2PetscLUT = -998
IA_petsc = 0
JA_petsc = 0
ASPAR_petsc = 0
oIA_petsc = 0
oJA_petsc = 0
oASPAR_petsc = 0
! to create IA_petsc JA_petsc we have to iterate over all nodes and
! their connected nodes and map to petsc order.
! the node numbers of the connected nodes in petsc order are not sorted.
! sort them with a simple bubble sort. yes, bubble sort is vey slow,
! but we have only a few numbers to sort (max 10 i assume).
J = 0
do IP_petsc = 1, nNodesWithoutInterfaceGhosts
IP = PLO2ALO(IP_petsc-1)+1
! fill with the largest numner petscInt can hold
toSort(:)%id = HUGE(0)
nToSort = 0
o_toSort(:)%id = HUGE(0)
o_nToSort = 0
! over all nodes in this row
do i = IA_P(IP) + 1, IA_P(IP+1)
! found a ghost node, treat them special
if(ALOold2ALO(JA(i)) .eq. -999) then
o_ntoSort = o_ntoSort + 1
! store the old position in ASPAR
o_toSort(o_nToSort)%userData = i
!> todo offdiagonal part with petsc global order? don't know why but it seems to work
o_toSort(o_nToSort)%id = AGO2PGO(iplg(JA(i))-1)
! not a ghost node
else
nToSort = nToSort + 1
! petsc local node number to sort for
toSort(nToSort)%id = ALO2PLO(JA_P(i))
! store the old position in ASPAR
toSort(nToSort)%userData = i
end if
end do
call bubbleSort(toSort, nToSort)
call bubbleSort(o_toSort, o_nToSort)
! write the sorted cols to the mappings
do i = 1, nToSort
J = J + 1
JA_petsc(J) = toSort(i)%id
CSR_App2PetscLUT(toSort(i)%userData) = J
end do
IA_petsc(IP_petsc+1) = IA_petsc(IP_petsc) + nToSort
do i = 1, o_nToSort
o_J = o_J + 1
oJA_petsc(o_J) = o_toSort(i)%id
o_CSR_App2PetscLUT(o_toSort(i)%userData) = o_J
end do
oIA_petsc(IP_petsc+1) = oIA_petsc(IP_petsc) + o_nToSort
end do
deallocate(toSort, o_toSort, stat=istat)
if(istat /= 0) then
write(DBG%FHNDL,*) __FILE__, " Line", __LINE__
stop 'wwm_petsc_parallel l.250'
endif
end SUBROUTINE
!> fill matrix, RHs, call solver
SUBROUTINE EIMPS_PETSC_PARALLEL(ISS, IDD)
USE DATAPOOL
use petscpool
use petscsys
use petscmat
! use petscvec
! use mpi
implicit none
integer, intent(in) :: ISS, IDD
integer :: I, J
integer :: IP, IPGL1, IE, POS
integer :: I1, I2, I3
integer :: POS_TRICK(3,2)
real(kind=8) :: DTK, TMP3
real(kind=8) :: LAMBDA(2), GTEMP1, GTEMP2, DELT, XIMP, DELFL, DELT5
real(kind=8) :: FL11, FL12, FL21, FL22, FL31, FL32, USFM, TEMP, FLHAB
real(kind=8) :: CRFS(3), K1, KM(3), K(3), TRIA03
real(kind=8) :: DELTAL(3,MNE)
real(kind=8) :: KP(3,MNE), NM(MNE)
real(kind=8) :: U(MNP), C(2,MNP)
real(kind=8) :: X(MNP)
real(kind=8) :: B(MNP)
real(kind=8) :: ASPAR(NNZ)
! solver timings
#ifdef TIMINGS
real(rkind) :: startTime, endTime
#endif
real, save :: solverTimeSum = 0
!
! Petsc stuff
!
PetscInt :: ncols
PetscInt :: eCol
PetscScalar :: eEntry
integer :: counter
KSPConvergedReason reason;
! solver iteration
PetscInt iteration
integer, save :: iterationSum = 0
call PetscLogStagePush(stageFill, petscErr);CHKERRQ(petscErr)
iteration = 0
if(ISS == 0 .and. IDD == 0) then
iterationSum = 0
solverTimeSum = 0
endif
POS_TRICK(1,1) = 2
POS_TRICK(1,2) = 3
POS_TRICK(2,1) = 3
POS_TRICK(2,2) = 1
POS_TRICK(3,1) = 1
POS_TRICK(3,2) = 2
CALL CADVXY(ISS,IDD,C)
!
! Calculate countour integral quantities ...
!
DO IE = 1, MNE
I1 = INE(1,IE)
I2 = INE(2,IE)
I3 = INE(3,IE)
LAMBDA(1) = ONESIXTH * (C(1,I1)+C(1,I2)+C(1,I3))
LAMBDA(2) = ONESIXTH * (C(2,I1)+C(2,I2)+C(2,I3))
K(1) = LAMBDA(1) * IEN(1,IE) + LAMBDA(2) * IEN(2,IE)
K(2) = LAMBDA(1) * IEN(3,IE) + LAMBDA(2) * IEN(4,IE)
K(3) = LAMBDA(1) * IEN(5,IE) + LAMBDA(2) * IEN(6,IE)
KP(1,IE) = MAX(0.0_rkind,K(1))
KP(2,IE) = MAX(0.0_rkind,K(2))
KP(3,IE) = MAX(0.0_rkind,K(3))
KM(1) = MIN(0.0_rkind,K(1))
KM(2) = MIN(0.0_rkind,K(2))
KM(3) = MIN(0.0_rkind,K(3))
FL11 = C(1,I2)*IEN(1,IE)+C(2,I2)*IEN(2,IE)
FL12 = C(1,I3)*IEN(1,IE)+C(2,I3)*IEN(2,IE)
FL21 = C(1,I3)*IEN(3,IE)+C(2,I3)*IEN(4,IE)
FL22 = C(1,I1)*IEN(3,IE)+C(2,I1)*IEN(4,IE)
FL31 = C(1,I1)*IEN(5,IE)+C(2,I1)*IEN(6,IE)
FL32 = C(1,I2)*IEN(5,IE)+C(2,I2)*IEN(6,IE)
CRFS(1) = - ONESIXTH * (TWO * FL31 &
& + FL32 &
& + FL21 &
& + TWO * FL22 )
CRFS(2) = - ONESIXTH * (TWO * FL32 &
& + TWO * FL11 &
& + FL12 &
& + FL31 )
CRFS(3) = - ONESIXTH * (TWO * FL12 &
& + TWO * FL21 &
& + FL22 &
& + FL11 )
DELTAL(:,IE) = CRFS(:)-KP(:,IE)
NM(IE) = 1.0_rkind/MIN(-THR,SUM(KM(:)))
END DO
U(:) = AC2(ISS,IDD,:)
J = 0 ! Counter ...
ASPAR = ZERO ! Mass matrix ...
B = ZERO ! Right hand side ...
!
! ... assembling the linear equation system ....
!
DO IP = 1, MNP
DO I = 1, CCON(IP)
J = J + 1
IE = IE_CELL(J)
POS = POS_CELL(J)
K1 = KP(POS,IE) ! Flux Jacobian
TRIA03 = ONETHIRD * TRIA(IE)
DTK = K1 * DT4A * IOBPD(IDD,IP) * IOBWB(IP) * IOBDP(IP)
TMP3 = DTK * NM(IE)
I1 = POSI(1,J) ! Position of the recent entry in the ASPAR matrix ... ASPAR is shown in fig. 42, p.122
I2 = POSI(2,J)
I3 = POSI(3,J)
ASPAR(I1) = TRIA03 + DTK - TMP3 * DELTAL(POS ,IE) + ASPAR(I1) ! Diagonal entry
ASPAR(I2) = - TMP3 * DELTAL(POS_TRICK(POS,1),IE) + ASPAR(I2) ! off diagonal entries ...
ASPAR(I3) = - TMP3 * DELTAL(POS_TRICK(POS,2),IE) + ASPAR(I3)
B(IP) = B(IP) + TRIA03 * U(IP)
END DO !I: loop over connected elements ...
END DO !IP
IF (LBCWA .OR. LBCSP) THEN
IF (LINHOM) THEN
DO IP = 1, IWBMNP
IPGL1 = IWBNDLC(IP)
ASPAR(I_DIAG(IPGL1)) = SI(IPGL1) ! Set boundary on the diagonal
B(IPGL1) = SI(IPGL1) * WBAC(ISS,IDD,IP)
END DO
ELSE
DO IP = 1, IWBMNP
IPGL1 = IWBNDLC(IP)
ASPAR(I_DIAG(IPGL1)) = SI(IPGL1) ! Set boundary on the diagonal
B(IPGL1) = SI(IPGL1) * WBAC(ISS,IDD,1)
END DO
ENDIF
END IF
IF (ICOMP .GE. 2 .AND. SMETHOD .GT. 0 .AND. LSOURCESWAM) THEN
DO IP = 1, MNP
IF (IOBWB(IP) .EQ. 1) THEN
!GTEMP1 = MAX((1.-DT4A*FL(IP,ID,IS)),1.)
!GTEMP2 = SL(IP,ID,IS)/GTEMP1/PI2/SPSIG(IS)
GTEMP1 = MAX((1.-DT4A*IMATDAA(IP,ISS,IDD)),1.)
GTEMP2 = IMATRAA(IP,ISS,IDD)/GTEMP1!/PI2/SPSIG(IS)
DELT = DT4S
XIMP = 1.0
DELT5 = XIMP*DELT
DELFL = COFRM4(ISS)*DELT
USFM = USNEW(IP)*MAX(FMEANWS(IP),FMEAN(IP))
TEMP = USFM*DELFL!/PI2/SPSIG(IS)
FLHAB = ABS(GTEMP2*DT4S)
FLHAB = MIN(FLHAB,TEMP)/DT4S
B(IP) = B(IP) + SIGN(FLHAB,GTEMP2) * DT4A * SI(IP) ! Add source term to the right hand side
ASPAR(I_DIAG(IP)) = ASPAR(I_DIAG(IP)) - SIGN(FLHAB,GTEMP2) * SI(IP)
!!B(IP) = B(IP) + GTEMP2 * DT4A * SI(IP) ! Add source term to the right hand side
!ASPAR(I_DIAG(IP)) = ASPAR(I_DIAG(IP)) - GTEMP2 * SI(IP)
!This is then for the shallow water physics take care about ISELECT
!ASPAR(I_DIAG(IP)) = ASPAR(I_DIAG(IP)) + IMATDAA(IP,IS,ID) * DT4A * SI(IP) ! Add source term to the diagonal
!B(IP) = B(IP) + IMATRAA(IP,IS,ID) * DT4A * SI(IP) ! Add source term to the right hand side
ENDIF
END DO
ELSE IF (ICOMP .GE. 2 .AND. SMETHOD .GT. 0 .AND. .NOT. LSOURCESWAM) THEN
DO IP = 1, MNP
IF (IOBWB(IP) .EQ. 1) THEN
ASPAR(I_DIAG(IP)) = ASPAR(I_DIAG(IP)) + IMATDAA(IP,ISS,IDD) * DT4A * SI(IP) ! Add source term to the diagonal
B(IP) = B(IP) + IMATRAA(IP,ISS,IDD) * DT4A * SI(IP) ! Add source term to the right hand side
ENDIF
END DO
ENDIF
! call checkAsparDiagonalAccuracy(ASPAR, IA, JA, ISS, IDD)
! fill the new matrix
ASPAR_petsc = 0
oASPAR_petsc = 0
counter = 1
ncols = 0
do i = 1, NP_RES
ncols = IA_P(i+1) - IA_P(i)
! this is a interface node (row). ignore it. just increase counter
if(ALOold2ALO(i) .eq. -999) then
counter = counter + ncols
cycle
end if
! insert col by col into matrix
do j = 1, ncols
if(CSR_App2PetscLUT(counter) == -999) then
oASPAR_petsc(o_CSR_App2PetscLUT(counter)) = ASPAR(counter)
else
ASPAR_petsc(CSR_App2PetscLUT(counter)) = ASPAR(counter)
endif
counter = counter + 1
end do
end do
call MatAssemblyBegin(matrix, MAT_FINAL_ASSEMBLY, petscErr)
CHKERRQ(petscErr)
call MatAssemblyEnd(matrix, MAT_FINAL_ASSEMBLY, petscErr)
CHKERRQ(petscErr)
!fill RHS vector
!iterate over all resident (and interface) nodes
!map it to petsc global ordering
!and insert the value from B into RHS vector
eEntry = 0;
call VecSet(myB, eEntry, petscErr);CHKERRQ(petscErr)
! call VecSet(myBAppOrder, eEntry, petscErr);CHKERRQ(petscErr)
do i= 1, np
! this is a interface node (row). ignore it. just increase counter
if(ALOold2ALO(i) .eq. -999) then
cycle
end if
! map to petsc global order
eCol = AGO2PGO(iplg(i) - 1 )
eEntry = B(i)
call VecSetValue(myB, eCol, eEntry, ADD_VALUES, petscErr)
CHKERRQ(petscErr)
end do
call VecAssemblyBegin(myB, petscErr);CHKERRQ(petscErr);
call VecAssemblyEnd(myB, petscErr);CHKERRQ(petscErr);
! Copy the old solution from AC2 to myX to make the solver faster
do i = 1, np
eCol = AGO2PGO(iplg(i)-1)
eEntry = AC2(ISS, IDD, i)
call VecSetValue(myX,eCol,eEntry,INSERT_VALUES,petscErr)
CHKERRQ(petscErr)
end do
call VecAssemblyBegin(myX, petscErr);CHKERRQ(petscErr);
call VecAssemblyEnd(myX, petscErr);CHKERRQ(petscErr);
! Solve
! To solve successive linear systems that have different preconditioner matrices (i.e., the matrix elements
! and/or the matrix data structure change), the user must call KSPSetOperators() and KSPSolve() for each
! solve.
if(samePreconditioner .eqv. .true.) call KSPSetOperators(Solver, matrix, matrix, SAME_PRECONDITIONER, petscErr);CHKERRQ(petscErr)
call PetscLogStagePop(petscErr);CHKERRQ(petscErr)
call PetscLogStagePush(stageSolve, petscErr);CHKERRQ(petscErr)
#ifdef TIMINGS
call WAV_MY_WTIME(startTime)
#endif
! Solve!
call KSPSolve(Solver, myB, myX, petscErr);CHKERRQ(petscErr);
#ifdef TIMINGS
call WAV_MY_WTIME(endTime)
#endif
call PetscLogStagePop(petscErr);CHKERRQ(petscErr)
call KSPGetConvergedReason(Solver, reason, petscErr);CHKERRQ(petscErr);
if (reason .LT. 0) then
!CALL WWM_ABORT('Failure to converge')
!write(stat%fhndl,*) 'Failure to converge'
endif
#ifdef PETSC_DEBUG
if(rank == 0) then
if(reason .LT. 0 ) then
write(DBG%FHNDL,*) "Failure to converge\n"
else
call KSPGetIterationNumber(Solver, iteration, petscErr)
CHKERRQ(petscErr)
! print only the mean number of iteration
iterationSum = iterationSum + iteration
solverTimeSum = solverTimeSum + (endTime - startTime)
if(ISS == MSC .and. IDD == MDC) then
write(DBG%FHNDL,*) "mean number of iterations", iterationSum / real((MSC*MDC))
print '("solver Time for all MSD MDC= ",f6.3," sec")', solverTimeSum
endif
endif
endif
#endif
X = 0.0_rkind
!get the solution back to fortran.
!iterate over all resident nodes (without interface and ghost nodes)
!map the solution from petsc local ordering back to app old local ordering
!(the app old ordering contains interface nodes)
call VecGetArrayF90(myX, myXtemp, petscErr); CHKERRQ(petscErr)
do i = 1, nNodesWithoutInterfaceGhosts
X(ipgl((PGO2AGO(PLO2PGO(i-1)))+1)%id) = myXtemp(i)
end do
call VecRestoreArrayF90(myX, myXtemp, petscErr)
CHKERRQ(petscErr);
!IF (SUM(X) .NE. SUM(X)) CALL WWM_ABORT('NaN in X')
! we have to fill the ghost and interface nodes with the solution from the other threads
! at least SUBROUTINE SOURCETERMS() make calculations on interface/ghost nodes which are
! normally set to 0, because they do net exist in petsc
call exchange_p2d(X)
AC2(ISS, IDD,:) = MAX(0.0_rkind,X)
END SUBROUTINE
!> cleanup memory. You never need to call this function by hand. It will automaticly called by PETSC_FINALIZE()
SUBROUTINE PETSC_FINALIZE_PARALLEL()
implicit none
! we deallocate only arrays who are declared in this file!
if(allocated(IA_petsc)) deallocate(IA_petsc)
if(allocated(JA_petsc)) deallocate(JA_petsc)
if(allocated(ASPAR_petsc)) deallocate(ASPAR_petsc)
if(allocated(oIA_petsc)) deallocate(oIA_petsc)
if(allocated(oJA_petsc)) deallocate(oJA_petsc)
if(allocated(oASPAR_petsc)) deallocate(oASPAR_petsc)
if(allocated(CSR_App2PetscLUT)) deallocate(CSR_App2PetscLUT)
if(allocated(o_CSR_App2PetscLUT)) deallocate(o_CSR_App2PetscLUT)
end SUBROUTINE
END MODULE
# endif
#endif