xref: /phasta/phSolver/compressible/au1mfg.f (revision 712d3df0b59ebebaaeaea358162c8d2c043c6e08)
1        subroutine Au1MFG (ypre,      y,    ac,     x,
2     &                     rmes,      res,       uBrg,
3     &                     BDiag,     iBC,       BC,
4     &                     iper,      ilwork,    shp,
5     &                     shgl,      shpb,
6     &                     shglb)
7c
8c----------------------------------------------------------------------
9c
10c This routine performs a matrix-vector product for the Matrix-Free
11c  Implicit/Iterative solver using a one-sided scheme.
12c
13c input:
14c  y      (nshg,ndof)           : Y-variables
15c  ypre   (nshg,nflow)           : preconditioned Y-variables
16c                                  (perturbed, no-scalars)
17c  x      (numnp,nsd)            : node coordinates
18c  rmes   (nshg,nflow)           : modified residual
19c  res    (nshg,nflow)           : residual
20c  uBrg   (nshg,nflow)           : Krylov space vector
21c  BDiag  (nshg,nflow,nflow)         : block-diagonal preconditioner
22c  iBC    (nshg)                : BC codes
23c  BC     (nshg,ndofBC)         : BC constraint parameters
24c  engBC  (nshg)                : energy for BC on density or pressure
25c  shp(b) (nshape,ngauss)        : element shape functions (boundary)
26c  shgl(b)(nsd,nshape,ngauss)    : local gradients of shape functions
27c
28c output:
29c  uBrg   (nshg,nflow)           : Krylov space vectors
30c
31c
32c Zdenek Johan, Winter 1991.  (Fortran 90)
33c----------------------------------------------------------------------
34c
35        include "common.h"
36        include "mpif.h"
37        include "auxmpi.h"
38c
39        dimension y(nshg,ndof),             ypre(nshg,nflow),
40     &            x(numnp,nsd),             ac(nshg,ndof),
41     &            rmes(nshg,nflow),          ytmp(nshg,nflow),
42     &            res(nshg,nflow),           uBrg(nshg,nflow),
43     &            BDiag(nshg,nflow,nflow),       iBC(nshg),
44     &            BC(nshg,ndofBC),          iper(nshg)
45c
46        dimension uBtmp(nshg,nflow),         tmpBC(nshg),
47     &            ilwork(nlwork)
48c
49
50        dimension shp(MAXTOP,maxsh,MAXQPT),
51     &            shgl(MAXTOP,nsd,maxsh,MAXQPT),
52     &            shpb(MAXTOP,maxsh,MAXQPT),
53     &            shglb(MAXTOP,nsd,maxsh,MAXQPT)
54
55c$$$        dimension shp(nshl,ngauss),  shgl(nsd,nshl,ngauss),
56c$$$     &            shpb(nshl,ngaussb),
57c$$$     &            shglb(nsd,nshl,ngaussb)
58c
59c.... calculate Rmod(y + eps u_i)
60c
61        uBtmp = zero
62c
63c        call yshuffle(ypre, 'new2old ')
64        uBrg = ypre + eGMRES * uBrg
65c
66        call i3LU (BDiag,  uBrg,  'backward')
67c
68        call yshuffle(uBrg, 'old2new ')
69c
70        call itrBC (uBrg, uBrg, iBC,  BC, iper, ilwork)
71c
72        call itrRes (uBrg,            y,
73     &               x,               shp,
74     &               shgl,            iBC,
75     &               BC,              shpb,
76     &               shglb,           uBtmp,
77     &               iper,            ilwork, ac)
78c
79        call i3LU (BDiag,  uBtmp,  'forward ')
80c
81c.... calculate ( Rmod(y + eps u_i) - Rmod(y) ) / eps
82c
83        uBrg = ( uBtmp - rmes ) / eGMRES
84c ... before returning lets put ypre back in the new format
85c         call yshuffle(ypre, 'old2new ')
86c
87c.... flop count
88c
89   !      flops = flops + 4*nflow*nshg
90c
91c.... end
92c
93        return
94        end
95
96
97
98
99