1 /* CalculiX - A 3-dimensional finite element program */
2 /* Copyright (C) 1998-2021 Guido Dhondt */
3
4 /* This program is free software; you can redistribute it and/or */
5 /* modify it under the terms of the GNU General Public License as */
6 /* published by the Free Software Foundation(version 2); */
7 /* */
8
9 /* This program is distributed in the hope that it will be useful, */
10 /* but WITHOUT ANY WARRANTY; without even the implied warranty of */
11 /* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the */
12 /* GNU General Public License for more details. */
13
14 /* You should have received a copy of the GNU General Public License */
15 /* along with this program; if not, write to the Free Software */
16 /* Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. */
17
18 #include <unistd.h>
19 #include <stdio.h>
20 #include <math.h>
21 #include <stdlib.h>
22 #include <pthread.h>
23 #include "CalculiX.h"
24 #ifdef SPOOLES
25 #include "spooles.h"
26 #endif
27 #ifdef SGI
28 #include "sgi.h"
29 #endif
30 #ifdef TAUCS
31 #include "tau.h"
32 #endif
33 #ifdef PARDISO
34 #include "pardiso.h"
35 #endif
36
37 static ITG num_cpus;
38
compfluid(double ** cop,ITG * nk,ITG ** ipkonfp,ITG * konf,char ** lakonfp,char ** sidefacep,ITG * ifreestream,ITG * nfreestream,ITG * isolidsurf,ITG * neighsolidsurf,ITG * nsolidsurf,ITG * nshcon,double * shcon,ITG * nrhcon,double * rhcon,double ** voldp,ITG * ntmat_,ITG * nodeboun,ITG * ndirboun,ITG * nboun,ITG * ipompc,ITG * nodempc,ITG * nmpc,ITG * ikmpc,ITG * ilmpc,ITG * ithermal,ITG * ikboun,ITG * ilboun,ITG * iturbulent,ITG * isolver,ITG * iexpl,double * ttime,double * time,double * dtime,ITG * nodeforc,ITG * ndirforc,double * xforc,ITG * nforc,ITG * nelemload,char * sideload,double * xload,ITG * nload,double * xbody,ITG * ipobody,ITG * nbody,ITG * ielmatf,char * matname,ITG * mi,ITG * ncmat_,double * physcon,ITG * istep,ITG * iinc,ITG * ibody,double * xloadold,double * xboun,double * coefmpc,ITG * nmethod,double * xforcold,double * xforcact,ITG * iamforc,ITG * iamload,double * xbodyold,double * xbodyact,double * t1old,double * t1,double * t1act,ITG * iamt1,double * amta,ITG * namta,ITG * nam,double * ampli,double * xbounold,double * xbounact,ITG * iamboun,ITG * itg,ITG * ntg,char * amname,double * t0,ITG ** nelemfacep,ITG * nface,double * cocon,ITG * ncocon,double * xloadact,double * tper,ITG * jmax,ITG * jout,char * set,ITG * nset,ITG * istartset,ITG * iendset,ITG * ialset,char * prset,char * prlab,ITG * nprint,double * trab,ITG * inotr,ITG * ntrans,char * filab,char * labmpc,double * sti,ITG * norien,double * orab,char * jobnamef,char * tieset,ITG * ntie,ITG * mcs,ITG * ics,double * cs,ITG * nkon,ITG * mpcfree,ITG * memmpc_,double * fmpc,ITG * nef,ITG ** inomatp,double * qfx,ITG * neifa,ITG * neiel,ITG * ielfa,ITG * ifaext,double * vfa,double * vel,ITG * ipnei,ITG * nflnei,ITG * nfaext,char * typeboun,ITG * neij,double * tincf,ITG * nactdoh,ITG * nactdohinv,ITG * ielorienf,char * jobnamec,ITG * ifatie,ITG * nstate_,double * xstate,char * orname,ITG * kon,double * ctrl,ITG * kode,double * velo,double * veloo,ITG * initial)39 void compfluid(double **cop,ITG *nk,ITG **ipkonfp,ITG *konf,char **lakonfp,
40 char **sidefacep,ITG *ifreestream,ITG *nfreestream,
41 ITG *isolidsurf,ITG *neighsolidsurf,ITG *nsolidsurf,
42 ITG *nshcon,double *shcon,ITG *nrhcon,double *rhcon,
43 double **voldp,ITG *ntmat_,ITG *nodeboun,ITG *ndirboun,
44 ITG *nboun,ITG *ipompc,ITG *nodempc,ITG *nmpc,ITG *ikmpc,
45 ITG *ilmpc,ITG *ithermal,ITG *ikboun,ITG *ilboun,
46 ITG *iturbulent,ITG *isolver,ITG *iexpl,double *ttime,
47 double *time,double *dtime,ITG *nodeforc,ITG *ndirforc,
48 double *xforc,ITG *nforc,ITG *nelemload,char *sideload,
49 double *xload,ITG *nload,double *xbody,ITG *ipobody,ITG *nbody,
50 ITG *ielmatf,char *matname,ITG *mi,ITG *ncmat_,double *physcon,
51 ITG *istep,ITG *iinc,ITG *ibody,double *xloadold,double *xboun,
52 double *coefmpc,ITG *nmethod,double *xforcold,double *xforcact,
53 ITG *iamforc,ITG *iamload,double *xbodyold,double *xbodyact,
54 double *t1old,double *t1,double *t1act,ITG *iamt1,double *amta,
55 ITG *namta,ITG *nam,double *ampli,double *xbounold,
56 double *xbounact,ITG *iamboun,ITG *itg,ITG *ntg,char *amname,
57 double *t0,ITG **nelemfacep,ITG *nface,double *cocon,
58 ITG *ncocon,double *xloadact,double *tper,ITG *jmax,ITG *jout,
59 char *set,ITG *nset,ITG *istartset,ITG *iendset,ITG *ialset,
60 char *prset,char *prlab,ITG *nprint,double *trab,ITG *inotr,
61 ITG *ntrans,char *filab,char *labmpc,double *sti,ITG *norien,
62 double *orab,char *jobnamef,char *tieset,ITG *ntie,ITG *mcs,
63 ITG *ics,double *cs,ITG *nkon,ITG *mpcfree,ITG *memmpc_,
64 double *fmpc,ITG *nef,ITG **inomatp,double *qfx,ITG *neifa,
65 ITG *neiel,ITG *ielfa,ITG *ifaext,double *vfa,double *vel,
66 ITG *ipnei,ITG *nflnei,ITG *nfaext,char *typeboun,ITG *neij,
67 double *tincf,ITG *nactdoh,ITG *nactdohinv,ITG *ielorienf,
68 char*jobnamec,ITG *ifatie,ITG *nstate_,double *xstate,
69 char *orname,ITG *kon,double *ctrl,ITG *kode,double *velo,
70 double *veloo,ITG *initial){
71
72 /* main computational fluid dynamics routine */
73
74 char cflag[1],*lakonf=NULL,fncvg[132]="",*lakon=NULL;
75
76 ITG *ipointer=NULL,*mast1=NULL,*irow=NULL,*icol=NULL,*jq=NULL,nzs=20000000,
77 compressible,*ifabou=NULL,*ja=NULL,*ikf=NULL,nfabou,im,iflag,*ipkon=NULL,
78 *ielprop=NULL,*ielmat=NULL,*ipkonf=NULL,last=0,icyclic,*iau6=NULL,
79 ifixdtimef=0,ithermalref,*integerglob=NULL,iincf,ipower=64,*konl=NULL,
80 iconvergence=0,i,*inum=NULL,j,k,ifreefa,isiz,*ipofano=NULL,*ifano=NULL,
81 *ia=NULL,*ielpropf=NULL,icent=0,isti=0,iqfx=0,nfield,ndim,iorienglob,
82 force=0,icfd=1,imach=0,ikappa=0,iatleastonepressurebc,iturb=0,*inoel=NULL,
83 *iponoel=NULL,icounter,ischeme=1,isimplec=0,iitf,iitg,iitp,*iam=NULL,
84 *jam=NULL,*iamorig=NULL,nz_num,*nestart=NULL,*ineighblock=NULL,
85 *neighblock=NULL,nneighblock,iit,iito,ip,ierrmax,ncfd,iitpt,*inlet=NULL,
86 *ipgradfa=NULL,*ipbount=NULL,*ipbounv1=NULL,*ipbounv2=NULL,
87 *ipbounv3=NULL,*ipbounp=NULL,*ibount=NULL,*ibounv1=NULL,*ibounv2=NULL,
88 *ibounv3=NULL,*ibounp=NULL,nbt,nbv1,nbv2,nbv3,nbp,nbpt,nbpv1,nbpv2,nbpv3,
89 nbpp,nkf,*iponofa=NULL,*inofa=NULL,symmetryflag=2,inputformat=4;
90
91 ITG nelt,isym,itol,itmax,iunit,lrgw,*igwk=NULL,ligw,ierr,*iwork=NULL,iter,
92 nsave,lenw,leniw;
93
94 double *umfa=NULL,reltime,*doubleglob=NULL,dtimefold,
95 *co=NULL,*vold=NULL,*coel=NULL,*cosa=NULL,*gradvel=NULL,*gradvfa=NULL,
96 *xxn=NULL,*xxi=NULL,*xle=NULL,*xlen=NULL,*xlet=NULL,timef,dtimef,
97 *cofa=NULL,*area=NULL,*xrlfa=NULL,reltimef,ttimef,*hcfa=NULL,*cvel=NULL,
98 *au=NULL,*ad=NULL,*b=NULL,*volume=NULL,*body=NULL,*dy=NULL,
99 *advfa=NULL,*ap=NULL,*bp=NULL,*xxj=NULL,*gradkel=NULL,*gradoel=NULL,
100 *cosb=NULL,hmin,tincfguess,*h=NULL,*gradelsh=NULL,
101 *hel=NULL,*hfa=NULL,*auv=NULL,*adv=NULL,*bv=NULL,*sel=NULL,*gamma=NULL,
102 *gradtfa=NULL,*gradtel=NULL,*umel=NULL,*cvfa=NULL,*gradpel=NULL,
103 *eei=NULL,*ener=NULL,*thicke=NULL,*eme=NULL,c[9],*gradkfa=NULL,
104 ptimef,*stn=NULL,*qfn=NULL,*hcel=NULL,*aua=NULL,a1,a2,a3,beta=0.,
105 *prop=NULL,*xxni=NULL,*xxnj=NULL,*xxicn=NULL,*xturb=NULL,
106 *xmach=NULL,*xkappa=NULL,*flux=NULL,*sc=NULL,*gradfash=NULL,
107 relnormt,relnormv,relnormp=0,relnormmax=1.e30,*temp=NULL,*yy=NULL,
108 *gradofa=NULL,betam=0.1,*gradpfa=NULL,*gradpcel=NULL,*gradpcfa=NULL,
109 *am=NULL,*f1=NULL,*of2=NULL,*xxna=NULL,*ale=NULL,*alet=NULL,
110 *ratio=NULL,dgmrestol=1.e-12,*bvcp=NULL,*bcp=NULL,*xbount=NULL,
111 *xbounv1=NULL,*xbounv2=NULL,*xbounv3=NULL,*xbounp=NULL,*adb=NULL,
112 *aub=NULL,sigma=0;
113
114 double tol,*rgwk=NULL,err,*sb=NULL,*sx=NULL,*rwork=NULL,*rf=NULL;
115
116 FILE *f3;
117
118 co=*cop;ipkonf=*ipkonfp;lakonf=*lakonfp;vold=*voldp;
119
120 #ifdef SGI
121 ITG token;
122 #endif
123
124 ncfd=(ITG)physcon[9];
125 // printf("ncfd=%d\n",ncfd);
126
127 strcpy(fncvg,jobnamec);
128 strcat(fncvg,".fcv");
129
130 if((f3=fopen(fncvg,"w"))==NULL){
131 // if((f3=fopen("fluidconvergence","w"))==NULL){
132 printf("*ERROR in compfluid: cannot open cvg file for writing...");
133 exit(0);
134 }
135 fprintf(f3,"temperature velocity pressure\n\n");
136
137 /* relative time at the end of the mechanical increment */
138
139 reltime=(*time)/(*tper);
140
141 /* open frd-file for fluids */
142
143 FORTRAN(openfilefluid,(jobnamef));
144
145 /* variables for multithreading procedure */
146
147 ITG sys_cpus;
148 char *env,*envloc,*envsys;
149
150 num_cpus=0;
151 sys_cpus=0;
152
153 /* explicit user declaration prevails */
154
155 envsys=getenv("NUMBER_OF_CPUS");
156 if(envsys){
157 sys_cpus=atoi(envsys);
158 if(sys_cpus<0) sys_cpus=0;
159 }
160
161 /* automatic detection of available number of processors */
162
163 if(sys_cpus==0){
164 sys_cpus = getSystemCPUs();
165 if(sys_cpus<1) sys_cpus=1;
166 }
167
168 /* local declaration prevails,if strictly positive */
169
170 envloc = getenv("CCX_NPROC_CFD");
171 if(envloc){
172 num_cpus=atoi(envloc);
173 if(num_cpus<0){
174 num_cpus=0;
175 }else if(num_cpus>sys_cpus){
176 num_cpus=sys_cpus;
177 }
178 }
179
180 /* else global declaration,if any,applies */
181
182 env = getenv("OMP_NUM_THREADS");
183 if(num_cpus==0){
184 if (env)
185 num_cpus = atoi(env);
186 if (num_cpus < 1) {
187 num_cpus=1;
188 }else if(num_cpus>sys_cpus){
189 num_cpus=sys_cpus;
190 }
191 }
192
193 // next line is to be inserted in a similar way for all other paralell parts
194
195 if(*nef<num_cpus) num_cpus=*nef;
196
197 printf(" Using up to %" ITGFORMAT " cpu(s) for CFD.\n",num_cpus);
198
199 // pthread_t tid[num_cpus];
200
201 /* *iexpl==0: structure:implicit,fluid:incompressible
202 *iexpl==1: structure:implicit,fluid:compressible
203 *iexpl==2: structure:explicit,fluid:incompressible
204 *iexpl==3: structure:explicit,fluid:compressible */
205
206 if((*iexpl==1)||(*iexpl==3)){
207 compressible=1;
208 }else{
209 compressible=0;
210 }
211
212 /* ischeme=1: ud
213 ischeme=2: modified smart */
214
215 ischeme=(ITG)floor(ctrl[47]);
216 if(compressible==1){
217 if(ischeme==1){
218 printf(" CFD scheme: upwind difference\n");
219 }else{
220 printf(" CFD scheme: modified smart\n");
221 }
222 }else{
223 printf(" CFD scheme: Gamma\n");
224 }
225
226 /* isimplec=0: simple
227 isimplec=1: simplec */
228
229 isimplec=(ITG)floor(ctrl[48]);
230 if(compressible==1){
231 if(isimplec==0){
232 printf(" algorithm: simple\n\n");
233 }else{
234 printf(" algorithm: simplec\n\n");
235 }
236 }
237
238 /* iitf: max number of transient iterations
239 iitg: max number of geometric iterations (extrapol_*.f)
240 iitp: max number of pressure iterations */
241
242 iitf=(ITG)floor(ctrl[49]);
243 iitg=(ITG)floor(ctrl[50]);
244 iitp=(ITG)floor(ctrl[51]);
245 iitpt=(ITG)floor(ctrl[52]);
246
247 /* if initial conditions are specified for the temperature,
248 it is assumed that the temperature is an unknown */
249
250 ithermalref=*ithermal;
251 if(*ithermal==1){
252 *ithermal=2;
253 }
254
255 /* determining the matrix structure */
256
257 NNEW(ipointer,ITG,*nef);
258 NNEW(mast1,ITG,nzs);
259 NNEW(irow,ITG,nzs);
260 NNEW(icol,ITG,*nef);
261 NNEW(jq,ITG,*nef+1);
262
263 mastructf(nk,konf,ipkonf,lakonf,nef,icol,jq,&mast1,&irow,
264 isolver,ipointer,&nzs,ipnei,neiel,mi);
265
266 SFREE(ipointer);SFREE(mast1);
267
268 NNEW(iau6,ITG,6**nef);
269 FORTRAN(create_iau6,(nef,ipnei,neiel,jq,irow,&nzs,iau6,lakonf));
270
271 if(compressible==0){
272 NNEW(ia,ITG,nzs+*nef);
273 NNEW(ja,ITG,*nef+1);
274 NNEW(aua,double,nzs+*nef);
275 FORTRAN(preconvert2slapcol,(irow,ia,jq,ja,&nzs,nef));
276 }
277
278 /* calculation geometric data */
279
280 NNEW(coel,double,3**nef);
281 NNEW(volume,double,*nef);
282 NNEW(h,double,*nflnei);
283 NNEW(sc,double,*nef);
284 DMEMSET(sc,0,*nef,1.);
285 NNEW(cosa,double,*nflnei);
286 NNEW(cosb,double,*nflnei);
287 NNEW(xxn,double,3**nflnei);
288 NNEW(xxna,double,3**nflnei);
289 NNEW(xxi,double,3**nflnei);
290 NNEW(xxj,double,3**nflnei);
291 NNEW(xxni,double,3**nflnei);
292 NNEW(xxicn,double,3**nflnei);
293 NNEW(xxnj,double,3**nflnei);
294 NNEW(xle,double,*nflnei);
295 NNEW(ale,double,*nflnei);
296 NNEW(xlen,double,*nflnei);
297 NNEW(xlet,double,*nflnei);
298 NNEW(alet,double,*nflnei);
299 NNEW(cofa,double,3**nface);
300 NNEW(area,double,*nface);
301 NNEW(xrlfa,double,3**nface);
302 NNEW(rf,double,3**nface);
303
304 /* closest distance to a node from a solid surface (dy)
305 distance from any element center to a solid surface (yy) */
306
307 if(*iturbulent>0){
308 NNEW(dy,double,*nsolidsurf);
309 if(*iturbulent>2){
310 NNEW(yy,double,*nef);
311 }
312 }
313
314 FORTRAN(geocfd,(nef,ipkonf,konf,lakonf,co,coel,cofa,nface,
315 ielfa,area,ipnei,neiel,xxn,xxi,xle,xlen,xlet,xrlfa,cosa,
316 volume,neifa,xxj,cosb,&hmin,ifatie,cs,tieset,&icyclic,c,
317 neij,physcon,isolidsurf,nsolidsurf,dy,xxni,xxnj,xxicn,
318 nflnei,iturbulent,rf,yy,vel,velo,veloo,xxna,ale,alet,h));
319
320 /* creating iam from neiel:
321 1) get rid of zero neighbors
322 2) get rid of cyclic symmetry neighbors
323 3) for parallel execution: get rid of neighbors belonging
324 to a different block */
325
326 NNEW(jam,ITG,*nef+1);
327 NNEW(iam,ITG,*nflnei+*nef);
328 NNEW(iamorig,ITG,*nflnei+*nef);
329 NNEW(nestart,ITG,num_cpus+1);
330 NNEW(ineighblock,ITG,num_cpus+1);
331 NNEW(neighblock,ITG,3**nflnei);
332
333 FORTRAN(createblock,(nef,ipnei,neiel,iam,jam,iamorig,nflnei,
334 &nz_num,&num_cpus,nestart,ineighblock,
335 neighblock,&icyclic,neifa,ifatie,&nneighblock));
336 RENEW(iam,ITG,nz_num);
337 RENEW(iamorig,ITG,nz_num);
338 RENEW(neighblock,ITG,3*nneighblock);
339 NNEW(am,double,nz_num);
340
341 /* storing pointers to the boundary conditions in ielfa */
342
343 NNEW(inlet,ITG,*nface);
344 NNEW(ifabou,ITG,7**nfaext);
345 FORTRAN(applyboun,(ifaext,nfaext,ielfa,ikboun,ilboun,
346 nboun,typeboun,nelemload,nload,sideload,isolidsurf,
347 nsolidsurf,ifabou,&nfabou,nface,nodeboun,ndirboun,ikmpc,
348 ilmpc,labmpc,nmpc,nactdohinv,&compressible,
349 &iatleastonepressurebc,ipkonf,kon,konf,inlet));
350 RENEW(ifabou,ITG,nfabou);
351
352 /* catalogueing the nodes for output purposes (interpolation at
353 the nodes */
354
355 NNEW(ipofano,ITG,*nk);
356 NNEW(ifano,ITG,2**nface*4);
357
358 FORTRAN(cataloguenodes,(ipofano,ifano,&ifreefa,ielfa,ifabou,ipkonf,
359 konf,lakonf,nface,nk));
360
361 RENEW(ifano,ITG,2*ifreefa);
362
363 /* determining the coefficients to get from the element center values:
364 - the nodal values (by extra/interpolation of the element center values)
365 - the facial values (by interpolation from the nodal values)
366 - the center gradient values (by use of the shape functions using the
367 nodal values)
368 - the facial gradient values (by use of the shape functions using the
369 nodal values) */
370
371 /* elem2node(coel,nef,&ncfd,&num_cpus,ipkonf,konf,lakonf,co,&hmin,ipnei,neiel,
372 xxn,nkon,nk,neifa,area,&ikf,&konl,&ratio,&gradelsh,&gradfash,nflnei,
373 &ipgradfa,&ipbount,&ipbounv1,&ipbounv2,&ipbounv3,&ipbounp,
374 &ibount,&ibounv1,&ibounv2,&ibounv3,&ibounp,&xbount,&xbounv1,
375 &xbounv2,&xbounv3,&xbounp,&nbt,&nbv1,&nbv2,&nbv3,&nbp,ithermal,
376 ipofano,ifano,&nbpt,&nbpv1,&nbpv2,&nbpv3,&nbpp,ielfa,ifabou,
377 &ifreefa,&nkf,&iponofa,&inofa,nface);*/
378
379 // SFREE(ikf);SFREE(konl);SFREE(ratio);
380 // SFREE(gradfash);SFREE(gradelsh);SFREE(ipgradfa);
381
382 /* material properties for athermal calculations
383 = calculation for which no initial thermal conditions
384 were defined */
385
386 NNEW(umfa,double,*nface);
387 NNEW(umel,double,*nef);
388
389 // if((*ithermal==0)||(*iturbulent>0)){
390 if(*ithermal==0){
391
392 /* athermal incompressible calculations */
393
394 /* calculating the dynamic viscosity at the element centers */
395
396 FORTRAN(calcumel,(nef,vel,shcon,nshcon,ielmatf,ntmat_,
397 ithermal,mi,umel));
398
399 }
400
401
402 if(*ithermal>0){
403 NNEW(hcfa,double,*nface);
404 NNEW(cvel,double,*nef);
405 NNEW(cvfa,double,*nface);
406 }
407
408 if(*nbody>0) NNEW(body,double,4**nef);
409
410 /* next section is for stationary calculations */
411
412 if(*nmethod==1){
413
414 /* boundary conditions at the end of the mechanical
415 increment */
416
417 FORTRAN(tempload,(xforcold,xforc,xforcact,iamforc,nforc,
418 xloadold,xload,xloadact,iamload,nload,ibody,xbody,nbody,
419 xbodyold,xbodyact,t1old,t1,t1act,iamt1,nk,amta,
420 namta,nam,ampli,time,&reltime,ttime,dtime,ithermal,nmethod,
421 xbounold,xboun,xbounact,iamboun,nboun,
422 nodeboun,ndirboun,nodeforc,ndirforc,istep,iinc,
423 co,vold,itg,ntg,amname,ikboun,ilboun,nelemload,sideload,mi,
424 ntrans,trab,inotr,vold,integerglob,doubleglob,tieset,istartset,
425 iendset,ialset,ntie,nmpc,ipompc,ikmpc,ilmpc,nodempc,coefmpc,
426 ipobody,iponoel,inoel,ipkon,kon,ielprop,prop,ielmat,
427 shcon,nshcon,rhcon,nrhcon,cocon,ncocon,ntmat_,lakon,set,nset));
428
429 /* body forces (gravity,centrifugal and Coriolis forces */
430
431 if(*nbody>0){
432 FORTRAN(inicalcbody,(nef,body,ipobody,ibody,xbody,coel,vel,lakonf,
433 nactdohinv,&icent));
434 }
435 }
436
437 /* extrapolating the velocity from the elements centers to the face
438 centers,thereby taking the boundary conditions into account */
439
440 NNEW(gradvel,double,9**nef);
441 NNEW(gradvfa,double,9**nface);
442
443 extrapol_velmain(nface,ielfa,xrlfa,vel,vfa,
444 ifabou,xbounact,ipnei,nef,&icyclic,c,ifatie,xxn,gradvel,
445 gradvfa,neifa,rf,area,volume,xle,xxi,xxj,xlet,
446 coefmpc,nmpc,labmpc,ipompc,nodempc,ifaext,nfaext,nactdoh,
447 &iitg,&num_cpus,&compressible,xxna,&ncfd,cofa);
448
449 /* extrapolation of the pressure at the element centers
450 to the face centers */
451
452 NNEW(gradpel,double,3**nef);
453 NNEW(gradpfa,double,3**nface);
454
455 extrapol_pelmain(nface,ielfa,xrlfa,vel,vfa,
456 ifabou,xbounact,nef,gradpel,gradpfa,neifa,rf,area,volume,
457 xle,xxi,&icyclic,xxn,ipnei,ifatie,
458 coefmpc,nmpc,labmpc,ipompc,nodempc,ifaext,nfaext,nactdoh,
459 &iitg,c,&num_cpus,&compressible,xxna,&ncfd);
460
461 /* generate fields for the pressure correction gradients */
462
463 NNEW(gradpcel,double,3**nef);
464 NNEW(gradpcfa,double,3**nface);
465
466 /* extrapolation of the temperature at the element centers
467 to the face centers */
468
469 if(*ithermal>0){
470
471 NNEW(gradtel,double,3**nef);
472 NNEW(gradtfa,double,3**nface);
473
474 extrapol_telmain(nface,ielfa,xrlfa,vel,vfa,
475 ifabou,xbounact,nef,gradtel,gradtfa,neifa,rf,area,volume,
476 xle,xxi,&icyclic,xxn,ipnei,ifatie,xload,xlet,xxj,
477 coefmpc,nmpc,labmpc,ipompc,nodempc,ifaext,nfaext,nactdoh,
478 &iitg,c,&num_cpus,&compressible,xxna,&ncfd);
479
480 /* calculating the heat conduction at the face centers */
481
482 FORTRAN(calchcfa,(nface,vfa,cocon,ncocon,ielmatf,ntmat_,
483 mi,ielfa,hcfa));
484
485 if(compressible==0){
486
487 /* calculating the specific heat at constant volume at the
488 face centers (secant value) */
489
490 FORTRAN(calccvfa,(nface,vfa,shcon,nshcon,ielmatf,ntmat_,
491 mi,ielfa,cvfa,physcon));
492 }else{
493
494 /* calculating the specific heat at constant volume at the
495 face centers (secant value) */
496
497 FORTRAN(calccvfacomp,(nface,vfa,shcon,nshcon,ielmatf,ntmat_,
498 mi,ielfa,cvfa,physcon));
499 }
500 }
501
502 // NNEW(flux,double,6**nef);
503 NNEW(flux,double,*nflnei);
504
505 if(compressible==0){
506
507 /* calculating the density at the element centers */
508
509 FORTRAN(calcrhoel,(nef,vel,rhcon,nrhcon,ielmatf,ntmat_,
510 ithermal,mi));
511
512 /* calculating the density at the face centers */
513
514 FORTRAN(calcrhofa,(nface,vfa,rhcon,nrhcon,ielmatf,ntmat_,
515 ithermal,mi,ielfa));
516
517 }else{
518
519 /* calculating the density at the element centers */
520
521 calcrhoelcompmain(nef,vel,shcon,ielmatf,ntmat_,mi,
522 &num_cpus);
523
524 /* calculating the density at the face centers */
525
526 if(ischeme==1){
527 hrr_udmain(nface,vfa,shcon,ielmatf,ntmat_,
528 mi,ielfa,ipnei,vel,nef,flux,
529 &num_cpus,xxi,xle,gradpel,gradtel,
530 neij);
531 }else{
532 hrr_mod_smartmain(nface,vfa,shcon,ielmatf,ntmat_,
533 mi,ielfa,ipnei,vel,nef,flux,
534 gradpel,gradtel,xxj,xlet,
535 &num_cpus);
536 }
537
538 }
539
540 /* calculating the initial mass flux */
541
542 FORTRAN(calcinitialflux,(area,vfa,xxna,ipnei,nef,neifa,lakonf,flux));
543
544 /* calculating the dynamic viscosity at the face centers */
545
546 FORTRAN(calcumfa,(nface,vfa,shcon,nshcon,ielmatf,ntmat_,
547 ithermal,mi,ielfa,umfa));
548
549 /* extrapolation of the turbulence variables at the element centers
550 to the face centers */
551
552 if(*iturbulent>0){
553
554 NNEW(gradkel,double,3**nef);
555 NNEW(gradkfa,double,3**nface);
556 NNEW(gradoel,double,3**nef);
557 NNEW(gradofa,double,3**nface);
558
559 DMEMSET(vel,7**nef,8**nef,1.);
560
561 extrapol_kelmain(nface,ielfa,xrlfa,vel,vfa,
562 ifabou,xbounact,nef,gradkel,gradkfa,neifa,rf,area,volume,
563 xle,xxi,&icyclic,xxn,ipnei,ifatie,xlet,xxj,
564 coefmpc,nmpc,labmpc,ipompc,nodempc,ifaext,nfaext,nactdoh,
565 umfa,physcon,&iitg,c,&num_cpus,&compressible,xxna,&ncfd,
566 inlet);
567
568 extrapol_oelmain(nface,ielfa,xrlfa,vel,vfa,
569 ifabou,xbounact,nef,gradoel,gradofa,neifa,rf,area,volume,
570 xle,xxi,&icyclic,xxn,ipnei,ifatie,xlet,xxj,
571 coefmpc,nmpc,labmpc,ipompc,nodempc,ifaext,nfaext,nactdoh,
572 umfa,physcon,&iitg,c,dy,&num_cpus,&compressible,xxna,&ncfd,
573 inlet);
574
575 if(*iturbulent>2){
576 NNEW(f1,double,*nef);
577 if(*iturbulent>3) NNEW(of2,double,*nef);
578 }
579
580 }
581
582 /* calculating the time increment */
583
584 FORTRAN(initincf,(nface,&hmin,vfa,umfa,cvfa,hcfa,ithermal,&dtimef,
585 &compressible));
586
587 /* start of the major loop */
588
589 NNEW(advfa,double,*nface);
590 NNEW(hfa,double,3**nface);
591
592 NNEW(ap,double,*nface);
593 NNEW(bp,double,*nface);
594
595 NNEW(au,double,*nflnei+*nef);
596 NNEW(ad,double,*nef);
597 NNEW(b,double,*nef);
598 NNEW(bcp,double,*nef);
599
600 NNEW(auv,double,*nflnei+*nef);
601
602 NNEW(bv,double,3**nef);
603 NNEW(bvcp,double,3**nef);
604 NNEW(hel,double,3**nef);
605 NNEW(sel,double,3**nef);
606
607 NNEW(rwork,double,*nef);
608
609 NNEW(inum,ITG,*nk);
610
611 // NNEW(velo,double,8**nef);
612 // NNEW(veloo,double,8**nef);
613
614 /* initializing velo and veloo */
615
616 if(*initial==1){
617 memcpy(&veloo[0],&vel[0],sizeof(double)*8**nef);
618 memcpy(&velo[0],&vel[0],sizeof(double)*8**nef);
619 }
620
621 /* check output requests */
622
623 if((strcmp1(&filab[1914],"MACH")==0)||
624 (strcmp1(&filab[3132],"PTF")==0)||
625 (strcmp1(&filab[3219],"TTF")==0)){
626 imach=1;
627 }
628
629 if((strcmp1(&filab[3132],"PTF")==0)||
630 (strcmp1(&filab[3219],"TTF")==0)){
631 ikappa=1;
632 }
633
634 if((strcmp1(&filab[2088],"TURB")==0)&&(*iturbulent>0)){
635 iturb=1;
636 }
637
638 for(i=0;i<*nprint;i++){
639 if(imach==0){
640 if((strcmp1(&prlab[6*i],"MACH")==0)||
641 (strcmp1(&prlab[6*i],"PTF")==0)||
642 (strcmp1(&prlab[6*i],"TTF")==0)){
643 imach=1;
644 }
645 }
646 if(ikappa==0){
647 if((strcmp1(&prlab[6*i],"PTF")==0)||
648 (strcmp1(&prlab[6*i],"TTF")==0)){
649 ikappa=1;
650 }
651 }
652 if(iturb==0){
653 if(strcmp1(&prlab[6*i],"TURB")==0){
654 iturb=1;
655 }
656 }
657 }
658
659 iincf=0;
660
661 /* if the user has specified a fixed fluid time increment,use it */
662
663 if(*tincf>0.){
664 dtimef=*tincf;
665 ifixdtimef=1;
666 }
667
668 printf("time increment for the CFD-calculations = %e\n\n",dtimef);
669
670 ttimef=*ttime;
671 timef=*time-*dtime;
672
673 if(compressible==0){
674 a1=1.5/dtimef;
675 a2=-2./dtimef;
676 a3=0.5/dtimef;
677 }else{
678 a1=1./dtimef;
679 a2=-a1;
680 a3=0.;
681 }
682
683 iito=iitf;
684
685 NNEW(temp,double,8**nef);
686 NNEW(gamma,double,*nface);
687
688 icounter=0;
689
690 do{
691
692 iincf++;
693
694 // printf("fluid increment = %d\n",iincf);
695
696 if((iincf/ipower)*ipower==iincf){
697 printf("fluid increment = %d\n",iincf);
698
699 /* reevaluating the time increment size
700 only for steady state compressible calculations*/
701
702 if((*nmethod==1)&&(ifixdtimef==0)){
703
704 if(*ithermal>0){
705 NNEW(hcel,double,*nef);
706 FORTRAN(calchcel,(vel,cocon,ncocon,ielmatf,ntmat_,mi,
707 hcel,nef));
708 }
709
710 dtimefold=dtimef;
711 FORTRAN(newtincf,(ithermal,&dtimef,&compressible,vel,
712 hcel,umel,cvel,h,sc,iturbulent,ipkonf,
713 nmethod,nef,lakonf,xxn,ipnei));
714
715 if(compressible==1){
716 a1=1./dtimef;
717 a2=-a1;
718 }else{
719 beta=dtimef/dtimefold;
720 a1=(2.+beta)/(1.+beta)/dtimef;
721 a2=-(1.+beta)/beta/dtimef;
722 a3=1./(beta*(1.+beta))/dtimef;
723 }
724 if(*ithermal>0) SFREE(hcel);
725 }
726 ipower*=2;
727 }
728
729 ierrmax=0;
730
731 timef+=dtimef;
732 if((*time<timef)&&(*nmethod==4)){
733 dtimefold=dtimef;
734 dtimef-=(timef-*time);
735 timef=*time;
736 last=1;
737 beta=dtimef/dtimefold;
738 a1=(2.+beta)/(1.+beta)/dtimef;
739 a2=-(1.+beta)/beta/dtimef;
740 a3=1./(beta*(1.+beta))/dtimef;
741 }
742
743 /* starting iterations till convergence of the fluid increment */
744
745 iit=0;
746
747 do{
748 iit++;
749
750 // printf(" iteration = %d\n",iit);
751
752 /* conditions for transient calculations */
753
754 if(*nmethod==4){
755
756 /* boundary conditions at end of fluid increment */
757
758 FORTRAN(tempload,(xforcold,xforc,xforcact,iamforc,nforc,
759 xloadold,xload,xloadact,iamload,nload,ibody,xbody,nbody,
760 xbodyold,xbodyact,t1old,t1,t1act,iamt1,nk,amta,
761 namta,nam,ampli,&timef,&reltimef,&ttimef,&dtimef,ithermal,nmethod,
762 xbounold,xboun,xbounact,iamboun,nboun,
763 nodeboun,ndirboun,nodeforc,ndirforc,istep,iinc,
764 co,vold,itg,ntg,amname,ikboun,ilboun,nelemload,sideload,mi,
765 ntrans,trab,inotr,vold,integerglob,doubleglob,tieset,istartset,
766 iendset,ialset,ntie,nmpc,ipompc,ikmpc,ilmpc,nodempc,coefmpc,
767 ipobody,iponoel,inoel,ipkon,kon,ielprop,prop,ielmat,
768 shcon,nshcon,rhcon,nrhcon,cocon,ncocon,ntmat_,lakon,set,nset));
769
770 /* body forces (gravity,centrifugal and Coriolis forces) */
771
772 if(*nbody>0){
773 FORTRAN(calcbody,(nef,body,ipobody,ibody,xbody,coel,vel,lakonf,
774 nactdohinv));
775 }
776
777 }else if(icent==1){
778
779 /* body forces (gravity,centrifugal and Coriolis forces;
780 only if centrifugal forces are active => the ensuing
781 Coriolis forces depend on the actual velocity) */
782
783 FORTRAN(calcbody,(nef,body,ipobody,ibody,xbody,coel,vel,lakonf,
784 nactdohinv));
785 }
786
787 /* updating of the material properties */
788
789 if(*ithermal>0){
790
791 if(compressible==0){
792
793 /* calculating material data
794 density (elements+faces)
795 heat capacity at constant volume (elements+faces)
796 dynamic viscosity (elements+faces)
797 heat conduction (faces) */
798
799 FORTRAN(materialdata_cfd,(nef,vel,shcon,nshcon,ielmatf,
800 ntmat_,mi,cvel,vfa,cocon,ncocon,physcon,cvfa,
801 ithermal,nface,umel,umfa,ielfa,hcfa,rhcon,nrhcon));
802
803 }else{
804
805 /* calculating material data
806 heat capacity at constant volume (elements+faces)
807 dynamic viscosity (elements+faces)
808 heat conduction (faces) */
809
810 materialdata_cfd_compmain(nef,vel,shcon,nshcon,ielmatf,
811 ntmat_,mi,cvel,vfa,cocon,ncocon,
812 physcon,cvfa,ithermal,nface,umel,
813 umfa,ielfa,hcfa,&num_cpus);
814 }
815
816 }
817
818 /* filling the lhs and rhs's for the balance of momentum
819 equations */
820
821 if(icounter==0){
822 DMEMSET(auv,0,*nflnei+*nef,0.);
823 DMEMSET(bv,0,3**nef,0.);
824 }
825
826 if(compressible==0){
827
828 /* calculate gamma (Ph.D. Thesis Jasak) */
829
830 calcgammavmain(nface,ielfa,vel,gradvel,gamma,xlet,xxj,ipnei,
831 &betam,nef,flux,&num_cpus);
832
833 mafillvmain(nef,ipnei,neifa,neiel,vfa,xxn,area,
834 auv,&auv[*nflnei],jq,irow,&nzs,bv,vel,cosa,umfa,
835 alet,ale,gradvfa,xxi,
836 body,volume,ielfa,lakonf,ifabou,nbody,
837 &dtimef,velo,veloo,sel,xrlfa,gamma,xxj,nactdohinv,&a1,
838 &a2,&a3,flux,&icyclic,c,ifatie,iau6,xxna,xxnj,
839 iturbulent,gradvel,of2,yy,umel,&ncfd,inlet,sc);
840
841 }else{
842
843 /* convection scheme */
844
845 if(ischeme==1){
846 hrv_udmain(nface,ielfa,vel,ipnei,nef,flux,vfa,&num_cpus,
847 xxi,xle,gradvel,neij);
848 }else{
849 hrv_mod_smartmain(nface,ielfa,vel,gradvel,xlet,xxj,ipnei,
850 nef,flux,vfa,&num_cpus);
851 }
852
853 // printf("mafillvcompmain\n");
854 mafillvcompmain(nef,ipnei,neifa,neiel,vfa,xxn,area,
855 auv,&auv[*nflnei],jq,irow,&nzs,bv,vel,cosa,umfa,
856 alet,ale,gradvfa,xxi,
857 body,volume,ielfa,lakonf,ifabou,nbody,
858 &dtimef,velo,veloo,sel,xrlfa,gamma,xxj,nactdohinv,&a1,
859 &a2,&a3,flux,&icyclic,c,ifatie,iau6,xxna,xxnj,
860 iturbulent,gradvel,of2,yy,umel,&ncfd,inlet,sc);
861 }
862
863 for(i=0;i<*nef;i++){rwork[i]=1./auv[*nflnei+i];}
864
865 /* underrelaxation of the velocity only for compressible
866 simple scheme */
867
868 if((compressible==1)&&(isimplec==0)){
869 memcpy(&temp[*nef],&vel[*nef],sizeof(double)*ncfd**nef);
870 }
871
872 /* reordering the lhs (getting rid of zeros) */
873
874 reorderlhsmain(auv,am,iamorig,&nz_num,&num_cpus);
875
876 /* modifying the rhs (taking common contributions
877 between the cpu-blocks into account) */
878
879 memcpy(&bvcp[0],&bv[0],sizeof(double)*3**nef);
880 for(i=0;i<ncfd;i++){
881 FORTRAN(reorderrhs,(auv,&bv[i**nef],&vel[(i+1)**nef],neighblock,&nneighblock));
882 }
883
884 /* calculation of the momentum residual */
885
886 calcresvfluidmain(nef,am,&bv[0],&auv[*nflnei],iam,jam,
887 &vel[*nef],&relnormv,nestart,&num_cpus,
888 &ncfd);
889
890 for(k=1;k<num_cpus+1;k++){
891 if(k>1){
892 memcpy(&bv[0],&bvcp[0],sizeof(double)*3**nef);
893 }
894 for(i=0;i<ncfd;i++){
895 if(k>1){
896 FORTRAN(reorderrhs,(auv,&bv[i**nef],&vel[(i+1)**nef],neighblock,
897 &nneighblock));
898 }
899 dgmresmain(nef,&bv[i**nef],&vel[(i+1)**nef],&nz_num,iam,jam,am,
900 &isym,&itol,&tol,&itmax,&iter,
901 &err,&ierr,&iunit,sb,sx,rgwk,&lrgw,igwk,
902 &ligw,rwork,iwork,nestart,&num_cpus,&dgmrestol);
903
904 if(ierr>1){
905 printf("*WARNING in compfluid: error message from predgmres (v_%d)=%d\n",i+1,ierr);
906 if(ierr>ierrmax)ierrmax=ierr;
907 }
908 }
909 }
910
911 /* for(i=0;i<ncfd;i++){
912 spooles(&auv[*nflnei],am,adb,aub,&sigma,&bv[i**nef],iam,jam,
913 nef,&nz_num,&symmetryflag,&inputformat,&nz_num);
914 memcpy(&vel[(i+1)**nef],&bv[i**nef],sizeof(double)**nef);
915 }*/
916
917 /* underrelaxation of the velocity only for compressible
918 simple scheme */
919
920 if((compressible==1)&&(isimplec==0)){
921 for(i=*nef;i<(ncfd+1)**nef;i++){vel[i]=0.8*vel[i]+0.2*temp[i];}
922 }
923
924 // printf("iitpt=%d,iitp=%d\n",iitpt,iitp);
925
926 for(j=0;j<iitpt;j++){
927 // printf("j=%d,iitpt=%d\n",j,iitpt);
928
929 /* generating vol/ad and v* at the face centers (advfa and hfa) */
930
931 if((compressible==0)||(isimplec==0)){
932 extrapolate_d_v_simplemain(nface,ielfa,xrlfa,&auv[*nflnei],advfa,hfa,
933 &icyclic,c,ifatie,vel,nef,volume,&num_cpus,
934 &ncfd);
935 }else{
936 extrapolate_d_v_simplecmain(nface,ielfa,xrlfa,&auv[*nflnei],advfa,hfa,
937 &icyclic,c,ifatie,vel,nef,volume,auv,
938 ipnei,&num_cpus,&ncfd);
939 }
940
941 /* calculate the mass flow based on the newly calculated
942 velocity */
943
944 calcfluxmain(area,vfa,xxna,ipnei,nef,neifa,flux,xxj,gradpfa,xlet,xle,vel,
945 advfa,ielfa,neiel,ifabou,hfa,&num_cpus);
946
947 /* calculating the lhs and rhs of the equation system to determine
948 p' (balance of mass) */
949
950 if(compressible==0){
951
952 /* incompressible media */
953
954 /* temporarily store the pressure in temp */
955
956 memcpy(&temp[4**nef],&vel[4**nef],sizeof(double)**nef);
957
958 /* first iteration: calculating both lhs and rhs */
959
960 DMEMSET(ad,0,*nef,0.);
961 DMEMSET(au,0,nzs,0.);
962 DMEMSET(b,0,*nef,0.);
963
964 mafillpmain(nef,lakonf,ipnei,neifa,neiel,vfa,area,
965 advfa,xlet,cosa,volume,au,ad,jq,irow,ap,
966 ielfa,ifabou,xle,b,xxn,nef,
967 &nzs,hfa,gradpel,bp,xxi,neij,xlen,cosb,
968 &iatleastonepressurebc,iau6,xxicn,flux);
969
970 FORTRAN(convert2slapcol,(au,ad,jq,&nzs,nef,aua));
971
972 nelt=nzs+*nef;
973 isym=1;
974
975 /* next line was changed from 10 to 3 on 22.12.2016 */
976
977 nsave=3;
978 itol=0;
979 tol=1.e-6;
980
981 /* next line was changed from 110 to 10 on 22.12.2016 */
982
983 itmax=10;
984 iunit=0;
985 lenw=131+17**nef+2*nelt;
986 NNEW(rgwk,double,lenw);
987 leniw=32+4**nef+2*nelt;
988 NNEW(igwk,ITG,leniw);
989
990 /* initial guess: 0 */
991
992 DMEMSET(vel,4**nef,5**nef,0.);
993
994 FORTRAN(dslugm,(nef,&b[0],&vel[4**nef],&nelt,ia,ja,aua,
995 &isym,&nsave,&itol,&tol,&itmax,&iter,
996 &err,&ierr,&iunit,rgwk,&lenw,igwk,&leniw));
997 SFREE(rgwk);SFREE(igwk);
998
999 /* non-orthogonal pressure correction */
1000
1001 /* calculate the p' gradient at the
1002 face centers */
1003
1004 for(ip=0;ip<iitp;ip++){
1005
1006 iflag=1;
1007 extrapol_dpelmain(nface,ielfa,xrlfa,vel,vfa,
1008 ifabou,xbounact,nef,gradpcel,gradpcfa,neifa,rf,
1009 area,volume,
1010 xle,xxi,&icyclic,xxn,ipnei,ifatie,
1011 coefmpc,nmpc,labmpc,ipompc,nodempc,ifaext,nfaext,
1012 nactdoh,&iflag,xxj,xlet,c,&num_cpus,xxna,&ncfd,
1013 &ip);
1014
1015 /* update the right hand side (taking skewness of
1016 elements into account) */
1017
1018 DMEMSET(b,0,*nef,0.);
1019 rhspmain(nef,lakonf,ipnei,neifa,neiel,vfa,area,
1020 advfa,xlet,cosa,volume,au,ad,jq,irow,ap,ielfa,ifabou,xle,
1021 b,xxn,nef,&nzs,hfa,gradpel,bp,xxi,neij,xlen,
1022 &iatleastonepressurebc,xxicn,flux,xxnj,gradpcfa,cosb);
1023
1024 /* calculate a better pressure correction p' */
1025
1026 nelt=nzs+*nef;
1027 isym=1;
1028 nsave=3;
1029 itol=0;
1030 tol=1.e-6;
1031 itmax=10;
1032 iunit=0;
1033 lenw=131+17**nef+2*nelt;
1034 NNEW(rgwk,double,lenw);
1035 leniw=32+4**nef+2*nelt;
1036 NNEW(igwk,ITG,leniw);
1037
1038 /* initial guess: 0 */
1039
1040 DMEMSET(vel,4**nef,5**nef,0.);
1041
1042 FORTRAN(dslugm,(nef,&b[0],&vel[4**nef],&nelt,ia,ja,aua,
1043 &isym,&nsave,&itol,&tol,&itmax,&iter,
1044 &err,&ierr,&iunit,rgwk,&lenw,igwk,&leniw));
1045 SFREE(rgwk);SFREE(igwk);
1046
1047 }// end loop iitp
1048
1049 /* calculate the p' gradient at the
1050 element centers */
1051
1052 /* 16.04.2020: should be iflag=1????? */
1053
1054
1055 iflag=0;
1056 extrapol_dpelmain(nface,ielfa,xrlfa,vel,vfa,
1057 ifabou,xbounact,nef,gradpcel,gradpcfa,neifa,rf,area,volume,
1058 xle,xxi,&icyclic,xxn,ipnei,ifatie,
1059 coefmpc,nmpc,labmpc,ipompc,nodempc,ifaext,nfaext,
1060 nactdoh,&iflag,xxj,xlet,c,&num_cpus,xxna,&ncfd,&ip);
1061
1062 /* correction of the velocity v* at the element centers due
1063 to the pressure change in order to get v** */
1064
1065 correctvelmain(&auv[*nflnei],nef,volume,gradpcel,vel,&num_cpus);
1066
1067 /* extrapolating the velocity from the elements centers to the face
1068 centers,thereby taking the boundary conditions into account */
1069
1070 extrapol_velmain(nface,ielfa,xrlfa,vel,vfa,
1071 ifabou,xbounact,ipnei,nef,&icyclic,c,ifatie,xxn,gradvel,
1072 gradvfa,neifa,rf,area,volume,xle,xxi,xxj,xlet,
1073 coefmpc,nmpc,labmpc,ipompc,nodempc,ifaext,nfaext,nactdoh,
1074 &iitg,&num_cpus,&compressible,xxna,&ncfd,cofa);
1075
1076 /* correcting the flux to get mf** */
1077
1078 correctfluxmain(nef,ipnei,neifa,neiel,flux,vfa,advfa,area,
1079 vel,alet,ielfa,ale,ifabou,&num_cpus,xxnj,
1080 gradpcfa);
1081
1082 /* correcting the pressure to get p* */
1083
1084 /* underrelaxation always except for the compressible simplec
1085 scheme */
1086
1087 for(i=4**nef;i<5**nef;i++){vel[i]=0.2*vel[i]+temp[i];}
1088
1089 extrapol_pelmain(nface,ielfa,xrlfa,vel,vfa,
1090 ifabou,xbounact,nef,gradpel,gradpfa,neifa,rf,area,volume,
1091 xle,xxi,&icyclic,xxn,ipnei,ifatie,
1092 coefmpc,nmpc,labmpc,ipompc,nodempc,ifaext,nfaext,nactdoh,
1093 &iitg,c,&num_cpus,&compressible,xxna,&ncfd);
1094
1095 }else{
1096
1097 /* compressible media */
1098
1099 /* temporarily store the pressure in temp */
1100
1101 memcpy(&temp[4**nef],&vel[4**nef],sizeof(double)**nef);
1102
1103 DMEMSET(au,0,*nflnei+*nef,0.);
1104 DMEMSET(b,0,*nef,0.);
1105
1106 // printf("mafillpcompmain\n");
1107 mafillpcompmain(nef,lakonf,ipnei,neifa,neiel,vfa,area,
1108 advfa,xlet,cosa,volume,au,&au[*nflnei],jq,irow,ap,
1109 ielfa,ifabou,xle,b,xxn,nef,
1110 &nzs,hfa,gradpel,bp,xxi,neij,xlen,cosb,
1111 ielmatf,mi,&a1,&a2,&a3,velo,veloo,&dtimef,shcon,
1112 ntmat_,vel,nactdohinv,xrlfa,flux,iau6,xxicn,
1113 gamma,inlet);
1114
1115 for(i=0;i<*nef;i++){rwork[i]=1./au[*nflnei+i];}
1116
1117 /* initial guess: 0 */
1118
1119 DMEMSET(vel,4**nef,5**nef,0.);
1120
1121 /* reordering the lhs (getting rid of zeros) */
1122
1123 reorderlhsmain(au,am,iamorig,&nz_num,&num_cpus);
1124
1125 /* calculation of the mass conservation residual */
1126
1127 FORTRAN(calcrespfluid,(nef,&b[0],&au[*nflnei],
1128 &temp[4**nef],&relnormp));
1129
1130 /* no change of rhs (reorderrhs) since initial guess is zero */
1131
1132 for(k=1;k<num_cpus+1;k++){
1133 if(k==1){
1134 memcpy(&bcp[0],&b[0],sizeof(double)**nef);
1135 }else{
1136 memcpy(&b[0],&bcp[0],sizeof(double)**nef);
1137 FORTRAN(reorderrhs,(au,&b[0],&vel[4**nef],neighblock,
1138 &nneighblock));
1139 }
1140 dgmresmain(nef,&b[0],&vel[4**nef],&nz_num,iam,jam,am,
1141 &isym,&itol,&tol,&itmax,&iter,
1142 &err,&ierr,&iunit,sb,sx,rgwk,&lrgw,igwk,
1143 &ligw,rwork,iwork,nestart,&num_cpus,&dgmrestol);
1144
1145 if(ierr>1){
1146 printf("*WARNING in compfluid: error message from predgmres (p)=%d\n",ierr);
1147 if(ierr>ierrmax)ierrmax=ierr;
1148 }
1149 }
1150
1151 /* spooles(&au[*nflnei],am,adb,aub,&sigma,&b[0],iam,jam,
1152 nef,&nz_num,&symmetryflag,&inputformat,&nz_num);
1153 memcpy(&vel[4**nef],&b[0],sizeof(double)**nef);*/
1154
1155 /* non-orthogonal pressure correction */
1156
1157 for(ip=0;ip<iitp;ip++){
1158
1159 /* calculate the p' gradient at the
1160 face centers */
1161
1162 iflag=1;
1163 extrapol_dpelmain(nface,ielfa,xrlfa,vel,vfa,
1164 ifabou,xbounact,nef,gradpcel,gradpcfa,neifa,rf,area,volume,
1165 xle,xxi,&icyclic,xxn,ipnei,ifatie,
1166 coefmpc,nmpc,labmpc,ipompc,nodempc,ifaext,nfaext,
1167 nactdoh,&iflag,xxj,xlet,c,&num_cpus,xxna,&ncfd,&ip);
1168
1169 /* update the right hand side (taking skewness of
1170 elements into account) */
1171
1172 DMEMSET(b,0,*nef,0.);
1173
1174 rhspcompmain(nef,lakonf,ipnei,neifa,neiel,vfa,area,
1175 advfa,xlet,cosa,volume,au,&au[*nflnei],jq,irow,ap,
1176 ielfa,ifabou,xle,b,xxn,nef,
1177 &nzs,hfa,gradpel,bp,xxi,neij,xlen,cosb,
1178 ielmatf,mi,&a1,&a2,&a3,velo,veloo,&dtimef,shcon,
1179 ntmat_,vel,nactdohinv,xrlfa,flux,iau6,xxicn,
1180 gamma,xxnj,gradpcfa);
1181
1182 for(i=0;i<*nef;i++){rwork[i]=1./au[*nflnei+i];}
1183
1184 /* reordering the lhs (getting rid of zeros) */
1185
1186 reorderlhsmain(au,am,iamorig,&nz_num,&num_cpus);
1187
1188 /* modifying the rhs (taking common contributions
1189 between the cpu-blocks into account) */
1190
1191 memcpy(&bcp[0],&b[0],sizeof(double)**nef);
1192 FORTRAN(reorderrhs,(au,&b[0],&vel[4**nef],neighblock,&nneighblock));
1193
1194 /* initial guess: 0 */
1195
1196 DMEMSET(vel,4**nef,5**nef,0.);
1197
1198 /* calculation of the mass conservation residual */
1199
1200 FORTRAN(calcrespfluid,(nef,&b[0],&au[*nflnei],
1201 &temp[4**nef],&relnormp));
1202
1203 for(k=1;k<num_cpus+1;k++){
1204 if(k>1){
1205 memcpy(&b[0],&bcp[0],sizeof(double)**nef);
1206 FORTRAN(reorderrhs,(au,&b[0],&vel[4**nef],neighblock,
1207 &nneighblock));
1208 }
1209 dgmresmain(nef,&b[0],&vel[4**nef],&nz_num,iam,jam,am,
1210 &isym,&itol,&tol,&itmax,&iter,
1211 &err,&ierr,&iunit,sb,sx,rgwk,&lrgw,igwk,
1212 &ligw,rwork,iwork,nestart,&num_cpus,&dgmrestol);
1213
1214 if(ierr>1){
1215 printf("*WARNING in compfluid: error message from predgmres (p)=%d\n",ierr);
1216 if(ierr>ierrmax)ierrmax=ierr;
1217 }
1218 }
1219
1220 /* spooles(&au[*nflnei],am,adb,aub,&sigma,&b[0],iam,jam,
1221 nef,&nz_num,&symmetryflag,&inputformat,&nz_num);
1222 memcpy(&vel[4**nef],&b[0],sizeof(double)**nef);*/
1223
1224 } // end loop iitp
1225
1226 /* calculate the p' gradient at the
1227 element centers */
1228
1229 /* 16.04.2020: should be iflag=1????? */
1230
1231 iflag=0;
1232 extrapol_dpelmain(nface,ielfa,xrlfa,vel,vfa,
1233 ifabou,xbounact,nef,gradpcel,gradpcfa,neifa,rf,area,volume,
1234 xle,xxi,&icyclic,xxn,ipnei,ifatie,
1235 coefmpc,nmpc,labmpc,ipompc,nodempc,ifaext,nfaext,
1236 nactdoh,&iflag,xxj,xlet,c,&num_cpus,xxna,&ncfd,&ip);
1237
1238 /* correction of the velocity v* at the element centers due
1239 to the pressure change in order to get v** */
1240
1241 // if((compressible==0)||(isimplec==0)){
1242 correctvelmain(&auv[*nflnei],nef,volume,gradpcel,vel,&num_cpus);
1243 /* }else{
1244 FORTRAN(correctvel_simplec,(&auv[*nflnei],nef,volume,gradpcel,vel,
1245 ipnei,auv));
1246 }*/
1247
1248 /* extrapolating the velocity from the elements centers to the face
1249 centers,thereby taking the boundary conditions into account */
1250
1251 extrapol_velmain(nface,ielfa,xrlfa,vel,vfa,
1252 ifabou,xbounact,ipnei,nef,&icyclic,c,ifatie,xxn,gradvel,
1253 gradvfa,neifa,rf,area,volume,xle,xxi,xxj,xlet,
1254 coefmpc,nmpc,labmpc,ipompc,nodempc,ifaext,nfaext,nactdoh,
1255 &iitg,&num_cpus,&compressible,xxna,&ncfd,cofa);
1256
1257 /* correcting the flux to get mf** */
1258
1259 correctfluxcompmain(nef,ipnei,neifa,neiel,flux,vfa,advfa,area,
1260 vel,alet,ielfa,ale,ifabou,ielmatf,mi,shcon,
1261 ntmat_,&num_cpus,xxnj,gradpcfa,inlet);
1262
1263 /* correcting the pressure to get p* */
1264
1265 /* underrelaxation always except for the compressible simplec
1266 scheme */
1267
1268 if((compressible==1)&&(isimplec==1)){
1269 for(i=4**nef;i<5**nef;i++){vel[i]=vel[i]+temp[i];}
1270 }else{
1271 for(i=4**nef;i<5**nef;i++){vel[i]=0.2*vel[i]+temp[i];}
1272 }
1273
1274 extrapol_pelmain(nface,ielfa,xrlfa,vel,vfa,
1275 ifabou,xbounact,nef,gradpel,gradpfa,neifa,rf,area,volume,
1276 xle,xxi,&icyclic,xxn,ipnei,ifatie,
1277 coefmpc,nmpc,labmpc,ipompc,nodempc,ifaext,nfaext,nactdoh,
1278 &iitg,c,&num_cpus,&compressible,xxna,&ncfd);
1279
1280 /* calculating the lhs and rhs of the energy equation */
1281
1282 DMEMSET(ad,0,*nef,0.);
1283 DMEMSET(au,0,*nflnei+*nef,0.);
1284 DMEMSET(b,0,*nef,0.);
1285
1286 /* underrelaxation of the temperature only for the
1287 compressible simple scheme */
1288
1289 if(isimplec==0){
1290 memcpy(&temp[0],&vel[0],sizeof(double)**nef);
1291 }
1292
1293 /* convective scheme */
1294
1295 // if(ischeme==1){
1296 hrt_udmain(nface,ielfa,vel,ipnei,nef,flux,vfa,&num_cpus,
1297 xxi,xle,gradtel,neij);
1298 /* }else{
1299 FORTRAN(hrt_mod_smart,(nface,ielfa,vel,gradtel,gamma,
1300 xlet,xxn,xxj,ipnei,&betam,nef,flux,vfa));
1301 }*/
1302
1303 // printf("mafilltcompmain\n");
1304 mafilltcompmain(nef,ipnei,neifa,neiel,vfa,xxn,area,
1305 au,&au[*nflnei],jq,irow,&nzs,b,vel,umel,alet,ale,gradtfa,xxi,
1306 body,volume,ielfa,lakonf,ifabou,nbody,nef,
1307 &dtimef,velo,veloo,cvfa,hcfa,cvel,gradvel,xload,gamma,
1308 xrlfa,xxj,nactdohinv,&a1,&a2,&a3,flux,iau6,xxni,xxnj,
1309 iturbulent,of2,sc);
1310
1311 for(i=0;i<*nef;i++){rwork[i]=1./au[*nflnei+i];}
1312
1313 /* reordering the lhs (getting rid of zeros) */
1314
1315 reorderlhsmain(au,am,iamorig,&nz_num,&num_cpus);
1316
1317 /* modifying the rhs (taking common contributions
1318 between the cpu-blocks into account) */
1319
1320 memcpy(&bcp[0],&b[0],sizeof(double)**nef);
1321 FORTRAN(reorderrhs,(au,&b[0],&vel[0],neighblock,&nneighblock));
1322
1323 /* calculation of the energy residual */
1324
1325 calcrestfluidmain(nef,am,&b[0],&au[*nflnei],iam,jam,
1326 &vel[0],&relnormt,nestart,&num_cpus);
1327
1328 for(k=1;k<num_cpus+1;k++){
1329 if(k>1){
1330 memcpy(&b[0],&bcp[0],sizeof(double)**nef);
1331 FORTRAN(reorderrhs,(au,&b[0],&vel[0],neighblock,&nneighblock));
1332 }
1333 dgmresmain(nef,&b[0],&vel[0],&nz_num,iam,jam,am,
1334 &isym,&itol,&tol,&itmax,&iter,
1335 &err,&ierr,&iunit,sb,sx,rgwk,&lrgw,igwk,
1336 &ligw,rwork,iwork,nestart,&num_cpus,&dgmrestol);
1337
1338 if(ierr>1){
1339 printf("*WARNING in compfluid: error message from predgmres (T)=%d\n",ierr);
1340 if(ierr>ierrmax)ierrmax=ierr;
1341 }
1342 }
1343
1344 /* spooles(&au[*nflnei],am,adb,aub,&sigma,&b[0],iam,jam,
1345 nef,&nz_num,&symmetryflag,&inputformat,&nz_num);
1346 memcpy(&vel[0],&b[0],sizeof(double)**nef);*/
1347
1348 /* underrelaxation of the temperature only for compressible
1349 simple scheme */
1350
1351 if(isimplec==0){
1352 for(i=0;i<*nef;i++){vel[i]=0.8*vel[i]+0.2*temp[i];}
1353 }
1354
1355 /* extrapolation of the temperature at the element centers
1356 to the face centers */
1357
1358 extrapol_telmain(nface,ielfa,xrlfa,vel,vfa,
1359 ifabou,xbounact,nef,gradtel,gradtfa,neifa,rf,area,volume,
1360 xle,xxi,&icyclic,xxn,ipnei,ifatie,xload,xlet,xxj,
1361 coefmpc,nmpc,labmpc,ipompc,nodempc,ifaext,nfaext,nactdoh,
1362 &iitg,c,&num_cpus,&compressible,xxna,&ncfd);
1363
1364 /* recalculating the density */
1365
1366 /* calculating the density at the element centers */
1367
1368 calcrhoelcompmain(nef,vel,shcon,ielmatf,ntmat_,mi,
1369 &num_cpus);
1370
1371 /* calculating the density at the face centers */
1372
1373 if(ischeme==1){
1374 hrr_udmain(nface,vfa,shcon,ielmatf,ntmat_,
1375 mi,ielfa,ipnei,vel,nef,flux,
1376 &num_cpus,xxi,xle,gradpel,
1377 gradtel,neij);
1378 }else{
1379 hrr_mod_smartmain(nface,vfa,shcon,ielmatf,ntmat_,
1380 mi,ielfa,ipnei,vel,nef,flux,
1381 gradpel,gradtel,xxj,xlet,
1382 &num_cpus);
1383 }
1384
1385 }
1386
1387 }
1388
1389 if((*ithermal>0)&&(compressible==0)){
1390
1391 /* calculating the lhs and rhs of the energy equation */
1392
1393 DMEMSET(ad,0,*nef,0.);
1394 DMEMSET(au,0,*nflnei+*nef,0.);
1395 DMEMSET(b,0,*nef,0.);
1396
1397 /* calculate gamma (Ph.D. Thesis Jasak) */
1398
1399 calcgammatmain(nface,ielfa,vel,gradtel,gamma,xlet,xxj,ipnei,
1400 &betam,nef,flux,&num_cpus);
1401
1402 mafilltmain(nef,ipnei,neifa,neiel,vfa,xxn,area,
1403 au,&au[*nflnei],jq,irow,&nzs,b,vel,umel,alet,ale,gradtfa,xxi,
1404 body,volume,ielfa,lakonf,ifabou,nbody,nef,
1405 &dtimef,velo,veloo,cvfa,hcfa,cvel,gradvel,xload,gamma,
1406 xrlfa,xxj,nactdohinv,&a1,&a2,&a3,flux,iau6,xxni,xxnj,
1407 iturbulent,of2,sc);
1408
1409 for(i=0;i<*nef;i++){rwork[i]=1./au[*nflnei+i];}
1410
1411 /* reordering the lhs (getting rid of zeros) */
1412
1413 reorderlhsmain(au,am,iamorig,&nz_num,&num_cpus);
1414
1415 /* modifying the rhs (taking common contributions
1416 between the cpu-blocks into account) */
1417
1418 FORTRAN(reorderrhs,(au,&b[0],&vel[0],neighblock,&nneighblock));
1419
1420 /* calculation of the energy residual */
1421
1422 calcrestfluidmain(nef,am,&b[0],&au[*nflnei],iam,jam,
1423 &vel[0],&relnormt,nestart,&num_cpus);
1424
1425 dgmresmain(nef,&b[0],&vel[0],&nz_num,iam,jam,am,
1426 &isym,&itol,&tol,&itmax,&iter,
1427 &err,&ierr,&iunit,sb,sx,rgwk,&lrgw,igwk,
1428 &ligw,rwork,iwork,nestart,&num_cpus,&dgmrestol);
1429
1430 if(ierr>1){
1431 printf("*WARNING in compfluid: error message from predgmres (T)=%d\n",ierr);
1432 if(ierr>ierrmax)ierrmax=ierr;
1433 }
1434
1435 /* extrapolation of the temperature at the element centers
1436 to the face centers */
1437
1438 extrapol_telmain(nface,ielfa,xrlfa,vel,vfa,
1439 ifabou,xbounact,nef,gradtel,gradtfa,neifa,rf,area,volume,
1440 xle,xxi,&icyclic,xxn,ipnei,ifatie,xload,xlet,xxj,
1441 coefmpc,nmpc,labmpc,ipompc,nodempc,ifaext,nfaext,nactdoh,
1442 &iitg,c,&num_cpus,&compressible,xxna,&ncfd);
1443
1444 /* recalculating the density */
1445
1446 /* calculating the density at the element centers */
1447
1448 FORTRAN(calcrhoel,(nef,vel,rhcon,nrhcon,ielmatf,ntmat_,
1449 ithermal,mi));
1450
1451 /* calculating the density at the face centers */
1452
1453 FORTRAN(calcrhofa,(nface,vfa,rhcon,nrhcon,ielmatf,ntmat_,
1454 ithermal,mi,ielfa));
1455
1456 }
1457
1458 if(*iturbulent>0){
1459
1460 /* calculating the lhs and rhs of the k-equation */
1461
1462 DMEMSET(au,0,*nflnei+*nef,0.);
1463 DMEMSET(b,0,*nef,0.);
1464
1465 if(compressible==0){
1466
1467 /* calculate gamma (Ph.D. Thesis Jasak) */
1468
1469 FORTRAN(calcgammak,(nface,ielfa,vel,gradkel,gamma,xlet,xxn,xxj,
1470 ipnei,&betam,nef,flux));
1471
1472 mafillkmain(nef,ipnei,neifa,neiel,vfa,xxn,area,
1473 au,&au[*nflnei],jq,irow,&nzs,b,vel,umfa,alet,ale,gradkfa,xxi,
1474 body,volume,ielfa,lakonf,ifabou,nbody,nef,
1475 &dtimef,velo,veloo,cvfa,hcfa,cvel,gradvel,xload,gamma,
1476 xrlfa,xxj,nactdohinv,&a1,&a2,&a3,flux,iau6,xxni,xxnj,
1477 iturbulent,f1,of2,yy,umel,gradkel,gradoel,sc);
1478 }else{
1479
1480 /* underrelaxation of the temperature only for the
1481 compressible simple scheme */
1482
1483 if(isimplec==0){
1484 memcpy(&temp[6**nef],&vel[6**nef],sizeof(double)**nef);
1485 }
1486
1487 /* convective scheme */
1488
1489 // if(ischeme==1){
1490 hrk_udmain(nface,ielfa,vel,ipnei,nef,flux,vfa,&num_cpus);
1491 /* }else{
1492 FORTRAN(hrk_mod_smart,(nface,ielfa,vel,gradtel,gamma,
1493 xlet,xxn,xxj,ipnei,&betam,nef,flux,vfa));
1494 }*/
1495
1496 mafillkcompmain(nef,ipnei,neifa,neiel,vfa,xxn,area,
1497 au,&au[*nflnei],jq,irow,&nzs,b,vel,umfa,alet,ale,gradkfa,xxi,
1498 body,volume,ielfa,lakonf,ifabou,nbody,nef,
1499 &dtimef,velo,veloo,cvfa,hcfa,cvel,gradvel,xload,
1500 xrlfa,xxj,nactdohinv,&a1,&a2,&a3,flux,iau6,xxni,xxnj,
1501 iturbulent,f1,of2,yy,umel,gradkel,gradoel,inlet,sc);
1502 }
1503
1504 for(i=0;i<*nef;i++){rwork[i]=1./au[*nflnei+i];}
1505
1506 /* reordering the lhs (getting rid of zeros) */
1507
1508 reorderlhsmain(au,am,iamorig,&nz_num,&num_cpus);
1509
1510 /* modifying the rhs (taking common contributions
1511 between the cpu-blocks into account) */
1512
1513 FORTRAN(reorderrhs,(au,&b[0],&vel[6**nef],neighblock,&nneighblock));
1514
1515 dgmresmain(nef,&b[0],&vel[6**nef],&nz_num,iam,jam,am,
1516 &isym,&itol,&tol,&itmax,&iter,
1517 &err,&ierr,&iunit,sb,sx,rgwk,&lrgw,igwk,
1518 &ligw,rwork,iwork,nestart,&num_cpus,&dgmrestol);
1519 if(ierr>1){
1520 printf("*WARNING in compfluid: error message from predgmres (k)=%d\n",ierr);
1521 if(ierr>ierrmax)ierrmax=ierr;
1522 }
1523
1524 /* underrelaxation of the temperature only for compressible
1525 simple scheme */
1526
1527 if((compressible==1)&&(isimplec==0)){
1528 for(i=6**nef;i<7**nef;i++){vel[i]=0.8*vel[i]+0.2*temp[i];}
1529 }
1530
1531 /* calculating the lhs and rhs of the omega-equation */
1532
1533 DMEMSET(au,0,*nflnei+*nef,0.);
1534 DMEMSET(b,0,*nef,0.);
1535
1536 if(compressible==0){
1537
1538 /* calculate gamma (Ph.D. Thesis Jasak) */
1539
1540 FORTRAN(calcgammao,(nface,ielfa,vel,gradoel,gamma,xlet,xxn,xxj,
1541 ipnei,&betam,nef,flux));
1542
1543 mafillomain(nef,ipnei,neifa,neiel,vfa,xxn,area,
1544 au,&au[*nflnei],jq,irow,&nzs,b,vel,umfa,alet,ale,gradofa,xxi,
1545 body,volume,ielfa,lakonf,ifabou,nbody,nef,
1546 &dtimef,velo,veloo,cvfa,hcfa,cvel,gradvel,xload,gamma,
1547 xrlfa,xxj,nactdohinv,&a1,&a2,&a3,flux,iau6,xxni,xxnj,
1548 iturbulent,f1,of2,gradkel,gradoel,sc);
1549 }else{
1550
1551 /* underrelaxation of the temperature only for the
1552 compressible simple scheme */
1553
1554 if(isimplec==0){
1555 memcpy(&temp[7**nef],&vel[7**nef],sizeof(double)**nef);
1556 }
1557
1558 /* convective scheme */
1559
1560 // if(ischeme==1){
1561 hro_udmain(nface,ielfa,vel,ipnei,nef,flux,vfa,&num_cpus);
1562 /* }else{
1563 FORTRAN(hro_mod_smart,(nface,ielfa,vel,gradtel,gamma,
1564 xlet,xxn,xxj,ipnei,&betam,nef,flux,vfa));
1565 }*/
1566
1567 mafillocompmain(nef,ipnei,neifa,neiel,vfa,xxn,area,
1568 au,&au[*nflnei],jq,irow,&nzs,b,vel,umfa,alet,ale,gradofa,xxi,
1569 body,volume,ielfa,lakonf,ifabou,nbody,nef,
1570 &dtimef,velo,veloo,cvfa,hcfa,cvel,gradvel,xload,
1571 xrlfa,xxj,nactdohinv,&a1,&a2,&a3,flux,iau6,xxni,xxnj,
1572 iturbulent,f1,of2,gradkel,gradoel,inlet,sc);
1573 }
1574
1575 for(i=0;i<*nef;i++){rwork[i]=1./au[*nflnei+i];}
1576
1577 /* reordering the lhs (getting rid of zeros) */
1578
1579 reorderlhsmain(au,am,iamorig,&nz_num,&num_cpus);
1580
1581 /* modifying the rhs (taking common contributions
1582 between the cpu-blocks into account) */
1583
1584 FORTRAN(reorderrhs,(au,&b[0],&vel[7**nef],neighblock,&nneighblock));
1585
1586 dgmresmain(nef,&b[0],&vel[7**nef],&nz_num,iam,jam,am,
1587 &isym,&itol,&tol,&itmax,&iter,
1588 &err,&ierr,&iunit,sb,sx,rgwk,&lrgw,igwk,
1589 &ligw,rwork,iwork,nestart,&num_cpus,&dgmrestol);
1590 if(ierr>1){
1591 printf("*WARNING in compfluid: error message from predgmres (om)=%d\n",ierr);
1592 if(ierr>ierrmax)ierrmax=ierr;
1593 }
1594
1595 /* underrelaxation of the temperature only for compressible
1596 simple scheme */
1597
1598 if((compressible==1)&&(isimplec==0)){
1599 for(i=7**nef;i<8**nef;i++){vel[i]=0.8*vel[i]+0.2*temp[i];}
1600 }
1601
1602 /* extrapolation of the turbulence variables at the element centers
1603 to the face centers */
1604 /* sum=0;
1605 for(i=6**nef;i<7**nef;i++){
1606 printf("i=%d,vel=%e\n",i+1,vel[i]);
1607 if(fabs(vel[i])>sum) sum=fabs(vel[i]);
1608 }
1609 printf("sum=%e\n\n",sum);*/
1610
1611 extrapol_kelmain(nface,ielfa,xrlfa,vel,vfa,
1612 ifabou,xbounact,nef,gradkel,gradkfa,neifa,rf,area,volume,
1613 xle,xxi,&icyclic,xxn,ipnei,ifatie,xlet,xxj,
1614 coefmpc,nmpc,labmpc,ipompc,nodempc,ifaext,nfaext,nactdoh,
1615 umfa,physcon,&iitg,c,&num_cpus,&compressible,xxna,&ncfd,
1616 inlet);
1617
1618 extrapol_oelmain(nface,ielfa,xrlfa,vel,vfa,
1619 ifabou,xbounact,nef,gradoel,gradofa,neifa,rf,area,volume,
1620 xle,xxi,&icyclic,xxn,ipnei,ifatie,xlet,xxj,
1621 coefmpc,nmpc,labmpc,ipompc,nodempc,ifaext,nfaext,nactdoh,
1622 umfa,physcon,&iitg,c,dy,&num_cpus,&compressible,xxna,
1623 &ncfd,inlet);
1624
1625 }
1626
1627 /* extrapolating the velocity from the elements centers to the face
1628 centers,thereby taking the boundary conditions into account */
1629
1630 /* extrapol_velmain(nface,ielfa,xrlfa,vel,vfa,
1631 ifabou,xbounact,ipnei,nef,&icyclic,c,ifatie,xxn,gradvel,
1632 gradvfa,neifa,rf,area,volume,xle,xxi,xxj,xlet,
1633 coefmpc,nmpc,labmpc,ipompc,nodempc,ifaext,nfaext,nactdoh,
1634 &iitg,&num_cpus,&compressible,xxna,&ncfd);*/
1635
1636 // FORTRAN(writevfa,(vfa,nface,nactdohinv,ielfa));
1637
1638 /* end subiterations */
1639
1640 relnormmax=relnormt;
1641 if(relnormv>relnormmax){relnormmax=relnormv;}
1642 if(relnormp>relnormmax){relnormmax=relnormp;}
1643
1644 // fprintf(f3,"%d %d %11.4e %11.4e %11.4e %11.4e\n",iincf,iit,dtimef,relnormt,relnormv,relnormp);
1645
1646 if(*nmethod==1){
1647
1648 /* steady state flow:
1649 calculate the velocity only once in each increment */
1650
1651 if(relnormmax<1.e-10) iconvergence=1;
1652 if(compressible==0){
1653
1654 /* check whether dtimef was changed in the last
1655 increment */
1656
1657 if(beta>0.){
1658 a1=1.5/dtimef;
1659 a2=-2./dtimef;
1660 a3=0.5/dtimef;
1661 beta=0.;
1662 }
1663 fprintf(f3,"%d %d %11.4e %11.4e %11.4e %11.4e\n",iincf,iit,dtimef,relnormt,relnormv,relnormp);
1664 break;
1665 }else if((relnormmax<1.e-3)||(iit==iitf)){
1666 fprintf(f3,"%d %d %11.4e %11.4e %11.4e %11.4e\n",iincf,iit,dtimef,relnormt,relnormv,relnormp);
1667 break;
1668 }
1669 }else{
1670
1671 /* transient flow:
1672 calculate the velocity repeatedly in each increment */
1673
1674 // if((relnormmax<1.e-8)||(iit==iitf)){
1675 if((relnormmax<1.e-5)||(iit==iitf)){
1676
1677 /* dynamic change of time increment for transient
1678 compressible flow */
1679
1680 if(compressible!=0){
1681
1682 /* compressible flow */
1683
1684 if((iito<3)&&(iit<3)){
1685
1686 /* fast convergence */
1687
1688 dtimef*=1.05;
1689 printf("increased time increment to %e\n",dtimef);
1690 a1=1./dtimef;
1691 a2=-a1;
1692 }else if(iit>iitf-1){
1693
1694 /* divergence */
1695
1696 timef-=dtimef;
1697 memcpy(&vel[0],&velo[0],sizeof(double)*8**nef);
1698 memcpy(&velo[0],&veloo[0],sizeof(double)*8**nef);
1699 dtimef*=0.25;
1700 printf("divergence: recalculated increment with reduced time increment %e\n",dtimef);
1701 a1=1./dtimef;
1702 a2=-a1;
1703 }else if((iito>10)&&(iit>10)){
1704
1705 /* slow convergence */
1706
1707 dtimef*=0.95;
1708 printf("decreased time increment to %e\n",dtimef);
1709 a1=1./dtimef;
1710 a2=-a1;
1711 }
1712 iito=iit;
1713 }
1714 fprintf(f3,"%d %d %11.4e %11.4e %11.4e %11.4e\n",iincf,iit,dtimef,relnormt,relnormv,relnormp);
1715 break;
1716 }
1717 }
1718
1719 }while(1);
1720
1721 if(((iincf/jout[1])*jout[1]==iincf)||(iconvergence==1)||
1722 (iincf==jmax[1])){
1723
1724 /* calculating the stress and the heat flow at the
1725 integration points, if requested */
1726
1727 if((strcmp1(&filab[3306],"SF ")==0)||
1728 (strcmp1(&filab[3480],"SVF ")==0))isti=1;
1729 if(strcmp1(&filab[3393],"HFLF")==0)iqfx=1;
1730 for(i=0;i<*nprint;i++){
1731 if(strcmp1(&prlab[6*i],"SVF")==0) isti=1;
1732 if(strcmp1(&prlab[6*i],"HFLF")==0)iqfx=1;
1733 }
1734
1735 /* calculating the heat conduction at the element centers */
1736
1737 if(iqfx==1){
1738 NNEW(hcel,double,*nef);
1739 FORTRAN(calchcel,(vel,cocon,ncocon,ielmatf,ntmat_,mi,
1740 hcel,nef));
1741 }
1742
1743 /* calculating the stress and/or the heat flux at the
1744 element centers */
1745
1746 if((isti==1)||(iqfx==1)){
1747 FORTRAN(calcstressheatflux,(sti,umel,gradvel,qfx,hcel,
1748 gradtel,nef,&isti,&iqfx,mi));
1749 if(iqfx==1)SFREE(hcel);
1750 }
1751
1752 /* extrapolating the stresses */
1753
1754 if((strcmp1(&filab[3306],"SF ")==0)||
1755 (strcmp1(&filab[3480],"SVF ")==0)){
1756 nfield=6;
1757 ndim=6;
1758 if((*norien>0)&&
1759 ((strcmp1(&filab[3311],"L")==0)||(strcmp1(&filab[3485],"L")==0))){
1760 iorienglob=1;
1761 }else{
1762 iorienglob=0;
1763 }
1764 strcpy1(&cflag[0],&filab[2962],1);
1765 NNEW(stn,double,6**nk);
1766 FORTRAN(extrapolate,(sti,stn,ipkonf,inum,konf,lakonf,
1767 &nfield,nk,nef,mi,&ndim,orab,ielorienf,co,&iorienglob,
1768 cflag,vold,&force,ielmatf,thicke,ielpropf,prop));
1769 }
1770
1771 /* extrapolating the heat flow */
1772
1773 if(strcmp1(&filab[3393],"HFLF")==0){
1774 nfield=3;
1775 ndim=3;
1776 if((*norien>0)&&(strcmp1(&filab[3398],"L")==0)){
1777 iorienglob=1;
1778 }else{
1779 iorienglob=0;
1780 }
1781 strcpy1(&cflag[0],&filab[3049],1);
1782 NNEW(qfn,double,3**nk);
1783 FORTRAN(extrapolate,(qfx,qfn,ipkonf,inum,konf,lakonf,
1784 &nfield,nk,nef,mi,&ndim,orab,ielorienf,co,&iorienglob,
1785 cflag,vold,&force,ielmatf,thicke,ielpropf,prop));
1786 }
1787
1788 /* extrapolating the facial values of the static temperature
1789 and/or the velocity and/or the static pressure to the nodes */
1790
1791 if(imach){NNEW(xmach,double,*nk);}
1792 if(ikappa){NNEW(xkappa,double,*nk);}
1793 if(iturb){NNEW(xturb,double,2**nk);}
1794
1795 FORTRAN(extrapolatefluid,(nk,ipofano,ifano,inum,vfa,vold,ielfa,
1796 ithermal,&imach,&ikappa,xmach,xkappa,shcon,
1797 nshcon,ntmat_,ielmatf,physcon,mi,&iturb,xturb,
1798 gradtfa,gradvfa,gradpfa,gradkfa,gradofa,co,
1799 cofa,ifabou));
1800
1801 /* storing the results in dat-format */
1802
1803 ptimef=ttimef+timef;
1804 FORTRAN(printoutfluid,(set,nset,istartset,iendset,ialset,nprint,
1805 prlab,prset,ipkonf,lakonf,sti,eei,
1806 xstate,ener,mi,nstate_,co,konf,qfx,
1807 &ptimef,trab,inotr,ntrans,orab,ielorienf,
1808 norien,vold,ielmatf,
1809 thicke,eme,xturb,physcon,nactdoh,
1810 ielpropf,prop,xkappa,xmach,ithermal,
1811 orname));
1812
1813 /* thermal flux and drag: storage in dat-format */
1814
1815 FORTRAN(printoutface,(co,rhcon,nrhcon,ntmat_,vold,shcon,nshcon,
1816 cocon,ncocon,&compressible,istartset,iendset,ipkonf,
1817 lakonf,konf,
1818 ialset,prset,&ptimef,nset,set,nprint,prlab,ielmatf,mi,
1819 ithermal,nactdoh,&icfd,time,stn));
1820
1821 /* storing the results in frd-format */
1822
1823 FORTRAN(frdfluid,(co,nk,konf,ipkonf,lakonf,nef,vold,kode,&timef,ielmatf,
1824 matname,filab,inum,ntrans,inotr,trab,mi,istep,
1825 stn,qfn,nactdohinv,xmach,xkappa,physcon,xturb,
1826 coel,vel,cofa,vfa,nface));
1827
1828 // FORTRAN(writevfa,(vfa,nface,nactdohinv,ielfa));
1829
1830 if((strcmp1(&filab[3306],"SF ")==0)||
1831 (strcmp1(&filab[3480],"SVF ")==0)){SFREE(stn);}
1832 if(strcmp1(&filab[3393],"HFLF")==0){SFREE(qfn);}
1833
1834 if(imach){SFREE(xmach);}
1835 if(ikappa){SFREE(xkappa);}
1836 if(iturb){SFREE(xturb);}
1837
1838 }
1839
1840 if(iincf==jmax[1]){
1841 printf("*INFO: maximum number of fluid increments reached\n\n");
1842 fclose(f3);
1843 break;
1844 // FORTRAN(stop,());
1845 }
1846 if(last==1){
1847 printf("*INFO: mechanical time increment reached: time=%e\n\n",*dtime);
1848 fclose(f3);
1849 // FORTRAN(stop,());
1850 break;
1851 }
1852 if(iconvergence==1){
1853 printf("*INFO: steady state reached\n\n");
1854 fclose(f3);
1855 // FORTRAN(stop,());
1856 break;
1857 }
1858
1859 memcpy(&veloo[0],&velo[0],sizeof(double)*8**nef);
1860 memcpy(&velo[0],&vel[0],sizeof(double)*8**nef);
1861
1862 }while(1);
1863
1864 FORTRAN(closefilefluid,());
1865
1866 SFREE(flux);
1867
1868 if(compressible==0){SFREE(ia);SFREE(ja);SFREE(aua);}
1869
1870 SFREE(irow);SFREE(icol);SFREE(jq);SFREE(iau6);
1871
1872 SFREE(coel);SFREE(cosa);SFREE(xxn);SFREE(xxi);SFREE(xle);SFREE(xlen);
1873 SFREE(xlet);SFREE(cofa);SFREE(area);SFREE(xrlfa);SFREE(volume);
1874 SFREE(cosb);SFREE(xxni);SFREE(xxnj);SFREE(xxicn);SFREE(xxj);
1875 SFREE(xxna);SFREE(rf);SFREE(ale);SFREE(alet);SFREE(inlet);SFREE(h);
1876 SFREE(sc);
1877 if(*iturbulent>0){
1878 SFREE(dy);
1879 if(*iturbulent>2) SFREE(yy);
1880 }
1881
1882 SFREE(ifabou);SFREE(umfa);SFREE(umel);
1883
1884 SFREE(gradvel);SFREE(gradvfa);SFREE(au);SFREE(ad);SFREE(b);SFREE(advfa);
1885 SFREE(ap);SFREE(bp);SFREE(gradpel);SFREE(rwork);SFREE(gradpfa);
1886 SFREE(gradpcel);SFREE(gradpcfa);SFREE(bvcp);SFREE(bcp);
1887 SFREE(hfa);SFREE(hel);SFREE(adv);SFREE(bv);SFREE(sel);
1888 SFREE(auv);
1889
1890 if(*ithermal>0){
1891 SFREE(gradtel);SFREE(gradtfa);SFREE(hcfa);SFREE(cvel);SFREE(cvfa);
1892 }
1893
1894 if(*iturbulent>0){
1895 SFREE(gradkel);SFREE(gradkfa);SFREE(gradoel);SFREE(gradofa);
1896 if(*iturbulent>2){
1897 SFREE(f1);
1898 if(*iturbulent>3) SFREE(of2);
1899 }
1900 }
1901
1902 SFREE(inum);
1903
1904 SFREE(ipofano);SFREE(ifano);
1905
1906 if(*nbody>0) SFREE(body);
1907
1908 *ithermal=ithermalref;
1909
1910 SFREE(temp);SFREE(gamma);
1911
1912 SFREE(iam);SFREE(jam);SFREE(iamorig);SFREE(am);
1913 SFREE(nestart);SFREE(ineighblock);SFREE(neighblock);
1914
1915 SFREE(iponofa);SFREE(inofa);
1916
1917 *cop=co;*ipkonfp=ipkonf;*lakonfp=lakonf;*voldp=vold;
1918
1919 return;
1920
1921 }
1922