1!
2!     CalculiX - A 3-dimensional finite element program
3!     Copyright (C) 1998-2021 Guido Dhondt
4!
5!     This program is free software; you can redistribute it and/or
6!     modify it under the terms of the GNU General Public License as
7!     published by the Free Software Foundation(version 2);
8!
9!
10!     This program is distributed in the hope that it will be useful,
11!     but WITHOUT ANY WARRANTY; without even the implied warranty of
12!     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13!     GNU General Public License for more details.
14!
15!     You should have received a copy of the GNU General Public License
16!     along with this program; if not, write to the Free Software
17!     Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
18!
19      subroutine cfds(inpc,textpart,nmethod,iperturb,isolver,
20     &     istep,istat,n,tinc,tper,tmin,tmax,idrct,ithermal,iline,ipol,
21     &     inl,ipoinp,inp,ipoinpc,alpha,ctrl,iexpl,tincf,ttime,physcon,
22     &     ier)
23!
24!     reading the input deck: *CFD
25!
26!     isolver=0: SPOOLES
27!     2: iterative solver with diagonal scaling
28!     3: iterative solver with Cholesky preconditioning
29!     4: sgi solver
30!     5: TAUCS
31!     7: pardiso
32!     8: pastix
33!
34!     iexpl=0:  structure:implicit, fluid:incompressible
35!     iexpl=1:  structure:implicit, fluid:compressible explicit
36!     iexpl=3:  structure:explicit, fluid:compressible implicit
37!
38      implicit none
39!
40      logical timereset,fem,shallowwater
41!
42      character*1 inpc(*)
43      character*20 solver
44      character*132 textpart(16)
45!
46      integer nmethod,iperturb(*),isolver,istep,istat,n,key,i,idrct,
47     &     ithermal(*),iline,ipol,inl,ipoinp(2,*),inp(3,*),ipoinpc(0:*),
48     &     iexpl,ier
49!
50      real*8 tinc,tper,tmin,tmax,alpha(*),ctrl(*),tincf,ttime,physcon(*)
51!
52      idrct=0
53      tmin=0.d0
54      tmax=0.d0
55      tincf=-1.d0
56      nmethod=4
57      timereset=.false.
58!
59!     ONLY CFD-CALCULATIONS WITH THE CBS-METHOD ARE ALLOWED
60!
61c      fem=.false.
62      fem=.true.
63      shallowwater=.false.
64!
65!     default: no turbulence model
66!
67      physcon(9)=0.5d0
68!
69!     default:
70!     FVM method: 2D/3D switch: default 3D
71!
72      physcon(10)=3.5d0
73!
74      if(iperturb(1).eq.0) then
75        iperturb(1)=2
76      elseif((iperturb(1).eq.1).and.(istep.gt.1)) then
77        write(*,*) '*ERROR reading *CFD: perturbation analysis is'
78        write(*,*) '       not provided in a *HEAT TRANSFER step.'
79        ier=1
80        return
81      endif
82!
83      if(istep.lt.1) then
84        write(*,*) '*ERROR reading *CFD: *HEAT TRANSFER can only '
85        write(*,*) '       be used within a STEP'
86        ier=1
87        return
88      endif
89!
90!     default solver
91!
92      solver='                    '
93      if(isolver.eq.0) then
94        solver(1:7)='SPOOLES'
95      elseif(isolver.eq.2) then
96        solver(1:16)='ITERATIVESCALING'
97      elseif(isolver.eq.3) then
98        solver(1:17)='ITERATIVECHOLESKY'
99      elseif(isolver.eq.4) then
100        solver(1:3)='SGI'
101      elseif(isolver.eq.5) then
102        solver(1:5)='TAUCS'
103      elseif(isolver.eq.7) then
104        solver(1:7)='PARDISO'
105      elseif(isolver.eq.8) then
106        solver(1:6)='PASTIX'
107      endif
108!
109      do i=2,n
110        if(textpart(i)(1:7).eq.'SOLVER=') then
111          read(textpart(i)(8:27),'(a20)') solver
112        elseif(textpart(i)(1:12).eq.'COMPRESSIBLE') then
113          iexpl=1
114        elseif(textpart(i)(1:12).eq.'IMPLICIT') then
115          iexpl=3
116        elseif(textpart(i)(1:12).eq.'SHALLOWWATER') then
117          shallowwater=.true.
118          iexpl=1
119        elseif(textpart(i)(1:11).eq.'STEADYSTATE') then
120          nmethod=1
121        elseif(textpart(i)(1:9).eq.'TIMERESET') then
122          timereset=.true.
123        elseif(textpart(i)(1:17).eq.'TOTALTIMEATSTART=') then
124          read(textpart(i)(18:37),'(f20.0)',iostat=istat) ttime
125        elseif(textpart(i)(1:16).eq.'TURBULENCEMODEL=') then
126!
127!     turbulence model
128!
129          if(textpart(i)(17:25).eq.'NONE') then
130            physcon(9)=0.5d0
131          elseif(textpart(i)(17:25).eq.'K-EPSILON') then
132            physcon(9)=1.5d0
133          elseif(textpart(i)(17:23).eq.'K-OMEGA') then
134            physcon(9)=2.5d0
135          elseif(textpart(i)(17:19).eq.'BSL') then
136            physcon(9)=3.5d0
137          elseif(textpart(i)(17:19).eq.'SST') then
138            physcon(9)=4.5d0
139          endif
140        elseif(textpart(i)(1:2).eq.'2D') then
141          physcon(10)=2.5d0
142        elseif(textpart(i)(1:9).eq.'SCHEME=UD') then
143          ctrl(48)=1.5d0
144        elseif(textpart(i)(1:15).eq.'SCHEME=MODSMART') then
145          ctrl(48)=2.5d0
146        elseif(textpart(i)(1:7).eq.'SIMPLEC') then
147          ctrl(49)=1.5d0
148        elseif(textpart(i)(1:3).eq.'FEM') then
149          fem=.true.
150        else
151          write(*,*)
152     &         '*WARNING reading *CFD: parameter not recognized:'
153          write(*,*) '         ',
154     &         textpart(i)(1:index(textpart(i),' ')-1)
155          call inputwarning(inpc,ipoinpc,iline,
156     &         "*CFD%")
157        endif
158      enddo
159      if(nmethod.eq.1) ctrl(27)=1.d30
160
161      if((ithermal(1).eq.0).and.(iexpl.eq.1).and.
162     &     (.not.shallowwater)) then
163        write(*,*) '*ERROR reading *CFD: please define initial '
164        write(*,*) '       conditions for the temperature'
165        ier=1
166        return
167      elseif(ithermal(1).gt.0) then
168        ithermal(1)=3
169      endif
170!
171      if(solver(1:7).eq.'SPOOLES') then
172        isolver=0
173      elseif(solver(1:16).eq.'ITERATIVESCALING') then
174        isolver=2
175      elseif(solver(1:17).eq.'ITERATIVECHOLESKY') then
176        isolver=3
177      elseif(solver(1:3).eq.'SGI') then
178        isolver=4
179      elseif(solver(1:5).eq.'TAUCS') then
180        isolver=5
181      elseif(solver(1:7).eq.'PARDISO') then
182        isolver=7
183      elseif(solver(1:6).eq.'PASTIX') then
184        isolver=8
185      else
186        write(*,*) '*WARNING reading *CFD: unknown solver;'
187        write(*,*) '         the default solver is used'
188      endif
189!
190!     the number of the turbulence model is increase by 20 for
191!     the CBS method for shallow water equations and by 10 for
192!     the CBS method for all other applications
193!
194      if(fem) then
195        physcon(9)=physcon(9)+10.d0
196      endif
197!
198      if(shallowwater) then
199        physcon(9)=physcon(9)+10.d0
200      endif
201!
202      call getnewline(inpc,textpart,istat,n,key,iline,ipol,inl,
203     &     ipoinp,inp,ipoinpc)
204      if((istat.lt.0).or.(key.eq.1)) then
205        if(nmethod.eq.4) then
206!
207!     transient cfd
208!
209          write(*,*) '*WARNING reading *CFD: a nonlinear geometric
210     &analysis is requested'
211          write(*,*) '         but no time increment nor step is speci
212     &fied'
213          write(*,*) '         the defaults (1,1) are used'
214          tinc=1.d0
215          tper=1.d0
216          tmin=1.d-5
217          tmax=1.d+30
218          tincf=-1.d0
219        elseif(nmethod.eq.1) then
220!
221!     steady state cfd: initial increment time and step
222!     time are set to large values since only the first
223!     increment is calculated
224!
225          tinc=1.d0
226          tper=1.d0
227          tmin=1.d-5
228          tmax=1.d+30
229          tincf=-1.d0
230        endif
231        if(timereset)ttime=ttime-tper
232        return
233      endif
234!
235      read(textpart(1)(1:20),'(f20.0)',iostat=istat) tinc
236      if(istat.gt.0) then
237        call inputerror(inpc,ipoinpc,iline,
238     &       "*CFD%",ier)
239        return
240      endif
241      read(textpart(2)(1:20),'(f20.0)',iostat=istat) tper
242      if(istat.gt.0) then
243        call inputerror(inpc,ipoinpc,iline,
244     &       "*CFD%",ier)
245        return
246      endif
247      read(textpart(3)(1:20),'(f20.0)',iostat=istat) tmin
248      if(istat.gt.0) then
249        call inputerror(inpc,ipoinpc,iline,
250     &       "*CFD%",ier)
251        return
252      endif
253      read(textpart(4)(1:20),'(f20.0)',iostat=istat) tmax
254      if(istat.gt.0) then
255        call inputerror(inpc,ipoinpc,iline,
256     &       "*CFD%",ier)
257        return
258      endif
259      read(textpart(5)(1:20),'(f20.0)',iostat=istat) tincf
260      if(istat.gt.0) then
261        call inputerror(inpc,ipoinpc,iline,
262     &       "*CFD%",ier)
263        return
264      endif
265!
266      if(tinc.le.0.d0) then
267        write(*,*) '*ERROR reading *CFD: initial increment size is
268     &negative'
269      endif
270      if(tper.le.0.d0) then
271        write(*,*) '*ERROR reading *CFD: step size is negative'
272      endif
273      if(tinc.gt.tper) then
274        write(*,*) '*ERROR reading *CFD: initial increment size exc
275     &eeds step size'
276      endif
277!
278      if(idrct.ne.1) then
279        if(dabs(tmin).lt.1.d-6*tper) then
280          tmin=min(tinc,1.d-6*tper)
281        endif
282        if(dabs(tmax).lt.1.d-10) then
283          tmax=1.d+30
284        endif
285        if(tinc.gt.dabs(tmax)) then
286          write(*,*) '*WARNING reading *CFD:'
287          write(*,*) '         the initial increment ',tinc
288          write(*,*) '         exceeds the maximum increment ',
289     &         tmax
290          write(*,*) '         the initial increment is reduced'
291          write(*,*) '         to the maximum value'
292          tinc=dabs(tmax)
293        endif
294      endif
295!
296      if(timereset)ttime=ttime-tper
297!
298      call getnewline(inpc,textpart,istat,n,key,iline,ipol,inl,
299     &     ipoinp,inp,ipoinpc)
300!
301      return
302      end
303