xref: /phasta/phSolver/common/qpbc.f (revision 595995161822a203c8467e0e4a253d7bd7d6df32)
1      subroutine qpbc( rmass, qres, iBC, iper, ilwork )
2c---------------------------------------------------------------------
3c
4c This routine satisfies the periodic boundary conditions
5c on the diffusive flux residual and mass matrix
6c
7c input:
8c     rmass   (nshg)              : mass matrix
9c     qres    (nshg,(nflow-1)*nsd) : diffusive flux vector
10c
11c output: modified qres and rmass
12c---------------------------------------------------------------------
13      include "common.h"
14
15      dimension rmass(nshg), qres(nshg,idflx),
16     &          iBC(nshg), iper(nshg),uv(nshg,2),
17     &          tmpvec(nshg,4), tmp(nshg)
18c
19      if(iabc==1) then   !are there any axisym bc's
20      do i=1,idflx/nsd
21         do j=1,2
22            istrt=j+(i-1)*(nflow-1)
23            uv(:,j)=qres(:,istrt)
24         enddo
25         call rotabc(uv, iBC, 'in ')
26         do j=1,2
27            istrt=j+(i-1)*(nflow-1)
28            qres(:,istrt)=uv(:,j)
29         enddo
30      enddo
31      endif
32c
33c
34c.... compute qi for node A, i.e., qres <-- qres/rmass
35c
36       if (numpe > 1) then
37          call commu (qres  , ilwork, idflx  , 'in ')
38          call commu (rmass , ilwork,  1            , 'in ')
39       endif
40c
41c  take care of periodic boundary conditions
42c  but not on surface tension terms in qres(:,10-12)
43c  that are used to compute normal vector
44c
45        idflow = (nflow-1)*nsd
46        do j= 1,nshg
47          if ((btest(iBC(j),10))) then
48            i = iper(j)
49            rmass(i) = rmass(i) + rmass(j)
50c            qres(i,:) = qres(i,:) + qres(j,:)
51            qres(i,1:idflow) = qres(i,1:idflow) + qres(j,1:idflow)
52          endif
53        enddo
54
55        do j= 1,nshg
56          if ((btest(iBC(j),10))) then
57            i = iper(j)
58            rmass(j) = rmass(i)
59c            qres(j,:) = qres(i,:)
60            qres(j,1:idflow) = qres(i,1:idflow)
61          endif
62        enddo
63c
64c.... invert the diagonal mass matrix and find q
65c
66        rmass = one/rmass
67
68       do i=1, idflx
69          qres(:,i) = rmass*qres(:,i)
70       enddo
71       if (isurf .eq. 1) then
72          idflow=(nflow-1)*nsd
73c
74c.... calculation of the unit normal vector
75c
76           tmp =  sqrt(qres(:,idflow+1)**2
77     &               + qres(:,idflow+2)**2
78     &               + qres(:,idflow+3)**2)
79           do i = 1, nshg
80             if (tmp(i) .lt. 0.0001) tmp(i) = 0.0001
81           end do
82           tmp = one/tmp
83
84          do i=1, nsd
85             qres(:,idflow+i) = tmp*qres(:,idflow+i)
86          enddo
87       endif
88
89       if(numpe > 1) then
90          call commu (qres, ilwork, idflx, 'out')
91       endif
92
93       if(iabc==1) then         !are there any axisym bc's
94c
95c       slave has masters value, for abc we need to rotate it
96c
97          do i=1,idflx/nsd
98             do j=1,2
99                istrt=j+(i-1)*(nflow-1)
100                uv(:,j)=qres(:,istrt)
101             enddo
102             call rotabc(uv, iBC, 'out')
103             do j=1,2
104                istrt=j+(i-1)*(nflow-1)
105                qres(:,istrt)=uv(:,j)
106             enddo
107          enddo
108       endif
109
110c
111c.... return
112c
113        return
114        end
115
116
117      subroutine qpbcSclr( rmass, qres, iBC, iper, ilwork )
118c---------------------------------------------------------------------
119c
120c This routine satisfies the periodic boundary conditions
121c on the diffusive flux residual and mass matrix
122c
123c input:
124c     rmass   (nshg)              : mass matrix
125c     qres    (nshg, nsd)         : diffusive flux vector
126c
127c output: modified qres and rmass
128c---------------------------------------------------------------------
129      include "common.h"
130
131      dimension rmass(nshg), qres(nshg,nsd),
132     &          iBC(nshg), iper(nshg)
133
134      if(iabc==1) !are there any axisym bc's
135     &         call rotabc(qres, iBC,  'in ')
136
137c
138c.... compute qi for node A, i.e., qres <-- qres/rmass
139c
140       if (numpe > 1) then
141          call commu (qres  , ilwork, nsd  , 'in ')
142          call commu (rmass , ilwork,  1   , 'in ')
143       endif
144
145c
146c  take care of periodic boundary conditions
147c
148        do j= 1,nshg
149          if (btest(iBC(j),10)) then
150            i = iper(j)
151            rmass(i) = rmass(i) + rmass(j)
152            qres(i,:) = qres(i,:) + qres(j,:)
153          endif
154        enddo
155
156        do j= 1,nshg
157          if (btest(iBC(j),10)) then
158            i = iper(j)
159            rmass(j) = rmass(i)
160            qres(j,:) = qres(i,:)
161          endif
162        enddo
163c
164c.... invert the diagonal mass matrix and find q
165c
166        rmass = one/rmass
167
168       do i=1, nsd
169          qres(:,i) = rmass*qres(:,i)
170       enddo
171
172       if(numpe > 1) then
173          call commu (qres, ilwork, nsd, 'out')
174       endif
175
176      if(iabc==1) !are there any axisym bc's
177     &         call rotabc(qres, iBC, 'out')
178
179c
180c.... return
181c
182        return
183        end
184
185