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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 751 - (hide annotations)
Mon Jun 26 01:46:34 2006 UTC (12 years, 11 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 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 147
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 jgs 150 /* Author: gross@access.edu.au */
33     /* Version: $Id$ */
34 jgs 147
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 jgs 150 char error_msg[LenErrorMsg_MAX];
49 jgs 147 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 jgs 150 Finley_resetError();
56 jgs 147
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 jgs 150 if (! Finley_noError()) return;
63 jgs 147
64     /* this function assumes NS=NN */
65     if (p.NN!=p.NS) {
66 jgs 150 Finley_setError(SYSTEM_ERROR,"__FILE__: for Finley_Assemble_PDE numNodes and numShapes have to be identical.");
67 jgs 147 return;
68     }
69     if (p.numDim!=p.numElementDim) {
70 jgs 150 Finley_setError(SYSTEM_ERROR,"__FILE__: Finley_Assemble_PDE accepts volume elements only.");
71 jgs 147 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 jgs 150 Finley_setError(TYPE_ERROR,"__FILE__: unexpected function space type for coefficient X");
83 jgs 147 }
84     if (! functionSpaceTypeEqual(funcspace,Y) ) {
85 jgs 150 Finley_setError(TYPE_ERROR,"__FILE__: unexpected function space type for coefficient Y");
86 jgs 147 }
87    
88     /* check if all function spaces are the same */
89    
90     if (! numSamplesEqual(X,p.numQuad,elements->numElements) ) {
91 jgs 150 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 jgs 147 }
94    
95     if (! numSamplesEqual(Y,p.numQuad,elements->numElements) ) {
96 jgs 150 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 jgs 147 }
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 jgs 150 sprintf(error_msg,"__FILE__: coefficient X, expected shape (%d,",dimensions[0]);
107     Finley_setError(TYPE_ERROR,error_msg);
108 jgs 147 }
109     }
110     if (!isEmpty(Y)) {
111     if (!isDataPointShapeEqual(Y,0,dimensions)) {
112 jgs 150 Finley_setError(TYPE_ERROR,"__FILE__: coefficient Y, rank 0 expected.");
113 jgs 147 }
114     }
115     } else {
116     if (!isEmpty(X)) {
117     dimensions[0]=p.numEqu;
118     dimensions[1]=p.numDim;
119     if (!isDataPointShapeEqual(X,2,dimensions)) {
120 jgs 150 sprintf(error_msg,"__FILE__: coefficient X, expected shape (%d,%d)",dimensions[0],dimensions[1]);
121     Finley_setError(TYPE_ERROR,error_msg);
122 jgs 147 }
123     }
124     if (!isEmpty(Y)) {
125     dimensions[0]=p.numEqu;
126     if (!isDataPointShapeEqual(Y,1,dimensions)) {
127 jgs 150 sprintf(error_msg,"__FILE__: coefficient Y, expected shape (%d,)",dimensions[0]);
128     Finley_setError(TYPE_ERROR,error_msg);
129 jgs 147 }
130     }
131     }
132    
133 jgs 150 if (Finley_noError()) {
134 jgs 147 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 bcumming 751 #ifndef PASO_MPI
202 jgs 147 Finley_Util_AddScatter(p.NN_row,index_row,p.numEqu,EM_F,getSampleData(F,0));
203 bcumming 751 #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 jgs 147 }
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 jgs 149 #ifdef Finley_TRACE
225 jgs 147 printf("timing: assemblage PDE Right hand Side: %.4e sec\n",Finley_timer()-time0);
226 jgs 149 #endif
227 jgs 147 }
228     }
229     /*
230     * $Log$
231 jgs 150 * 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 jgs 149 * 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 jgs 147 * 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 jgs 150 * 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 jgs 149 * 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 jgs 147 * 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