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