xref: /phasta/phSolver/common/proces.f (revision ee150812e03b04b642c5c2bd3124e8eb86e98922)
1        subroutine proces
2c
3c----------------------------------------------------------------------
4c
5c This subroutine generates the problem data and calls the solution
6c  driver.
7c
8c
9c Zdenek Johan, Winter 1991.  (Fortran 90)
10c----------------------------------------------------------------------
11c
12        use readarrays          ! used to access x, iper, ilwork
13        use turbsa          ! used to access d2wall
14        use dtnmod
15        use periodicity
16        use pvsQbi
17        include "common.h"
18        include "mpif.h"
19c
20c arrays in the following 2 lines are now dimensioned in readnblk
21c        dimension x(numnp,nsd)
22c        dimension iper(nshg), ilwork(nlwork)
23c
24        dimension y(nshg,ndof),
25     &            iBC(nshg),
26     &            BC(nshg,ndofBC),
27     &            ac(nshg,ndof)
28c
29c.... shape function declarations
30c
31        dimension shp(MAXTOP,maxsh,MAXQPT),
32     &            shgl(MAXTOP,nsd,maxsh,MAXQPT),
33     &            shpb(MAXTOP,maxsh,MAXQPT),
34     &            shglb(MAXTOP,nsd,maxsh,MAXQPT)
35c
36c  stuff for dynamic model s.w.avg and wall model
37c
38c        dimension ifath(numnp),    velbar(nfath,nflow),
39c     &            nsons(nfath)
40
41        dimension velbar(nfath,nflow)
42c
43c stuff to interpolate profiles at inlet
44c
45        real*8 bcinterp(100,ndof+1),interp_mask(ndof)
46        logical exlog
47
48c        if ((irscale .ge. 0) .and. (myrank.eq.master)) then
49c           call setSPEBC(numnp, point2nfath, nsonmax)
50c        endif
51c
52c.... generate the geometry and boundary conditions data
53c
54        call gendat (y,              ac,             point2x,
55     &               iBC,            BC,
56     &               point2iper,     point2ilwork,   shp,
57     &               shgl,           shpb,           shglb,
58     &               point2ifath,    velbar,         point2nsons )
59        call setper(nshg)
60        call perprep(iBC,point2iper,nshg)
61        if (iLES/10 .eq. 1) then
62        call keeplhsG ! Allocating space for the mass (Gram) matrix to be
63                      ! used for projection filtering and reconstruction
64                      ! of the strain rate tensor.
65
66        call setrls   ! Allocating space for the resolved Leonard stresses
67c                         See bardmc.f
68        endif
69c
70c.... time averaged statistics
71c
72        if (ioform .eq. 2) then
73           call initStats(point2x, iBC, point2iper, point2ilwork)
74        endif
75c
76c.... RANS turbulence model
77c
78        if (iRANS .lt. 0) then
79           call initTurb( point2x )
80        endif
81c
82c.... p vs. Q boundary
83c
84           call initNABI( point2x, shpb )
85c
86c.... check for execution mode
87c
88        if (iexec .eq. 0) then
89           lstep = 0
90           call restar ('out ',  y  ,ac)
91           return
92        endif
93c
94c.... initialize AutoSponge
95c
96        if(matflg(5,1).ge.4) then ! cool case (sponge)
97           call initSponge( y,point2x)
98        endif
99c
100c
101c.... adjust BC's to interpolate from file
102c
103
104        inquire(file="inlet.dat",exist=exlog)
105        if(exlog) then
106           open (unit=654,file="inlet.dat",status="old")
107           read(654,*) ninterp,ibcmatch,(interp_mask(j),j=1,ndof)
108           do i=1,ninterp
109              read(654,*) (bcinterp(i,j),j=1,ndof+1) ! distance to wall+
110                        ! ndof but note order of BC's
111                        ! p,t,u,v,w,scalars. Also note that must be in
112                        ! increasing distance from the wall order.
113
114           enddo
115           do i=1,nshg  ! only correct for linears at this time
116              if(mod(iBC(i),1024).eq.ibcmatch) then
117                 iupper=0
118                 do j=2,ninterp
119                    if(bcinterp(j,1).gt.d2wall(i)) then !bound found
120                       xi=(d2wall(i)-bcinterp(j-1,1))/
121     &                    (bcinterp(j,1)-bcinterp(j-1,1))
122                       iupper=j
123                       exit
124                    endif
125                 enddo
126                 if(iupper.eq.0) then
127                    write(*,*) 'failure in finterp, ynint=',d2wall(i)
128                    stop
129                 endif
130                 if(interp_mask(1).ne.zero) then
131                    BC(i,1)=(xi*bcinterp(iupper,2)
132     &                +(one-xi)*bcinterp(iupper-1,2))*interp_mask(1)
133                 endif
134                 if(interp_mask(2).ne.zero) then
135                    BC(i,2)=(xi*bcinterp(iupper,3)
136     &                +(one-xi)*bcinterp(iupper-1,3))*interp_mask(2)
137                 endif
138                 if(interp_mask(3).ne.zero) then
139                    BC(i,3)=(xi*bcinterp(iupper,4)
140     &                +(one-xi)*bcinterp(iupper-1,4))*interp_mask(3)
141                 endif
142                 if(interp_masK(4).ne.zero) then
143                    BC(i,4)=(xi*bcinterp(iupper,5)
144     &                +(one-xi)*bcinterp(iupper-1,5))*interp_mask(4)
145                 endif
146                 if(interp_mask(5).ne.zero) then
147                    BC(i,5)=(xi*bcinterp(iupper,6)
148     &                +(one-xi)*bcinterp(iupper-1,6))*interp_mask(5)
149                 endif
150                 if(interp_mask(6).ne.zero) then
151                    BC(i,7)=(xi*bcinterp(iupper,7)
152     &                +(one-xi)*bcinterp(iupper-1,7))*interp_mask(6)
153                 endif
154                 if(interp_mask(7).ne.zero) then
155                    BC(i,8)=(xi*bcinterp(iupper,8)
156     &                +(one-xi)*bcinterp(iupper-1,8))*interp_mask(7)
157                 endif
158              endif
159           enddo
160        endif
161c$$$$$$$$$$$$$$$$$$$$
162
163c
164c
165c.... call the semi-discrete predictor multi-corrector iterative driver
166c
167        call itrdrv (y,              ac,
168     &               uold,           point2x,
169     &               iBC,            BC,
170     &               point2iper,     point2ilwork,   shp,
171     &               shgl,           shpb,           shglb,
172     &               point2ifath,    velbar,         point2nsons )
173c
174c.... return
175c
176c
177c.... stop CPU-timer
178c
179CAD        call timer ('End     ')
180c
181c.... close echo file
182c
183
184!MR CHANGE
185      if (numpe > 1) call MPI_BARRIER(MPI_COMM_WORLD, ierr)
186      if(myrank.eq.0)  then
187          write(*,*) 'process - before closing iecho'
188      endif
189!MR CHANGE END
190
191        close (iecho)
192
193!MR CHANGE
194      if (numpe > 1) call MPI_BARRIER(MPI_COMM_WORLD, ierr)
195      if(myrank.eq.0)  then
196          write(*,*) 'process - after closing iecho'
197      endif
198!MR CHANGE END
199
200
201c
202c.... end of the program
203c
204CAD        write(6,*) 'Life: ', second(0) - ttim(100)
205        deallocate(point2iper)
206        if(numpe.gt.1) then
207          call Dctypes(point2ilwork(1))
208          deallocate(point2ilwork)
209        endif
210        deallocate(point2x)
211        deallocate(point2nsons)
212        deallocate(point2ifath)
213        deallocate(uold)
214        deallocate(wnrm)
215        deallocate(otwn)
216        call finalizeDtN
217        call clearper
218        call finalizeNABI
219
220        return
221        end
222
223
224