Newer
Older
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
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! |
! ||===\\ |
! || \\ |
! || || //==\\ || || //==|| ||/==\\ |
! || || || || || || || || || || |
! || // || || || || || || || |
! ||===// \\==// \\==\\ \\==\\ || |
! |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! |
! BUILD_SYSTEM_WSMP & IFIND_LOC OCT. 2008 |
! |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
subroutine build_system_wsmp (n,n_iproc,n_iproc_st,nz_loc,iproc_col,ia,ja, &
ldb,nrhs,avals,b,params,osolve,ndof,mat,vo, &
istep,iter,iter_nl,threadinfo,weightel_glob)
!------------------------------------------------------------------------------|
!(((((((((((((((( Purpose of the routine ))))))))))))))))))))))))))))))))))))))
!------------------------------------------------------------------------------|
! This routine is used to build the information on the system of equations
! and the elemental matrices that will be sent to the solver.
! It first uses the make_cut routine to build the fe matrices and the rhs
! vector taking into account cut cells
! Note that parallelism of the solver also reauires parallelism of the
! matrix building phase. In other words not all elements are built in all
! processors. However because the problem is divided in columns, that is in
! degrees of freedom, it is necessary that some elements are computed by more
! than one processor.
! In the current version, build_system is called (with ndof=3) to form the viscous
! stiffness matrix and (with ndof=1) to form the temperature stiffness matrix.
!------------------------------------------------------------------------------|
!(((((((((((((((( declaration of the subroutine arguments ))))))))))))))))))))
!------------------------------------------------------------------------------|
use threads
use definitions
implicit none
integer n,nz_loc,n_iproc,n_iproc_st
logical iproc_col(n)
integer ia(n_iproc+1)
integer ja(nz_loc)
integer ldb,nrhs
double precision avals(nz_loc)
double precision b(ldb,nrhs)
type(parameters) params
type(octreesolve) osolve
integer ndof
type (material) mat(0:params%nmat)
type (void) vo
integer istep
integer iter
integer iter_nl
type (thread) threadinfo
double precision weightel_glob(osolve%nleaves)
!------------------------------------------------------------------------------|
!(((((((((((((((( declaration of the subroutine internal variables )))))))))))))
!------------------------------------------------------------------------------|
integer, external :: ifind_loc_wsmp
double precision,dimension(:,:),allocatable::lsf_el
double precision ael(8*ndof,8*ndof),bel(8*ndof),patch1(5,5),patch2(3,3)
double precision,dimension(:),allocatable::penalty_loc,penalty_glob,eviscosity,weightel
double precision r0,s0,t0,rst,amin
integer i,j,k,itodo,idof,i1,i2,iii,jjj,err,jj,i1loc,k1,k2,ii
integer ileaves,icut,iloc,idof1,idof2,nelem_iproc,nelem_min,nelem_max
integer inode(5)
integer levelmax,level
integer iproc,nproc,ierr
double precision rhs(osolve%nnode*ndof)
character*72 :: shift
logical, dimension(:),allocatable :: is_plastic,vbounded,vbounded2
integer, dimension(:), allocatable :: kfix
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
INCLUDE 'mpif.h'
call mpi_comm_size (mpi_comm_world,nproc,ierr)
call mpi_comm_rank (mpi_comm_world,iproc,ierr)
data patch1 / 1., 1., 1., 1.,-4., &
1., 1., 1., 1.,-4., &
1., 1., 1., 1.,-4., &
1., 1., 1., 1.,-4., &
-4.,-4.,-4.,-4.,16./
data patch2 / 1.,-2., 1., &
-2., 4.,-2., &
1.,-2., 1./
if (ndof==3) then
allocate(kfix(osolve%nnode*3))
kfix=osolve%kfix
elseif (ndof==1) then
allocate(kfix(osolve%nnode))
kfix=osolve%kfixt
end if
allocate (is_plastic(osolve%nleaves),stat=threadinfo%err) ; call heap (threadinfo,'is_plastic','build_system',size(is_plastic),'bool',+1)
allocate (vbounded(osolve%nleaves),stat=threadinfo%err) ; call heap (threadinfo,'vbounded','build_system',size(vbounded),'bool',+1)
allocate (vbounded2(osolve%nleaves),stat=threadinfo%err) ; call heap (threadinfo,'vbounded','build_system',size(vbounded2),'bool',+1)
allocate (eviscosity(osolve%nleaves),stat=threadinfo%err) ; call heap (threadinfo,'eviscosity','build_system',size(eviscosity),'dp',+1)
allocate (weightel(osolve%nleaves),stat=threadinfo%err) ; call heap (threadinfo,'weightel','build_system',size(weightel),'dp',+1)
shift=' '
is_plastic=.false.
vbounded=.false.
vbounded2=.false.
eviscosity=0.d0
avals=0.d0
b=0.d0
nelem_iproc=0
weightel=0.d0
allocate (lsf_el(params%mpe,osolve%nlsf),stat=err) ; if (err.ne.0) call stop_run ('Error alloc lsf_el in build_system$')
amin=1.d6
do ileaves=1,osolve%nleaves
! first test if in the void
if (vo%leaf(ileaves).eq.0) then
! then test if this leaf touches the colums given to the current proc
itodo=0
do k=1,params%mpe
do idof=1,ndof
i1=(vo%ftr(osolve%icon(k,ileaves))-1)*ndof+idof
if (iproc_col(i1)) itodo=1
enddo
enddo
! what follows is only performed for a subset of elements
if (itodo.eq.1) then
nelem_iproc=nelem_iproc+1
level=0
levelmax=params%levelcut
do k=1,params%mpe
do i=1,osolve%nlsf
lsf_el(k,i)=osolve%lsf(osolve%icon(k,ileaves),i)
enddo
enddo
r0=-1.d0
s0=-1.d0
t0=-1.d0
rst=2.d0
! make_cut is used to build the elemental matrix taking into
! account the cut cells
call make_cut (level,levelmax,ndof,ael,bel,osolve%icon(1,ileaves), &
osolve%x,osolve%y,osolve%z,kfix,mat, &
osolve%u,osolve%v,osolve%w,osolve%temp, &
osolve%spressure(ileaves),osolve%strain,is_plastic(ileaves),&
osolve%nnode,rhs,lsf_el,osolve%nlsf,r0,s0,t0,rst,icut, &
ileaves,eviscosity(ileaves),vbounded(ileaves),params,threadinfo, &
weightel(ileaves))
! transfer elemental matrix and rhs into wsmpmatrix structure
! different for symmetric or nonsymmetric matrices
do k1=1,params%mpe
do idof1=1,ndof
i1=(vo%ftr(osolve%icon(k1,ileaves))-1)*ndof+idof1
if (iproc_col(i1)) then
iii=(k1-1)*ndof+idof1
i1loc=i1-n_iproc_st+1
do k2=1,params%mpe
do idof2=1,ndof
i2=(vo%ftr(osolve%icon(k2,ileaves))-1)*ndof+idof2
jjj=(k2-1)*ndof+idof2
if (i2.ge.i1.or.ndof.eq.1) then
iloc=ifind_loc_wsmp(i1loc,i2,ja,ia,n_iproc,nz_loc)
avals(iloc)=avals(iloc)+ael(jjj,iii)
if (iii.eq.jjj) amin=min(amin,abs(ael(jjj,iii)))
endif
enddo
enddo
b(i1loc,1)=b(i1loc,1)+bel(iii)
endif
enddo
enddo
endif
endif
enddo
!print*,'check diag',iproc,amin
deallocate(lsf_el)
deallocate(kfix)
call mpi_reduce(nelem_iproc,nelem_min,1,mpi_integer,mpi_min,0,mpi_comm_world,ierr)
call mpi_reduce(nelem_iproc,nelem_max,1,mpi_integer,mpi_max,0,mpi_comm_world,ierr)
if (iproc.eq.0) write (*,'(a,2i12)') shift//'nelem per proc (min/max)',nelem_min,nelem_max
if (ndof==3) then
osolve%eviscosity=0.d0
call mpi_allreduce (eviscosity,osolve%eviscosity,osolve%nleaves,mpi_double_precision,mpi_max,mpi_comm_world,ierr)
if (iproc.eq.0) write (*,'(a,2es13.5)') shift//'viscosity range ', minval(osolve%eviscosity,osolve%eviscosity>0.d0), maxval(osolve%eviscosity)
osolve%is_plastic=.false.
call mpi_allreduce (is_plastic,osolve%is_plastic,osolve%nleaves,mpi_logical,mpi_lor,mpi_comm_world,ierr)
call mpi_allreduce (vbounded,vbounded2,osolve%nleaves,mpi_logical,mpi_lor,mpi_comm_world,ierr)
if (iproc.eq.0) write (*,'(a,i10)') shift//'viscosity capped in', count(vbounded2)
call mpi_allreduce (weightel,weightel_glob,osolve%nleaves,mpi_double_precision,mpi_min,mpi_comm_world,ierr)
end if
call heap (threadinfo,'is_plastic','build_system',size(is_plastic),'bool',-1) ; deallocate (is_plastic)
call heap (threadinfo,'vbounded','build_system',size(vbounded),'bool',-1) ; deallocate (vbounded)
call heap (threadinfo,'vbounded2','build_system',size(vbounded2),'bool',-1) ; deallocate (vbounded2)
call heap (threadinfo,'eviscosity','build_system',size(eviscosity),'dp',-1) ; deallocate (eviscosity)
call heap (threadinfo,'weightel','build_system',size(weightel),'dp',-1) ; deallocate (weightel)
!=====================
!=====[bad faces]=====
!=====================
if (osolve%nface.eq.0) return
allocate (penalty_loc(osolve%nface),stat=threadinfo%err) ; call heap (threadinfo,'penalty_loc','build_system',size(penalty_loc),'dp',+1)
allocate (penalty_glob(osolve%nface),stat=threadinfo%err) ; call heap (threadinfo,'penalty_glob','build_system',size(penalty_glob),'dp',+1)
if (params%penalty.gt.0.d0) then
! first case assumes that a unique, uniform value is used
! for the penalty parameter
penalty_glob=params%penalty
else
! second case assumes that the penalty parameter is scaled according to
! the value of the stiffness matrix
! first go around the values of the stiffness matrix that are going to be
! touched by the linear constraints
! and deduce from it the value of a local penalty parameter (penalty_loc)
! penalty_loc=0.d0
! do i=1,osolve%nface
! if (vo%face(i).eq.0) then
! ! First the 4 constraints along the edges
! penalty_loc(i)=0.d0
! do j=1,4
! inode(1)=vo%ftr(osolve%iface(j,i))
! inode(2)=vo%ftr(osolve%iface(j+4,i))
! inode(3)=vo%ftr(osolve%iface(1+mod(j,4),i))
! do idof=1,ndof
! do ii=1,3
! i1=(inode(ii)-1)*ndof+idof
! if (wsmp%iproc_col(i1)) then
! do jj=1,3
! i2=(inode(jj)-1)*ndof+idof
! if (i1.ge.i2.or.ndof.eq.1) then
! iloc=ifind_loc_wsmp(i1,i2,ja,ia,vo%nnode*ndof,nz_loc)
! penalty_loc(i)=max(penalty_loc(i),abs(wsmp%avals(iloc)))
! endif
! enddo
! endif
! enddo
! enddo
! enddo
! ! then the other constraints across the four corners and the centre
! inode(1)=vo%ftr(osolve%iface(1,i))
! inode(2)=vo%ftr(osolve%iface(2,i))
! inode(3)=vo%ftr(osolve%iface(3,i))
! inode(4)=vo%ftr(osolve%iface(4,i))
! inode(5)=vo%ftr(osolve%iface(9,i))
! do k=1,ndof
! do ii=1,5
! i1=(inode(ii)-1)*ndof+k
! if (wsmp%iproc_col(i1)) then
! do jj=1,5
! i2=(inode(jj)-1)*ndof+k
! if (i1.ge.i2.or.ndof.eq.1) then
! iloc=ifind_loc_wsmp(i1,i2,ja,ia,vo%nnode*ndof,nz_loc)
! penalty_loc(i)=max(penalty_loc(i),abs(wsmp%avals(iloc)))
! endif
! enddo
! endif
! enddo
! enddo
! endif
! enddo
! applies an additional factor 10^5 to the resulting local penalty parameter
! penalty_loc=penalty_loc*1.d5
! call mpi_allreduce (penalty_loc,penalty_glob,osolve%nface,mpi_double_precision,mpi_max,mpi_comm_world,ierr)
endif
! now, in both cases, applies the resulting constraints by multiplying the
! stiffness matrix by a pattern (patch1 et patch2) scaled by the penalty parameter
do i=1,osolve%nface
if (vo%face(i).eq.0) then
! First the 4 constraints along the edges
do j=1,4
inode(1)=vo%ftr(osolve%iface(j,i))
inode(2)=vo%ftr(osolve%iface(j+4,i))
inode(3)=vo%ftr(osolve%iface(1+mod(j,4),i))
do k=1,ndof
do ii=1,3
i1=(inode(ii)-1)*ndof+k
i1loc=i1-n_iproc_st+1
if (iproc_col(i1)) then
do jj=1,3
i2=(inode(jj)-1)*ndof+k
if (i2.ge.i1.or.ndof.eq.1) then
iloc=ifind_loc_wsmp(i1loc,i2,ja,ia,n_iproc,nz_loc)
avals(iloc)=avals(iloc)+patch2(ii,jj)*penalty_glob(i)
endif
enddo
endif
enddo
enddo
enddo
! then the other constraints across the four corners and the centre
inode(1)=vo%ftr(osolve%iface(1,i))
inode(2)=vo%ftr(osolve%iface(2,i))
inode(3)=vo%ftr(osolve%iface(3,i))
inode(4)=vo%ftr(osolve%iface(4,i))
inode(5)=vo%ftr(osolve%iface(9,i))
do idof=1,ndof
do ii=1,5
i1=(inode(ii)-1)*ndof+idof
i1loc=i1-n_iproc_st+1
if (iproc_col(i1)) then
do jj=1,5
i2=(inode(jj)-1)*ndof+idof
if (i2.ge.i1.or.ndof.eq.1) then
iloc=ifind_loc_wsmp(i1loc,i2,ja,ia,n_iproc,nz_loc)
avals(iloc)=avals(iloc)+patch1(ii,jj)*penalty_glob(i)
endif
enddo
endif
enddo
enddo
endif
enddo
call heap (threadinfo,'penalty_glob','build_system',size(penalty_glob),'dp',-1) ; deallocate (penalty_glob)
call heap (threadinfo,'penalty_loc','build_system',size(penalty_loc),'dp',-1) ; deallocate (penalty_loc)
end subroutine build_system_wsmp
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
! function used to navigate through the global matrix
! i and j are the global indices in the global matrix
! irn and jcn are needed by MUMPS to
! locate the nonzero numbers of the global stiffness matrix
! n is the total number of dofs
! nz is the total number of nonzero numbers in the global matrix
! there is an internal check that could be removed now that this
! routine and the storage algorithm developed have been fully tested
! but we will leave it in for the time being
integer function ifind_loc_wsmp (i1,i2,ja,ia,n_iproc,nz_loc)
implicit none
integer i1,i2
integer ia(n_iproc+1),ja(nz_loc)
integer n_iproc,nz_loc
integer istart,k
istart=ia(i1)
1 continue
if (ja(istart).eq.i2) then
ifind_loc_wsmp=istart
return
endif
istart=istart+1
if (istart.le.nz_loc) goto 1
!print*,'istart=',istart
!print*,'nz',nz
!print*,'i',i
!print*,'idg(i)',idg(i)
!print*,'j',j
!print*,(ja(k),k=idg(i),idg(i+1)-1)
STOP 'error in ifind_loc_wsmp'
end function ifind_loc_wsmp
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|