source: issm/trunk/src/c/parallel/ProcessResults.cpp@ 1037

Last change on this file since 1037 was 1037, checked in by Mathieu Morlighem, 16 years ago

Added core routine in control parallel
File size: 12.8 KB
Line 
1/*!\file: ProcessResults.cpp
2 * \brief: go through results dataset, and for each result, process it for easier retrieval
3 * by the Matlab side. This usually means splitting the velocities from the g-size nodeset
4 * to the grid set (ug->vx,vy,vz), same for pressure (p_g->pressure), etc ... It also implies
5 * departitioning of the results.
6 * This phase is necessary prior to outputting the results on disk.
7 */
8
9#ifdef HAVE_CONFIG_H
10 #include "config.h"
11#else
12#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
13#endif
14
15#undef __FUNCT__
16#define __FUNCT__ "ProcessResults"
17
18#include "../DataSet/DataSet.h"
19#include "../objects/objects.h"
20#include "../EnumDefinitions/EnumDefinitions.h"
21#include "../shared/shared.h"
22
23void ProcessResults(DataSet** presults,FemModel* fems,int analysis_type){
24
25 int i,n;
26 Result* result=NULL;
27 Result* newresult=NULL;
28
29 /*input: */
30 DataSet* results=NULL;
31
32 /*output: */
33 DataSet* newresults=NULL;
34
35 /*fem diagnostic models: */
36 FemModel* fem_dh=NULL;
37 FemModel* fem_dv=NULL;
38 FemModel* fem_dhu=NULL;
39 FemModel* fem_ds=NULL;
40 FemModel* fem_sl=NULL;
41
42 /*fem thermal models: */
43 FemModel* fem_t=NULL;
44
45 /*fem prognostic models: */
46 FemModel* fem_p=NULL;
47
48 /*fem control models: */
49 FemModel* fem_c=NULL;
50
51 /*some parameters*/
52 int ishutter;
53 int ismacayealpattyn;
54 int isstokes;
55 int dim;
56
57 /*intermediary: */
58 Vec u_g=NULL;
59 double* u_g_serial=NULL;
60 double* vx=NULL;
61 double* vy=NULL;
62 double* vz=NULL;
63 double* vel=NULL;
64 Vec p_g=NULL;
65 double* p_g_serial=NULL;
66 double* pressure=NULL;
67 double* partition=NULL;
68 double yts;
69
70 double* param_g=NULL;
71 double* parameter=NULL;
72
73 Vec s_g=NULL;
74 double* s_g_serial=NULL;
75 double* surface=NULL;
76
77 Vec b_g=NULL;
78 double* b_g_serial=NULL;
79 double* bed=NULL;
80
81 Vec h_g=NULL;
82 double* h_g_serial=NULL;
83 double* thickness=NULL;
84
85 Vec t_g=NULL;
86 double* t_g_serial=NULL;
87 double* temperature=NULL;
88
89 Vec m_g=NULL;
90 double* m_g_serial=NULL;
91 double* melting=NULL;
92
93 int numberofnodes;
94
95 /*recover input results: */
96 results=*presults;
97
98 /*Initialize new results: */
99 newresults=new DataSet(ResultsEnum());
100
101 /*Recover femmodels first: */
102 if(analysis_type==DiagnosticAnalysisEnum()){
103
104 fem_dh=fems+0;
105 fem_dv=fems+1;
106 fem_ds=fems+2;
107 fem_dhu=fems+3;
108 fem_sl=fems+4;
109
110 /*some flags needed: */
111 fem_dh->parameters->FindParam((void*)&dim,"dim");
112 fem_dhu->parameters->FindParam((void*)&ishutter,"ishutter");
113 fem_ds->parameters->FindParam((void*)&isstokes,"isstokes");
114 fem_dh->parameters->FindParam((void*)&ismacayealpattyn,"ismacayealpattyn");
115 }
116
117 if(analysis_type==TransientAnalysisEnum()){
118
119 fem_dh=fems+0;
120 fem_dv=fems+1;
121 fem_ds=fems+2;
122 fem_dhu=fems+3;
123 fem_sl=fems+4;
124 fem_p=fems+5;
125
126 /*some flags needed: */
127 fem_dh->parameters->FindParam((void*)&dim,"dim");
128 fem_dhu->parameters->FindParam((void*)&ishutter,"ishutter");
129 fem_ds->parameters->FindParam((void*)&isstokes,"isstokes");
130 fem_dh->parameters->FindParam((void*)&ismacayealpattyn,"ismacayealpattyn");
131
132 if (dim==3){
133 fem_t=fems+6;
134 }
135 }
136
137 if(analysis_type==PrognosticAnalysisEnum()){
138 fem_p=fems+0;
139 }
140
141 if(analysis_type==ThermalAnalysisEnum()){
142 fem_t=fems+0;
143 }
144
145 if(analysis_type==ControlAnalysisEnum()){
146 fem_dh=fems+0; //load u_g
147 fem_c=fems+0; //load param_g
148
149 /*some flags needed: */
150 fem_dh->parameters->FindParam((void*)&dim,"dim");
151 fem_dh->parameters->FindParam((void*)&ishutter,"ishutter");
152 fem_dh->parameters->FindParam((void*)&isstokes,"isstokes");
153 fem_dh->parameters->FindParam((void*)&ismacayealpattyn,"ismacayealpattyn");
154 }
155
156 for(n=0;n<results->Size();n++){
157 result=(Result*)results->GetObjectByOffset(n);
158
159 if(strcmp(result->GetFieldName(),"u_g")==0){
160
161 /*Ok, are we dealing with velocities coming from MacAyeal, Pattyin, Hutter, on 2,3 dofs or
162 *Stokes on 4 dofs: */
163 result->GetField(&u_g);
164 VecToMPISerial(&u_g_serial,u_g);
165
166 //2d results -> 2 dofs per node
167 if (dim==2){
168 /*ok, 2 dofs, on number of nodes: */
169 if(ismacayealpattyn){
170 fem_dh->parameters->FindParam((void*)&numberofnodes,"numberofnodes");
171 VecToMPISerial(&partition,fem_dh->partition);
172 fem_dh->parameters->FindParam((void*)&yts,"yts");
173 }
174 else{
175 fem_dhu->parameters->FindParam((void*)&numberofnodes,"numberofnodes");
176 VecToMPISerial(&partition,fem_dhu->partition);
177 fem_dhu->parameters->FindParam((void*)&yts,"yts");
178 }
179 vx=(double*)xmalloc(numberofnodes*sizeof(double));
180 vy=(double*)xmalloc(numberofnodes*sizeof(double));
181 vz=(double*)xmalloc(numberofnodes*sizeof(double));
182 vel=(double*)xmalloc(numberofnodes*sizeof(double));
183
184 for(i=0;i<numberofnodes;i++){
185 vx[i]=u_g_serial[2*(int)partition[i]+0]*yts;
186 vy[i]=u_g_serial[2*(int)partition[i]+1]*yts;
187 vz[i]=0;
188 vel[i]=sqrt(pow(vx[i],2)+pow(vy[i],2)+pow(vz[i],2));
189 }
190 }
191 //3d results -> 3 or 4 (stokes) dofs per node
192 else{
193 if(!isstokes){
194 /*ok, 3 dofs, on number of nodes: */
195 if(ismacayealpattyn){
196 fem_dh->parameters->FindParam((void*)&numberofnodes,"numberofnodes");
197 VecToMPISerial(&partition,fem_dh->partition);
198 fem_dh->parameters->FindParam((void*)&yts,"yts");
199 }
200 else{
201 fem_dhu->parameters->FindParam((void*)&numberofnodes,"numberofnodes");
202 VecToMPISerial(&partition,fem_dhu->partition);
203 fem_dhu->parameters->FindParam((void*)&yts,"yts");
204 }
205 vx=(double*)xmalloc(numberofnodes*sizeof(double));
206 vy=(double*)xmalloc(numberofnodes*sizeof(double));
207 vz=(double*)xmalloc(numberofnodes*sizeof(double));
208 vel=(double*)xmalloc(numberofnodes*sizeof(double));
209
210 for(i=0;i<numberofnodes;i++){
211 vx[i]=u_g_serial[3*(int)partition[i]+0]*yts;
212 vy[i]=u_g_serial[3*(int)partition[i]+1]*yts;
213 vz[i]=u_g_serial[3*(int)partition[i]+2]*yts;
214 vel[i]=sqrt(pow(vx[i],2)+pow(vy[i],2)+pow(vz[i],2));
215 }
216 }
217 else{
218 /* 4 dofs on number of nodes. discard pressure: */
219 fem_ds->parameters->FindParam((void*)&numberofnodes,"numberofnodes");
220 VecToMPISerial(&partition,fem_ds->partition);
221 fem_ds->parameters->FindParam((void*)&yts,"yts");
222 vx=(double*)xmalloc(numberofnodes*sizeof(double));
223 vy=(double*)xmalloc(numberofnodes*sizeof(double));
224 vz=(double*)xmalloc(numberofnodes*sizeof(double));
225 vel=(double*)xmalloc(numberofnodes*sizeof(double));
226 for(i=0;i<numberofnodes;i++){
227 vx[i]=u_g_serial[4*(int)partition[i]+0]*yts;
228 vy[i]=u_g_serial[4*(int)partition[i]+1]*yts;
229 vz[i]=u_g_serial[4*(int)partition[i]+2]*yts;
230 vel[i]=sqrt(pow(vx[i],2)+pow(vy[i],2)+pow(vz[i],2));
231 }
232 }
233 }
234
235 /*Ok, add vx,vy and vz to newresults: */
236 newresult=new Result(newresults->Size()+1,result->GetTime(),result->GetStep(),"vx",vx,numberofnodes);
237 newresults->AddObject(newresult);
238
239 newresult=new Result(newresults->Size()+1,result->GetTime(),result->GetStep(),"vy",vy,numberofnodes);
240 newresults->AddObject(newresult);
241
242 newresult=new Result(newresults->Size()+1,result->GetTime(),result->GetStep(),"vz",vz,numberofnodes);
243 newresults->AddObject(newresult);
244
245 newresult=new Result(newresults->Size()+1,result->GetTime(),result->GetStep(),"vel",vel,numberofnodes);
246 newresults->AddObject(newresult);
247
248 /*do some cleanup: */
249 xfree((void**)&u_g_serial);
250 xfree((void**)&partition);
251 }
252 else if(strcmp(result->GetFieldName(),"p_g")==0){
253 /*easy, p_g is of size numberofnodes, on 1 dof, just repartition: */
254 result->GetField(&p_g);
255 VecToMPISerial(&p_g_serial,p_g);
256
257 if(!isstokes){
258 if(ismacayealpattyn){
259 fem_dh->parameters->FindParam((void*)&numberofnodes,"numberofnodes");
260 VecToMPISerial(&partition,fem_dh->partition);
261 }
262 else{
263 fem_dhu->parameters->FindParam((void*)&numberofnodes,"numberofnodes");
264 VecToMPISerial(&partition,fem_dhu->partition);
265 }
266 }
267 else{
268 fem_ds->parameters->FindParam((void*)&numberofnodes,"numberofnodes");
269 VecToMPISerial(&partition,fem_ds->partition);
270 }
271
272 pressure=(double*)xmalloc(numberofnodes*sizeof(double));
273
274 for(i=0;i<numberofnodes;i++){
275 pressure[i]=p_g_serial[(int)partition[i]];
276 }
277
278 /*Ok, add pressure,vy and vz to newresults: */
279 newresult=new Result(newresults->Size()+1,result->GetTime(),result->GetStep(),"pressure",pressure,numberofnodes);
280 newresults->AddObject(newresult);
281
282 /*do some cleanup: */
283 xfree((void**)&p_g_serial);
284 xfree((void**)&partition);
285 }
286 else if(strcmp(result->GetFieldName(),"t_g")==0){
287 /*easy, t_g is of size numberofnodes, on 1 dof, just repartition: */
288 result->GetField(&t_g);
289 VecToMPISerial(&t_g_serial,t_g);
290 fem_t->parameters->FindParam((void*)&numberofnodes,"numberofnodes");
291 VecToMPISerial(&partition,fem_t->partition);
292
293 temperature=(double*)xmalloc(numberofnodes*sizeof(double));
294
295 for(i=0;i<numberofnodes;i++){
296 temperature[i]=t_g_serial[(int)partition[i]];
297 }
298
299 /*Ok, add pressure to newresults: */
300 newresult=new Result(newresults->Size()+1,result->GetTime(),result->GetStep(),"temperature",temperature,numberofnodes);
301 newresults->AddObject(newresult);
302
303 /*do some cleanup: */
304 xfree((void**)&t_g_serial);
305 xfree((void**)&partition);
306 }
307 else if(strcmp(result->GetFieldName(),"m_g")==0){
308 /*easy, m_g is of size numberofnodes, on 1 dof, just repartition: */
309 result->GetField(&m_g);
310 VecToMPISerial(&m_g_serial,m_g);
311 fem_t->parameters->FindParam((void*)&numberofnodes,"numberofnodes");
312 fem_t->parameters->FindParam((void*)&yts,"yts");
313 VecToMPISerial(&partition,fem_t->partition);
314
315 melting=(double*)xmalloc(numberofnodes*sizeof(double));
316
317 for(i=0;i<numberofnodes;i++){
318 melting[i]=m_g_serial[(int)partition[i]]*yts;
319 }
320
321 /*Ok, add pressure to newresults: */
322 newresult=new Result(newresults->Size()+1,result->GetTime(),result->GetStep(),"melting",melting,numberofnodes);
323 newresults->AddObject(newresult);
324
325 /*do some cleanup: */
326 xfree((void**)&m_g_serial);
327 xfree((void**)&partition);
328 }
329 else if(strcmp(result->GetFieldName(),"h_g")==0){
330 /*easy, h_g is of size numberofnodes, on 1 dof, just repartition: */
331 result->GetField(&h_g);
332 VecToMPISerial(&h_g_serial,h_g);
333 fem_p->parameters->FindParam((void*)&numberofnodes,"numberofnodes");
334 VecToMPISerial(&partition,fem_p->partition);
335
336 thickness=(double*)xmalloc(numberofnodes*sizeof(double));
337
338 for(i=0;i<numberofnodes;i++){
339 thickness[i]=h_g_serial[(int)partition[i]];
340 }
341
342 /*Ok, add pressure to newresults: */
343 newresult=new Result(newresults->Size()+1,result->GetTime(),result->GetStep(),"thickness",thickness,numberofnodes);
344 newresults->AddObject(newresult);
345
346 /*do some cleanup: */
347 xfree((void**)&h_g_serial);
348 xfree((void**)&partition);
349 }
350 else if(strcmp(result->GetFieldName(),"s_g")==0){
351 /*easy, s_g is of size numberofnodes, on 1 dof, just repartition: */
352 result->GetField(&s_g);
353 VecToMPISerial(&s_g_serial,s_g);
354 fem_p->parameters->FindParam((void*)&numberofnodes,"numberofnodes");
355 VecToMPISerial(&partition,fem_p->partition);
356
357 surface=(double*)xmalloc(numberofnodes*sizeof(double));
358
359 for(i=0;i<numberofnodes;i++){
360 surface[i]=s_g_serial[(int)partition[i]];
361 }
362
363 /*Ok, add pressure to newresults: */
364 newresult=new Result(newresults->Size()+1,result->GetTime(),result->GetStep(),"surface",surface,numberofnodes);
365 newresults->AddObject(newresult);
366
367 /*do some cleanup: */
368 xfree((void**)&s_g_serial);
369 xfree((void**)&partition);
370 }
371 else if(strcmp(result->GetFieldName(),"b_g")==0){
372 /*easy, b_g is of size numberofnodes, on 1 dof, just repartition: */
373 result->GetField(&b_g);
374 VecToMPISerial(&b_g_serial,b_g);
375 fem_p->parameters->FindParam((void*)&numberofnodes,"numberofnodes");
376 VecToMPISerial(&partition,fem_p->partition);
377
378 bed=(double*)xmalloc(numberofnodes*sizeof(double));
379
380 for(i=0;i<numberofnodes;i++){
381 bed[i]=b_g_serial[(int)partition[i]];
382 }
383
384 /*Ok, add pressure to newresults: */
385 newresult=new Result(newresults->Size()+1,result->GetTime(),result->GetStep(),"bed",bed,numberofnodes);
386 newresults->AddObject(newresult);
387
388 /*do some cleanup: */
389 xfree((void**)&b_g_serial);
390 xfree((void**)&partition);
391 }
392 else if(strcmp(result->GetFieldName(),"param_g")==0){
393 /*easy, param_g is of size numberofnodes, on 1 dof, just repartition: */
394 result->GetField(&param_g);
395 fem_dh->parameters->FindParam((void*)&numberofnodes,"numberofnodes");
396 VecToMPISerial(&partition,fem_dh->partition);
397
398 parameter=(double*)xmalloc(numberofnodes*sizeof(double));
399
400 for(i=0;i<numberofnodes;i++){
401 parameter[i]=param_g[2*(int)partition[i]];
402 }
403
404 /*Ok, add parameter to newresults: */
405 newresult=new Result(newresults->Size()+1,result->GetTime(),result->GetStep(),"parameter",parameter,numberofnodes);
406 newresults->AddObject(newresult);
407
408 /*do some cleanup: */
409 xfree((void**)&partition);
410 }
411 else{
412 /*Just copy the result into the new results dataset: */
413 newresults->AddObject(result->copy());
414 }
415 }
416
417 /*Delete results: */
418 delete results;
419
420
421 /*Assign output pointers:*/
422 *presults=newresults;
423}
Note: See TracBrowser for help on using the repository browser.