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 char *lakonf1;
26
27 static ITG num_cpus,*nef1,*ipnei1,*neifa1,*neiel1,*jq1,*irow1,*nzs1,*ielfa1,
28 *ifabou1,*nbody1,*nactdohinv1,*icyclic1,*ifatie1,*iau61,*iturbulent1,
29 *ncfd1,*inlet1;
30
31 static double *auv1,*adv1,*bv1,*vfa1,*xxn1,*area1,*vel1,
32 *cosa1,*umfa1,*alet1,*ale1,*gradvfa1,*xxi1,*body1,*volume1,*dtimef1,
33 *velo1,*veloo1,*sel1,*xrlfa1,*gamma1,*xxj1,*a11,*a21,*a31,*flux1,
34 *c1,*xxna1,*xxnj1,*gradvel1,*of21,*yy1,*umel1,*sc1;
35
mafillvmain(ITG * nef,ITG * ipnei,ITG * neifa,ITG * neiel,double * vfa,double * xxn,double * area,double * auv,double * adv,ITG * jq,ITG * irow,ITG * nzs,double * bv,double * vel,double * cosa,double * umfa,double * alet,double * ale,double * gradvfa,double * xxi,double * body,double * volume,ITG * ielfa,char * lakonf,ITG * ifabou,ITG * nbody,double * dtimef,double * velo,double * veloo,double * sel,double * xrlfa,double * gamma,double * xxj,ITG * nactdohinv,double * a1,double * a2,double * a3,double * flux,ITG * icyclic,double * c,ITG * ifatie,ITG * iau6,double * xxna,double * xxnj,ITG * iturbulent,double * gradvel,double * of2,double * yy,double * umel,ITG * ncfd,ITG * inlet,double * sc)36 void mafillvmain(ITG *nef,ITG *ipnei,ITG *neifa,ITG *neiel,
37 double *vfa,double *xxn,double *area,double *auv,double *adv,
38 ITG *jq,ITG *irow,ITG *nzs,double *bv,double *vel,double *cosa,
39 double *umfa,double *alet,double *ale,double *gradvfa,
40 double *xxi,double *body,double *volume,
41 ITG *ielfa,char *lakonf,ITG *ifabou,ITG *nbody,
42 double *dtimef,double *velo,double *veloo,
43 double *sel,double *xrlfa,double *gamma,double *xxj,
44 ITG *nactdohinv,double *a1,double *a2,double *a3,
45 double *flux,ITG *icyclic,double *c,ITG *ifatie,ITG *iau6,
46 double *xxna,double *xxnj,ITG *iturbulent,double *gradvel,
47 double *of2,double *yy,double *umel,ITG *ncfd,ITG *inlet,
48 double *sc){
49
50 ITG i,j;
51
52 /* variables for multithreading procedure */
53
54 ITG sys_cpus,*ithread=NULL;
55 char *env,*envloc,*envsys;
56
57 num_cpus = 0;
58 sys_cpus=0;
59
60 /* explicit user declaration prevails */
61
62 envsys=getenv("NUMBER_OF_CPUS");
63 if(envsys){
64 sys_cpus=atoi(envsys);
65 if(sys_cpus<0) sys_cpus=0;
66 }
67
68 /* automatic detection of available number of processors */
69
70 if(sys_cpus==0){
71 sys_cpus = getSystemCPUs();
72 if(sys_cpus<1) sys_cpus=1;
73 }
74
75 /* local declaration prevails, if strictly positive */
76
77 envloc = getenv("CCX_NPROC_CFD");
78 if(envloc){
79 num_cpus=atoi(envloc);
80 if(num_cpus<0){
81 num_cpus=0;
82 }else if(num_cpus>sys_cpus){
83 num_cpus=sys_cpus;
84 }
85
86 }
87
88 /* else global declaration, if any, applies */
89
90 env = getenv("OMP_NUM_THREADS");
91 if(num_cpus==0){
92 if (env)
93 num_cpus = atoi(env);
94 if (num_cpus < 1) {
95 num_cpus=1;
96 }else if(num_cpus>sys_cpus){
97 num_cpus=sys_cpus;
98 }
99 }
100
101 // next line is to be inserted in a similar way for all other paralell parts
102
103 if(*nef<num_cpus) num_cpus=*nef;
104
105 pthread_t tid[num_cpus];
106
107 /* calculating the stiffness and/or mass matrix
108 (symmetric part) */
109
110 nef1=nef;ipnei1=ipnei;neifa1=neifa;neiel1=neiel;vfa1=vfa;xxn1=xxn;
111 area1=area;jq1=jq;irow1=irow;nzs1=nzs;vel1=vel;cosa1=cosa;umfa1=umfa;
112 alet1=alet;ale1=ale;gradvfa1=gradvfa;xxi1=xxi;body1=body;volume1=volume;
113 ielfa1=ielfa;lakonf1=lakonf;ifabou1=ifabou;nbody1=nbody;
114 dtimef1=dtimef;velo1=velo;veloo1=veloo;sel1=sel;xrlfa1=xrlfa;
115 gamma1=gamma;xxj1=xxj;nactdohinv1=nactdohinv;a11=a1;a21=a2;a31=a3;
116 flux1=flux;icyclic1=icyclic;c1=c;ifatie1=ifatie;iau61=iau6;
117 adv1=adv;auv1=auv;bv1=bv;xxna1=xxna;xxnj1=xxnj;iturbulent1=iturbulent;
118 gradvel1=gradvel;of21=of2;yy1=yy;umel1=umel;ncfd1=ncfd;inlet1=inlet;
119 sc1=sc;
120
121 /* create threads and wait */
122
123 NNEW(ithread,ITG,num_cpus);
124 for(i=0; i<num_cpus; i++) {
125 ithread[i]=i;
126 pthread_create(&tid[i], NULL, (void *)mafillvmt, (void *)&ithread[i]);
127 }
128 for(i=0; i<num_cpus; i++) pthread_join(tid[i], NULL);
129
130 SFREE(ithread);
131
132 return;
133
134 }
135
136 /* subroutine for multithreading of mafillv */
137
mafillvmt(ITG * i)138 void *mafillvmt(ITG *i){
139
140 ITG nefa,nefb,nefdelta;
141
142 // ceil -> floor
143
144 nefdelta=(ITG)floor(*nef1/(double)num_cpus);
145 nefa=*i*nefdelta+1;
146 nefb=(*i+1)*nefdelta;
147 // next line! -> all parallel sections
148 if((*i==num_cpus-1)&&(nefb<*nef1)) nefb=*nef1;
149
150 FORTRAN(mafillv,(nef1,ipnei1,neifa1,neiel1,vfa1,xxn1,area1,
151 auv1,adv1,jq1,irow1,nzs1,bv1,
152 vel1,cosa1,umfa1,alet1,ale1,gradvfa1,xxi1,
153 body1,volume1,ielfa1,lakonf1,ifabou1,nbody1,
154 dtimef1,velo1,veloo1,sel1,xrlfa1,gamma1,xxj1,nactdohinv1,
155 a11,a21,a31,flux1,&nefa,&nefb,icyclic1,c1,ifatie1,iau61,
156 xxna1,xxnj1,iturbulent1,gradvel1,of21,yy1,umel1,
157 ncfd1,inlet1,sc1));
158
159 return NULL;
160 }
161