/[escript]/trunk-mpi-branch/finley/src/Assemble_PDE.c
ViewVC logotype

Contents of /trunk-mpi-branch/finley/src/Assemble_PDE.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1196 - (show annotations)
Fri Jun 15 03:45:48 2007 UTC (11 years, 10 months ago) by ksteube
File MIME type: text/plain
File size: 17363 byte(s)
Use of PAPI on solver is now enabled with papi_instrument_solver=1 in scons/ess_options.py.
Can instrument other blocks of code with blockpapi.c.
Added interval timers to grad, integrate and Assemble_PDE.

1 /*
2 ************************************************************
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 */
12
13
14 /**************************************************************/
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 /* Author: gross@access.edu.au */
43 /* Version: $Id$ */
44
45 /**************************************************************/
46
47 #include "Assemble.h"
48 #include "Util.h"
49 #include "escript/blocktimer.h"
50 #ifdef _OPENMP
51 #include <omp.h>
52 #endif
53
54
55 /**************************************************************/
56
57 void Finley_Assemble_PDE(Finley_NodeFile* nodes,Finley_ElementFile* elements,Paso_SystemMatrix* S, escriptDataC* F,
58 escriptDataC* A, escriptDataC* B, escriptDataC* C, escriptDataC* D, escriptDataC* X, escriptDataC* Y ) {
59
60 bool_t reducedIntegrationOrder=FALSE;
61 char error_msg[LenErrorMsg_MAX];
62 Assemble_Parameters p;
63 double time0;
64 dim_t dimensions[ESCRIPT_MAX_DATA_RANK];
65 type_t funcspace;
66 double blocktimer_start = blocktimer_time();
67
68 Finley_resetError();
69
70 {
71 int iam, numCPUs;
72 #ifdef Paso_MPI
73 iam = elements->elementDistribution->MPIInfo->rank;
74 numCPUs = elements->elementDistribution->MPIInfo->size;
75 #endif
76 }
77
78 if (nodes==NULL || elements==NULL) return;
79 if (S==NULL && isEmpty(F)) return;
80
81 if (isEmpty(F) && !isEmpty(X) && !isEmpty(F)) {
82 Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: right hand side coefficients are non-zero bat no right hand side vector given.");
83 }
84
85 if (S==NULL && !isEmpty(A) && !isEmpty(B) && !isEmpty(C) && !isEmpty(D)) {
86 Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: coefficients are non-zero but no matrix is given.");
87 }
88
89 /* get the functionspace for this assemblage call */
90 funcspace=UNKNOWN;
91 updateFunctionSpaceType(funcspace,A);
92 updateFunctionSpaceType(funcspace,B);
93 updateFunctionSpaceType(funcspace,C);
94 updateFunctionSpaceType(funcspace,D);
95 updateFunctionSpaceType(funcspace,X);
96 updateFunctionSpaceType(funcspace,Y);
97 if (funcspace==UNKNOWN) return; /* all data are empty */
98
99 /* check if all function spaces are the same */
100 if (! functionSpaceTypeEqual(funcspace,A) ) {
101 Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: unexpected function space type for coefficient A");
102 }
103 if (! functionSpaceTypeEqual(funcspace,B) ) {
104 Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: unexpected function space type for coefficient B");
105 }
106 if (! functionSpaceTypeEqual(funcspace,C) ) {
107 Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: unexpected function space type for coefficient C");
108 }
109 if (! functionSpaceTypeEqual(funcspace,D) ) {
110 Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: unexpected function space type for coefficient D");
111 }
112 if (! functionSpaceTypeEqual(funcspace,X) ) {
113 Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: unexpected function space type for coefficient X");
114 }
115 if (! functionSpaceTypeEqual(funcspace,Y) ) {
116 Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: unexpected function space type for coefficient Y");
117 }
118 if (! Finley_noError()) return;
119
120 /* check if all function spaces are the same */
121 if (funcspace==FINLEY_ELEMENTS) {
122 reducedIntegrationOrder=FALSE;
123 } else if (funcspace==FINLEY_FACE_ELEMENTS) {
124 reducedIntegrationOrder=FALSE;
125 } else if (funcspace==FINLEY_CONTACT_ELEMENTS_1) {
126 reducedIntegrationOrder=FALSE;
127 } else if (funcspace==FINLEY_CONTACT_ELEMENTS_2) {
128 reducedIntegrationOrder=FALSE;
129 } else if (funcspace==FINLEY_REDUCED_ELEMENTS) {
130 reducedIntegrationOrder=TRUE;
131 } else if (funcspace==FINLEY_REDUCED_FACE_ELEMENTS) {
132 reducedIntegrationOrder=TRUE;
133 } else if (funcspace==FINLEY_REDUCED_CONTACT_ELEMENTS_1) {
134 reducedIntegrationOrder=TRUE;
135 } else if (funcspace==FINLEY_REDUCED_CONTACT_ELEMENTS_2) {
136 reducedIntegrationOrder=TRUE;
137 } else {
138 Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: assemblage failed because of illegal function space.");
139 }
140 if (! Finley_noError()) return;
141
142 /* set all parameters in p*/
143 Assemble_getAssembleParameters(nodes,elements,S,F, reducedIntegrationOrder, &p);
144 if (! Finley_noError()) return;
145
146 /* check if all function spaces are the same */
147
148 if (! numSamplesEqual(A,p.numQuad,elements->numElements) ) {
149 sprintf(error_msg,"Finley_Assemble_PDE: sample points of coefficient A don't match (%d,%d)",p.numQuad,elements->numElements);
150 Finley_setError(TYPE_ERROR,error_msg);
151 }
152
153 if (! numSamplesEqual(B,p.numQuad,elements->numElements) ) {
154 sprintf(error_msg,"Finley_Assemble_PDE: sample points of coefficient B don't match (%d,%d)",p.numQuad,elements->numElements);
155 Finley_setError(TYPE_ERROR,error_msg);
156 }
157
158 if (! numSamplesEqual(C,p.numQuad,elements->numElements) ) {
159 sprintf(error_msg,"Finley_Assemble_PDE: sample points of coefficient C don't match (%d,%d)",p.numQuad,elements->numElements);
160 Finley_setError(TYPE_ERROR,error_msg);
161 }
162
163 if (! numSamplesEqual(D,p.numQuad,elements->numElements) ) {
164 sprintf(error_msg,"Finley_Assemble_PDE: sample points of coefficient D don't match (%d,%d)",p.numQuad,elements->numElements);
165 Finley_setError(TYPE_ERROR,error_msg);
166 }
167
168 if (! numSamplesEqual(X,p.numQuad,elements->numElements) ) {
169 sprintf(error_msg,"Finley_Assemble_PDE: sample points of coefficient X don't match (%d,%d)",p.numQuad,elements->numElements);
170 Finley_setError(TYPE_ERROR,error_msg);
171 }
172
173 if (! numSamplesEqual(Y,p.numQuad,elements->numElements) ) {
174 sprintf(error_msg,"Finley_Assemble_PDE: sample points of coefficient Y don't match (%d,%d)",p.numQuad,elements->numElements);
175 Finley_setError(TYPE_ERROR,error_msg);
176 }
177
178 /* check the dimensions: */
179
180 if (p.numEqu==1 && p.numComp==1) {
181 if (!isEmpty(A)) {
182 dimensions[0]=p.numDim;
183 dimensions[1]=p.numDim;
184 if (!isDataPointShapeEqual(A,2,dimensions)) {
185 sprintf(error_msg,"Finley_Assemble_PDE: coefficient A: illegal shape, expected shape (%d,%d)",dimensions[0],dimensions[1]);
186 Finley_setError(TYPE_ERROR,error_msg);
187 }
188 }
189 if (!isEmpty(B)) {
190 dimensions[0]=p.numDim;
191 if (!isDataPointShapeEqual(B,1,dimensions)) {
192 sprintf(error_msg,"Finley_Assemble_PDE: coefficient B: illegal shape (%d,)",dimensions[0]);
193 Finley_setError(TYPE_ERROR,error_msg);
194 }
195 }
196 if (!isEmpty(C)) {
197 dimensions[0]=p.numDim;
198 if (!isDataPointShapeEqual(C,1,dimensions)) {
199 sprintf(error_msg,"Finley_Assemble_PDE: coefficient C, expected shape (%d,)",dimensions[0]);
200 Finley_setError(TYPE_ERROR,error_msg);
201 }
202 }
203 if (!isEmpty(D)) {
204 if (!isDataPointShapeEqual(D,0,dimensions)) {
205 Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: coefficient D, rank 0 expected.");
206 }
207 }
208 if (!isEmpty(X)) {
209 dimensions[0]=p.numDim;
210 if (!isDataPointShapeEqual(X,1,dimensions)) {
211 sprintf(error_msg,"Finley_Assemble_PDE: coefficient X, expected shape (%d,",dimensions[0]);
212 Finley_setError(TYPE_ERROR,error_msg);
213 }
214 }
215 if (!isEmpty(Y)) {
216 if (!isDataPointShapeEqual(Y,0,dimensions)) {
217 Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: coefficient Y, rank 0 expected.");
218 }
219 }
220 } else {
221 if (!isEmpty(A)) {
222 dimensions[0]=p.numEqu;
223 dimensions[1]=p.numDim;
224 dimensions[2]=p.numComp;
225 dimensions[3]=p.numDim;
226 if (!isDataPointShapeEqual(A,4,dimensions)) {
227 sprintf(error_msg,"Finley_Assemble_PDE: coefficient A, expected shape (%d,%d,%d,%d)",dimensions[0],dimensions[1],dimensions[2],dimensions[3]);
228 Finley_setError(TYPE_ERROR,error_msg);
229 }
230 }
231 if (!isEmpty(B)) {
232 dimensions[0]=p.numEqu;
233 dimensions[1]=p.numDim;
234 dimensions[2]=p.numComp;
235 if (!isDataPointShapeEqual(B,3,dimensions)) {
236 sprintf(error_msg,"Finley_Assemble_PDE: coefficient B, expected shape (%d,%d,%d)",dimensions[0],dimensions[1],dimensions[2]);
237 Finley_setError(TYPE_ERROR,error_msg);
238 }
239 }
240 if (!isEmpty(C)) {
241 dimensions[0]=p.numEqu;
242 dimensions[1]=p.numComp;
243 dimensions[2]=p.numDim;
244 if (!isDataPointShapeEqual(C,3,dimensions)) {
245 sprintf(error_msg,"Finley_Assemble_PDE: coefficient C, expected shape (%d,%d,%d)",dimensions[0],dimensions[1],dimensions[2]);
246 Finley_setError(TYPE_ERROR,error_msg);
247 }
248 }
249 if (!isEmpty(D)) {
250 dimensions[0]=p.numEqu;
251 dimensions[1]=p.numComp;
252 if (!isDataPointShapeEqual(D,2,dimensions)) {
253 sprintf(error_msg,"Finley_Assemble_PDE: coefficient D, expected shape (%d,%d)",dimensions[0],dimensions[1]);
254 Finley_setError(TYPE_ERROR,error_msg);
255 }
256 }
257 if (!isEmpty(X)) {
258 dimensions[0]=p.numEqu;
259 dimensions[1]=p.numDim;
260 if (!isDataPointShapeEqual(X,2,dimensions)) {
261 sprintf(error_msg,"Finley_Assemble_PDE: coefficient X, expected shape (%d,%d)",dimensions[0],dimensions[1]);
262 Finley_setError(TYPE_ERROR,error_msg);
263 }
264 }
265 if (!isEmpty(Y)) {
266 dimensions[0]=p.numEqu;
267 if (!isDataPointShapeEqual(Y,1,dimensions)) {
268 sprintf(error_msg,"Finley_Assemble_PDE: coefficient Y, expected shape (%d,)",dimensions[0]);
269 Finley_setError(TYPE_ERROR,error_msg);
270 }
271 }
272 }
273 if (Finley_noError()) {
274 time0=Finley_timer();
275 if (p.numEqu == p. numComp) {
276 if (p.numEqu > 1) {
277 /* system of PDESs */
278 if (p.numDim==3) {
279 if (p.row_NS == p.col_NS && p.row_NS == p.row_NN && p.col_NS == p.col_NN ) {
280 Finley_Assemble_PDE_System2_3D(p,elements,S,F,A,B,C,D,X,Y);
281 } else if ( p.row_NS == p.col_NS && 2*p.row_NS == p.row_NN && 2*p.col_NS == p.col_NN ) {
282 if ( !isEmpty(A) || !isEmpty(B) || !isEmpty(C) || !isEmpty(X) ) {
283 Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: Contact elements require A, B, C and X to be empty.");
284 } else {
285 Finley_Assemble_PDE_System2_C(p,elements,S,F,D,Y);
286 }
287 } else {
288 Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE supports numShape=NumNodes or 2*numShape=NumNodes only.");
289 }
290 } else if (p.numDim==2) {
291 if ((p.row_NS == p.col_NS) && (p.row_NS == p.row_NN) && (p.col_NS == p.col_NN )) {
292 Finley_Assemble_PDE_System2_2D(p,elements,S,F,A,B,C,D,X,Y);
293 } else if ( p.row_NS == p.col_NS && 2*p.row_NS == p.row_NN && 2*p.col_NS == p.col_NN ) {
294 if ( !isEmpty(A) || !isEmpty(B) || !isEmpty(C) || !isEmpty(X) ) {
295 Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: Contact elements require A, B, C and X to be empty.");
296 } else {
297 Finley_Assemble_PDE_System2_C(p,elements,S,F,D,Y);
298 }
299 } else {
300 Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE supports numShape=NumNodes or 2*numShape=NumNodes only.");
301 }
302 } else if (p.numDim==2) {
303 if (p.row_NS == p.col_NS && p.row_NS == p.row_NN && p.col_NS == p.col_NN ) {
304 Finley_Assemble_PDE_System2_1D(p,elements,S,F,A,B,C,D,X,Y);
305 } else if ( p.row_NS == p.col_NS && 2*p.row_NS == p.row_NN && 2*p.col_NS == p.col_NN ) {
306 if ( !isEmpty(A) || !isEmpty(B) || !isEmpty(C) || !isEmpty(X) ) {
307 Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: Contact elements require A, B, C and X to be empty.");
308 } else {
309 Finley_Assemble_PDE_System2_C(p,elements,S,F,D,Y);
310 }
311 } else {
312 Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE supports numShape=NumNodes or 2*numShape=NumNodes only.");
313 }
314 } else {
315 Finley_setError(VALUE_ERROR,"Finley_Assemble_PDE supports spatial dimensions 1,2,3 only.");
316 }
317 } else {
318 /* single PDES */
319 if (p.numDim==3) {
320 if (p.row_NS == p.col_NS && p.row_NS == p.row_NN && p.col_NS == p.col_NN ) {
321 Finley_Assemble_PDE_Single2_3D(p,elements,S,F,A,B,C,D,X,Y);
322 } else if ( p.row_NS == p.col_NS && 2*p.row_NS == p.row_NN && 2*p.col_NS == p.col_NN ) {
323 if ( !isEmpty(A) || !isEmpty(B) || !isEmpty(C) || !isEmpty(X) ) {
324 Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: Contact elements require A, B, C and X to be empty.");
325 } else {
326 Finley_Assemble_PDE_Single2_C(p,elements,S,F,D,Y);
327 }
328 } else {
329 Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE supports numShape=NumNodes or 2*numShape=NumNodes only.");
330 }
331 } else if (p.numDim==2) {
332 if ((p.row_NS == p.col_NS) && (p.row_NS == p.row_NN) && (p.col_NS == p.col_NN )) {
333 Finley_Assemble_PDE_Single2_2D(p,elements,S,F,A,B,C,D,X,Y);
334 } else if ( p.row_NS == p.col_NS && 2*p.row_NS == p.row_NN && 2*p.col_NS == p.col_NN ) {
335 if ( !isEmpty(A) || !isEmpty(B) || !isEmpty(C) || !isEmpty(X) ) {
336 Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: Contact elements require A, B, C and X to be empty.");
337 } else {
338 Finley_Assemble_PDE_Single2_C(p,elements,S,F,D,Y);
339 }
340 } else {
341 Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE supports numShape=NumNodes or 2*numShape=NumNodes only.");
342 }
343 } else if (p.numDim==2) {
344 if (p.row_NS == p.col_NS && p.row_NS == p.row_NN && p.col_NS == p.col_NN ) {
345 Finley_Assemble_PDE_Single2_1D(p,elements,S,F,A,B,C,D,X,Y);
346 } else if ( p.row_NS == p.col_NS && 2*p.row_NS == p.row_NN && 2*p.col_NS == p.col_NN ) {
347 if ( !isEmpty(A) || !isEmpty(B) || !isEmpty(C) || !isEmpty(X) ) {
348 Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: Contact elements require A, B, C and X to be empty.");
349 } else {
350 Finley_Assemble_PDE_Single2_C(p,elements,S,F,D,Y);
351 }
352 } else {
353 Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE supports numShape=NumNodes or 2*numShape=NumNodes only.");
354 }
355 } else {
356 Finley_setError(VALUE_ERROR,"Finley_Assemble_PDE supports spatial dimensions 1,2,3 only.");
357 }
358 }
359 } else {
360 Finley_setError(VALUE_ERROR,"Finley_Assemble_PDE requires number of equations == number of solutions .");
361 }
362 #ifdef Finley_TRACE
363 printf("timing: assemblage PDE: %.4e sec\n",Finley_timer()-time0);
364 #endif
365 }
366 blocktimer_increment("Finley_Assemble_PDE()", blocktimer_start);
367 }
368 /*
369 * $Log$
370 * Revision 1.8 2005/09/15 03:44:21 jgs
371 * Merge of development branch dev-02 back to main trunk on 2005-09-15
372 *
373 * Revision 1.7 2005/09/01 03:31:35 jgs
374 * Merge of development branch dev-02 back to main trunk on 2005-09-01
375 *
376 * Revision 1.6 2005/08/12 01:45:42 jgs
377 * erge of development branch dev-02 back to main trunk on 2005-08-12
378 *
379 * Revision 1.5.2.3 2005/09/07 06:26:17 gross
380 * the solver from finley are put into the standalone package paso now
381 *
382 * Revision 1.5.2.2 2005/08/24 02:02:18 gross
383 * timing output switched off. solver output can be swiched through getSolution(verbose=True) now.
384 *
385 * Revision 1.5.2.1 2005/08/03 08:54:27 gross
386 * contact element assemblage was called with wrong element table pointer
387 *
388 * Revision 1.5 2005/07/08 04:07:46 jgs
389 * Merge of development branch back to main trunk on 2005-07-08
390 *
391 * Revision 1.4 2004/12/15 07:08:32 jgs
392 * *** empty log message ***
393 * Revision 1.1.1.1.2.2 2005/06/29 02:34:47 gross
394 * some changes towards 64 integers in finley
395 *
396 * Revision 1.1.1.1.2.1 2004/11/24 01:37:12 gross
397 * some changes dealing with the integer overflow in memory allocation. Finley solves 4M unknowns now
398 *
399 *
400 *
401 */

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26