/[escript]/trunk/finley/src/Assemble_LumpedSystem.c
ViewVC logotype

Annotation of /trunk/finley/src/Assemble_LumpedSystem.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 751 - (hide annotations)
Mon Jun 26 01:46:34 2006 UTC (12 years, 6 months ago) by bcumming
Original Path: trunk/finley/src/Assemble_PDE.c
File MIME type: text/plain
File size: 17136 byte(s)
Changes relating to the MPI version of escript
The standard OpenMP version of escript is unchanged

- updated data types (Finley_Mesh, Finley_NodeFile, etc) to store meshes
  over multiple MPI processes.
- added CommBuffer code in Paso for communication of Data associated
  with distributed meshes
- updates in Finley and Escript to support distributed data and operations
  on distributed data (such as interpolation).
- construction of RHS in MPI, so that simple explicit schemes (such as
  /docs/examples/wave.py without IO and the Locator) can run in MPI.
- updated mesh generation for first order line, rectangle and brick
  meshes and second order line meshes in MPI.        
- small changes to trunk/SConstruct and trunk/scons/ess_options.py to
  build the MPI version, these changes are turned off by default.

1 jgs 150 /*
2 elspeth 616 ************************************************************
3     * Copyright 2006 by ACcESS MNRF *
4     * *
5     * http://www.access.edu.au *
6     * Primary Business: Queensland, Australia *
7     * Licensed under the Open Software License version 3.0 *
8     * http://www.opensource.org/licenses/osl-3.0.php *
9     * *
10     ************************************************************
11 jgs 150 */
12 jgs 82
13 jgs 150
14 jgs 82 /**************************************************************/
15    
16     /* assembles the system of numEq PDEs into the stiffness matrix S and right hand side F */
17    
18     /* -div(A*grad u)-div(B*u)+C*grad u + D*u= -div X + Y */
19    
20     /* -(A_{k,i,m,j} u_m,j)_i-(B_{k,i,m} u_m)_i+C_{k,m,j} u_m,j-D_{k,m} u_m = -(X_{k,i})_i + Y_k */
21    
22     /* u has numComp components. */
23    
24     /* Shape of the coefficients: */
25    
26     /* A = numEqu x numDim x numComp x numDim */
27     /* B = numDim x numEqu x numComp */
28     /* C = numEqu x numDim x numComp */
29     /* D = numEqu x numComp */
30     /* X = numEqu x numDim */
31     /* Y = numEqu */
32    
33     /* The coefficients A,B,C,D,X and Y have to be defined on the integartion points or not present (=NULL). */
34    
35     /* S and F have to be initialized before the routine is called. S or F can be NULL. In this case the left or */
36     /* the right hand side of the PDE is not processed. */
37    
38     /* The routine does not consider any boundary conditions. */
39    
40     /**************************************************************/
41    
42 jgs 150 /* Author: gross@access.edu.au */
43     /* Version: $Id$ */
44 jgs 82
45     /**************************************************************/
46    
47     #include "Assemble.h"
48     #include "Util.h"
49     #ifdef _OPENMP
50     #include <omp.h>
51     #endif
52    
53    
54     /**************************************************************/
55    
56 jgs 150 void Finley_Assemble_PDE(Finley_NodeFile* nodes,Finley_ElementFile* elements,Paso_SystemMatrix* S, escriptDataC* F,
57 jgs 82 escriptDataC* A, escriptDataC* B, escriptDataC* C, escriptDataC* D, escriptDataC* X, escriptDataC* Y ) {
58    
59 jgs 150 char error_msg[LenErrorMsg_MAX];
60 jgs 82 double *EM_S=NULL,*EM_F=NULL,*V=NULL,*dVdv=NULL,*dSdV=NULL,*Vol=NULL,*dvdV=NULL;
61     double time0;
62 jgs 123 dim_t dimensions[ESCRIPT_MAX_DATA_RANK],e,q;
63 jgs 82 Assemble_Parameters p;
64 jgs 123 index_t *index_row=NULL,*index_col=NULL,color;
65 jgs 150 Finley_resetError();
66 jgs 82
67     if (nodes==NULL || elements==NULL) return;
68     if (S==NULL && isEmpty(F)) return;
69    
70     /* set all parameters in p*/
71     Assemble_getAssembleParameters(nodes,elements,S,F,&p);
72 jgs 150 if (! Finley_noError()) return;
73 jgs 82
74     /* this function assumes NS=NN */
75     if (p.NN!=p.NS) {
76 gross 532 Finley_setError(SYSTEM_ERROR,"Finley_Assemble_PDE: for Finley_Assemble_PDE numNodes and numShapes have to be identical.");
77 jgs 82 return;
78     }
79     if (p.numDim!=p.numElementDim) {
80 gross 532 Finley_setError(SYSTEM_ERROR,"Finley_Assemble_PDE: Finley_Assemble_PDE accepts volume elements only.");
81 jgs 82 return;
82     }
83     /* get a functionspace */
84 jgs 123 type_t funcspace=UNKNOWN;
85 jgs 82 updateFunctionSpaceType(funcspace,A);
86     updateFunctionSpaceType(funcspace,B);
87     updateFunctionSpaceType(funcspace,C);
88     updateFunctionSpaceType(funcspace,D);
89     updateFunctionSpaceType(funcspace,X);
90     updateFunctionSpaceType(funcspace,Y);
91     if (funcspace==UNKNOWN) return; /* all data are empty */
92    
93     /* check if all function spaces are the same */
94    
95     if (! functionSpaceTypeEqual(funcspace,A) ) {
96 gross 532 Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: unexpected function space type for coefficient A");
97 jgs 82 }
98     if (! functionSpaceTypeEqual(funcspace,B) ) {
99 gross 532 Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: unexpected function space type for coefficient B");
100 jgs 82 }
101     if (! functionSpaceTypeEqual(funcspace,C) ) {
102 gross 532 Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: unexpected function space type for coefficient C");
103 jgs 82 }
104     if (! functionSpaceTypeEqual(funcspace,D) ) {
105 gross 532 Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: unexpected function space type for coefficient D");
106 jgs 82 }
107     if (! functionSpaceTypeEqual(funcspace,X) ) {
108 gross 532 Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: unexpected function space type for coefficient X");
109 jgs 82 }
110     if (! functionSpaceTypeEqual(funcspace,Y) ) {
111 gross 532 Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: unexpected function space type for coefficient Y");
112 jgs 82 }
113    
114     /* check if all function spaces are the same */
115    
116     if (! numSamplesEqual(A,p.numQuad,elements->numElements) ) {
117 gross 532 sprintf(error_msg,"Finley_Assemble_PDE: sample points of coefficient A don't match (%d,%d)",p.numQuad,elements->numElements);
118 jgs 150 Finley_setError(TYPE_ERROR,error_msg);
119 jgs 82 }
120    
121     if (! numSamplesEqual(B,p.numQuad,elements->numElements) ) {
122 gross 532 sprintf(error_msg,"Finley_Assemble_PDE: sample points of coefficient B don't match (%d,%d)",p.numQuad,elements->numElements);
123 jgs 150 Finley_setError(TYPE_ERROR,error_msg);
124 jgs 82 }
125    
126     if (! numSamplesEqual(C,p.numQuad,elements->numElements) ) {
127 gross 532 sprintf(error_msg,"Finley_Assemble_PDE: sample points of coefficient C don't match (%d,%d)",p.numQuad,elements->numElements);
128 jgs 150 Finley_setError(TYPE_ERROR,error_msg);
129 jgs 82 }
130    
131     if (! numSamplesEqual(D,p.numQuad,elements->numElements) ) {
132 gross 532 sprintf(error_msg,"Finley_Assemble_PDE: sample points of coefficient D don't match (%d,%d)",p.numQuad,elements->numElements);
133 jgs 150 Finley_setError(TYPE_ERROR,error_msg);
134 jgs 82 }
135    
136     if (! numSamplesEqual(X,p.numQuad,elements->numElements) ) {
137 gross 532 sprintf(error_msg,"Finley_Assemble_PDE: sample points of coefficient X don't match (%d,%d)",p.numQuad,elements->numElements);
138 jgs 150 Finley_setError(TYPE_ERROR,error_msg);
139 jgs 82 }
140    
141     if (! numSamplesEqual(Y,p.numQuad,elements->numElements) ) {
142 gross 532 sprintf(error_msg,"Finley_Assemble_PDE: sample points of coefficient Y don't match (%d,%d)",p.numQuad,elements->numElements);
143 jgs 150 Finley_setError(TYPE_ERROR,error_msg);
144 jgs 82 }
145    
146     /* check the dimensions: */
147    
148     if (p.numEqu==1 && p.numComp==1) {
149     if (!isEmpty(A)) {
150     dimensions[0]=p.numDim;
151     dimensions[1]=p.numDim;
152     if (!isDataPointShapeEqual(A,2,dimensions)) {
153 gross 532 sprintf(error_msg,"Finley_Assemble_PDE: coefficient A: illegal shape, expected shape (%d,%d)",dimensions[0],dimensions[1]);
154 jgs 150 Finley_setError(TYPE_ERROR,error_msg);
155 jgs 82 }
156     }
157     if (!isEmpty(B)) {
158     dimensions[0]=p.numDim;
159     if (!isDataPointShapeEqual(B,1,dimensions)) {
160 gross 532 sprintf(error_msg,"Finley_Assemble_PDE: coefficient B: illegal shape (%d,)",dimensions[0]);
161 jgs 150 Finley_setError(TYPE_ERROR,error_msg);
162 jgs 82 }
163     }
164     if (!isEmpty(C)) {
165     dimensions[0]=p.numDim;
166     if (!isDataPointShapeEqual(C,1,dimensions)) {
167 gross 532 sprintf(error_msg,"Finley_Assemble_PDE: coefficient C, expected shape (%d,)",dimensions[0]);
168 jgs 150 Finley_setError(TYPE_ERROR,error_msg);
169 jgs 82 }
170     }
171     if (!isEmpty(D)) {
172     if (!isDataPointShapeEqual(D,0,dimensions)) {
173 gross 532 Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: coefficient D, rank 0 expected.");
174 jgs 82 }
175     }
176     if (!isEmpty(X)) {
177     dimensions[0]=p.numDim;
178     if (!isDataPointShapeEqual(X,1,dimensions)) {
179 gross 532 sprintf(error_msg,"Finley_Assemble_PDE: coefficient X, expected shape (%d,",dimensions[0]);
180 jgs 150 Finley_setError(TYPE_ERROR,error_msg);
181 jgs 82 }
182     }
183     if (!isEmpty(Y)) {
184     if (!isDataPointShapeEqual(Y,0,dimensions)) {
185 gross 532 Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: coefficient Y, rank 0 expected.");
186 jgs 82 }
187     }
188     } else {
189     if (!isEmpty(A)) {
190     dimensions[0]=p.numEqu;
191     dimensions[1]=p.numDim;
192     dimensions[2]=p.numComp;
193     dimensions[3]=p.numDim;
194     if (!isDataPointShapeEqual(A,4,dimensions)) {
195 gross 532 sprintf(error_msg,"Finley_Assemble_PDE: coefficient A, expected shape (%d,%d,%d,%d)",dimensions[0],dimensions[1],dimensions[2],dimensions[3]);
196 jgs 150 Finley_setError(TYPE_ERROR,error_msg);
197 jgs 82 }
198     }
199     if (!isEmpty(B)) {
200     dimensions[0]=p.numEqu;
201     dimensions[1]=p.numDim;
202     dimensions[2]=p.numComp;
203     if (!isDataPointShapeEqual(B,3,dimensions)) {
204 gross 532 sprintf(error_msg,"Finley_Assemble_PDE: coefficient B, expected shape (%d,%d,%d)",dimensions[0],dimensions[1],dimensions[2]);
205 jgs 150 Finley_setError(TYPE_ERROR,error_msg);
206 jgs 82 }
207     }
208     if (!isEmpty(C)) {
209     dimensions[0]=p.numEqu;
210     dimensions[1]=p.numComp;
211     dimensions[2]=p.numDim;
212     if (!isDataPointShapeEqual(C,3,dimensions)) {
213 gross 532 sprintf(error_msg,"Finley_Assemble_PDE: coefficient C, expected shape (%d,%d,%d)",dimensions[0],dimensions[1],dimensions[2]);
214 jgs 150 Finley_setError(TYPE_ERROR,error_msg);
215 jgs 82 }
216     }
217     if (!isEmpty(D)) {
218     dimensions[0]=p.numEqu;
219     dimensions[1]=p.numComp;
220     if (!isDataPointShapeEqual(D,2,dimensions)) {
221 gross 532 sprintf(error_msg,"Finley_Assemble_PDE: coefficient D, expected shape (%d,%d)",dimensions[0],dimensions[1]);
222 jgs 150 Finley_setError(TYPE_ERROR,error_msg);
223 jgs 82 }
224     }
225     if (!isEmpty(X)) {
226     dimensions[0]=p.numEqu;
227     dimensions[1]=p.numDim;
228     if (!isDataPointShapeEqual(X,2,dimensions)) {
229 gross 532 sprintf(error_msg,"Finley_Assemble_PDE: coefficient X, expected shape (%d,%d)",dimensions[0],dimensions[1]);
230 jgs 150 Finley_setError(TYPE_ERROR,error_msg);
231 jgs 82 }
232     }
233     if (!isEmpty(Y)) {
234     dimensions[0]=p.numEqu;
235     if (!isDataPointShapeEqual(Y,1,dimensions)) {
236 gross 532 sprintf(error_msg,"Finley_Assemble_PDE: coefficient Y, expected shape (%d,)",dimensions[0]);
237 jgs 150 Finley_setError(TYPE_ERROR,error_msg);
238 jgs 82 }
239     }
240     }
241    
242 jgs 150 if (Finley_noError()) {
243 jgs 82 time0=Finley_timer();
244     #pragma omp parallel private(index_col,index_row,EM_S,EM_F,V,dVdv,dSdV,Vol,dvdV,color,q) \
245     firstprivate(elements,nodes,S,F,A,B,C,D,X,Y)
246     {
247     EM_S=EM_F=V=dVdv=dSdV=Vol=dvdV=NULL;
248     index_row=index_col=NULL;
249    
250     /* allocate work arrays: */
251    
252 jgs 102 EM_S=(double*) THREAD_MEMALLOC(p.NN_row*p.NN_col*p.numEqu*p.numComp,double);
253     EM_F=(double*) THREAD_MEMALLOC(p.NN_row*p.numEqu,double);
254     V=(double*) THREAD_MEMALLOC(p.NN*p.numDim,double);
255     dVdv=(double*) THREAD_MEMALLOC(p.numDim*p.numDim*p.numQuad,double);
256     dvdV=(double*) THREAD_MEMALLOC(p.numDim*p.numDim*p.numQuad,double);
257     dSdV=(double*) THREAD_MEMALLOC(p.NS_row*p.numQuad*p.numDim,double);
258     Vol=(double*) THREAD_MEMALLOC(p.numQuad,double);
259 jgs 123 index_col=(index_t*) THREAD_MEMALLOC(p.NN_col,index_t);
260     index_row=(index_t*) THREAD_MEMALLOC(p.NN_row,index_t);
261 jgs 82
262     if (! (Finley_checkPtr(EM_S) || Finley_checkPtr(EM_F) || Finley_checkPtr(V) || Finley_checkPtr(index_col) ||
263     Finley_checkPtr(index_row) || Finley_checkPtr(dVdv) || Finley_checkPtr(dSdV) || Finley_checkPtr(Vol) )) {
264    
265     /* open loop over all colors: */
266 bcumming 751 #ifndef PASO_MPI
267 jgs 123 for (color=elements->minColor;color<=elements->maxColor;color++) {
268 jgs 82 /* open loop over all elements: */
269     #pragma omp for private(e) schedule(static)
270     for(e=0;e<elements->numElements;e++){
271     if (elements->Color[e]==color) {
272 bcumming 751 #else
273     for(e=0;e<elements->numElements;e++){
274     #endif
275 gross 532 //============================
276 jgs 82 for (q=0;q<p.NN_row;q++) index_row[q]=p.label_row[elements->Nodes[INDEX2(p.row_node[q],e,p.NN)]];
277     /* gather V-coordinates of nodes into V: */
278     Finley_Util_Gather_double(p.NN,&(elements->Nodes[INDEX2(0,e,p.NN)]),p.numDim,nodes->Coordinates,V);
279     /* calculate dVdv(i,j,q)=V(i,k)*DSDv(k,j,q) */
280     Finley_Util_SmallMatMult(p.numDim,p.numDim*p.numQuad,dVdv,p.NS,V,p.referenceElement->dSdv);
281     /* dvdV=invert(dVdv) inplace */
282     Finley_Util_InvertSmallMat(p.numQuad,p.numDim,dVdv,dvdV,Vol);
283     /* calculate dSdV=DSDv*DvDV */
284     Finley_Util_SmallMatSetMult(p.numQuad,p.NS_row,p.numDim,dSdV,p.numDim,p.referenceElement_row->dSdv,dvdV);
285     /* scale volume: */
286     for (q=0;q<p.numQuad;q++) Vol[q]=ABS(Vol[q]*p.referenceElement->QuadWeights[q]);
287 gross 532 //============================
288 jgs 82
289     /* integration for the stiffness matrix: */
290     /* in order to optimze the number of operations the case of constants coefficience needs a bit more work */
291     /* to make use of some symmetry. */
292    
293     if (S!=NULL) {
294     for (q=0;q<p.NN_row*p.NN_col*p.numEqu*p.numComp;q++) EM_S[q]=0;
295     if (p.numEqu==1 && p.numComp==1) {
296     Finley_Assemble_PDEMatrix_Single2(p.NS_row,p.numDim,p.numQuad,
297     p.referenceElement_row->S,dSdV,Vol,p.NN_row,EM_S,
298     getSampleData(A,e),isExpanded(A),
299     getSampleData(B,e),isExpanded(B),
300     getSampleData(C,e),isExpanded(C),
301     getSampleData(D,e),isExpanded(D));
302     } else {
303     Finley_Assemble_PDEMatrix_System2(p.NS_row,p.numDim,p.numQuad,p.numEqu,p.numComp,
304     p.referenceElement_row->S,dSdV,Vol,p.NN_row,EM_S,
305     getSampleData(A,e),isExpanded(A),
306     getSampleData(B,e),isExpanded(B),
307     getSampleData(C,e),isExpanded(C),
308     getSampleData(D,e),isExpanded(D));
309     }
310     for (q=0;q<p.NN_col;q++) index_col[q]=p.label_col[elements->Nodes[INDEX2(p.col_node[q],e,p.NN)]];
311 jgs 150 Finley_Assemble_addToSystemMatrix(S,p.NN_row,index_row,p.numEqu,p.NN_col,index_col,p.numComp,EM_S);
312 jgs 82 }
313     if (!isEmpty(F)) {
314     for (q=0;q<p.NN_row*p.numEqu;q++) EM_F[q]=0;
315     if (p.numEqu==1) {
316     Finley_Assemble_RHSMatrix_Single(p.NS_row,p.numDim,p.numQuad,
317     p.referenceElement_row->S,dSdV,Vol,p.NN_row,EM_F,
318     getSampleData(X,e),isExpanded(X),
319     getSampleData(Y,e),isExpanded(Y));
320     } else {
321     Finley_Assemble_RHSMatrix_System(p.NS_row,p.numDim,p.numQuad,
322     p.numEqu,p.referenceElement_row->S,dSdV,Vol,p.NN_row,EM_F,
323     getSampleData(X,e),isExpanded(X),
324     getSampleData(Y,e),isExpanded(Y));
325     }
326     /* add */
327 bcumming 751 #ifndef PASO_MPI
328 jgs 82 Finley_Util_AddScatter(p.NN_row,index_row,p.numEqu,EM_F,getSampleData(F,0));
329 bcumming 751 #else
330     Finley_Util_AddScatter_upperBound(p.NN_row,index_row,p.numEqu,EM_F,getSampleData(F,0),p.degreeOfFreedomUpperBound);
331     #endif
332 jgs 82 }
333     }
334     }
335 bcumming 751 #ifndef PASO_MPI
336 jgs 82 }
337     }
338 bcumming 751 #endif
339 jgs 82 /* clean up */
340     THREAD_MEMFREE(EM_S);
341     THREAD_MEMFREE(EM_F);
342     THREAD_MEMFREE(V);
343     THREAD_MEMFREE(dVdv);
344     THREAD_MEMFREE(dvdV);
345     THREAD_MEMFREE(dSdV);
346     THREAD_MEMFREE(Vol);
347     THREAD_MEMFREE(index_col);
348     THREAD_MEMFREE(index_row);
349     }
350 jgs 149 #ifdef Finley_TRACE
351 jgs 82 printf("timing: assemblage PDE: %.4e sec\n",Finley_timer()-time0);
352 jgs 149 #endif
353 jgs 82 }
354     }
355     /*
356     * $Log$
357 jgs 150 * Revision 1.8 2005/09/15 03:44:21 jgs
358     * Merge of development branch dev-02 back to main trunk on 2005-09-15
359     *
360 jgs 149 * Revision 1.7 2005/09/01 03:31:35 jgs
361     * Merge of development branch dev-02 back to main trunk on 2005-09-01
362     *
363 jgs 147 * Revision 1.6 2005/08/12 01:45:42 jgs
364     * erge of development branch dev-02 back to main trunk on 2005-08-12
365     *
366 jgs 150 * Revision 1.5.2.3 2005/09/07 06:26:17 gross
367     * the solver from finley are put into the standalone package paso now
368     *
369 jgs 149 * Revision 1.5.2.2 2005/08/24 02:02:18 gross
370     * timing output switched off. solver output can be swiched through getSolution(verbose=True) now.
371     *
372 jgs 147 * Revision 1.5.2.1 2005/08/03 08:54:27 gross
373     * contact element assemblage was called with wrong element table pointer
374     *
375 jgs 123 * Revision 1.5 2005/07/08 04:07:46 jgs
376     * Merge of development branch back to main trunk on 2005-07-08
377     *
378 jgs 102 * Revision 1.4 2004/12/15 07:08:32 jgs
379 jgs 97 * *** empty log message ***
380 jgs 123 * Revision 1.1.1.1.2.2 2005/06/29 02:34:47 gross
381     * some changes towards 64 integers in finley
382 jgs 82 *
383 jgs 123 * Revision 1.1.1.1.2.1 2004/11/24 01:37:12 gross
384     * some changes dealing with the integer overflow in memory allocation. Finley solves 4M unknowns now
385 jgs 97 *
386 jgs 82 *
387 jgs 123 *
388 jgs 82 */

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision

  ViewVC Help
Powered by ViewVC 1.1.26