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

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

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