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 #ifdef PASTIX
37 #include "pastix.h"
38 #endif
39
40 static char *sideload1,*covered1=NULL;
41
42 static ITG *kontri1,*nloadtr1,*idist=NULL,*ntrit1,*mi1,*jqrad1,
43 *irowrad1,*nzsrad1,num_cpus,*ntri1,*ntr1,ng1;
44
45 static double *vold1,*co1,*pmid1,*e11,*e21,*e31,*adview=NULL,*auview=NULL,*dist=NULL,
46 *area1,sidemean1;
47
radflowload(ITG * itg,ITG * ieg,ITG * ntg,ITG * ntr,double * adrad,double * aurad,double * bcr,ITG * ipivr,double * ac,double * bc,ITG * nload,char * sideload,ITG * nelemload,double * xloadact,char * lakon,ITG * ipiv,ITG * ntmat_,double * vold,double * shcon,ITG * nshcon,ITG * ipkon,ITG * kon,double * co,ITG * kontri,ITG * ntri,ITG * nloadtr,double * tarea,double * tenv,double * physcon,double * erad,double ** adviewp,double ** auviewp,ITG * nflow,ITG * ikboun,double * xbounact,ITG * nboun,ITG * ithermal,ITG * iinc,ITG * iit,double * cs,ITG * mcs,ITG * inocs,ITG * ntrit,ITG * nk,double * fenv,ITG * istep,double * dtime,double * ttime,double * time,ITG * ilboun,ITG * ikforc,ITG * ilforc,double * xforcact,ITG * nforc,double * cam,ITG * ielmat,ITG * nteq,double * prop,ITG * ielprop,ITG * nactdog,ITG * nacteq,ITG * nodeboun,ITG * ndirboun,ITG * network,double * rhcon,ITG * nrhcon,ITG * ipobody,ITG * ibody,double * xbodyact,ITG * nbody,ITG * iviewfile,char * jobnamef,double * ctrl,double * xloadold,double * reltime,ITG * nmethod,char * set,ITG * mi,ITG * istartset,ITG * iendset,ITG * ialset,ITG * nset,ITG * ineighe,ITG * nmpc,ITG * nodempc,ITG * ipompc,double * coefmpc,char * labmpc,ITG * iemchange,ITG * nam,ITG * iamload,ITG * jqrad,ITG * irowrad,ITG * nzsrad,ITG * icolrad,ITG * ne,ITG * iaxial,double * qa,double * cocon,ITG * ncocon,ITG * iponoel,ITG * inoel,ITG * nprop,char * amname,ITG * namta,double * amta)48 void radflowload(ITG *itg,ITG *ieg,ITG *ntg,ITG *ntr,double *adrad,
49 double *aurad,double *bcr,ITG *ipivr,
50 double *ac,double *bc,ITG *nload,char *sideload,
51 ITG *nelemload,double *xloadact,char *lakon,ITG *ipiv,
52 ITG *ntmat_,double *vold,double *shcon,
53 ITG *nshcon,ITG *ipkon,ITG *kon,double *co,
54 ITG *kontri,
55 ITG *ntri,ITG *nloadtr,double *tarea,double *tenv,
56 double *physcon,double *erad,double **adviewp,
57 double **auviewp,
58 ITG *nflow,ITG *ikboun,
59 double *xbounact,ITG *nboun,ITG *ithermal,
60 ITG *iinc,ITG *iit,double *cs, ITG *mcs, ITG *inocs,
61 ITG *ntrit,ITG *nk, double *fenv,ITG *istep,double *dtime,
62 double *ttime,double *time,ITG *ilboun,ITG *ikforc,
63 ITG *ilforc,double *xforcact,ITG *nforc,double *cam,
64 ITG *ielmat,ITG *nteq,double *prop,ITG *ielprop,ITG *nactdog,
65 ITG *nacteq,ITG *nodeboun,ITG *ndirboun,
66 ITG *network, double *rhcon, ITG *nrhcon, ITG *ipobody,
67 ITG *ibody, double *xbodyact, ITG *nbody,ITG *iviewfile,
68 char *jobnamef, double *ctrl, double *xloadold,
69 double *reltime, ITG *nmethod, char *set, ITG *mi,
70 ITG * istartset,ITG* iendset,ITG *ialset,ITG *nset,
71 ITG *ineighe, ITG *nmpc, ITG *nodempc,ITG *ipompc,
72 double *coefmpc,char *labmpc, ITG *iemchange,ITG *nam,
73 ITG *iamload,ITG *jqrad,ITG *irowrad,ITG *nzsrad,
74 ITG *icolrad,ITG *ne,ITG *iaxial,double *qa,
75 double *cocon,ITG *ncocon,ITG *iponoel,ITG *inoel,
76 ITG *nprop,char *amname,ITG *namta,double *amta){
77
78 /* network=0: no network
79 network=1: purely thermal (presence of "Dx"- and/or of "D " network elements; declared
80 by the user to be purely thermal (on the *STEP card); simultaneous solution)
81 network=2: purely thermal (alternating solution; this becomes a simultaneous solution in
82 the absence of "Dx"-elements)
83 network=3: general case (temperatures, fluxes and pressures unknown)
84 network=4: purely aerodynamic, i.e. only fluxes and pressures unknown
85
86 "D "-elements (D followed by a blank) alone do not trigger the alternating solution
87 (are not counted in envtemp.f as true network elements) */
88
89 ITG nhrs=1,info=0,i,j,iin=0,icntrl,icutb=0,iin_abs=0,mt=mi[1]+1,im,
90 symmetryflag=2,inputformat=1,node,ichannel,*ithread=NULL,iplausi;
91
92 static ITG ifactorization=0;
93
94 double camt[2],camf[2],camp[2],qat,qaf,ramt,ramf,ramp,
95 cam1t=0.,cam1f=0.,cam1p=0.,sidemean,qa0,qau,ea,*prop_store=NULL,
96 cam2t=0.,cam2f=0.,cam2p=0.,dtheta=1.,*v=NULL,cama[2],cam1a=0.,
97 cam2a=0.,vamt=0.,vamf=0.,vamp=0.,vama=0.,cam0t=0.,cam0f=0.,
98 cam0p=0.,cam0a=0.,sigma=0.,*adbrad=NULL,*aubrad=NULL,*q=NULL,
99 *area=NULL,*pmid=NULL,*e1=NULL,*e2=NULL,*e3=NULL,
100 qamt,qamf,qamtold,qamfold;
101
102 adview=*adviewp;auview=*auviewp;
103
104 qa0=ctrl[20];qau=ctrl[21];ea=ctrl[23];
105
106 /* check whether there are any gas temperature nodes; this check should
107 NOT be done on nteq, since also for zero equations the temperature
108 of the gas nodes with boundary conditions must be stored in v
109 (in initialgas) */
110
111 NNEW(v,double,mt**nk);
112
113 /* gas networks */
114
115 if(*ntg!=0) {
116 #ifdef COMPANY
117 NNEW(prop_store,double,*nprop);
118 memcpy(&prop_store[0],&prop[0],sizeof(double)**nprop);
119 FORTRAN(propertynet,(ieg,nflow,prop,ielprop,lakon,&iin,
120 prop_store,ttime,time,nam,amname,namta,amta));
121 FORTRAN(checkinputvaluesnet,(ieg,nflow,prop,ielprop,lakon));
122 #else
123 if(*iit==-1) FORTRAN(checkinputvaluesnet,(ieg,nflow,prop,ielprop,lakon));
124 #endif
125 icntrl=0;
126 while(icntrl==0) {
127
128 if(iin==0){
129
130 /* since in nonlingeo.c radflowload.c is called before
131 results.c at the start of a new increment vold does
132 not contain the spc's at the end of the increment yet;
133 for this purpose voldwithspc is locally introduced
134 in radflowload.c */
135
136 memcpy(&v[0],&vold[0],sizeof(double)*mt**nk);
137
138 /* resetting ineighe to 0 for renewed call of
139 radflowload (not for cut-backs without
140 leaving radflowload */
141
142 /* for a cut-back iin is reset to 0, iin_abs is not */
143
144 if(iin_abs==0) DMEMSET(ineighe,0,*ntg,0);
145
146 /* initialization pressurized flow
147 (no free surface: gas networks or
148 water networks with fully wetted perimeter*/
149
150 FORTRAN(initialnet,(itg,ieg,ntg,ac,bc,lakon,v,
151 ipkon,kon,nflow,
152 ikboun,nboun,prop,ielprop,nactdog,ndirboun,
153 nodeboun,xbounact,ielmat,ntmat_,shcon,nshcon,
154 physcon,ipiv,nteq,rhcon,nrhcon,ipobody,ibody,
155 xbodyact,co,nbody,network,&iin_abs,vold,set,
156 istep,iit,mi,ineighe,ilboun,&ichannel,iaxial,
157 nmpc,labmpc,ipompc,nodempc,coefmpc,ttime,time,
158 iponoel,inoel));
159
160 /* initialization for channels with free surface */
161
162 if(ichannel==1){
163 FORTRAN(initialchannel,(itg,ieg,ntg,ac,bc,lakon,v,
164 ipkon,kon,nflow,
165 ikboun,nboun,prop,ielprop,nactdog,ndirboun,
166 nodeboun,xbounact,ielmat,ntmat_,shcon,nshcon,
167 physcon,ipiv,nteq,rhcon,nrhcon,ipobody,ibody,
168 xbodyact,co,nbody,network,&iin_abs,vold,set,
169 istep,iit,mi,ineighe,ilboun,ttime,time,iaxial));
170 }
171
172 /* storing the residual in the rhs vector */
173
174 FORTRAN(resultnet,(itg,ieg,ntg,bc,nload,sideload,
175 nelemload,xloadact,
176 lakon,ntmat_,v,shcon,nshcon,ipkon,kon,co,nflow,
177 iinc,istep,dtime,ttime,time,
178 ikforc,ilforc,xforcact,
179 nforc,ielmat,nteq,prop,ielprop,nactdog,nacteq,&iin,
180 physcon,camt,camf,camp,rhcon,nrhcon,ipobody,
181 ibody,xbodyact,nbody,&dtheta,vold,xloadold,
182 reltime,nmethod,set,mi,ineighe,cama,&vamt,
183 &vamf,&vamp,&vama,nmpc,nodempc,ipompc,coefmpc,
184 labmpc,iaxial,&qat,&qaf,&ramt,&ramf,&ramp,
185 cocon,ncocon,iponoel,inoel,&iplausi));
186
187 /* iniializing qamt and qamf (mean typical energy flow
188 and mass flow */
189
190 if(qau>1.e-10){qamt=qau;}
191 else if(qa0>1.e-10){qamt=qa0;}
192 else if(qat>1.e-10){qamt=qat;}
193 else {qamt=1.e-2;}
194
195 if(qau>1.e-10){qamf=qau;}
196 else if(qa0>1.e-10){qamf=qa0;}
197 else if(qaf>1.e-10){qamf=qaf;}
198 else {qamf=1.e-2;}
199 }
200
201 iin++;
202 iin_abs++;
203 printf(" gas iteration %" ITGFORMAT " \n \n",iin);
204
205 /* store actual values of typical energy flow and
206 mass flow */
207
208 qamtold=qamt;
209 qamfold=qamf;
210
211 /* filling the lhs matrix */
212
213 FORTRAN(mafillnet,(itg,ieg,ntg,ac,nload,sideload,
214 nelemload,xloadact,lakon,ntmat_,v,
215 shcon,nshcon,ipkon,kon,co,nflow,iinc,
216 istep,dtime,ttime,time,
217 ielmat,nteq,prop,ielprop,nactdog,nacteq,
218 physcon,rhcon,nrhcon,ipobody,ibody,xbodyact,
219 nbody,vold,xloadold,reltime,nmethod,set,mi,
220 nmpc,nodempc,ipompc,coefmpc,labmpc,iaxial,
221 cocon,ncocon,iponoel,inoel));
222
223 /* solving the system of equations */
224
225 if(*nteq>0){
226 FORTRAN(dgesv,(nteq,&nhrs,ac,nteq,ipiv,bc,nteq,&info));
227 }
228
229 /*spooles(ac,au,adb,aub,&sigma,bc,icol,irow,nteq,nteq,
230 &symmetryflag,&inputformat);*/
231
232 if (info!=0) {
233 printf(" *WARNING in radflowload: singular matrix\n");
234
235 FORTRAN(mafillnet,(itg,ieg,ntg,ac,nload,sideload,
236 nelemload,xloadact,lakon,ntmat_,v,
237 shcon,nshcon,ipkon,kon,co,nflow,iinc,
238 istep,dtime,ttime,time,
239 ielmat,nteq,prop,ielprop,nactdog,nacteq,
240 physcon,rhcon,nrhcon,ipobody,ibody,xbodyact,
241 nbody,vold,xloadold,reltime,nmethod,set,mi,
242 nmpc,nodempc,ipompc,coefmpc,labmpc,iaxial,
243 cocon,ncocon,iponoel,inoel));
244
245 FORTRAN(equationcheck,(ac,nteq,nactdog,itg,ntg,nacteq,network));
246
247 iin=0;
248
249 }
250 else {
251
252 /* storing the residual in the rhs vector */
253
254 FORTRAN(resultnet,(itg,ieg,ntg,bc,nload,sideload,nelemload,
255 xloadact,lakon,ntmat_,v,shcon,nshcon,ipkon,kon,co,
256 nflow,iinc,istep,dtime,ttime,time,ikforc,ilforc,xforcact,
257 nforc,ielmat,nteq,prop,ielprop,nactdog,nacteq,
258 &iin,physcon,camt,camf,camp,rhcon,nrhcon,ipobody,
259 ibody,xbodyact,nbody,&dtheta,vold,xloadold,
260 reltime,nmethod,set,mi,ineighe,cama,&vamt,
261 &vamf,&vamp,&vama,nmpc,nodempc,ipompc,coefmpc,labmpc,
262 iaxial,&qat,&qaf,&ramt,&ramf,&ramp,cocon,ncocon,iponoel,
263 inoel,&iplausi));
264
265 /* updating the mean typical energy flow and mass flow */
266
267 if(qau<1.e-10){
268 if(qat>ea*qamt){qamt=(qamtold*iin+qat)/(iin+1);}
269 else {qamt=qamtold;}
270 if(qaf>ea*qamf){qamf=(qamfold*iin+qaf)/(iin+1);}
271 else {qamf=qamfold;}
272 }
273
274 /* printing the largest corrections */
275
276 if(*network!=4){
277 cam2t=cam1t;
278 cam1t=cam0t;
279 cam0t=camt[0];
280 printf
281 (" mean typical energy flow since start of network iterations= %e\n",qamt);
282 printf
283 (" largest energy flow residual in present network iteration= %e\n",ramt);
284 printf
285 (" largest change of gas temperature since start of network iterations= %e\n",vamt);
286 if((ITG)camt[1]==0){
287 printf
288 (" largest correction to gas temperature in present network iteration= %e\n\n",
289 fabs(camt[0]));
290 }else{
291 printf
292 (" largest correction to gas temperature in present network iteration= %e in node %" ITGFORMAT "\n\n",
293 fabs(camt[0]),(ITG)camt[1]);
294 }
295 }
296
297 if(*network>2){
298 cam2f=cam1f;
299 cam1f=cam0f;
300 cam0f=camf[0];
301 printf
302 (" mean typical mass flow since start of network iterations= %e\n",qamf);
303 printf
304 (" largest mass flow residual in present network iteration= %e\n",ramf);
305 printf(" largest change of gas massflow since start of network iterations= %e\n",vamf);
306 if((ITG)camf[1]==0){
307 printf(" largest correction to gas massflow in present network iteration= %e\n\n",
308 fabs(camf[0]));
309 }else{
310 printf(" largest correction to gas massflow in present network iteration= %e in node %" ITGFORMAT "\n\n",
311 fabs(camf[0]),(ITG)camf[1]);
312 }
313
314 cam2p=cam1p;
315 cam1p=cam0p;
316 cam0p=camp[0];
317 printf
318 (" largest element equation residual in present network iteration= %e\n",ramp);
319 printf(" largest change of gas pressure since start of network iterations= %e\n",vamp);
320 if((ITG)camp[1]==0){
321 printf(" largest correction to gas pressure in present network iteration= %e\n\n",
322 fabs(camp[0]));
323 }else{
324 printf(" largest correction to gas pressure in present network iteration= %e in node %" ITGFORMAT "\n\n",
325 fabs(camp[0]),(ITG)camp[1]);
326 }
327
328 cam2a=cam1a;
329 cam1a=cam0a;
330 cam0a=cama[0];
331 printf(" largest change of geometry since start of network iterations= %e\n",vama);
332 if((ITG)cama[1]==0){
333 printf(" largest correction to geometry in present network iteration= %e\n",
334 fabs(cama[0]));
335 }else{
336 printf(" largest correction to geometry in present network iteration= %e in node %" ITGFORMAT "\n",
337 fabs(cama[0]),(ITG)cama[1]);
338 }
339 }
340 }
341
342 printf("\n");
343
344 /* for purely thermal calculations no iterations are
345 deemed necessary */
346
347 if(*network<=2) {icntrl=1;}
348 else {
349
350 /* check the convergence */
351
352 checkconvnet(&icutb,&iin,
353 &cam1t,&cam1f,&cam1p,&cam2t,&cam2f,&cam2p,&cam0t,&cam0f,
354 &cam0p,&icntrl,&dtheta,ctrl,&cam1a,&cam2a,&cam0a,
355 &vamt,&vamf,&vamp,&vama,qa,&qamt,&qamf,&ramt,&ramf,&ramp,
356 &iplausi,&ichannel,iaxial);
357 }
358 }
359
360 /* storing network output as boundary conditions for
361 the structure */
362
363 FORTRAN(flowresult,(ntg,itg,cam,vold,v,nload,sideload,
364 nelemload,xloadact,nactdog,network,mi,ne,ipkon,lakon,kon));
365
366 /* extra output for hydraulic jump (fluid channels) */
367
368 #ifdef NETWORKOUT
369 if(*network>2){
370 FORTRAN(flowoutput,(itg,ieg,ntg,nteq,bc,lakon,ntmat_,
371 v,shcon,nshcon,ipkon,kon,co,nflow, dtime,ttime,time,
372 ielmat,prop,ielprop,nactdog,nacteq,&iin,physcon,
373 camt,camf,camp,rhcon,nrhcon,
374 vold,jobnamef,set,istartset,iendset,ialset,nset,
375 mi,iaxial,istep,iit));
376 }
377 #endif
378 #ifdef COMPANY
379 memcpy(&prop[0],&prop_store[0],sizeof(double)**nprop);
380 SFREE(prop_store);
381 #endif
382 }
383
384 /* radiation */
385
386 if(*ntr>0){
387
388 /* variables for multithreading procedure */
389
390 ITG sys_cpus;
391 char *env,*envloc,*envsys;
392
393 num_cpus = 0;
394 sys_cpus=0;
395
396 /* explicit user declaration prevails */
397
398 envsys=getenv("NUMBER_OF_CPUS");
399 if(envsys){
400 sys_cpus=atoi(envsys);
401 if(sys_cpus<0) sys_cpus=0;
402 }
403
404 /* automatic detection of available number of processors */
405
406 if(sys_cpus==0){
407 sys_cpus = getSystemCPUs();
408 if(sys_cpus<1) sys_cpus=1;
409 }
410
411 /* local declaration prevails, if strictly positive */
412
413 envloc = getenv("CCX_NPROC_VIEWFACTOR");
414 if(envloc){
415 num_cpus=atoi(envloc);
416 if(num_cpus<0){
417 num_cpus=0;
418 }else if(num_cpus>sys_cpus){
419 num_cpus=sys_cpus;
420 }
421
422 }
423
424 /* else global declaration, if any, applies */
425
426 env = getenv("OMP_NUM_THREADS");
427 if(num_cpus==0){
428 if (env)
429 num_cpus = atoi(env);
430 if (num_cpus < 1) {
431 num_cpus=1;
432 }else if(num_cpus>sys_cpus){
433 num_cpus=sys_cpus;
434 }
435 }
436
437 // next line is to be inserted in a similar way for all other paralell parts
438
439 if(*ntr<num_cpus) num_cpus=*ntr;
440
441 pthread_t tid[num_cpus];
442
443 /* update of vold for the boundary conditions (only in the first
444 iteration of each increment and if not already done in
445 initialnet.f) */
446
447 /* in principle ok, does lead for thermomech.inp to much more
448 iterations, results the same; */
449
450 /* if(*iit<1){
451 if(*ntg==0){
452 for(i=0;i<*nboun;i++){
453 vold[mt*(nodeboun[i]-1)+ndirboun[i]]=xbounact[i];
454 }
455 }
456 }*/
457
458 /*the default sink temperature is updated at the start of each
459 increment */
460
461 for(i=0;i<*ntr;i++){
462 node=nelemload[2*nloadtr[i]-1];
463 if(node!=0){
464 tenv[i]=vold[mt*(node-1)]-physcon[0];
465 }else if(*iit<=0){
466 tenv[i]=xloadact[2*nloadtr[i]-1]-physcon[0];
467 }
468 }
469
470 /* for pure thermal steps the viewfactors have to be
471 calculated only once, for thermo-mechanical steps
472 (ithermal=3) they are recalculated in each iteration
473 unless they are read from file */
474
475 if(((*ithermal==3)&&(*iviewfile>=0))||(*iit==-1)){
476 if(*iviewfile<0){
477
478 /* reading viewfactors from file */
479
480 FORTRAN(readview,(ntr,adview,auview,fenv,nzsrad,ithermal,
481 jobnamef));
482
483 }else{
484
485 /* determining geometric data to calculate the viewfactors */
486
487 NNEW(area,double,*ntrit);
488 NNEW(pmid,double,3**ntrit);
489 NNEW(e1,double,3**ntrit);
490 NNEW(e2,double,3**ntrit);
491 NNEW(e3,double,4**ntrit);
492
493 FORTRAN(geomview,(vold,co,pmid,e1,e2,e3,kontri,area,
494 cs,mcs,inocs,ntrit,nk,mi,&sidemean));
495
496 RENEW(adview,double,num_cpus**ntr);
497 RENEW(auview,double,num_cpus*2**nzsrad);
498
499 NNEW(dist,double,num_cpus**ntrit);
500 NNEW(idist,ITG,num_cpus**ntrit);
501
502 DMEMSET(adview,0,num_cpus**ntr,0.);
503 DMEMSET(auview,0,num_cpus*2**nzsrad,0.);
504
505 sideload1=sideload;vold1=vold;co1=co;pmid1=pmid;
506 e11=e1;e21=e2;e31=e3;kontri1=kontri;ntr1=ntr;
507 nloadtr1=nloadtr;area1=area;ntri1=ntri;
508 ntrit1=ntrit;mi1=mi;jqrad1=jqrad;irowrad1=irowrad;
509 nzsrad1=nzsrad;sidemean1=sidemean;
510
511 /* size of the square mesh used to detect
512 the visibility of a triangle; the denser
513 the mesh,the more accurate the results */
514
515 ng1=1280;
516 // ng1=2560;
517 NNEW(covered1,char,num_cpus*ng1*ng1);
518
519 /* calculating the viewfactors */
520
521 printf(" Using up to %" ITGFORMAT " cpu(s) for the viewfactor calculation.\n\n", num_cpus);
522
523 /* create threads and wait */
524
525 NNEW(ithread,ITG,num_cpus);
526 for(i=0; i<num_cpus; i++) {
527 ithread[i]=i;
528 pthread_create(&tid[i], NULL, (void *)calcviewmt, (void *)&ithread[i]);
529 }
530 for(i=0; i<num_cpus; i++) pthread_join(tid[i], NULL);
531
532 for(i=0;i<*ntr;i++){
533 for(j=1;j<num_cpus;j++){
534 adview[i]+=adview[i+j**ntr];
535 }
536 }
537 RENEW(adview,double,*ntr);
538
539 for(i=0;i<2**nzsrad;i++){
540 for(j=1;j<num_cpus;j++){
541 auview[i]+=auview[i+j*2**nzsrad];
542 }
543 }
544 RENEW(auview,double,2**nzsrad);
545
546 SFREE(dist);SFREE(idist);SFREE(e1);SFREE(e2);SFREE(e3);
547 SFREE(pmid);SFREE(ithread);SFREE(covered1);
548
549 /* postprocessing the viewfactors */
550
551 FORTRAN(postview,(ntr,sideload,nelemload,kontri,ntri,nloadtr,
552 tenv,adview,auview,area,fenv,jqrad,irowrad,
553 nzsrad));
554
555 SFREE(area);
556
557 if(*iviewfile>=2){
558
559 /* writing viewfactors to file */
560
561 FORTRAN(writeview,(ntr,adview,auview,fenv,nzsrad,
562 jobnamef));
563 }
564
565 if(*iviewfile==3){
566
567 /* calculation of viewfactors only */
568
569 FORTRAN(stop,());
570 }
571
572 }
573 }
574
575 /* assembling the radiation matrix */
576
577 FORTRAN(radmatrix,(ntr,adrad,aurad,bcr,sideload,nelemload,
578 xloadact,lakon,vold,ipkon,kon,co,nloadtr,tarea,tenv,physcon,
579 erad,adview,auview,ithermal,iinc,iit,fenv,istep,dtime,ttime,
580 time,iviewfile,xloadold,reltime,nmethod,mi,iemchange,nam,
581 iamload,jqrad,irowrad,nzsrad));
582
583 /* factoring the system of equations */
584
585 /* the left hand side of the radiation matrix has probably
586 changed if
587 - the viewfactors were updated
588 - a new step was started and NO CHANGE is not active
589 - the emissivity coefficients were changed
590 - a new increment was started in a stationary calculation
591 (since the emissivity coefficients are ramped)
592 in that case the LU decomposition has to be repeated
593 (i.e. call of dgesv) */
594
595 if(((*ithermal==3)&&(*iviewfile>=0))||
596 ((*iit==-1)&&(*iviewfile!=-2))||
597 (*iemchange==1)||
598 ((*iit==0)&&(abs(*nmethod)==1))){
599
600 #if defined(PASTIX)
601 ITG symmetry = 1;
602 ITG inputformat = 1;
603 pastix_factor_main_as(adrad,aurad,adbrad,aubrad,&sigma,icolrad,
604 irowrad,ntr,nzsrad, &symmetry, &inputformat, jqrad, nzsrad);
605 ifactorization=1;
606 #elif defined(PARDISO)
607 if(ifactorization==1) pardiso_cleanup_as(ntr,&symmetryflag);
608 pardiso_factor_as(adrad,aurad,adbrad,aubrad,&sigma,icolrad,
609 irowrad,ntr,nzsrad,jqrad);
610 ifactorization=1;
611 #elif defined(SPOOLES)
612 if(ifactorization==1) spooles_cleanup_rad();
613 spooles_factor_rad(adrad,aurad,adbrad,aubrad,&sigma,
614 icolrad,irowrad,ntr,nzsrad,
615 &symmetryflag,&inputformat);
616 ifactorization=1;
617 #else
618 printf("*ERROR in radflowload: the SPOOLES library is not linked\n\n");
619 FORTRAN(stop,());
620 #endif
621
622 }
623
624 /* solving the system of equations */
625
626 #if defined(PASTIX)
627 ITG symmetry = 1;
628 ITG nrhs = 1;
629 pastix_solve_as(bcr,ntr,&symmetry,&nrhs);
630 #elif defined(PARDISO)
631 pardiso_solve_as(bcr,ntr);
632
633 #elif defined(SPOOLES)
634 spooles_solve_rad(bcr,ntr);
635 #endif
636
637 if (info!=0){
638 printf("*ERROR IN RADFLOWLOAD: SINGULAR MATRIX*\n");}
639
640 else{
641 NNEW(q,double,*ntr);
642 FORTRAN(radresult,(ntr,xloadact,bcr,nloadtr,tarea,
643 tenv,physcon,erad,auview,fenv,
644 irowrad,jqrad,nzsrad,q));
645 SFREE(q);
646 }
647
648 }
649
650 SFREE(v);
651
652 *adviewp=adview;*auviewp=auview;
653
654 return;
655
656 }
657
658 /* subroutine for multithreading of calcview */
659
calcviewmt(ITG * i)660 void *calcviewmt(ITG *i){
661
662 ITG indexad,indexau,indexdi,ntria,ntrib,nedelta,indexcovered;
663
664 indexad=*i**ntr1;
665 indexau=*i*2**nzsrad1;
666 indexdi=*i**ntrit1;
667 indexcovered=*i*ng1*ng1;
668
669 nedelta=(ITG)ceil(*ntri1/(double)num_cpus);
670 ntria=*i*nedelta+1;
671 ntrib=(*i+1)*nedelta;
672 if(ntrib>*ntri1) ntrib=*ntri1;
673
674 // printf("i=%" ITGFORMAT ",ntria=%" ITGFORMAT ",ntrib=%" ITGFORMAT "\n",i,ntria,ntrib);
675 // printf("indexad=%" ITGFORMAT ",indexau=%" ITGFORMAT ",indexdi=%" ITGFORMAT "\n",indexad,indexau,indexdi);
676
677 FORTRAN(calcview,(sideload1,vold1,co1,pmid1,e11,e21,e31,
678 kontri1,nloadtr1,&adview[indexad],
679 &auview[indexau],&dist[indexdi],&idist[indexdi],area1,
680 ntrit1,mi1,jqrad1,irowrad1,nzsrad1,&sidemean1,
681 &ntria,&ntrib,&covered1[indexcovered],&ng1));
682
683 return NULL;
684 }
685
686