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
25 static ITG *ielfa1,*icyclic1,*ifatie1,*ifabou1,*ipnei1,*nef1,*num_cpus1,
26 *neifa1,*nface1,*ncfd1,*inlet1;
27
28 static double *xrlfa1,*vfap1,*vel1,*c1,*xbounact1,*xxn1,*area1,*volume1,
29 *gradkel1,*gradkfa1,*xxi1,*vfa1,*rf1,*xle1,*xlet1,*xxj1,*umfa1,
30 constant1,*xxna1;
31
extrapol_kelmain(ITG * nface,ITG * ielfa,double * xrlfa,double * vel,double * vfa,ITG * ifabou,double * xbounact,ITG * nef,double * gradkel,double * gradkfa,ITG * neifa,double * rf,double * area,double * volume,double * xle,double * xxi,ITG * icyclic,double * xxn,ITG * ipnei,ITG * ifatie,double * xlet,double * xxj,double * coefmpc,ITG * nmpc,char * labmpc,ITG * ipompc,ITG * nodempc,ITG * ifaext,ITG * nfaext,ITG * nactdoh,double * umfa,double * physcon,ITG * iitg,double * c,ITG * num_cpus,ITG * compressible,double * xxna,ITG * ncfd,ITG * inlet)32 void extrapol_kelmain(ITG *nface,ITG *ielfa,double *xrlfa,double *vel,
33 double *vfa,ITG *ifabou,double *xbounact,ITG *nef,
34 double *gradkel,double *gradkfa,ITG *neifa,double *rf,
35 double *area,double *volume,double *xle,double *xxi,
36 ITG *icyclic,double *xxn,ITG *ipnei,ITG *ifatie,
37 double *xlet,double *xxj,double *coefmpc,
38 ITG *nmpc,char *labmpc,ITG *ipompc,ITG *nodempc,ITG *ifaext,
39 ITG *nfaext,ITG *nactdoh,double *umfa,double *physcon,ITG *iitg,
40 double *c,ITG *num_cpus,ITG *compressible,double *xxna,
41 ITG *ncfd,ITG *inlet){
42
43 ITG i,is,ie,m;
44
45 double *vfap=NULL,constant;
46
47 /* variables for multithreading procedure */
48
49 ITG *ithread=NULL;
50
51 num_cpus1=num_cpus;
52
53 pthread_t tid[*num_cpus];
54
55 NNEW(vfap,double,8**nface);
56 constant=1.7393e-3*physcon[4]/(physcon[6]*physcon[7]);
57
58 /* inter/extrapolation of k at the center of the elements
59 to the center of the faces */
60
61 ielfa1=ielfa;xrlfa1=xrlfa;vfap1=vfap;vel1=vel;ifabou1=ifabou;
62 nef1=nef;nface1=nface;constant1=constant;umfa1=umfa;inlet1=inlet;
63
64 /* create threads and wait */
65
66 NNEW(ithread,ITG,*num_cpus);
67 for(i=0; i<*num_cpus; i++) {
68 ithread[i]=i;
69 pthread_create(&tid[i], NULL, (void *)extrapol_kel1mt, (void *)&ithread[i]);
70 }
71 for(i=0; i<*num_cpus; i++) pthread_join(tid[i], NULL);
72
73 SFREE(ithread);
74
75 /* multiple point constraints */
76
77 if(*nmpc>0){
78 is=6;
79 ie=6;
80 FORTRAN(applympc,(nface,ielfa,&is,&ie,ifabou,ipompc,vfap,coefmpc,
81 nodempc,ipnei,neifa,labmpc,xbounact,nactdoh,
82 ifaext,nfaext));
83 }
84
85 /* calculate the gradient of the turbulent kinetic energy at the center of
86 the elements */
87
88 ipnei1=ipnei;neifa1=neifa;vfa1=vfap;area1=area;xxna1=xxna;volume1=volume;
89 gradkel1=gradkel;nef1=nef;ncfd1=ncfd;
90
91 /* create threads and wait */
92
93 NNEW(ithread,ITG,*num_cpus);
94 for(i=0; i<*num_cpus; i++) {
95 ithread[i]=i;
96 pthread_create(&tid[i], NULL, (void *)extrapol_kel2mt, (void *)&ithread[i]);
97 }
98 for(i=0; i<*num_cpus; i++) pthread_join(tid[i], NULL);
99
100 SFREE(ithread);
101
102 /* interpolate/extrapolate the turbulent kinetic energy gradient from the
103 center of the elements to the center of the faces */
104
105 ielfa1=ielfa;xrlfa1=xrlfa;icyclic1=icyclic;ifatie1=ifatie;gradkfa1=gradkfa;
106 gradkel1=gradkel;c1=c;ipnei1=ipnei;xxi1=xxi;nface1=nface;ncfd1=ncfd;
107
108 /* create threads and wait */
109
110 NNEW(ithread,ITG,*num_cpus);
111 for(i=0; i<*num_cpus; i++) {
112 ithread[i]=i;
113 pthread_create(&tid[i], NULL, (void *)extrapol_kel3mt, (void *)&ithread[i]);
114 }
115 for(i=0; i<*num_cpus; i++) pthread_join(tid[i], NULL);
116
117 SFREE(ithread);
118
119 /* correct the facial turbulent kinetic energy gradients:
120 Moukalled et al. p 289 */
121
122 ielfa1=ielfa;ipnei1=ipnei;vel1=vel;xlet1=xlet;gradkfa1=gradkfa;xxj1=xxj;
123 nef1=nef;nface1=nface;ncfd1=ncfd;
124
125 /* create threads and wait */
126
127 NNEW(ithread,ITG,*num_cpus);
128 for(i=0; i<*num_cpus; i++) {
129 ithread[i]=i;
130 pthread_create(&tid[i], NULL, (void *)extrapol_kel5mt, (void *)&ithread[i]);
131 }
132 for(i=0; i<*num_cpus; i++) pthread_join(tid[i], NULL);
133
134 SFREE(ithread);
135
136 /* inter/extrapolation of k at the center of the elements
137 to the center of the faces: taking the skewness of the
138 elements into account (using the gradient at the center
139 of the faces)
140 Moukalled et al. p 279 */
141
142 ielfa1=ielfa;vfa1=vfa;vfap1=vfap;gradkfa1=gradkfa;rf1=rf;ifabou1=ifabou;
143 ipnei1=ipnei;vel1=vel;xxi1=xxi;xle1=xle;nef1=nef;nface1=nface;inlet1=inlet;
144
145 /* create threads and wait */
146
147 NNEW(ithread,ITG,*num_cpus);
148 for(i=0; i<*num_cpus; i++) {
149 ithread[i]=i;
150 pthread_create(&tid[i], NULL, (void *)extrapol_kel4mt, (void *)&ithread[i]);
151 }
152 for(i=0; i<*num_cpus; i++) pthread_join(tid[i], NULL);
153
154 SFREE(ithread);
155
156 /* multiple point constraints */
157
158 if(*nmpc>0){
159 is=6;
160 ie=6;
161 FORTRAN(applympc,(nface,ielfa,&is,&ie,ifabou,ipompc,vfa,coefmpc,
162 nodempc,ipnei,neifa,labmpc,xbounact,nactdoh,
163 ifaext,nfaext));
164 }
165
166 /* correction loops */
167
168 for(m=0;m<*iitg;m++){
169
170 /* calculate the gradient of the turbulent kinetic energies at the center of
171 the elements */
172
173 ipnei1=ipnei;neifa1=neifa;vfa1=vfa;area1=area;xxna1=xxna;volume1=volume;
174 gradkel1=gradkel;nef1=nef;ncfd1=ncfd;
175
176 /* create threads and wait */
177
178 NNEW(ithread,ITG,*num_cpus);
179 for(i=0; i<*num_cpus; i++) {
180 ithread[i]=i;
181 pthread_create(&tid[i], NULL, (void *)extrapol_kel2mt, (void *)&ithread[i]);
182 }
183 for(i=0; i<*num_cpus; i++) pthread_join(tid[i], NULL);
184
185 SFREE(ithread);
186
187 /* interpolate/extrapolate the turbulent kinetic energy gradient from the
188 center of the elements to the center of the faces */
189
190 ielfa1=ielfa;xrlfa1=xrlfa;icyclic1=icyclic;ifatie1=ifatie;gradkfa1=gradkfa;
191 gradkel1=gradkel;c1=c;ipnei1=ipnei;xxi1=xxi;nface1=nface;ncfd1=ncfd;
192
193 /* create threads and wait */
194
195 NNEW(ithread,ITG,*num_cpus);
196 for(i=0; i<*num_cpus; i++) {
197 ithread[i]=i;
198 pthread_create(&tid[i], NULL, (void *)extrapol_kel3mt, (void *)&ithread[i]);
199 }
200 for(i=0; i<*num_cpus; i++) pthread_join(tid[i], NULL);
201
202 SFREE(ithread);
203
204 /* correct the facial turbulent kinetic energy gradients:
205 Moukalled et al. p 289 */
206
207 ielfa1=ielfa;ipnei1=ipnei;vel1=vel;xlet1=xlet;gradkfa1=gradkfa;xxj1=xxj;
208 nef1=nef;nface1=nface;ncfd1=ncfd;
209
210 /* create threads and wait */
211
212 NNEW(ithread,ITG,*num_cpus);
213 for(i=0; i<*num_cpus; i++) {
214 ithread[i]=i;
215 pthread_create(&tid[i], NULL, (void *)extrapol_kel5mt, (void *)&ithread[i]);
216 }
217 for(i=0; i<*num_cpus; i++) pthread_join(tid[i], NULL);
218
219 SFREE(ithread);
220
221
222 /* inter/extrapolation of k at the center of the elements
223 to the center of the faces: taking the skewness of the
224 elements into account (using the gradient at the center
225 of the faces)
226 Moukalled et al. p 279 */
227
228 ielfa1=ielfa;vfa1=vfa;vfap1=vfap;gradkfa1=gradkfa;rf1=rf;ifabou1=ifabou;
229 ipnei1=ipnei;vel1=vel;xxi1=xxi;xle1=xle;nef1=nef;nface1=nface;inlet1=inlet;
230
231 /* create threads and wait */
232
233 NNEW(ithread,ITG,*num_cpus);
234 for(i=0; i<*num_cpus; i++) {
235 ithread[i]=i;
236 pthread_create(&tid[i], NULL, (void *)extrapol_kel4mt, (void *)&ithread[i]);
237 }
238 for(i=0; i<*num_cpus; i++) pthread_join(tid[i], NULL);
239
240 SFREE(ithread);
241
242 /* multiple point constraints */
243
244 if(*nmpc>0){
245 is=6;
246 ie=6;
247 FORTRAN(applympc,(nface,ielfa,&is,&ie,ifabou,ipompc,vfa,coefmpc,
248 nodempc,ipnei,neifa,labmpc,xbounact,nactdoh,
249 ifaext,nfaext));
250 }
251 }
252
253 SFREE(vfap);
254
255 return;
256
257 }
258
259 /* subroutine for multithreading of extrapol_kel1 */
260
extrapol_kel1mt(ITG * i)261 void *extrapol_kel1mt(ITG *i){
262
263 ITG nfacea,nfaceb,nfacedelta;
264
265 nfacedelta=(ITG)floor(*nface1/(double)(*num_cpus1));
266 nfacea=*i*nfacedelta+1;
267 nfaceb=(*i+1)*nfacedelta;
268 if((*i==*num_cpus1-1)&&(nfaceb<*nface1)) nfaceb=*nface1;
269
270 FORTRAN(extrapol_kel1,(ielfa1,xrlfa1,vfap1,vel1,
271 ifabou1,nef1,umfa1,&constant1,
272 inlet1,&nfacea,&nfaceb));
273
274 return NULL;
275 }
276
277 /* subroutine for multithreading of extrapol_kel2 */
278
extrapol_kel2mt(ITG * i)279 void *extrapol_kel2mt(ITG *i){
280
281 ITG nefa,nefb,nefdelta;
282
283 nefdelta=(ITG)floor(*nef1/(double)(*num_cpus1));
284 nefa=*i*nefdelta+1;
285 nefb=(*i+1)*nefdelta;
286 if((*i==*num_cpus1-1)&&(nefb<*nef1)) nefb=*nef1;
287
288 FORTRAN(extrapol_kel2,(ipnei1,neifa1,vfa1,area1,xxna1,volume1,gradkel1,
289 &nefa,&nefb,ncfd1));
290
291 return NULL;
292 }
293
294 /* subroutine for multithreading of extrapol_kel3 */
295
extrapol_kel3mt(ITG * i)296 void *extrapol_kel3mt(ITG *i){
297
298 ITG nfacea,nfaceb,nfacedelta;
299
300 nfacedelta=(ITG)floor(*nface1/(double)(*num_cpus1));
301 nfacea=*i*nfacedelta+1;
302 nfaceb=(*i+1)*nfacedelta;
303 if((*i==*num_cpus1-1)&&(nfaceb<*nface1)) nfaceb=*nface1;
304
305 FORTRAN(extrapol_kel3,(ielfa1,xrlfa1,icyclic1,ifatie1,gradkfa1,
306 gradkel1,c1,ipnei1,xxi1,
307 &nfacea,&nfaceb,ncfd1));
308
309 return NULL;
310 }
311
312 /* subroutine for multithreading of extrapol_kel4 */
313
extrapol_kel4mt(ITG * i)314 void *extrapol_kel4mt(ITG *i){
315
316 ITG nfacea,nfaceb,nfacedelta;
317
318 nfacedelta=(ITG)floor(*nface1/(double)(*num_cpus1));
319 nfacea=*i*nfacedelta+1;
320 nfaceb=(*i+1)*nfacedelta;
321 if((*i==*num_cpus1-1)&&(nfaceb<*nface1)) nfaceb=*nface1;
322
323 FORTRAN(extrapol_kel4,(ielfa1,vfa1,vfap1,gradkfa1,rf1,ifabou1,
324 ipnei1,vel1,xxi1,xle1,nef1,inlet1,&nfacea,&nfaceb));
325
326 return NULL;
327 }
328
329 /* subroutine for multithreading of extrapol_kel5 */
330
extrapol_kel5mt(ITG * i)331 void *extrapol_kel5mt(ITG *i){
332
333 ITG nfacea,nfaceb,nfacedelta;
334
335 nfacedelta=(ITG)floor(*nface1/(double)(*num_cpus1));
336 nfacea=*i*nfacedelta+1;
337 nfaceb=(*i+1)*nfacedelta;
338 if((*i==*num_cpus1-1)&&(nfaceb<*nface1)) nfaceb=*nface1;
339
340 FORTRAN(extrapol_kel5,(ielfa1,ipnei1,vel1,xlet1,gradkfa1,xxj1,
341 nef1,&nfacea,&nfaceb,ncfd1));
342
343 return NULL;
344 }
345