Rocstar  1.0
Rocstar multiphysics simulation application
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
createM.f90
Go to the documentation of this file.
1 !*********************************************************************
2 !* Illinois Open Source License *
3 !* *
4 !* University of Illinois/NCSA *
5 !* Open Source License *
6 !* *
7 !* Copyright@2008, University of Illinois. All rights reserved. *
8 !* *
9 !* Developed by: *
10 !* *
11 !* Center for Simulation of Advanced Rockets *
12 !* *
13 !* University of Illinois *
14 !* *
15 !* www.csar.uiuc.edu *
16 !* *
17 !* Permission is hereby granted, free of charge, to any person *
18 !* obtaining a copy of this software and associated documentation *
19 !* files (the "Software"), to deal with the Software without *
20 !* restriction, including without limitation the rights to use, *
21 !* copy, modify, merge, publish, distribute, sublicense, and/or *
22 !* sell copies of the Software, and to permit persons to whom the *
23 !* Software is furnished to do so, subject to the following *
24 !* conditions: *
25 !* *
26 !* *
27 !* @ Redistributions of source code must retain the above copyright *
28 !* notice, this list of conditions and the following disclaimers. *
29 !* *
30 !* @ Redistributions in binary form must reproduce the above *
31 !* copyright notice, this list of conditions and the following *
32 !* disclaimers in the documentation and/or other materials *
33 !* provided with the distribution. *
34 !* *
35 !* @ Neither the names of the Center for Simulation of Advanced *
36 !* Rockets, the University of Illinois, nor the names of its *
37 !* contributors may be used to endorse or promote products derived *
38 !* from this Software without specific prior written permission. *
39 !* *
40 !* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, *
41 !* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES *
42 !* OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND *
43 !* NONINFRINGEMENT. IN NO EVENT SHALL THE CONTRIBUTORS OR *
44 !* COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER *
45 !* LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, *
46 !* ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE *
47 !* USE OR OTHER DEALINGS WITH THE SOFTWARE. *
48 !*********************************************************************
49 !* Please acknowledge The University of Illinois Center for *
50 !* Simulation of Advanced Rockets in works and publications *
51 !* resulting from this software or its derivatives. *
52 !*********************************************************************
53 SUBROUTINE createm(global)
54 
55  USE implicit_global
56  USE comp_row_global
57  USE rocstar_rocfraccomm
58 
59  IMPLICIT NONE
60 
61  include 'mpif.h'
62 
63  TYPE(rocfrac_global) :: global
64 
65  INTEGER :: i, j, m, p, inode, jnode, counter, idof, jdof
66 
67  REAL(kind=wp) :: tempkval
68  INTEGER, ALLOCATABLE, DIMENSION(:) :: amounttoreceive, amounttosend
69 
70  REAL(kind=wp), ALLOCATABLE, DIMENSION(:) :: bufsnd
71 
72  REAL(kind=wp),ALLOCATABLE,DIMENSION(:) :: tempmg, mg
73  REAL(kind=wp),ALLOCATABLE,DIMENSION(:) :: aval_m_temp
74  INTEGER,ALLOCATABLE,DIMENSION(:) :: rp_m_temp, cval_m_temp
75 
76  REAL(kind=wp), DIMENSION(1:3,1:8) :: ri = reshape( &
77  (/-0.577350269189626,-0.577350269189626,-0.577350269189626, &
78  0.577350269189626,-0.577350269189626,-0.577350269189626, &
79  0.577350269189626, 0.577350269189626,-0.577350269189626, &
80  -0.577350269189626, 0.577350269189626,-0.577350269189626, &
81  -0.577350269189626,-0.577350269189626, 0.577350269189626, &
82  0.577350269189626,-0.577350269189626, 0.577350269189626, &
83  0.577350269189626, 0.577350269189626, 0.577350269189626, &
84  -0.577350269189626, 0.577350269189626, 0.577350269189626/),(/3,8/) )
85 
86 
87 
88  ! Construct the global mass matrix
89  IF(myid==0) print*,'CONSTRUCTING THE MASS MATRIX'
90 
91  ! Construct the mass matrix as a consistent matrix
92  CALL implicit_v3d8_mass_consistent(global%NumElVol,global%NumNp,global%NumMatVol,global%MeshCoor,global%ElConnVol,global%MatIdVol,ri,global%rho,global%ElConnVol)
93  nnz_m = nnz_temp
94  ALLOCATE(cval_m(1:nnz_m))
95  ALLOCATE(aval_m(1:nnz_m))
96  ALLOCATE(rp_m(1:3*gnumnp+1))
97  rp_m = rp_temp
98  cval_m = cval_temp
99  aval_m = aval_temp
100  DEALLOCATE(rp_temp)
101  DEALLOCATE(cval_temp)
102  DEALLOCATE(aval_temp)
103 
104  ! Communicate mass to the other procs
105  IF (nprocs > 1) THEN
106  CALL mpi_barrier(rocstar_communicator,ierr)
107  ALLOCATE(req_rcv(1:nprocs))
108  ALLOCATE(req_snd(1:nprocs))
109  ALLOCATE(stat_snd(1:mpi_status_size,1:nprocs))
110  ALLOCATE(stat_rcv(1:mpi_status_size,1:nprocs))
111  DO i = 1, nprocs
112  req_rcv(i) = 0
113  req_snd(i) = 0
114  DO j = 1, mpi_status_size
115  stat_rcv(1,i) = 0
116  stat_snd(j,i) = 0
117  ENDDO
118  ENDDO
119  ALLOCATE(amounttoreceive(1:numcommprocsfrom1))
120  DO i = 1, numcommprocsfrom1
121  CALL mpi_irecv(amounttoreceive(i),1, &
122  mpi_integer,commprocsfrom1(i),10,rocstar_communicator, &
123  req_rcv(i),ierr)
124  ENDDO
125  ALLOCATE(amounttosend(1:numcommprocs1))
126  DO i = 1, numcommprocs1
127  counter = 0
128  DO j = 1, numcommnodes1(i)
129  DO p = 1, 3
130  DO m = 1, 3*gnumnp
131  tempkval = 0.0
132  inode = local2global(commnodes1(i,j))
133  jnode = int((m-0.5)/3)+1
134  idof = p
135  jdof = m - 3*jnode + 3
136  CALL comp_row_getval(3*gnumnp,3*gnumnp,nnz_m,1,rp_m,cval_m,aval_m, &
137  inode*3-3+idof,jnode*3-3+jdof,tempkval)
138  IF (tempkval /= 0.0) THEN
139  counter = counter + 1
140  !print*,myid,Global2Local(inode)*3-3+idof,Global2Local(jnode)*3-3+jdof,tempKval
141  ENDIF
142  ENDDO
143  ENDDO
144  !print*,myid,' sending stuff about node ',Local2Global(CommNodes1(i,j)),' to ',CommProcs1(i),' : ',counter
145  ENDDO
146  amounttosend(i) = counter
147  !CALL MPI_ISEND(AmountToSend(i),1,MPI_INTEGER, &
148  ! CommProcs1(i),10,ROCSTAR_COMMUNICATOR,req_snd(i),ierr)
149  CALL mpi_send(amounttosend(i),1,mpi_integer, &
150  commprocs1(i),10,rocstar_communicator,ierr)
151  ENDDO
152  !CALL MPI_WAITALL(NumCommProcsFrom1,req_rcv,stat_rcv,ierr)
153  !CALL MPI_WAITALL(NumCommProcs1,req_snd,stat_snd,ierr)
154  CALL mpi_waitall(nprocs,req_rcv,stat_rcv,ierr)
155  !CALL MPI_WAITALL(nprocs,req_snd,stat_snd,ierr)
156 
157  ! Communicate parts of the M matrix to other procs
158  CALL mpi_barrier(rocstar_communicator,ierr)
159  ALLOCATE(frmproc(1:numcommprocsfrom1))
160  DO i = 1, numcommprocsfrom1
161  ALLOCATE(frmproc(i)%rcvbuf(1:3*amounttoreceive(i)))
162  CALL mpi_irecv(frmproc(i)%rcvbuf(1),3*amounttoreceive(i), &
163  mpi_double_precision,commprocsfrom1(i),10,rocstar_communicator, &
164  req_rcv(i),ierr)
165  ENDDO
166  DO i = 1, numcommprocs1
167  ALLOCATE(bufsnd(1:3*amounttosend(i)))
168  bufsnd(:) = 0.0
169  counter = 0
170  DO j = 1, numcommnodes1(i)
171  DO p = 1, 3
172  DO m = 1, 3*gnumnp
173  tempkval = 0.0
174  inode = local2global(commnodes1(i,j))
175  jnode = int((m-0.5)/3)+1
176  idof = p
177  jdof = m - 3*jnode + 3
178  CALL comp_row_getval(3*gnumnp,3*gnumnp,nnz_m,1,rp_m,cval_m,aval_m, &
179  inode*3-3+idof,jnode*3-3+jdof,tempkval)
180  IF (tempkval /= 0.0) THEN
181  counter = counter + 1
182  bufsnd(3*counter-2) = local2global(commnodes1(i,j))*3-3+p
183  bufsnd(3*counter-1) = m
184  bufsnd(3*counter) = tempkval
185  ENDIF
186  ENDDO
187  ENDDO
188  ENDDO
189  !print*,myid,' sending to ',CommProcs1(i),' : ',bufsnd(:)
190  !CALL MPI_ISEND(bufsnd,3*AmountToSend(i),MPI_DOUBLE_PRECISION, &
191  ! CommProcs1(i),10,ROCSTAR_COMMUNICATOR,req_snd(i),ierr)
192  CALL mpi_send(bufsnd,3*amounttosend(i),mpi_double_precision, &
193  commprocs1(i),10,rocstar_communicator,ierr)
194  DEALLOCATE(bufsnd)
195  ENDDO
196  !CALL MPI_WAITALL(NumCommProcsFrom1,req_rcv,stat_rcv,ierr)
197  !CALL MPI_WAITALL(NumCommProcs1,req_snd,stat_snd,ierr)
198  CALL mpi_waitall(nprocs,req_rcv,stat_rcv,ierr)
199  !CALL MPI_WAITALL(nprocs,req_snd,stat_snd,ierr)
200  DEALLOCATE(req_rcv)
201  DEALLOCATE(req_snd)
202  DEALLOCATE(stat_snd)
203  DEALLOCATE(stat_rcv)
204  ENDIF
205 
206  IF (nprocs > 1) THEN
207  DO i = 1, numcommprocsfrom1
208  !print*,myid,' received from ',CommProcsFrom1(i), ' : ',frmproc(i)%rcvbuf(:)
209  DO j = 1, amounttoreceive(i)
210  CALL comp_row_addval(3*gnumnp,3*gnumnp,nnz_m,1,rp_m,cval_m,aval_m, &
211  int(frmproc(i)%rcvbuf(3*j-2)),int(frmproc(i)%rcvbuf(3*j-1)),frmproc(i)%rcvbuf(3*j))
212  IF (nnz_temp /= nnz_m) THEN
213  DEALLOCATE(cval_m,aval_m)
214  ALLOCATE(cval_m(nnz_temp),aval_m(nnz_temp))
215  nnz_m = nnz_temp
216  rp_m = rp_temp
217  cval_m = cval_temp
218  ENDIF
219  aval_m = aval_temp
220  DEALLOCATE(rp_temp,cval_temp,aval_temp)
221  ENDDO
222  ENDDO
223  DEALLOCATE(frmproc)
224  DEALLOCATE(amounttoreceive)
225  DEALLOCATE(amounttosend)
226  ENDIF
227 
228  ! Resize the matrix to the size needed on this proc
229  CALL comp_row_resize(3*gnumnp,3*gnumnp,nnz_m,1,rp_m,cval_m,aval_m,nstart_km,nrows_km)
230  DEALLOCATE(rp_m,cval_m,aval_m)
231  ALLOCATE(rp_m(1:nrows_km+1))
232  ALLOCATE(cval_m(1:nnz_temp))
233  ALLOCATE(aval_m(1:nnz_temp))
234  nnz_m = nnz_temp
235  rp_m = rp_temp
236  cval_m = cval_temp
237  aval_m = aval_temp
238  DEALLOCATE(rp_temp,cval_temp,aval_temp)
239 
240 
241 
242 
243  ! Construct the mass matrix as lumped
244  ALLOCATE(tempmg(1:global%NumNp))
245  ALLOCATE(mg(1:3*lnumnp))
246  ALLOCATE(cval_m_temp(1:3*lnumnp))
247  ALLOCATE(aval_m_temp(1:3*lnumnp))
248  ALLOCATE(rp_m_temp(1:3*lnumnp+1))
249  tempmg(1:global%NumNp) = 0.0
250  CALL implicit_v3d8_mass(global%NumElVol,global%NumNp,global%NumMatVol,global%MeshCoor,global%ElConnVol,global%MatIdVol,ri,global%rho,tempmg,1,global%NumElVol)
251  DO i = 1, global%NumNp
252  tempmg(i) = 1.0 / tempmg(i) ! Get the mass matrix, not its inverse
253  ENDDO
254 
255  IF (nprocs > 1) THEN
256  CALL mpi_barrier(rocstar_communicator,ierr)
257  !ALLOCATE(req_rcv(1:NumCommProcsFrom1))
258  !ALLOCATE(req_snd(1:NumCommProcs1))
259  !ALLOCATE(stat_snd(1:MPI_STATUS_SIZE,1:NumCommProcs1))
260  !ALLOCATE(stat_rcv(1:MPI_STATUS_SIZE,1:NumCommProcsFrom1))
261  ALLOCATE(req_rcv(1:nprocs))
262  ALLOCATE(req_snd(1:nprocs))
263  ALLOCATE(stat_snd(1:mpi_status_size,1:nprocs))
264  ALLOCATE(stat_rcv(1:mpi_status_size,1:nprocs))
265  DO i = 1, nprocs
266  req_rcv(i) = 0
267  req_snd(i) = 0
268  DO j = 1, mpi_status_size
269  stat_rcv(1,i) = 0
270  stat_snd(j,i) = 0
271  ENDDO
272  ENDDO
273  ALLOCATE(frmproc(1:numcommprocsfrom1))
274  DO i = 1, numcommprocsfrom1
275  ALLOCATE(frmproc(i)%rcvbuf(1:2*numcommnodesfrom1(i)))
276  CALL mpi_irecv(frmproc(i)%rcvbuf(1),2*numcommnodesfrom1(i), &
277  mpi_double_precision,commprocsfrom1(i),10,rocstar_communicator, &
278  req_rcv(i),ierr)
279  ENDDO
280  DO i = 1, numcommprocs1
281  ALLOCATE(bufsnd(2*numcommnodes1(i)))
282  DO j = 1, numcommnodes1(i)
283  bufsnd(2*j-1) = local2global(commnodes1(i,j))
284  bufsnd(2*j) = tempmg(commnodes1(i,j))
285  ENDDO
286  !CALL MPI_ISEND(bufsnd,2*NumCommNodes1(i),MPI_DOUBLE_PRECISION, &
287  ! CommProcs1(i),10,ROCSTAR_COMMUNICATOR,req_snd(i),ierr)
288  CALL mpi_send(bufsnd,2*numcommnodes1(i),mpi_double_precision, &
289  commprocs1(i),10,rocstar_communicator,ierr)
290  DEALLOCATE(bufsnd)
291  ENDDO
292  !CALL MPI_WAITALL(NumCommProcsFrom1,req_rcv,stat_rcv,ierr)
293  !CALL MPI_WAITALL(NumCommProcs1,req_snd,stat_snd,ierr)
294  CALL mpi_waitall(nprocs,req_rcv,stat_rcv,ierr)
295  !CALL MPI_WAITALL(nprocs,req_snd,stat_snd,ierr)
296  DEALLOCATE(req_rcv)
297  DEALLOCATE(req_snd)
298  DEALLOCATE(stat_snd)
299  DEALLOCATE(stat_rcv)
300 
301  ! Add the communicated mass to the current mass
302  DO i = 1, numcommprocsfrom1
303  DO j = 1, numcommnodesfrom1(i)
304  tempmg(global2local(int(frmproc(i)%rcvbuf(2*j-1)))) = &
305  tempmg(global2local(int(frmproc(i)%rcvbuf(2*j-1)))) + frmproc(i)%rcvbuf(2*j)
306  ENDDO
307  ENDDO
308 
309  DEALLOCATE(frmproc)
310  ENDIF
311 
312  ! Put the mass into compressed row format
313  counter = 0
314  DO m = 1, gnumnp
315  DO i = 1, global%NumNp
316  IF (local2global(i) == m) THEN
317  IF (nodeproc(i) == myid) THEN
318  DO j = 1, 3
319  counter = counter + 1
320  mg(counter) = tempmg(i)
321  ENDDO
322  ENDIF
323  ENDIF
324  ENDDO
325  ENDDO
326  DO i = 1, 3*lnumnp
327  rp_m_temp(i) = i-1
328  cval_m_temp(i) = i-1 + (nstart_km - 1)
329  aval_m_temp(i) = mg(i)
330  ENDDO
331  rp_m_temp(3*lnumnp+1) = 3*lnumnp
332  DEALLOCATE(tempmg)
333  DEALLOCATE(mg)
334 
335  ! Average the consistent and lumped mass matrices to get a higher-order mass matrix
336  CALL comp_row_add(3*lnumnp,3*gnumnp,nrows_km,nrows_km,nnz_m,3*lnumnp,1, &
337  rp_m,cval_m,0.5*aval_m,1,rp_m_temp,cval_m_temp,0.5*aval_m_temp)
338  DEALLOCATE(rp_m,cval_m,aval_m)
339  DEALLOCATE(rp_m_temp,cval_m_temp,aval_m_temp)
340  ALLOCATE(rp_m(nrows_km+1),cval_m(nnz_temp),aval_m(nnz_temp))
341  nnz_m = nnz_temp
342  rp_m = rp_temp
343  cval_m = cval_temp
344  aval_m = aval_temp
345  DEALLOCATE(rp_temp,cval_temp,aval_temp)
346 
347 END SUBROUTINE createm
348 
rocstar_rocfraccomm
Definition: ROCSTAR_RocFracComm.f90:53
Aff_transformation_2::m
FT m(int i, int j) const
Definition: Aff_transformation_2.h:228
comp_row_global
Definition: comp_row_global.f90:26
inode
**********************************************************************Rocstar Simulation Suite Illinois Rocstar LLC All rights reserved ****Illinois Rocstar LLC IL **www illinoisrocstar com **sales illinoisrocstar com WITHOUT WARRANTY OF ANY **EXPRESS OR INCLUDING BUT NOT LIMITED TO THE WARRANTIES **OF FITNESS FOR A PARTICULAR PURPOSE AND **NONINFRINGEMENT IN NO EVENT SHALL THE CONTRIBUTORS OR **COPYRIGHT HOLDERS BE LIABLE FOR ANY DAMAGES OR OTHER WHETHER IN AN ACTION OF TORT OR **Arising OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE **USE OR OTHER DEALINGS WITH THE SOFTWARE **********************************************************************INTERFACE SUBROUTINE inode
Definition: RFLO_CalcGradDummySymm.h:32
comp_row_add
subroutine comp_row_add(ndim, gndim, nrows1, nrows2, nnz1, nnz2, nstart1, rp1, cval1, aval1, nstart2, rp2, cval2, aval2)
Definition: comp_row_support.f90:89
implicit_global
Definition: implicit_global.f90:76
createm
subroutine createm(global)
Definition: createM.f90:53
comp_row_getval
subroutine comp_row_getval(ndim, nrows, nnz, nstart, rp, cval, aval, ipos, jpos, val)
Definition: comp_row_support.f90:511
comp_row_resize
subroutine comp_row_resize(ndim, nrows, nnz, nstart, rp, cval, aval, newnstart, newnrows)
Definition: comp_row_support.f90:592
implicit_v3d8_mass
subroutine implicit_v3d8_mass(NumEl, NumNP, NumMat, coor, nodes, MatType, ri, rho, xm, iElStart, iElEnd)
Definition: implicit_v3d8_mass.f90:82
i
blockLoc i
Definition: read.cpp:79
rocstar_rocfrac::rocfrac_global
Definition: ROCSTAR_RocFrac.f90:80
Aff_transformation_repS2::print
virtual std::ostream & print(std::ostream &os) const
Definition: Aff_transformationS2.h:277
j
j indices j
Definition: Indexing.h:6
implicit_v3d8_mass_consistent
subroutine implicit_v3d8_mass_consistent(NumEl, NumNP, NumMat, coor, nodes, MatType, ri, rho, ElConnVol)
Definition: implicit_v3d8_mass_consistent.f90:83
jnode
**********************************************************************Rocstar Simulation Suite Illinois Rocstar LLC All rights reserved ****Illinois Rocstar LLC IL **www illinoisrocstar com **sales illinoisrocstar com WITHOUT WARRANTY OF ANY **EXPRESS OR INCLUDING BUT NOT LIMITED TO THE WARRANTIES **OF FITNESS FOR A PARTICULAR PURPOSE AND **NONINFRINGEMENT IN NO EVENT SHALL THE CONTRIBUTORS OR **COPYRIGHT HOLDERS BE LIABLE FOR ANY DAMAGES OR OTHER WHETHER IN AN ACTION OF TORT OR **Arising OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE **USE OR OTHER DEALINGS WITH THE SOFTWARE **********************************************************************INTERFACE SUBROUTINE jnode
Definition: RFLO_CalcGradDummySymm.h:32
comp_row_addval
subroutine comp_row_addval(ndim, nrows, nnz, nstart, rp, cval, aval, ipos, jpos, val)
Definition: comp_row_support.f90:410