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

Contents of /trunk/finley/src/Assemble_PDE_RHS.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 751 - (show annotations)
Mon Jun 26 01:46:34 2006 UTC (13 years, 3 months ago) by bcumming
File MIME type: text/plain
File size: 10507 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 /*
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 /* assembles the system of numEq PDEs into side F */
16
17 /* -div X + Y */
18
19 /* Shape of the coefficients: */
20
21 /* X = numEqu x numDim */
22 /* Y = numEqu */
23
24 /* The coefficients X and Y have to be defined on the integartion points or not present (=NULL). */
25
26 /* F has to be initialized before the routine is called. F can be NULL. */
27
28 /* The routine does not consider any boundary conditions. */
29
30 /**************************************************************/
31
32 /* Author: gross@access.edu.au */
33 /* Version: $Id$ */
34
35 /**************************************************************/
36
37 #include "Assemble.h"
38 #include "Util.h"
39 #ifdef _OPENMP
40 #include <omp.h>
41 #endif
42
43
44 /**************************************************************/
45
46 void Finley_Assemble_PDE_RHS(Finley_NodeFile* nodes,Finley_ElementFile* elements,escriptDataC* F,escriptDataC* X, escriptDataC* Y ) {
47
48 char error_msg[LenErrorMsg_MAX];
49 double *EM_F=NULL,*V=NULL,*dVdv=NULL,*dSdV=NULL,*Vol=NULL,*dvdV=NULL;
50 double time0;
51 dim_t dimensions[ESCRIPT_MAX_DATA_RANK],e,q;
52 Assemble_Parameters p;
53 index_t *index_row=NULL,color;
54 bool_t assemble;
55 Finley_resetError();
56
57 if (nodes==NULL || elements==NULL) return;
58 if (isEmpty(F) || (isEmpty(Y) && isEmpty(X)) ) return;
59
60 /* set all parameters in p*/
61 Assemble_getAssembleParameters(nodes,elements,NULL,F,&p);
62 if (! Finley_noError()) return;
63
64 /* this function assumes NS=NN */
65 if (p.NN!=p.NS) {
66 Finley_setError(SYSTEM_ERROR,"__FILE__: for Finley_Assemble_PDE numNodes and numShapes have to be identical.");
67 return;
68 }
69 if (p.numDim!=p.numElementDim) {
70 Finley_setError(SYSTEM_ERROR,"__FILE__: Finley_Assemble_PDE accepts volume elements only.");
71 return;
72 }
73 /* get a functionspace */
74 type_t funcspace=UNKNOWN;
75 updateFunctionSpaceType(funcspace,X);
76 updateFunctionSpaceType(funcspace,Y);
77 if (funcspace==UNKNOWN) return; /* all data are empty */
78
79 /* check if all function spaces are the same */
80
81 if (! functionSpaceTypeEqual(funcspace,X) ) {
82 Finley_setError(TYPE_ERROR,"__FILE__: unexpected function space type for coefficient X");
83 }
84 if (! functionSpaceTypeEqual(funcspace,Y) ) {
85 Finley_setError(TYPE_ERROR,"__FILE__: unexpected function space type for coefficient Y");
86 }
87
88 /* check if all function spaces are the same */
89
90 if (! numSamplesEqual(X,p.numQuad,elements->numElements) ) {
91 sprintf(error_msg,"__FILE__: sample points of coefficient X don't match (%d,%d)",p.numQuad,elements->numElements);
92 Finley_setError(TYPE_ERROR,error_msg);
93 }
94
95 if (! numSamplesEqual(Y,p.numQuad,elements->numElements) ) {
96 sprintf(error_msg,"__FILE__: sample points of coefficient Y don't match (%d,%d)",p.numQuad,elements->numElements);
97 Finley_setError(TYPE_ERROR,error_msg);
98 }
99
100 /* check the dimensions: */
101
102 if (p.numEqu==1 && p.numComp==1) {
103 if (!isEmpty(X)) {
104 dimensions[0]=p.numDim;
105 if (!isDataPointShapeEqual(X,1,dimensions)) {
106 sprintf(error_msg,"__FILE__: coefficient X, expected shape (%d,",dimensions[0]);
107 Finley_setError(TYPE_ERROR,error_msg);
108 }
109 }
110 if (!isEmpty(Y)) {
111 if (!isDataPointShapeEqual(Y,0,dimensions)) {
112 Finley_setError(TYPE_ERROR,"__FILE__: coefficient Y, rank 0 expected.");
113 }
114 }
115 } else {
116 if (!isEmpty(X)) {
117 dimensions[0]=p.numEqu;
118 dimensions[1]=p.numDim;
119 if (!isDataPointShapeEqual(X,2,dimensions)) {
120 sprintf(error_msg,"__FILE__: coefficient X, expected shape (%d,%d)",dimensions[0],dimensions[1]);
121 Finley_setError(TYPE_ERROR,error_msg);
122 }
123 }
124 if (!isEmpty(Y)) {
125 dimensions[0]=p.numEqu;
126 if (!isDataPointShapeEqual(Y,1,dimensions)) {
127 sprintf(error_msg,"__FILE__: coefficient Y, expected shape (%d,)",dimensions[0]);
128 Finley_setError(TYPE_ERROR,error_msg);
129 }
130 }
131 }
132
133 if (Finley_noError()) {
134 time0=Finley_timer();
135 #pragma omp parallel private(index_row,EM_F,V,dVdv,dSdV,Vol,dvdV,color,q,assemble)
136 {
137 EM_F=V=dVdv=dSdV=Vol=dvdV=NULL;
138 index_row=NULL;
139
140 /* allocate work arrays: */
141
142 EM_F=(double*) THREAD_MEMALLOC(p.NN_row*p.numEqu,double);
143 V=(double*) THREAD_MEMALLOC(p.NN*p.numDim,double);
144 dVdv=(double*) THREAD_MEMALLOC(p.numDim*p.numDim*p.numQuad,double);
145 dvdV=(double*) THREAD_MEMALLOC(p.numDim*p.numDim*p.numQuad,double);
146 dSdV=(double*) THREAD_MEMALLOC(p.NS_row*p.numQuad*p.numDim,double);
147 Vol=(double*) THREAD_MEMALLOC(p.numQuad,double);
148 index_row=(index_t*) THREAD_MEMALLOC(p.NN_row,index_t);
149
150 if (! (Finley_checkPtr(EM_F) || Finley_checkPtr(V) ||
151 Finley_checkPtr(index_row) || Finley_checkPtr(dVdv) || Finley_checkPtr(dSdV) || Finley_checkPtr(Vol) )) {
152
153 /* open loop over all colors: */
154 for (color=elements->minColor;color<=elements->maxColor;color++) {
155 /* open loop over all elements: */
156 #pragma omp for private(e) schedule(dynamic)
157 for(e=0;e<elements->numElements;e++){
158 if (elements->Color[e]==color) {
159 /* check if element has to be processed */
160 if (isEmpty(X)) {
161 assemble=FALSE;
162 } else {
163 if (isExpanded(X)) {
164 assemble=Finley_Util_anyNonZeroDouble(p.numEqu*p.numDim*p.numQuad,getSampleData(X,e));
165 } else {
166 assemble=Finley_Util_anyNonZeroDouble(p.numEqu*p.numDim,getSampleData(X,e));
167 }
168 }
169 if (!assemble && !isEmpty(Y)) {
170 if (isExpanded(Y)) {
171 assemble=Finley_Util_anyNonZeroDouble(p.numEqu*p.numQuad,getSampleData(Y,e));
172 } else {
173 assemble=Finley_Util_anyNonZeroDouble(p.numEqu,getSampleData(Y,e));
174 }
175 }
176 if (assemble) {
177 for (q=0;q<p.NN_row;q++) index_row[q]=p.label_row[elements->Nodes[INDEX2(p.row_node[q],e,p.NN)]];
178 /* gather V-coordinates of nodes into V: */
179 Finley_Util_Gather_double(p.NN,&(elements->Nodes[INDEX2(0,e,p.NN)]),p.numDim,nodes->Coordinates,V);
180 /* calculate dVdv(i,j,q)=V(i,k)*DSDv(k,j,q) */
181 Finley_Util_SmallMatMult(p.numDim,p.numDim*p.numQuad,dVdv,p.NS,V,p.referenceElement->dSdv);
182 /* dvdV=invert(dVdv) inplace */
183 Finley_Util_InvertSmallMat(p.numQuad,p.numDim,dVdv,dvdV,Vol);
184 /* calculate dSdV=DSDv*DvDV */
185 Finley_Util_SmallMatSetMult(p.numQuad,p.NS_row,p.numDim,dSdV,p.numDim,p.referenceElement_row->dSdv,dvdV);
186 /* scale volume: */
187 for (q=0;q<p.numQuad;q++) Vol[q]=ABS(Vol[q]*p.referenceElement->QuadWeights[q]);
188 for (q=0;q<p.NN_row*p.numEqu;q++) EM_F[q]=0;
189 if (p.numEqu==1) {
190 Finley_Assemble_RHSMatrix_Single(p.NS_row,p.numDim,p.numQuad,
191 p.referenceElement_row->S,dSdV,Vol,p.NN_row,EM_F,
192 getSampleData(X,e),isExpanded(X),
193 getSampleData(Y,e),isExpanded(Y));
194 } else {
195 Finley_Assemble_RHSMatrix_System(p.NS_row,p.numDim,p.numQuad,
196 p.numEqu,p.referenceElement_row->S,dSdV,Vol,p.NN_row,EM_F,
197 getSampleData(X,e),isExpanded(X),
198 getSampleData(Y,e),isExpanded(Y));
199 }
200 /* add */
201 #ifndef PASO_MPI
202 Finley_Util_AddScatter(p.NN_row,index_row,p.numEqu,EM_F,getSampleData(F,0));
203 #else
204 /* add influence to local degrees of freedom */
205 /* it might be more efficient to use a local buffer for copying into using Finley_Util_AddScatter()
206 then doing a memcpy at the end */
207 printf( "Adding element %d's contribution to the RHS vector...\n", e );
208 Finley_Util_AddScatter_upperBound(p.NN_row,index_row,p.numEqu,EM_F,getSampleData(F,0),p.degreeOfFreedomUpperBound);
209 #endif
210 }
211 }
212 }
213 }
214 }
215 /* clean up */
216 THREAD_MEMFREE(EM_F);
217 THREAD_MEMFREE(V);
218 THREAD_MEMFREE(dVdv);
219 THREAD_MEMFREE(dvdV);
220 THREAD_MEMFREE(dSdV);
221 THREAD_MEMFREE(Vol);
222 THREAD_MEMFREE(index_row);
223 }
224 #ifdef Finley_TRACE
225 printf("timing: assemblage PDE Right hand Side: %.4e sec\n",Finley_timer()-time0);
226 #endif
227 }
228 }
229 /*
230 * $Log$
231 * Revision 1.4 2005/09/15 03:44:21 jgs
232 * Merge of development branch dev-02 back to main trunk on 2005-09-15
233 *
234 * Revision 1.3 2005/09/01 03:31:35 jgs
235 * Merge of development branch dev-02 back to main trunk on 2005-09-01
236 *
237 * Revision 1.2 2005/08/12 01:45:42 jgs
238 * erge of development branch dev-02 back to main trunk on 2005-08-12
239 *
240 * Revision 1.1.2.3 2005/09/07 06:26:17 gross
241 * the solver from finley are put into the standalone package paso now
242 *
243 * Revision 1.1.2.2 2005/08/24 02:02:18 gross
244 * timing output switched off. solver output can be swiched through getSolution(verbose=True) now.
245 *
246 * Revision 1.1.2.1 2005/08/04 22:41:11 gross
247 * some extra routines for finley that might speed-up RHS assembling in some cases (not actived right now)
248 *
249 *
250 */

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26