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

Last change on this file since 3373 was 3373, checked in by Mathieu Morlighem, 15 years ago

Added parallel version of pragnostic2 (to be completed)
File size: 16.7 KB
RevLine 
[643]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#include "../DataSet/DataSet.h"
16#include "../objects/objects.h"
17#include "../EnumDefinitions/EnumDefinitions.h"
18#include "../shared/shared.h"
19
[2112]20void ProcessResults(DataSet** pnewresults, DataSet* results,Model* model,int analysis_type){
[643]21
22 int i,n;
23 Result* result=NULL;
24 Result* newresult=NULL;
[659]25
[643]26 /*output: */
27 DataSet* newresults=NULL;
28
[667]29 /*fem diagnostic models: */
[643]30 FemModel* fem_dh=NULL;
31 FemModel* fem_dv=NULL;
32 FemModel* fem_dhu=NULL;
33 FemModel* fem_ds=NULL;
34 FemModel* fem_sl=NULL;
35
[694]36 /*fem thermal models: */
37 FemModel* fem_t=NULL;
38
[667]39 /*fem prognostic models: */
40 FemModel* fem_p=NULL;
41
[758]42 /*fem control models: */
43 FemModel* fem_c=NULL;
44
[694]45 /*some parameters*/
[643]46 int ishutter;
47 int ismacayealpattyn;
48 int isstokes;
[675]49 int dim;
[643]50
51 /*intermediary: */
52 Vec u_g=NULL;
53 double* u_g_serial=NULL;
54 double* vx=NULL;
55 double* vy=NULL;
56 double* vz=NULL;
[659]57 double* vel=NULL;
[643]58 Vec p_g=NULL;
59 double* p_g_serial=NULL;
60 double* pressure=NULL;
61 double* partition=NULL;
62 double yts;
63
[758]64 double* param_g=NULL;
65 double* parameter=NULL;
[1805]66 Vec riftproperties=NULL;
67 double* riftproperties_serial=NULL;
68 int numrifts=0;
[758]69
[853]70 Vec s_g=NULL;
71 double* s_g_serial=NULL;
72 double* surface=NULL;
73
[3086]74 Vec sx_g=NULL;
75 double* sx_g_serial=NULL;
76 double* slopex=NULL;
77
78 Vec sy_g=NULL;
79 double* sy_g_serial=NULL;
80 double* slopey=NULL;
81
[853]82 Vec b_g=NULL;
83 double* b_g_serial=NULL;
84 double* bed=NULL;
85
[667]86 Vec h_g=NULL;
87 double* h_g_serial=NULL;
88 double* thickness=NULL;
89
[694]90 Vec t_g=NULL;
91 double* t_g_serial=NULL;
92 double* temperature=NULL;
93
94 Vec m_g=NULL;
95 double* m_g_serial=NULL;
96 double* melting=NULL;
97
[2892]98 Vec grad_g=NULL;
99 double* grad_g_serial=NULL;
100 double* gradient=NULL;
101
[2722]102 Vec v_g=NULL;
103 double* v_g_serial=NULL;
104
[643]105 int numberofnodes;
106
107 /*Initialize new results: */
108 newresults=new DataSet(ResultsEnum());
109
[1881]110 /*some flags needed: */
111 model->FindParam(&dim,"dim");
112 model->FindParam(&ishutter,"ishutter");
113 model->FindParam(&isstokes,"isstokes");
114 model->FindParam(&ismacayealpattyn,"ismacayealpattyn");
[643]115
[1881]116 /*Recover femmodels first: */
117 fem_dh=model->GetFormulation(DiagnosticAnalysisEnum(),HorizAnalysisEnum());
118 fem_c=fem_dh;
119 fem_dv=model->GetFormulation(DiagnosticAnalysisEnum(),VertAnalysisEnum());
120 fem_ds=model->GetFormulation(DiagnosticAnalysisEnum(),StokesAnalysisEnum());
121 fem_dhu=model->GetFormulation(DiagnosticAnalysisEnum(),HutterAnalysisEnum());
[3086]122 fem_sl=model->GetFormulation(SlopecomputeAnalysisEnum());
[2714]123 if(analysis_type==PrognosticAnalysisEnum()){
124 fem_p=model->GetFormulation(PrognosticAnalysisEnum());
125 }
[3373]126 if(analysis_type==Prognostic2AnalysisEnum()){
127 fem_p=model->GetFormulation(Prognostic2AnalysisEnum());
128 }
[2718]129 if(analysis_type==TransientAnalysisEnum()){
130 fem_p=model->GetFormulation(PrognosticAnalysisEnum());
131 }
[2714]132 if(analysis_type==BalancedthicknessAnalysisEnum()){
133 fem_p=model->GetFormulation(BalancedthicknessAnalysisEnum());
134 }
[2722]135 if(analysis_type==BalancedvelocitiesAnalysisEnum()){
136 fem_p=model->GetFormulation(BalancedvelocitiesAnalysisEnum());
137 }
[1881]138 fem_t=model->GetFormulation(ThermalAnalysisEnum());
[928]139
[643]140 for(n=0;n<results->Size();n++){
141 result=(Result*)results->GetObjectByOffset(n);
[853]142
[643]143 if(strcmp(result->GetFieldName(),"u_g")==0){
[853]144
[675]145 /*Ok, are we dealing with velocities coming from MacAyeal, Pattyin, Hutter, on 2,3 dofs or
[643]146 *Stokes on 4 dofs: */
147 result->GetField(&u_g);
148 VecToMPISerial(&u_g_serial,u_g);
149
[675]150 //2d results -> 2 dofs per node
151 if (dim==2){
[643]152 /*ok, 2 dofs, on number of nodes: */
153 if(ismacayealpattyn){
[2333]154 fem_dh->parameters->FindParam(&numberofnodes,"numberofnodes");
[2316]155 VecToMPISerial(&partition,fem_dh->partition->vector);
[2333]156 fem_dh->parameters->FindParam(&yts,"yts");
[643]157 }
158 else{
[2333]159 fem_dhu->parameters->FindParam(&numberofnodes,"numberofnodes");
[2316]160 VecToMPISerial(&partition,fem_dhu->partition->vector);
[2333]161 fem_dhu->parameters->FindParam(&yts,"yts");
[643]162 }
163 vx=(double*)xmalloc(numberofnodes*sizeof(double));
164 vy=(double*)xmalloc(numberofnodes*sizeof(double));
165 vz=(double*)xmalloc(numberofnodes*sizeof(double));
[659]166 vel=(double*)xmalloc(numberofnodes*sizeof(double));
[643]167
168 for(i=0;i<numberofnodes;i++){
169 vx[i]=u_g_serial[2*(int)partition[i]+0]*yts;
170 vy[i]=u_g_serial[2*(int)partition[i]+1]*yts;
171 vz[i]=0;
[659]172 vel[i]=sqrt(pow(vx[i],2)+pow(vy[i],2)+pow(vz[i],2));
[643]173 }
[675]174 }
175 //3d results -> 3 or 4 (stokes) dofs per node
[643]176 else{
[675]177 if(!isstokes){
178 /*ok, 3 dofs, on number of nodes: */
179 if(ismacayealpattyn){
[2333]180 fem_dh->parameters->FindParam(&numberofnodes,"numberofnodes");
[2316]181 VecToMPISerial(&partition,fem_dh->partition->vector);
[2333]182 fem_dh->parameters->FindParam(&yts,"yts");
[675]183 }
184 else{
[2333]185 fem_dhu->parameters->FindParam(&numberofnodes,"numberofnodes");
[2316]186 VecToMPISerial(&partition,fem_dhu->partition->vector);
[2333]187 fem_dhu->parameters->FindParam(&yts,"yts");
[675]188 }
189 vx=(double*)xmalloc(numberofnodes*sizeof(double));
190 vy=(double*)xmalloc(numberofnodes*sizeof(double));
191 vz=(double*)xmalloc(numberofnodes*sizeof(double));
192 vel=(double*)xmalloc(numberofnodes*sizeof(double));
193
194 for(i=0;i<numberofnodes;i++){
195 vx[i]=u_g_serial[3*(int)partition[i]+0]*yts;
196 vy[i]=u_g_serial[3*(int)partition[i]+1]*yts;
197 vz[i]=u_g_serial[3*(int)partition[i]+2]*yts;
198 vel[i]=sqrt(pow(vx[i],2)+pow(vy[i],2)+pow(vz[i],2));
199 }
[643]200 }
[675]201 else{
202 /* 4 dofs on number of nodes. discard pressure: */
[2333]203 fem_ds->parameters->FindParam(&numberofnodes,"numberofnodes");
[2316]204 VecToMPISerial(&partition,fem_ds->partition->vector);
[2333]205 fem_ds->parameters->FindParam(&yts,"yts");
[675]206 vx=(double*)xmalloc(numberofnodes*sizeof(double));
207 vy=(double*)xmalloc(numberofnodes*sizeof(double));
208 vz=(double*)xmalloc(numberofnodes*sizeof(double));
209 vel=(double*)xmalloc(numberofnodes*sizeof(double));
210 for(i=0;i<numberofnodes;i++){
211 vx[i]=u_g_serial[4*(int)partition[i]+0]*yts;
212 vy[i]=u_g_serial[4*(int)partition[i]+1]*yts;
213 vz[i]=u_g_serial[4*(int)partition[i]+2]*yts;
214 vel[i]=sqrt(pow(vx[i],2)+pow(vy[i],2)+pow(vz[i],2));
215 }
216 }
[643]217 }
218
219 /*Ok, add vx,vy and vz to newresults: */
220 newresult=new Result(newresults->Size()+1,result->GetTime(),result->GetStep(),"vx",vx,numberofnodes);
221 newresults->AddObject(newresult);
222
223 newresult=new Result(newresults->Size()+1,result->GetTime(),result->GetStep(),"vy",vy,numberofnodes);
224 newresults->AddObject(newresult);
225
226 newresult=new Result(newresults->Size()+1,result->GetTime(),result->GetStep(),"vz",vz,numberofnodes);
227 newresults->AddObject(newresult);
228
[659]229 newresult=new Result(newresults->Size()+1,result->GetTime(),result->GetStep(),"vel",vel,numberofnodes);
230 newresults->AddObject(newresult);
231
[643]232 /*do some cleanup: */
233 xfree((void**)&u_g_serial);
234 xfree((void**)&partition);
[1964]235 xfree((void**)&vx);
236 xfree((void**)&vy);
237 xfree((void**)&vz);
238 xfree((void**)&vel);
[1943]239 VecFree(&u_g);
[643]240 }
241 else if(strcmp(result->GetFieldName(),"p_g")==0){
242 /*easy, p_g is of size numberofnodes, on 1 dof, just repartition: */
243 result->GetField(&p_g);
244 VecToMPISerial(&p_g_serial,p_g);
245
246 if(!isstokes){
247 if(ismacayealpattyn){
[2333]248 fem_dh->parameters->FindParam(&numberofnodes,"numberofnodes");
[2316]249 VecToMPISerial(&partition,fem_dh->partition->vector);
[643]250 }
251 else{
[2333]252 fem_dhu->parameters->FindParam(&numberofnodes,"numberofnodes");
[2316]253 VecToMPISerial(&partition,fem_dhu->partition->vector);
[643]254 }
255 }
256 else{
[2333]257 fem_ds->parameters->FindParam(&numberofnodes,"numberofnodes");
[2316]258 VecToMPISerial(&partition,fem_ds->partition->vector);
[643]259 }
260
261 pressure=(double*)xmalloc(numberofnodes*sizeof(double));
262
263 for(i=0;i<numberofnodes;i++){
264 pressure[i]=p_g_serial[(int)partition[i]];
265 }
266
267 /*Ok, add pressure,vy and vz to newresults: */
268 newresult=new Result(newresults->Size()+1,result->GetTime(),result->GetStep(),"pressure",pressure,numberofnodes);
269 newresults->AddObject(newresult);
270
271 /*do some cleanup: */
272 xfree((void**)&p_g_serial);
273 xfree((void**)&partition);
[1964]274 xfree((void**)&pressure);
[1943]275 VecFree(&p_g);
[643]276 }
[694]277 else if(strcmp(result->GetFieldName(),"t_g")==0){
278 /*easy, t_g is of size numberofnodes, on 1 dof, just repartition: */
279 result->GetField(&t_g);
280 VecToMPISerial(&t_g_serial,t_g);
[2333]281 fem_t->parameters->FindParam(&numberofnodes,"numberofnodes");
[2316]282 VecToMPISerial(&partition,fem_t->partition->vector);
[694]283
284 temperature=(double*)xmalloc(numberofnodes*sizeof(double));
285
286 for(i=0;i<numberofnodes;i++){
287 temperature[i]=t_g_serial[(int)partition[i]];
288 }
289
290 /*Ok, add pressure to newresults: */
291 newresult=new Result(newresults->Size()+1,result->GetTime(),result->GetStep(),"temperature",temperature,numberofnodes);
292 newresults->AddObject(newresult);
293
294 /*do some cleanup: */
295 xfree((void**)&t_g_serial);
296 xfree((void**)&partition);
[1964]297 xfree((void**)&temperature);
[1943]298 VecFree(&t_g);
[694]299 }
[2892]300 else if(strcmp(result->GetFieldName(),"grad_g")==0){
301
302 /*easy, grad_g is of size numberofnodes, on 1 dof, just repartition: */
303 result->GetField(&grad_g);
304 VecToMPISerial(&grad_g_serial,grad_g);
305 fem_c->parameters->FindParam(&numberofnodes,"numberofnodes");
306 VecToMPISerial(&partition,fem_c->partition->vector);
307
308 gradient=(double*)xmalloc(numberofnodes*sizeof(double));
309
310 for(i=0;i<numberofnodes;i++){
311 gradient[i]=grad_g_serial[(int)partition[i]];
312 }
313
314 /*Ok, add gradient to newresults: */
315 newresult=new Result(newresults->Size()+1,result->GetTime(),result->GetStep(),"gradient",gradient,numberofnodes);
316 newresults->AddObject(newresult);
317
318 /*do some cleanup: */
319 xfree((void**)&grad_g_serial);
320 xfree((void**)&partition);
321 xfree((void**)&gradient);
322 VecFree(&grad_g);
323 }
[694]324 else if(strcmp(result->GetFieldName(),"m_g")==0){
325 /*easy, m_g is of size numberofnodes, on 1 dof, just repartition: */
326 result->GetField(&m_g);
327 VecToMPISerial(&m_g_serial,m_g);
[2333]328 fem_t->parameters->FindParam(&numberofnodes,"numberofnodes");
329 fem_t->parameters->FindParam(&yts,"yts");
[2316]330 VecToMPISerial(&partition,fem_t->partition->vector);
[694]331
332 melting=(double*)xmalloc(numberofnodes*sizeof(double));
333
334 for(i=0;i<numberofnodes;i++){
[695]335 melting[i]=m_g_serial[(int)partition[i]]*yts;
[694]336 }
337
338 /*Ok, add pressure to newresults: */
339 newresult=new Result(newresults->Size()+1,result->GetTime(),result->GetStep(),"melting",melting,numberofnodes);
340 newresults->AddObject(newresult);
341
342 /*do some cleanup: */
343 xfree((void**)&m_g_serial);
344 xfree((void**)&partition);
[1964]345 xfree((void**)&melting);
[1943]346 VecFree(&m_g);
[694]347 }
[667]348 else if(strcmp(result->GetFieldName(),"h_g")==0){
349 /*easy, h_g is of size numberofnodes, on 1 dof, just repartition: */
350 result->GetField(&h_g);
351 VecToMPISerial(&h_g_serial,h_g);
[2333]352 fem_p->parameters->FindParam(&numberofnodes,"numberofnodes");
[2316]353 VecToMPISerial(&partition,fem_p->partition->vector);
[667]354
355 thickness=(double*)xmalloc(numberofnodes*sizeof(double));
356
357 for(i=0;i<numberofnodes;i++){
358 thickness[i]=h_g_serial[(int)partition[i]];
359 }
360
361 /*Ok, add pressure to newresults: */
362 newresult=new Result(newresults->Size()+1,result->GetTime(),result->GetStep(),"thickness",thickness,numberofnodes);
363 newresults->AddObject(newresult);
364
365 /*do some cleanup: */
366 xfree((void**)&h_g_serial);
[1964]367 xfree((void**)&thickness);
[667]368 xfree((void**)&partition);
[1943]369 VecFree(&h_g);
[667]370 }
[2722]371 else if(strcmp(result->GetFieldName(),"v_g")==0){
372 /*easy, v_g is of size numberofnodes, on 1 dof, just repartition: */
373 result->GetField(&v_g);
374 VecToMPISerial(&v_g_serial,v_g);
375 fem_p->parameters->FindParam(&numberofnodes,"numberofnodes");
376 VecToMPISerial(&partition,fem_p->partition->vector);
377
378 vel=(double*)xmalloc(numberofnodes*sizeof(double));
379
380 for(i=0;i<numberofnodes;i++){
381 vel[i]=v_g_serial[(int)partition[i]];
382 }
383
384 /*Ok, add pressure to newresults: */
385 newresult=new Result(newresults->Size()+1,result->GetTime(),result->GetStep(),"vel",vel,numberofnodes);
386 newresults->AddObject(newresult);
387
388 /*do some cleanup: */
389 xfree((void**)&v_g_serial);
390 xfree((void**)&vel);
391 xfree((void**)&partition);
392 VecFree(&v_g);
393 }
[853]394 else if(strcmp(result->GetFieldName(),"s_g")==0){
395 /*easy, s_g is of size numberofnodes, on 1 dof, just repartition: */
396 result->GetField(&s_g);
397 VecToMPISerial(&s_g_serial,s_g);
[2333]398 fem_p->parameters->FindParam(&numberofnodes,"numberofnodes");
[2316]399 VecToMPISerial(&partition,fem_p->partition->vector);
[853]400
401 surface=(double*)xmalloc(numberofnodes*sizeof(double));
402
403 for(i=0;i<numberofnodes;i++){
404 surface[i]=s_g_serial[(int)partition[i]];
405 }
406
407 /*Ok, add pressure to newresults: */
408 newresult=new Result(newresults->Size()+1,result->GetTime(),result->GetStep(),"surface",surface,numberofnodes);
409 newresults->AddObject(newresult);
410
411 /*do some cleanup: */
412 xfree((void**)&s_g_serial);
413 xfree((void**)&partition);
[1964]414 xfree((void**)&surface);
[1943]415 VecFree(&s_g);
[853]416 }
[3086]417 else if(strcmp(result->GetFieldName(),"sx_g")==0){
418 /*easy, s_g is of size numberofnodes, on 1 dof, just repartition: */
419 result->GetField(&sx_g);
420 VecToMPISerial(&sx_g_serial,sx_g);
421 fem_sl->parameters->FindParam(&numberofnodes,"numberofnodes");
422 VecToMPISerial(&partition,fem_sl->partition->vector);
423
424 slopex=(double*)xmalloc(numberofnodes*sizeof(double));
425
426 for(i=0;i<numberofnodes;i++){
427 slopex[i]=sx_g_serial[(int)partition[i]];
428 }
429
430 /*Ok, add pressure to newresults: */
431 newresult=new Result(newresults->Size()+1,result->GetTime(),result->GetStep(),"slopex",slopex,numberofnodes);
432 newresults->AddObject(newresult);
433
434 /*do some cleanup: */
435 xfree((void**)&sx_g_serial);
436 xfree((void**)&partition);
437 xfree((void**)&slopex);
438 VecFree(&sx_g);
439 }
440 else if(strcmp(result->GetFieldName(),"sy_g")==0){
441 /*easy, s_g is of size numberofnodes, on 1 dof, just repartition: */
442 result->GetField(&sy_g);
443 VecToMPISerial(&sy_g_serial,sy_g);
444 fem_sl->parameters->FindParam(&numberofnodes,"numberofnodes");
445 VecToMPISerial(&partition,fem_sl->partition->vector);
446
447 slopey=(double*)xmalloc(numberofnodes*sizeof(double));
448
449 for(i=0;i<numberofnodes;i++){
450 slopey[i]=sy_g_serial[(int)partition[i]];
451 }
452
453 /*Ok, add pressure to newresults: */
454 newresult=new Result(newresults->Size()+1,result->GetTime(),result->GetStep(),"slopey",slopey,numberofnodes);
455 newresults->AddObject(newresult);
456
457 /*do some cleanup: */
458 xfree((void**)&sy_g_serial);
459 xfree((void**)&partition);
460 xfree((void**)&slopey);
461 VecFree(&sy_g);
462 }
[853]463 else if(strcmp(result->GetFieldName(),"b_g")==0){
464 /*easy, b_g is of size numberofnodes, on 1 dof, just repartition: */
465 result->GetField(&b_g);
466 VecToMPISerial(&b_g_serial,b_g);
[2333]467 fem_p->parameters->FindParam(&numberofnodes,"numberofnodes");
[2316]468 VecToMPISerial(&partition,fem_p->partition->vector);
[853]469
470 bed=(double*)xmalloc(numberofnodes*sizeof(double));
471
472 for(i=0;i<numberofnodes;i++){
473 bed[i]=b_g_serial[(int)partition[i]];
474 }
475
476 /*Ok, add pressure to newresults: */
477 newresult=new Result(newresults->Size()+1,result->GetTime(),result->GetStep(),"bed",bed,numberofnodes);
478 newresults->AddObject(newresult);
479
480 /*do some cleanup: */
481 xfree((void**)&b_g_serial);
482 xfree((void**)&partition);
[1964]483 xfree((void**)&bed);
[1943]484 VecFree(&b_g);
[853]485 }
[669]486 else if(strcmp(result->GetFieldName(),"param_g")==0){
487 /*easy, param_g is of size numberofnodes, on 1 dof, just repartition: */
[758]488 result->GetField(&param_g);
[2333]489 fem_dh->parameters->FindParam(&numberofnodes,"numberofnodes");
[2316]490 VecToMPISerial(&partition,fem_dh->partition->vector);
[669]491
492 parameter=(double*)xmalloc(numberofnodes*sizeof(double));
493
494 for(i=0;i<numberofnodes;i++){
[1046]495 parameter[i]=param_g[(int)partition[i]];
[669]496 }
497
498 /*Ok, add parameter to newresults: */
499 newresult=new Result(newresults->Size()+1,result->GetTime(),result->GetStep(),"parameter",parameter,numberofnodes);
500 newresults->AddObject(newresult);
501
502 /*do some cleanup: */
503 xfree((void**)&partition);
[1943]504 xfree((void**)&param_g);
[1964]505 xfree((void**)&parameter);
[669]506 }
[1805]507 else if(strcmp(result->GetFieldName(),"riftproperties")==0){
508 result->GetField(&riftproperties);
[2333]509 fem_dh->parameters->FindParam(&numrifts,"numrifts");
[1805]510 VecToMPISerial(&riftproperties_serial,riftproperties);
511
512 /*Ok, add parameter to newresults: */
513 newresult=new Result(newresults->Size()+1,result->GetTime(),result->GetStep(),"riftproperties",riftproperties_serial,numrifts);
514 newresults->AddObject(newresult);
[1943]515 xfree((void**)&riftproperties);
[1805]516
517 }
[669]518 else{
519 /*Just copy the result into the new results dataset: */
[758]520 newresults->AddObject(result->copy());
[669]521 }
[643]522 }
523
524 /*Assign output pointers:*/
[2112]525 *pnewresults=newresults;
[643]526}
Note: See TracBrowser for help on using the repository browser.