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