root/src/compfluid.c

/* [<][>][^][v][top][bottom][index][help] */

DEFINITIONS

This source file includes following definitions.
  1. compfluid

   1 /*     CalculiX - A 3-dimensional finite element program                 */
   2 /*              Copyright (C) 1998-2015 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 
  39 void compfluid(double **cop, ITG *nk, ITG **ipkonfp, ITG **konp, char **lakonfp,
  40     char **sidefacep, ITG *ifreestream, 
  41     ITG *nfreestream, ITG *isolidsurf, ITG *neighsolidsurf,
  42     ITG *nsolidsurf, ITG **iponoelp, ITG **inoelp, ITG *nshcon, double *shcon,
  43     ITG *nrhcon, double *rhcon, double **voldp, ITG *ntmat_,ITG *nodeboun, 
  44     ITG *ndirboun, ITG *nboun, ITG **ipompcp,ITG **nodempcp, ITG *nmpc,
  45     ITG **ikmpcp, ITG **ilmpcp, ITG *ithermal, ITG *ikboun, ITG *ilboun,
  46     ITG *iturbulent, ITG *isolver, ITG *iexpl, double *vcontu, double *ttime,
  47     double *time, double *dtime, ITG *nodeforc,ITG *ndirforc,double *xforc,
  48     ITG *nforc, ITG *nelemload, char *sideload, double *xload,ITG *nload,
  49     double *xbody,ITG *ipobody,ITG *nbody, ITG *ielmatf, char *matname,
  50     ITG *mi, ITG *ncmat_, double *physcon, ITG *istep, ITG *iinc,
  51     ITG *ibody, double *xloadold, double *xboun,
  52     double **coefmpcp, 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, double *xbounact,
  56     ITG *iamboun, ITG *itg, ITG *ntg, char *amname, double *t0, 
  57     ITG **nelemfacep,
  58     ITG *nface, double *cocon, ITG *ncocon, double *xloadact, double *tper,
  59     ITG *jmax, ITG *jout, char *set, ITG *nset, ITG *istartset,
  60     ITG *iendset, ITG *ialset, char *prset, char *prlab, ITG *nprint,
  61     double *trab, ITG *inotr, ITG *ntrans, char *filab, char **labmpcp, 
  62     double *sti, ITG *norien, double *orab, char *jobnamef,char *tieset,
  63     ITG *ntie, ITG *mcs, ITG *ics, double *cs, ITG *nkon, ITG *mpcfree,
  64     ITG *memmpc_,double **fmpcp,ITG *nef,ITG **inomatp,double *qfx,
  65     ITG *neifa,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 
  69     /* main computational fluid dynamics routine */
  70   
  71   char cflag[1],*labmpc=NULL,*lakonf=NULL,*sideface=NULL;
  72 
  73   char matvec[7]="MATVEC",msolve[7]="MSOLVE";
  74 
  75   ITG *ipointer=NULL,*mast1=NULL,*irow=NULL,*icol=NULL,*jq=NULL,
  76       nzs=20000000,neq,kode,compressible,*ifabou=NULL,*ja=NULL,
  77       *nodempc=NULL,*ipompc=NULL,*ikmpc=NULL,*ilmpc=NULL,nfabou,im,
  78       *ipkonf=NULL,*kon=NULL,*nelemface=NULL,*inoel=NULL,last=0,
  79       *iponoel=NULL,*inomat=NULL,ithermalref,*integerglob=NULL,iit,
  80       iconvergence=0,symmetryflag,inputformat,i,*inum=NULL,iitf,ifreefa,
  81       *iponofa=NULL,*inofa=NULL,is,ie,*ia=NULL,nstate_,*ielpropf=NULL,
  82       icent=0,isti=0,iqfx=0,nfield,ndim,iorienglob,force=0,icfdout=1;
  83 
  84   ITG nelt,isym,itol,itmax,iunit,lrgw,*igwk=NULL,ligw,ierr,*iwork=NULL,iter,
  85       nsave,lenw,leniw;
  86 
  87   double *coefmpc=NULL,*fmpc=NULL,*umfa=NULL,reltime,*doubleglob=NULL,
  88       *co=NULL,*vold=NULL,*coel=NULL,*cosa=NULL,*gradvel=NULL,*gradvfa=NULL,
  89       *xxn=NULL,*xxi=NULL,*xle=NULL,*xlen=NULL,*xlet=NULL,timef,dtimef,
  90       *cofa=NULL,*area=NULL,*xrlfa=NULL,reltimef,ttimef,*hcfa=NULL,*cpel=NULL,
  91       *au=NULL,*ad=NULL,*b=NULL,*volume=NULL,*body=NULL,sigma=0.,betam,
  92       *adb=NULL,*aub=NULL,*advfa=NULL,*ap=NULL,*bp=NULL,*xxj=NULL,
  93       *v=NULL,*velo=NULL,*veloo=NULL,*gammat=NULL,*cosb=NULL,dmin,tincfguess,
  94       *hel=NULL,*hfa=NULL,*auv=NULL,*adv=NULL,*bv=NULL,*sel=NULL,*gamma=NULL,
  95       *gradtfa=NULL,*gradtel=NULL,*umel=NULL,*cpfa=NULL,*gradpel=NULL,
  96       *fn=NULL,*eei=NULL,*xstate=NULL,*ener=NULL,*thicke=NULL,*eme=NULL,
  97       ptimef,*stn=NULL,*qfn=NULL,*hcel=NULL,*aua=NULL,a1,a2,a3,beta,
  98       *prop=NULL;
  99 
 100   double tol,*rgwk=NULL,err,*sb=NULL,*sx=NULL,*rwork=NULL;
 101 
 102   nodempc=*nodempcp;ipompc=*ipompcp;ikmpc=*ikmpcp;ilmpc=*ilmpcp;
 103   coefmpc=*coefmpcp;labmpc=*labmpcp;fmpc=*fmpcp;co=*cop;
 104   ipkonf=*ipkonfp;lakonf=*lakonfp;kon=*konp;
 105   nelemface=*nelemfacep;sideface=*sidefacep;inoel=*inoelp;
 106   iponoel=*iponoelp;vold=*voldp;inomat=*inomatp;
 107 
 108 #ifdef SGI
 109   ITG token;
 110 #endif
 111 
 112   /* relative time at the end of the mechanical increment */
 113 
 114   reltime=(*time)/(*tper);
 115 
 116   /* open frd-file for fluids */
 117 
 118   FORTRAN(openfilefluid,(jobnamef));
 119 
 120   /* variables for multithreading procedure */
 121 
 122   ITG sys_cpus;
 123   char *env,*envloc,*envsys;
 124       
 125   num_cpus = 0;
 126   sys_cpus=0;
 127   
 128   /* explicit user declaration prevails */
 129   
 130   envsys=getenv("NUMBER_OF_CPUS");
 131   if(envsys){
 132       sys_cpus=atoi(envsys);
 133       if(sys_cpus<0) sys_cpus=0;
 134   }
 135   
 136   /* automatic detection of available number of processors */
 137   
 138   if(sys_cpus==0){
 139       sys_cpus = getSystemCPUs();
 140       if(sys_cpus<1) sys_cpus=1;
 141   }
 142   
 143   /* local declaration prevails, if strictly positive */
 144   
 145   envloc = getenv("CCX_NPROC_CFD");
 146   if(envloc){
 147       num_cpus=atoi(envloc);
 148       if(num_cpus<0){
 149           num_cpus=0;
 150       }else if(num_cpus>sys_cpus){
 151           num_cpus=sys_cpus;
 152       }
 153   }
 154   
 155   /* else global declaration, if any, applies */
 156   
 157   env = getenv("OMP_NUM_THREADS");
 158   if(num_cpus==0){
 159       if (env)
 160           num_cpus = atoi(env);
 161       if (num_cpus < 1) {
 162           num_cpus=1;
 163       }else if(num_cpus>sys_cpus){
 164           num_cpus=sys_cpus;
 165       }
 166   }
 167   
 168 // next line is to be inserted in a similar way for all other paralell parts
 169   
 170   if(*nef<num_cpus) num_cpus=*nef;
 171   
 172   printf(" Using up to %" ITGFORMAT " cpu(s) for CFD.\n", num_cpus);
 173   
 174   pthread_t tid[num_cpus];
 175 
 176   
 177   kode=0;
 178   
 179   /*  *iexpl==0:  structure:implicit, fluid:incompressible
 180       *iexpl==1:  structure:implicit, fluid:compressible
 181       *iexpl==2:  structure:explicit, fluid:incompressible
 182       *iexpl==3:  structure:explicit, fluid:compressible */
 183 
 184   if((*iexpl==1)||(*iexpl==3)){
 185       compressible=1;
 186   }else{
 187       compressible=0;
 188   }
 189 
 190   /* if initial conditions are specified for the temperature, 
 191      it is assumed that the temperature is an unknown */
 192 
 193   ithermalref=*ithermal;
 194   if(*ithermal==1){
 195     *ithermal=2;
 196   }
 197 
 198   /* determining the matrix structure */
 199   
 200   NNEW(ipointer,ITG,3**nk);
 201   NNEW(mast1,ITG,nzs);
 202   NNEW(irow,ITG,nzs);
 203   NNEW(ia,ITG,nzs);
 204   NNEW(icol,ITG,*nef);
 205   NNEW(jq,ITG,*nef+1);
 206   NNEW(ja,ITG,*nef+1);
 207 //  NNEW(nactdoh,ITG,*nef);
 208 
 209   mastructf(nk,kon,ipkonf,lakonf,nef,icol,jq,&mast1,&irow,
 210             isolver,&neq,ipointer,&nzs,ipnei,neiel,mi);
 211 
 212 //  printf("Unterschied start\n");
 213 //  for(i=0;i<*ne;i++){if(i+1!=nactdoh[i]){printf("Unterschied i=%d,nactdoh[i]=%d\n",i+1,nactdoh[i]);}}
 214 //  printf("Unterschied end\n");
 215 
 216   SFREE(ipointer);SFREE(mast1);
 217  
 218   /* calculation geometric data */
 219 
 220   NNEW(coel,double,3**nef);
 221   NNEW(volume,double,*nef);
 222   NNEW(cosa,double,*nflnei);
 223   NNEW(cosb,double,*nflnei);
 224   NNEW(xxn,double,3**nflnei);
 225   NNEW(xxi,double,3**nflnei);
 226   NNEW(xxj,double,3**nflnei);
 227   NNEW(xle,double,*nflnei);
 228   NNEW(xlen,double,*nflnei);
 229   NNEW(xlet,double,*nflnei);
 230   NNEW(cofa,double,3**nface);
 231   NNEW(area,double,*nface);
 232   NNEW(xrlfa,double,3**nface);
 233 
 234   FORTRAN(initialcfd,(nef,ipkonf,kon,lakonf,co,coel,cofa,nface,
 235           ielfa,area,ipnei,neiel,xxn,xxi,xle,xlen,xlet,xrlfa,cosa,
 236           volume,neifa,xxj,cosb,vel,&dmin));
 237 
 238   /* storing pointers to the boundary conditions in ielfa */
 239 
 240   NNEW(ifabou,ITG,7**nfaext);
 241   FORTRAN(applyboun,(ifaext,nfaext,ielfa,ikboun,ilboun,
 242        nboun,typeboun,nelemload,nload,sideload,isolidsurf,nsolidsurf,
 243        ifabou,&nfabou,nface,nodeboun,ndirboun,ikmpc,ilmpc,labmpc,nmpc,
 244        nactdohinv));
 245   RENEW(ifabou,ITG,nfabou);
 246 
 247   /* catalogueing the nodes for output purposes (interpolation at
 248      the nodes */
 249   
 250   NNEW(iponofa,ITG,*nk);
 251   NNEW(inofa,ITG,2**nface*4);
 252 
 253   FORTRAN(cataloguenodes,(iponofa,inofa,&ifreefa,ielfa,ifabou,ipkonf,
 254                           kon,lakonf,nface,nk));
 255 
 256   RENEW(inofa,ITG,2*ifreefa);
 257 
 258   /* material properties for athermal calculations 
 259      = calculation for which no initial thermal conditions
 260      were defined */
 261 
 262   NNEW(umfa,double,*nface);
 263   NNEW(umel,double,*nef);
 264       
 265   /* calculating the density at the element centers */
 266   
 267   FORTRAN(calcrhoel,(nef,vel,rhcon,nrhcon,ielmatf,ntmat_,
 268                      ithermal,mi));
 269   
 270   /* calculating the density at the face centers */
 271   
 272   FORTRAN(calcrhofa,(nface,vfa,rhcon,nrhcon,ielmatf,ntmat_,
 273                      ithermal,mi,ielfa));
 274   
 275   /* calculating the dynamic viscosity at the face centers */
 276   
 277   FORTRAN(calcumfa,(nface,vfa,shcon,nshcon,ielmatf,ntmat_,
 278                     ithermal,mi,ielfa,umfa));
 279   
 280   /* calculating the dynamic viscosity at the element centers */
 281   
 282   FORTRAN(calcumel,(nef,vel,shcon,nshcon,ielmatf,ntmat_,
 283                     ithermal,mi,umel));
 284   
 285   if(*ithermal!=0){
 286       NNEW(hcfa,double,*nface);
 287       NNEW(cpel,double,*nef);
 288       NNEW(cpfa,double,*nface);
 289   }
 290       
 291 
 292   if(*nbody>0) NNEW(body,double,4**nef);
 293 
 294   /* v is a auxiliary field: set to zero for the calls to
 295      tempload */
 296 
 297   NNEW(v,double,5**nk);
 298 
 299   /* next section is for stationary calculations */
 300   
 301   if(*nmethod==1){
 302       
 303       /* boundary conditions at the end of the mechanical
 304          increment */
 305       
 306       FORTRAN(tempload,(xforcold,xforc,xforcact,iamforc,nforc,
 307              xloadold,xload,xloadact,iamload,nload,ibody,xbody,nbody,
 308              xbodyold,xbodyact,t1old,t1,t1act,iamt1,nk,amta,
 309              namta,nam,ampli,time,&reltime,ttime,dtime,ithermal,nmethod,
 310              xbounold,xboun,xbounact,iamboun,nboun,
 311              nodeboun,ndirboun,nodeforc,ndirforc,istep,iinc,
 312              co,v,itg,ntg,amname,ikboun,ilboun,nelemload,sideload,mi,
 313              ntrans,trab,inotr,vold,integerglob,doubleglob,tieset,istartset,
 314              iendset,ialset,ntie,nmpc,ipompc,ikmpc,ilmpc,nodempc,coefmpc));
 315 
 316       /* body forces (gravity, centrifugal and Coriolis forces */
 317 
 318       if(*nbody>0){
 319           FORTRAN(inicalcbody,(nef,body,ipobody,ibody,xbody,coel,vel,lakonf,
 320                             nactdohinv,&icent));
 321       }
 322   }
 323 
 324   /* extrapolating the velocity from the elements centers to the face
 325      centers, thereby taking the boundary conditions into account */
 326 
 327   FORTRAN(extrapolate_vel,(nface,ielfa,xrlfa,vel,vfa,
 328                            ifabou,xbounact,ipnei,nef));
 329 
 330   /* applying MPC's to the faces */
 331 
 332   if(*nmpc>0){
 333       is=1;ie=3;
 334       FORTRAN(applympc,(nface,ielfa,&is,&ie,ifabou,ipompc,vfa,coefmpc,
 335                         nodempc,ipnei,neifa,labmpc,xbounact,nactdoh));
 336   }
 337 
 338   /* extrapolation of the pressure at the element centers
 339      to the face centers */
 340   
 341   FORTRAN(extrapolate_pel,(nface,ielfa,xrlfa,vel,vfa,ifabou,
 342                            xbounact,nef));
 343 
 344   /* applying MPC's to the faces */
 345 
 346   if(*nmpc>0){
 347       is=4;ie=4;
 348       FORTRAN(applympc,(nface,ielfa,&is,&ie,ifabou,ipompc,vfa,coefmpc,
 349                         nodempc,ipnei,neifa,labmpc,xbounact,nactdoh));
 350   }
 351 
 352   /* extrapolation of the temperature at the element centers
 353      to the face centers */
 354 
 355   if(*ithermal>0){
 356       FORTRAN(extrapolate_tel,(nface,ielfa,xrlfa,vel,vfa,
 357                                ifabou,xbounact,ipnei,nef));
 358 
 359       /* applying MPC's to the faces */
 360       
 361       if(*nmpc>0){
 362           is=0;ie=0;
 363           FORTRAN(applympc,(nface,ielfa,&is,&ie,ifabou,ipompc,vfa,coefmpc,
 364                             nodempc,ipnei,neifa,labmpc,xbounact,nactdoh));
 365       }
 366   }
 367 
 368   /* calculating the maximum velocity (for the determination of the
 369      time increment) */
 370 
 371   FORTRAN(calcguesstincf,(nface,&dmin,vfa,umfa,&tincfguess));
 372 
 373   /* start of the major loop */
 374 
 375   NNEW(gradvel,double,9**nef);
 376   NNEW(gradvfa,double,9**nface);
 377 
 378   NNEW(advfa,double,*nface);
 379   NNEW(hfa,double,3**nface);
 380 
 381   NNEW(ap,double,*nface);
 382   NNEW(bp,double,*nface);
 383 
 384   NNEW(au,double,2*nzs+neq);
 385   NNEW(aua,double,nzs+neq);
 386 //  NNEW(au,double,2*nzs);
 387   NNEW(ad,double,neq);
 388   NNEW(b,double,neq);
 389 
 390   NNEW(auv,double,2*nzs+neq);
 391 //  NNEW(auv,double,2*nzs);
 392   NNEW(adv,double,neq);
 393   NNEW(bv,double,3*neq);
 394   NNEW(hel,double,3*neq);
 395   NNEW(sel,double,3*neq);
 396 
 397   NNEW(rwork,double,neq);
 398 
 399   NNEW(gradpel,double,3**nef);
 400 
 401   if(*ithermal>0){
 402       NNEW(gradtel,double,3**nef);
 403       NNEW(gradtfa,double,3**nface);
 404   }
 405 
 406   NNEW(inum,ITG,*nk);
 407 
 408   NNEW(velo,double,6**nef);
 409   NNEW(veloo,double,6**nef);
 410 
 411   /* initializing velo and veloo */
 412 
 413   memcpy(&veloo[0],&vel[0],sizeof(double)*6**nef);
 414   memcpy(&velo[0],&vel[0],sizeof(double)*6**nef);
 415 
 416   iit=0;
 417 
 418   a1=1.5;
 419   a2=-2.;
 420   a3=0.5;
 421 
 422   if(*tincf<=0.) *tincf=tincfguess;
 423   printf("time increment for the CFD-calculations = %e\n\n",*tincf);
 424 
 425   ttimef=*ttime;
 426   timef=*time-*dtime;
 427   dtimef=*tincf;
 428 
 429   do{
 430 
 431       iit++;
 432 
 433       printf("iteration = %d\n",iit);
 434 
 435       timef+=dtimef;
 436       if((*time<timef)&&(*nmethod==4)){
 437           dtimef-=(timef-*time);
 438           timef=*time;
 439           last=1;
 440           beta=dtimef/(*tincf);
 441           a1=(2.+beta)/(1.+beta);
 442           a2=-(1.+beta)/beta;
 443           a3=1./(beta*(1.+beta));
 444       }
 445 
 446       /* conditions for transient calculations */
 447      
 448       if(*nmethod==4){
 449 
 450           /* boundary conditions at end of fluid increment */
 451 
 452           FORTRAN(tempload,(xforcold,xforc,xforcact,iamforc,nforc,
 453              xloadold,xload,xloadact,iamload,nload,ibody,xbody,nbody,
 454              xbodyold,xbodyact,t1old,t1,t1act,iamt1,nk,amta,
 455              namta,nam,ampli,&timef,&reltimef,&ttimef,&dtimef,ithermal,nmethod,
 456              xbounold,xboun,xbounact,iamboun,nboun,
 457              nodeboun,ndirboun,nodeforc,ndirforc,istep,iinc,
 458              co,v,itg,ntg,amname,ikboun,ilboun,nelemload,sideload,mi,
 459              ntrans,trab,inotr,vold,integerglob,doubleglob,tieset,istartset,
 460              iendset,ialset,ntie,nmpc,ipompc,ikmpc,ilmpc,nodempc,coefmpc));
 461 
 462           /* body forces (gravity, centrifugal and Coriolis forces) */
 463       
 464           if(*nbody>0){
 465               FORTRAN(calcbody,(nef,body,ipobody,ibody,xbody,coel,vel,lakonf,
 466                                 nactdohinv));
 467           }
 468 
 469       }else if(icent==1){
 470           
 471           /* body forces (gravity, centrifugal and Coriolis forces;
 472              only if centrifugal forces are active => the ensuing
 473              Coriolis forces depend on the actual velocity) */
 474           
 475           FORTRAN(calcbody,(nef,body,ipobody,ibody,xbody,coel,vel,lakonf,
 476                                 nactdohinv));
 477       }
 478 
 479       /* updating of the material properties */
 480 
 481       if(*ithermal>0){
 482           
 483           /* calculating the density at the element centers */
 484           
 485           FORTRAN(calcrhoel,(nef,vel,rhcon,nrhcon,ielmatf,ntmat_,
 486                              ithermal,mi));
 487           
 488           /* calculating the density at the face centers */
 489           
 490           FORTRAN(calcrhofa,(nface,vfa,rhcon,nrhcon,ielmatf,ntmat_,
 491                              ithermal,mi,ielfa));
 492           
 493           /* calculating the dynamic viscosity at the face centers */
 494           
 495           FORTRAN(calcumfa,(nface,vfa,shcon,nshcon,ielmatf,ntmat_,
 496                              ithermal,mi,ielfa,umfa));
 497 
 498           /* calculating the dynamic viscosity at the element centers */
 499           
 500           FORTRAN(calcumel,(nef,vel,shcon,nshcon,ielmatf,ntmat_,
 501                             ithermal,mi,umel));
 502           
 503           /* calculating the heat conduction at the face centers */
 504           
 505           FORTRAN(calchcfa,(nface,vfa,cocon,ncocon,ielmatf,ntmat_,
 506                             mi,ielfa,hcfa));
 507 
 508           /* calculating the specific heat at constant pressure at the 
 509              element centers (secant value) */
 510           
 511           FORTRAN(calccpel,(nef,vel,shcon,nshcon,ielmatf,ntmat_,
 512                             mi,cpel,physcon));
 513 
 514           /* calculating the specific heat at constant pressure at the 
 515              face centers (secant value) */
 516           
 517           FORTRAN(calccpfa,(nface,vfa,shcon,nshcon,ielmatf,ntmat_,
 518                             mi,ielfa,cpfa,physcon));
 519 
 520       }
 521           
 522       /* calculating the gradient of the velocity at the element
 523          centers */
 524       
 525       DMEMSET(gradvel,0,9**nef,0.);
 526       FORTRAN(calcgradvel,(nef,lakonf,ipnei,vfa,area,xxn,gradvel,neifa,
 527                          volume));
 528       
 529       /* extrapolating the gradient of the velocity from the element
 530          centers to the face centers */
 531       
 532       FORTRAN(extrapolate_gradvel,(nface,ielfa,xrlfa,gradvel,gradvfa));
 533       
 534       /* calculate gamma (Ph.D. Thesis Jasak) */
 535 
 536       /*       betam=0.1;
 537       NNEW(gamma,double,*nface);
 538       FORTRAN(calcgamma,(nface,ielfa,vel,gradvfa,gamma,xlet,xxn,xxj,
 539       ipnei,&betam,nef));*/
 540 
 541       /* filling the lhs and rhs's for the balance of momentum
 542          equations */
 543       
 544       DMEMSET(adv,0,neq,0.);
 545       DMEMSET(auv,0,2*nzs,0.);
 546       DMEMSET(bv,0,3*neq,0.);
 547       
 548       FORTRAN(mafillv,(nef,ipnei,neifa,neiel,vfa,xxn,area,
 549                        auv,adv,jq,irow,&nzs,bv,vel,cosa,umfa,xlet,xle,gradvfa,xxi,
 550                        body,volume,&compressible,ielfa,lakonf,ifabou,nbody,&neq,
 551                        &dtimef,velo,veloo,sel,xrlfa,gamma,xxj,nactdohinv,&a1,
 552                        &a2,&a3));
 553 //      SFREE(gamma);
 554       
 555       /* LU decomposition (asymmetric system) */
 556       
 557 /*      inputformat=1;
 558       symmetryflag=2;
 559       
 560       if(*isolver==0){
 561 #ifdef SPOOLES
 562           spooles_factor(adv,auv,adb,aub,&sigma,icol,irow,&neq,&nzs,
 563                          &symmetryflag,&inputformat,&nzs);
 564 #else
 565           printf("*ERROR in compfluid: the SPOOLES library is not linked\n\n");
 566           FORTRAN(stop,());
 567 #endif
 568       }
 569       else if(*isolver==7){
 570 #ifdef PARDISO
 571           pardiso_factor(adv,auv,adb,aub,&sigma,icol,irow,&neq,&nzs,
 572                          &symmetryflag,&inputformat,jq,&nzs);
 573 #else
 574           printf("*ERROR in compfluid: the PARDISO library is not linked\n\n");
 575           FORTRAN(stop,());
 576 #endif
 577 }*/
 578       
 579       /* solving the system of equations (for x, y and z
 580          separately */
 581       
 582 /*      for(i=0;i<3;i++){
 583           
 584           if(*isolver==0){
 585 #ifdef SPOOLES
 586               spooles_solve(&bv[i*neq],&neq);
 587 #endif
 588           }
 589           else if(*isolver==7){
 590 #ifdef PARDISO
 591               pardiso_solve(&bv[i*neq],&neq,&symmetryflag);
 592 #endif
 593           }
 594           }*/
 595       
 596       /* free memory */
 597       
 598 /*      if(*isolver==0){
 599 #ifdef SPOOLES
 600           spooles_cleanup();
 601 #endif
 602       }
 603       else if(*isolver==7){
 604 #ifdef PARDISO
 605           pardiso_cleanup(&neq,&symmetryflag);
 606 #endif
 607 }*/
 608 
 609       memcpy(&auv[2*nzs],adv,sizeof(double)*neq);
 610       nelt=2*nzs+neq;
 611       isym=0;
 612       itol=0;
 613       tol=0.;
 614       itmax=0;
 615       iunit=0;
 616       lrgw=131+16*neq;
 617       NNEW(rgwk,double,lrgw);
 618       NNEW(igwk,ITG,20);
 619       igwk[0]=10;
 620       igwk[1]=10;
 621       igwk[2]=0;
 622       igwk[3]=1;
 623       igwk[4]=10;
 624       ligw=20;
 625       for(i=0;i<neq;i++){rwork[i]=1./adv[i];}
 626       FORTRAN(predgmres,(&neq,&bv[0],&vel[neq],&nelt,irow,jq,auv,
 627                       &isym,&itol,&tol,&itmax,&iter,
 628                       &err,&ierr,&iunit,sb,sx,rgwk,&lrgw,igwk,
 629                       &ligw,rwork,iwork));
 630 //            printf(" err1=%e\n",err);
 631 //      memcpy(&bv[0],&vel[neq],sizeof(double)*neq);
 632       FORTRAN(predgmres,(&neq,&bv[neq],&vel[2*neq],&nelt,irow,jq,auv,
 633                       &isym,&itol,&tol,&itmax,&iter,
 634                       &err,&ierr,&iunit,sb,sx,rgwk,&lrgw,igwk,
 635                       &ligw,rwork,iwork));
 636 //            printf(" err2=%e\n",err);
 637 //      memcpy(&bv[neq],&vel[2*neq],sizeof(double)*neq);
 638       FORTRAN(predgmres,(&neq,&bv[2*neq],&vel[3*neq],&nelt,irow,jq,auv,
 639                       &isym,&itol,&tol,&itmax,&iter,
 640                       &err,&ierr,&iunit,sb,sx,rgwk,&lrgw,igwk,
 641                       &ligw,rwork,iwork));
 642 //            printf(" err3=%e\n",err);
 643 //      memcpy(&bv[2*neq],&vel[3*neq],sizeof(double)*neq);
 644       SFREE(rgwk);SFREE(igwk);
 645       
 646       /* storing the solution into vel */
 647       
 648 //      FORTRAN(calcvel,(ne,nactdoh,vel,bv,&neq,ne));
 649 
 650       /* calculating the pressure gradient at the element
 651          centers */
 652 
 653       DMEMSET(gradpel,0,3**nef,0.);
 654       FORTRAN(calcgradpel,(nef,lakonf,ipnei,vfa,area,xxn,gradpel,neifa,
 655                                volume));
 656       
 657       for(iitf=0;iitf<2;iitf++){
 658 
 659           memcpy(&hel[0],&sel[0],sizeof(double)*(3*neq));
 660           
 661           /* completing hel with the neighboring velocity contributions */
 662           
 663           FORTRAN(complete_hel,(&neq,&vel[neq],hel,adv,auv,jq,irow,&nzs));
 664 //          FORTRAN(complete_hel,(&neq,bv,hel,adv,auv,jq,irow,&nzs));
 665           
 666           /* generating ad and h at the face centers (advfa and hfa) */
 667           
 668           FORTRAN(extrapolate_ad_h,(nface,ielfa,xrlfa,adv,advfa,hel,hfa));
 669           
 670           /* calculating the lhs and rhs of the equation system to determine
 671              p (balance of mass) */
 672           
 673           DMEMSET(b,0,neq,0.);
 674 
 675           if(iitf==0){
 676 
 677               /* first iteration: calculating both lhs and rhs */
 678 
 679               DMEMSET(ad,0,neq,0.);
 680               DMEMSET(au,0,nzs,0.);
 681               
 682               FORTRAN(mafillp,(nef,lakonf,ipnei,neifa,neiel,vfa,area,
 683                                advfa,xlet,cosa,volume,au,ad,jq,irow,ap,
 684                                ielfa,ifabou,xle,b,xxn,&compressible,&neq,
 685                                &nzs,hfa,gradpel,bp,xxi,neij,xlen,cosb));
 686 
 687               FORTRAN(convert2slapcol,(au,ad,irow,ia,jq,ja,&nzs,&neq,aua));
 688               
 689               /* LU decomposition of the p system (symmetric system) */
 690               
 691 /*            inputformat=0;
 692               symmetryflag=0;
 693               
 694               if(*isolver==0){
 695 #ifdef SPOOLES
 696                   spooles_factor(ad,au,adb,aub,&sigma,icol,irow,&neq,&nzs,
 697                                  &symmetryflag,&inputformat,&nzs);
 698 #else
 699                   printf("*ERROR in compfluid: the SPOOLES library is not linked\n\n");
 700                   FORTRAN(stop,());
 701 #endif
 702               }
 703               else if(*isolver==7){
 704 #ifdef PARDISO
 705                   pardiso_factor(ad,au,adb,aub,&sigma,icol,irow,&neq,&nzs,
 706                                  &symmetryflag,&inputformat,jq,&nzs);
 707 #else
 708                   printf("*ERROR in compfluid: the PARDISO library is not linked\n\n");
 709                   FORTRAN(stop,());
 710 #endif
 711 }*/
 712           }else{
 713 
 714               /* calculating the pressure gradient at the element
 715                  centers */
 716 
 717               DMEMSET(gradpel,0,3**nef,0.);
 718               FORTRAN(calcgradpel,(nef,lakonf,ipnei,vfa,area,xxn,gradpel,neifa,
 719               volume));
 720 
 721               /* second, third.. iteration: calculate the rhs only */
 722               
 723               FORTRAN(rhsp,(nef,lakonf,nactdoh,ipnei,neifa,neiel,vfa,area,
 724                       advfa,xlet,cosa,volume,au,ad,jq,irow,ap,ielfa,ifabou,xle,
 725                       b,xxn,&compressible,&neq,&nzs,hfa,bp,neij,xxi,
 726                       gradpel,xlen));
 727 
 728           }
 729           
 730           /* solving the system of equations for p  */
 731           
 732 /*        if(*isolver==0){
 733 #ifdef SPOOLES
 734               spooles_solve(b,&neq);
 735 #endif
 736           }
 737           else if(*isolver==7){
 738 #ifdef PARDISO
 739               pardiso_solve(b,&neq,&symmetryflag);
 740 #endif
 741 }*/
 742 
 743           /* dslugm; attention: convert2slapcol changes au! */
 744 
 745 //        FORTRAN(convert2slapcol,(au,ad,irow,ia,jq,ja,&nzs,&neq,aua));
 746 
 747           nelt=nzs+neq;
 748           isym=1;
 749           nsave=10;
 750           itol=0;
 751           tol=0.;
 752           itmax=110;
 753           iunit=0;
 754           lenw=131+17*neq+2*nelt;
 755           NNEW(rgwk,double,lenw);
 756           leniw=32+4*neq+2*nelt;
 757           NNEW(igwk,ITG,leniw);
 758           
 759           FORTRAN(dslugm,(&neq,&b[0],&vel[4*neq],&nelt,ia,ja,aua,
 760                           &isym,&nsave,&itol,&tol,&itmax,&iter,
 761                           &err,&ierr,&iunit,rgwk,&lenw,igwk,&leniw));
 762           SFREE(rgwk);SFREE(igwk);
 763 
 764           if(ierr!=0){
 765               printf("       error message from dslugm: ierr = %d\n",ierr);
 766               printf(" err=%e\n",err);
 767               }
 768           
 769           /* storing the solution p into vel(4,*) */
 770           
 771 //        FORTRAN(calcpel,(nef,nactdoh,vel,b,nef));
 772           
 773           /* extrapolation of the pressure at the element centers
 774              to the face centers */
 775           
 776           FORTRAN(extrapolate_pel,(nface,ielfa,xrlfa,vel,vfa,ifabou,
 777                                    xbounact,nef));
 778 
 779           /* applying MPC's to the faces */
 780           
 781           if(*nmpc>0){
 782               is=4;ie=4;
 783               FORTRAN(applympc,(nface,ielfa,&is,&ie,ifabou,ipompc,vfa,coefmpc,
 784                                 nodempc,ipnei,neifa,labmpc,xbounact,nactdoh));
 785           }
 786           
 787           /* correction of the velocity at the element centers due
 788              to the pressure change */
 789 
 790           FORTRAN(correctvel,(hel,adv,vfa,ipnei,area,&vel[neq],xxn,neifa,
 791                               lakonf,nef,&neq));
 792 
 793       }
 794       
 795       /* free memory */
 796       
 797       /*   if(*isolver==0){
 798 #ifdef SPOOLES
 799           spooles_cleanup();
 800 #endif
 801       }
 802       else if(*isolver==7){
 803 #ifdef PARDISO
 804           pardiso_cleanup(&neq,&symmetryflag);
 805 #endif
 806 }*/
 807       
 808       if(*ithermal>0){
 809 
 810           /* adding the velocity correction at the face centers
 811              due to the balance of mass =>
 812              the resulting mass flux is correct,
 813              the face velocity vectors do not have to be correct
 814              only needed if extra equations (temperature,
 815              turbulence have to be solved)  */
 816 
 817           FORTRAN(correctvfa,(nface,ielfa,area,vfa,ap,bp,xxn,
 818                               ifabou,ipnei,nef,neifa,hfa,vel,xbounact));
 819 
 820           /* calculating the temperature gradient at the element
 821              centers */
 822 
 823           DMEMSET(gradtel,0,3**nef,0.);
 824           FORTRAN(calcgradtel,(nef,lakonf,ipnei,vfa,area,xxn,gradtel,neifa,
 825                                volume));
 826 
 827           /* extrapolating the temperature gradient from the element
 828              centers to the face centers */
 829 
 830           FORTRAN(extrapolate_gradtel,(nface,ielfa,xrlfa,gradtel,gradtfa));
 831       
 832       /* calculate gammat (Ph.D. Thesis Jasak) */
 833 
 834 //        betam=0.1;
 835 //        NNEW(gammat,double,*nface);
 836 //        FORTRAN(calcgammat,(nface,ielfa,vel,gradtfa,gammat,xlet,xxn,xxj,
 837 //                            ipnei,&betam,ne));
 838 
 839           /* calculating the lhs and rhs of the energy equation */
 840 
 841           DMEMSET(ad,0,neq,0.);
 842           DMEMSET(au,0,2*nzs,0.);
 843           DMEMSET(b,0,neq,0.);
 844 
 845           FORTRAN(mafillt,(nef,ipnei,neifa,neiel,vfa,xxn,area,
 846                au,ad,jq,irow,&nzs,b,vel,umel,xlet,xle,gradtfa,xxi,
 847                body,volume,&compressible,ielfa,lakonf,ifabou,nbody,&neq,
 848                &dtimef,velo,veloo,cpfa,hcfa,cpel,gradvel,xload,gammat,
 849                xrlfa,xxj,nactdohinv,&a1,&a2,&a3));
 850 
 851 //        SFREE(gammat);
 852           
 853           /* solving the asymmetric system of equations */
 854 
 855           /* inputformat=1;
 856           symmetryflag=2;
 857           
 858           if(*isolver==0){
 859 #ifdef SPOOLES
 860               spooles(ad,au,adb,aub,&sigma,b,icol,irow,&neq,&nzs,
 861                       &symmetryflag,&inputformat,&nzs);
 862 #else
 863               printf("*ERROR in compfluid: the SPOOLES library is not linked\n\n");
 864               FORTRAN(stop,());
 865 #endif
 866           }
 867           else if(*isolver==7){
 868 #ifdef PARDISO
 869               pardiso_main(ad,au,adb,aub,&sigma,b,icol,irow,&neq,&nzs,
 870                            &symmetryflag,&inputformat,jq,&nzs);
 871 #else
 872               printf("*ERROR in compfluid: the PARDISO library is not linked\n\n");
 873               FORTRAN(stop,());
 874 #endif
 875 }*/
 876       
 877       /* free memory */
 878       
 879       /*    if(*isolver==0){
 880 #ifdef SPOOLES
 881           spooles_cleanup();
 882 #endif
 883       }
 884       else if(*isolver==7){
 885 #ifdef PARDISO
 886           pardiso_cleanup(&neq,&symmetryflag);
 887 #endif
 888 }*/
 889 
 890           memcpy(&au[2*nzs],ad,sizeof(double)*neq);
 891           nelt=2*nzs+neq;
 892           isym=0;
 893           itol=0;
 894           tol=0.;
 895           itmax=0;
 896           iunit=0;
 897           lrgw=131+16*neq;
 898           NNEW(rgwk,double,lrgw);
 899           NNEW(igwk,ITG,20);
 900           igwk[0]=10;
 901           igwk[1]=10;
 902           igwk[2]=0;
 903           igwk[3]=1;
 904           igwk[4]=10;
 905           ligw=20;
 906           for(i=0;i<neq;i++){rwork[i]=1./ad[i];}
 907           FORTRAN(predgmres,(&neq,&b[0],&vel[0],&nelt,irow,jq,au,
 908                              &isym,&itol,&tol,&itmax,&iter,
 909                              &err,&ierr,&iunit,sb,sx,rgwk,&lrgw,igwk,
 910                              &ligw,rwork,iwork));
 911           SFREE(rgwk);SFREE(igwk);
 912           if(ierr>0){
 913               printf("*WARNING in compfluid: error message from predgmres=%e\n",err);
 914           }
 915           
 916           /* storing the solution t into vel(0,*) */
 917           
 918 //        FORTRAN(calctel,(ne,nactdoh,vel,b,ne));
 919 
 920           /* extrapolation of the temperature at the element centers
 921              to the face centers */
 922 
 923           FORTRAN(extrapolate_tel,(nface,ielfa,xrlfa,vel,vfa,
 924                                    ifabou,xbounact,ipnei,nef));
 925 
 926           /* applying MPC's to the faces */
 927           
 928           if(*nmpc>0){
 929               is=0;ie=0;
 930               FORTRAN(applympc,(nface,ielfa,&is,&ie,ifabou,ipompc,vfa,coefmpc,
 931                                 nodempc,ipnei,neifa,labmpc,xbounact,nactdoh));
 932           }
 933 
 934       }
 935       
 936       /* storing the solution into vel */
 937       
 938 //      FORTRAN(calcvel,(ne,nactdoh,vel,bv,&neq,ne));
 939 
 940       /* extrapolating the velocity from the elements centers to the face
 941          centers, thereby taking the boundary conditions into account */
 942       
 943       FORTRAN(extrapolate_vel,(nface,ielfa,xrlfa,vel,vfa,
 944                                ifabou,xbounact,ipnei,nef));
 945 
 946       /* applying MPC's to the faces */
 947       
 948       if(*nmpc>0){
 949           is=1;ie=3;
 950           FORTRAN(applympc,(nface,ielfa,&is,&ie,ifabou,ipompc,vfa,coefmpc,
 951                             nodempc,ipnei,neifa,labmpc,xbounact,nactdoh));
 952       }
 953 
 954 //      FORTRAN(writevfa,(vfa,nface,nactdohinv,ielfa));
 955       
 956       
 957       if(((iit/jout[1])*jout[1]==iit)||(iconvergence==1)||
 958          (iit==jmax[1])){
 959 
 960           /* calculating the stress and the heat flow at the
 961              integration points, if requested */
 962 
 963           if((strcmp1(&filab[3306],"SF  ")==0)||
 964              (strcmp1(&filab[3480],"SVF ")==0))isti=1;
 965           if(strcmp1(&filab[3393],"HFLF")==0)iqfx=1;
 966           for(i=0;i<*nprint;i++){
 967               if(strcmp1(&prlab[6*i],"SVF")==0) isti=1;
 968               if(strcmp1(&prlab[6*i],"HFLF")==0)iqfx=1;
 969           }
 970 
 971           /* calculating the heat conduction at the element centers */
 972 
 973           if(iqfx==1){
 974               NNEW(hcel,double,*nef);
 975               FORTRAN(calchcel,(vel,cocon,ncocon,ielmatf,ntmat_,mi,
 976                                hcel,nef));
 977           }
 978 
 979           /* calculating the stress and/or the heat flux at the
 980              element centers */
 981 
 982           if((isti==1)||(iqfx==1)){
 983               FORTRAN(calcstressheatflux,(sti,umel,gradvel,qfx,hcel,
 984                                           gradtel,nef,&isti,&iqfx,mi));
 985               if(iqfx==1)SFREE(hcel);
 986           }
 987  
 988           /* extrapolating the stresses */
 989 
 990           if((strcmp1(&filab[3306],"SF  ")==0)||
 991              (strcmp1(&filab[3480],"SVF ")==0)){
 992               nfield=6;
 993               ndim=6;
 994               if((*norien>0)&&
 995                  ((strcmp1(&filab[3311],"L")==0)||(strcmp1(&filab[3485],"L")==0))){
 996                   iorienglob=1;
 997               }else{
 998                   iorienglob=0;
 999               }
1000               strcpy1(&cflag[0],&filab[2962],1);
1001               NNEW(stn,double,6**nk);
1002               FORTRAN(extrapolate,(sti,stn,ipkonf,inum,kon,lakonf,
1003                       &nfield,nk,nef,mi,&ndim,orab,ielorienf,co,&iorienglob,
1004                       cflag,nelemload,nload,nodeboun,nboun,ndirboun,
1005                       vold,ithermal,&force,&icfdout,ielmatf,thicke,filab));
1006           }
1007 
1008           /* extrapolating the heat flow */
1009 
1010           
1011           if(strcmp1(&filab[3393],"HFLF")==0){
1012               nfield=3;
1013               ndim=3;
1014               if((*norien>0)&&(strcmp1(&filab[3398],"L")==0)){
1015                   iorienglob=1;
1016               }else{
1017                   iorienglob=0;
1018               }
1019               strcpy1(&cflag[0],&filab[3049],1);
1020               NNEW(qfn,double,3**nk);
1021               FORTRAN(extrapolate,(qfx,qfn,ipkonf,inum,kon,lakonf,
1022                       &nfield,nk,nef,mi,&ndim,orab,ielorienf,co,&iorienglob,
1023                       cflag,nelemload,nload,nodeboun,nboun,ndirboun,
1024                       vold,ithermal,&force,&icfdout,ielmatf,thicke,filab));
1025           }
1026           
1027           /* extrapolating the facial values of the static temperature 
1028              and/or the velocity and/or the static pressure to the nodes */
1029           
1030           FORTRAN(extrapolatefluid,(nk,iponofa,inofa,inum,vfa,vold,ielfa,
1031                       ithermal));
1032 
1033           /* storing the results in dat-format */
1034 
1035           ptimef=ttimef+timef;
1036           FORTRAN(printoutfluid,(set,nset,istartset,iendset,ialset,nprint,
1037                                  prlab,prset,v,t1,fn,ipkonf,lakonf,sti,eei,
1038                                  xstate,ener,mi,&nstate_,ithermal,co,kon,qfx,
1039                                  &ptimef,trab,inotr,ntrans,orab,ielorienf,
1040                                  norien,nk,nef,inum,filab,vold,ielmatf,
1041                                  thicke,eme,vcontu,physcon,nactdoh,
1042                                  ielpropf,prop));
1043           
1044           /* storing the results in frd-format */
1045           
1046           FORTRAN(frdfluid,(co,nk,kon,ipkonf,lakonf,nef,vold,&kode,&timef,ielmatf,
1047                             matname,filab,inum,ntrans,inotr,trab,mi,istep,
1048                             stn,qfn,nactdohinv));
1049 
1050           if((strcmp1(&filab[3306],"SF  ")==0)||
1051              (strcmp1(&filab[3480],"SVF ")==0)){SFREE(stn);}
1052           if(strcmp1(&filab[3393],"HFLF")==0){SFREE(qfn);}
1053 
1054       }
1055       
1056       if(iit==jmax[1]){
1057           printf("*INFO: maximum number of fluid increments reached\n\n");
1058           FORTRAN(stop,());
1059       }
1060       if(last==1){FORTRAN(stop,());}
1061       
1062       memcpy(&veloo[0],&velo[0],sizeof(double)*6**nef);
1063       memcpy(&velo[0],&vel[0],sizeof(double)*6**nef);
1064       
1065   }while(1);
1066   
1067   FORTRAN(closefilefluid,());
1068 
1069   SFREE(irow);SFREE(ia);SFREE(icol);SFREE(jq);SFREE(ja);
1070 //SFREE(nactdoh);
1071   
1072   SFREE(coel);SFREE(cosa);SFREE(xxn);SFREE(xxi);SFREE(xle);SFREE(xlen);
1073   SFREE(xlet);SFREE(cofa);SFREE(area);SFREE(xrlfa);SFREE(volume);
1074   SFREE(cosb);SFREE(xxj);
1075 
1076   SFREE(ifabou);SFREE(umfa);SFREE(umel);
1077 
1078   SFREE(gradvel);SFREE(gradvfa);SFREE(au);SFREE(ad);SFREE(b);SFREE(advfa);
1079   SFREE(ap);SFREE(bp);SFREE(gradpel);SFREE(rwork);SFREE(aua);
1080   SFREE(hfa);SFREE(hel);SFREE(auv);SFREE(adv);SFREE(bv);SFREE(sel);
1081 
1082   if(*ithermal>0){
1083       SFREE(gradtel);SFREE(gradtfa);SFREE(hcfa);SFREE(cpel);SFREE(cpfa);
1084   }
1085 
1086   SFREE(inum);SFREE(v);SFREE(velo);SFREE(veloo);
1087 
1088   SFREE(iponofa);SFREE(inofa);
1089 
1090   if(*nbody>0) SFREE(body);
1091 
1092   *ithermal=ithermalref;
1093 
1094   return;
1095   
1096 } 

/* [<][>][^][v][top][bottom][index][help] */