xref: /phasta/phSolver/compressible/itrres.f (revision 712d3df0b59ebebaaeaea358162c8d2c043c6e08)
1        subroutine ItrRes (yp,        yc,        x,
2     &                     shp,       shgl,
3     &                     iBC,       BC,        shpb,
4     &                     shglb,     rmes,      iper,
5     &                     ilwork,    ac)
6c
7c----------------------------------------------------------------------
8c
9c This routine calculates the modified residual vector.
10c
11c
12c Zdenek Johan, Winter 1991.  (Fortran 90)
13c----------------------------------------------------------------------
14c
15        use pointer_data
16c
17        include "common.h"
18        include "mpif.h"
19        include "auxmpi.h"
20c
21        dimension yp(nshg,nflow),             yc(nshg,ndof),
22     &            x(numnp,nsd),             ac(nshg,ndof),
23     &            iBC(nshg),                BC(nshg,ndofBC),
24     &            rmes(nshg,nflow),          ilwork(nlwork),
25     &            iper(nshg)
26c
27        dimension shp(MAXTOP,maxsh,MAXQPT),
28     &            shgl(MAXTOP,nsd,maxsh,MAXQPT),
29     &            shpb(MAXTOP,maxsh,MAXQPT),
30     &            shglb(MAXTOP,nsd,maxsh,MAXQPT)
31
32        real*8, allocatable :: tmpshp(:,:), tmpshgl(:,:,:)
33        real*8, allocatable :: tmpshpb(:,:), tmpshglb(:,:,:)
34
35        ttim(81) = ttim(81) - secs(0.0)
36
37c
38c.... -------------------->   interior elements   <--------------------
39c
40        jump  = 0
41        ires  = 2
42        iprec = 0
43c
44c.... loop over the element-blocks
45c
46        do iblk = 1, nelblk
47c
48c.... set up the parameters
49c$$$c
50c$$$          iel    = lcblk(1,iblk)
51c$$$          nenl   = lcblk(5,iblk)
52c$$$          mattyp = lcblk(7,iblk)
53c$$$          ndofl  = lcblk(8,iblk)
54c$$$          npro   = lcblk(1,iblk+1) - iel
55c
56          nenl   = lcblk(5,iblk)   ! no. of vertices per element
57          iel    = lcblk(1,iblk)
58          lelCat = lcblk(2,iblk)
59          lcsyst = lcblk(3,iblk)
60          iorder = lcblk(4,iblk)
61          nenl   = lcblk(5,iblk)   ! no. of vertices per element
62          nshl   = lcblk(10,iblk)
63          mattyp = lcblk(7,iblk)
64          ndofl  = lcblk(8,iblk)
65          nsymdl = lcblk(9,iblk)
66          npro   = lcblk(1,iblk+1) - iel
67          ngauss = nint(lcsyst)
68c
69c
70c.... compute and assemble the residuals and the preconditioner
71c
72          allocate (tmpshp(nshl, MAXQPT))
73          allocate (tmpshgl(nsd,nshl,MAXQPT))
74
75          tmpshp(1:nshl,:) = shp(lcsyst,1:nshl,:)
76          tmpshgl(:,1:nshl,:) = shgl(lcsyst,:,1:nshl,:)
77
78          call AsIRes (yp,                      yc,
79     &                 x,                       mxmudmi(iblk)%p,
80     &                 tmpshp,                  tmpshgl,
81     &                 mien(iblk)%p,            mmat(iblk)%p,
82     &                 rmes,                    ac)
83
84          deallocate (tmpshp)
85          deallocate (tmpshgl)
86
87c
88c.... end of interior element loop
89c
90        enddo
91c
92c.... -------------------->   boundary elements   <--------------------
93c
94        if (Navier .eq. 1 .and. Jactyp.ne.0) then
95c
96c.... loop over the elements
97c
98        do iblk = 1, nelblb
99c
100c.... set up the parameters
101c$$$c
102c$$$          iel    = lcblkb(1,iblk)
103c$$$          nenl   = lcblkb(5,iblk)
104c$$$          nenbl  = lcblkb(6,iblk)
105c$$$          mattyp = lcblkb(7,iblk)
106c$$$          ndofl  = lcblkb(8,iblk)
107c$$$          npro   = lcblkb(1,iblk+1) - iel
108c$$$c
109c
110          iel    = lcblkb(1,iblk)
111          lelCat = lcblkb(2,iblk)
112          lcsyst = lcblkb(3,iblk)
113          iorder = lcblkb(4,iblk)
114          nenl   = lcblkb(5,iblk)  ! no. of vertices per element
115          nenbl  = lcblkb(6,iblk)  ! no. of vertices per bdry. face
116          mattyp = lcblkb(7,iblk)
117          ndofl  = lcblkb(8,iblk)
118          nshl   = lcblkb(9,iblk)
119          nshlb  = lcblkb(10,iblk)
120          npro   = lcblkb(1,iblk+1) - iel
121          if(lcsyst.eq.3) lcsyst=nenbl
122          ngaussb = nintb(lcsyst)
123c
124          allocate (tmpshpb(nshl,MAXQPT))
125          allocate (tmpshglb(nsd,nshl,MAXQPT))
126
127          tmpshpb(1:nshl,:) = shpb(lcsyst,1:nshl,:)
128          tmpshglb(:,1:nshl,:) = shglb(lcsyst,:,1:nshl,:)
129
130c
131c.... compute and assemble the residuals
132c
133
134          call AsBRes (yp,   yc,              x,
135     &                 tmpshpb,           tmpshglb,
136     &                 mienb(iblk)%p,     mmatb(iblk)%p,
137     &                 miBCB(iblk)%p,     mBCB(iblk)%p,
138     &                 rmes)
139c
140
141          deallocate (tmpshpb)
142          deallocate (tmpshglb)
143c.... end of boundary element loop
144c
145        enddo
146
147        endif
148
149        ttim(81) = ttim(81) + secs(0.0)
150
151c
152c.... ---------------------->   communications  <-----------------------
153c
154        if((iabc==1)) !are there any axisym bc's
155     &       call rotabc(rmes(1,2), iBC, 'in ')
156c
157        if (numpe > 1) then
158           call commu (rmes, ilwork, nflow, 'in ')
159        endif
160
161c
162c.... ---------------------->   post processing  <----------------------
163c
164c.... satisfy the BCs on the modified residual
165c
166        call bc3Res (yc,  iBC,  BC,  rmes, iper, ilwork)
167c
168c.... return
169c
170        return
171        end
172