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

Last change on this file since 2316 was 2316, checked in by Eric.Larour, 15 years ago

Started switching to new DofVec framework.
partition, tpartition and now yg are dofvecs.
File size: 13.0 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
[643]93 int numberofnodes;
94
95 /*Initialize new results: */
96 newresults=new DataSet(ResultsEnum());
97
[1881]98 /*some flags needed: */
99 model->FindParam(&dim,"dim");
100 model->FindParam(&ishutter,"ishutter");
101 model->FindParam(&isstokes,"isstokes");
102 model->FindParam(&ismacayealpattyn,"ismacayealpattyn");
[643]103
[1881]104 /*Recover femmodels first: */
105 fem_dh=model->GetFormulation(DiagnosticAnalysisEnum(),HorizAnalysisEnum());
106 fem_c=fem_dh;
107 fem_dv=model->GetFormulation(DiagnosticAnalysisEnum(),VertAnalysisEnum());
108 fem_ds=model->GetFormulation(DiagnosticAnalysisEnum(),StokesAnalysisEnum());
109 fem_dhu=model->GetFormulation(DiagnosticAnalysisEnum(),HutterAnalysisEnum());
110 fem_sl=model->GetFormulation(SlopeComputeAnalysisEnum());
111 fem_p=model->GetFormulation(PrognosticAnalysisEnum());
112 fem_t=model->GetFormulation(ThermalAnalysisEnum());
[928]113
[643]114 for(n=0;n<results->Size();n++){
115 result=(Result*)results->GetObjectByOffset(n);
[853]116
[643]117 if(strcmp(result->GetFieldName(),"u_g")==0){
[853]118
[675]119 /*Ok, are we dealing with velocities coming from MacAyeal, Pattyin, Hutter, on 2,3 dofs or
[643]120 *Stokes on 4 dofs: */
121 result->GetField(&u_g);
122 VecToMPISerial(&u_g_serial,u_g);
123
[675]124 //2d results -> 2 dofs per node
125 if (dim==2){
[643]126 /*ok, 2 dofs, on number of nodes: */
127 if(ismacayealpattyn){
128 fem_dh->parameters->FindParam((void*)&numberofnodes,"numberofnodes");
[2316]129 VecToMPISerial(&partition,fem_dh->partition->vector);
[643]130 fem_dh->parameters->FindParam((void*)&yts,"yts");
131 }
132 else{
133 fem_dhu->parameters->FindParam((void*)&numberofnodes,"numberofnodes");
[2316]134 VecToMPISerial(&partition,fem_dhu->partition->vector);
[643]135 fem_dhu->parameters->FindParam((void*)&yts,"yts");
136 }
137 vx=(double*)xmalloc(numberofnodes*sizeof(double));
138 vy=(double*)xmalloc(numberofnodes*sizeof(double));
139 vz=(double*)xmalloc(numberofnodes*sizeof(double));
[659]140 vel=(double*)xmalloc(numberofnodes*sizeof(double));
[643]141
142 for(i=0;i<numberofnodes;i++){
143 vx[i]=u_g_serial[2*(int)partition[i]+0]*yts;
144 vy[i]=u_g_serial[2*(int)partition[i]+1]*yts;
145 vz[i]=0;
[659]146 vel[i]=sqrt(pow(vx[i],2)+pow(vy[i],2)+pow(vz[i],2));
[643]147 }
[675]148 }
149 //3d results -> 3 or 4 (stokes) dofs per node
[643]150 else{
[675]151 if(!isstokes){
152 /*ok, 3 dofs, on number of nodes: */
153 if(ismacayealpattyn){
154 fem_dh->parameters->FindParam((void*)&numberofnodes,"numberofnodes");
[2316]155 VecToMPISerial(&partition,fem_dh->partition->vector);
[675]156 fem_dh->parameters->FindParam((void*)&yts,"yts");
157 }
158 else{
159 fem_dhu->parameters->FindParam((void*)&numberofnodes,"numberofnodes");
[2316]160 VecToMPISerial(&partition,fem_dhu->partition->vector);
[675]161 fem_dhu->parameters->FindParam((void*)&yts,"yts");
162 }
163 vx=(double*)xmalloc(numberofnodes*sizeof(double));
164 vy=(double*)xmalloc(numberofnodes*sizeof(double));
165 vz=(double*)xmalloc(numberofnodes*sizeof(double));
166 vel=(double*)xmalloc(numberofnodes*sizeof(double));
167
168 for(i=0;i<numberofnodes;i++){
169 vx[i]=u_g_serial[3*(int)partition[i]+0]*yts;
170 vy[i]=u_g_serial[3*(int)partition[i]+1]*yts;
171 vz[i]=u_g_serial[3*(int)partition[i]+2]*yts;
172 vel[i]=sqrt(pow(vx[i],2)+pow(vy[i],2)+pow(vz[i],2));
173 }
[643]174 }
[675]175 else{
176 /* 4 dofs on number of nodes. discard pressure: */
177 fem_ds->parameters->FindParam((void*)&numberofnodes,"numberofnodes");
[2316]178 VecToMPISerial(&partition,fem_ds->partition->vector);
[675]179 fem_ds->parameters->FindParam((void*)&yts,"yts");
180 vx=(double*)xmalloc(numberofnodes*sizeof(double));
181 vy=(double*)xmalloc(numberofnodes*sizeof(double));
182 vz=(double*)xmalloc(numberofnodes*sizeof(double));
183 vel=(double*)xmalloc(numberofnodes*sizeof(double));
184 for(i=0;i<numberofnodes;i++){
185 vx[i]=u_g_serial[4*(int)partition[i]+0]*yts;
186 vy[i]=u_g_serial[4*(int)partition[i]+1]*yts;
187 vz[i]=u_g_serial[4*(int)partition[i]+2]*yts;
188 vel[i]=sqrt(pow(vx[i],2)+pow(vy[i],2)+pow(vz[i],2));
189 }
190 }
[643]191 }
192
193 /*Ok, add vx,vy and vz to newresults: */
194 newresult=new Result(newresults->Size()+1,result->GetTime(),result->GetStep(),"vx",vx,numberofnodes);
195 newresults->AddObject(newresult);
196
197 newresult=new Result(newresults->Size()+1,result->GetTime(),result->GetStep(),"vy",vy,numberofnodes);
198 newresults->AddObject(newresult);
199
200 newresult=new Result(newresults->Size()+1,result->GetTime(),result->GetStep(),"vz",vz,numberofnodes);
201 newresults->AddObject(newresult);
202
[659]203 newresult=new Result(newresults->Size()+1,result->GetTime(),result->GetStep(),"vel",vel,numberofnodes);
204 newresults->AddObject(newresult);
205
[643]206 /*do some cleanup: */
207 xfree((void**)&u_g_serial);
208 xfree((void**)&partition);
[1964]209 xfree((void**)&vx);
210 xfree((void**)&vy);
211 xfree((void**)&vz);
212 xfree((void**)&vel);
[1943]213 VecFree(&u_g);
[643]214 }
215 else if(strcmp(result->GetFieldName(),"p_g")==0){
216 /*easy, p_g is of size numberofnodes, on 1 dof, just repartition: */
217 result->GetField(&p_g);
218 VecToMPISerial(&p_g_serial,p_g);
219
220 if(!isstokes){
221 if(ismacayealpattyn){
222 fem_dh->parameters->FindParam((void*)&numberofnodes,"numberofnodes");
[2316]223 VecToMPISerial(&partition,fem_dh->partition->vector);
[643]224 }
225 else{
226 fem_dhu->parameters->FindParam((void*)&numberofnodes,"numberofnodes");
[2316]227 VecToMPISerial(&partition,fem_dhu->partition->vector);
[643]228 }
229 }
230 else{
231 fem_ds->parameters->FindParam((void*)&numberofnodes,"numberofnodes");
[2316]232 VecToMPISerial(&partition,fem_ds->partition->vector);
[643]233 }
234
235 pressure=(double*)xmalloc(numberofnodes*sizeof(double));
236
237 for(i=0;i<numberofnodes;i++){
238 pressure[i]=p_g_serial[(int)partition[i]];
239 }
240
241 /*Ok, add pressure,vy and vz to newresults: */
242 newresult=new Result(newresults->Size()+1,result->GetTime(),result->GetStep(),"pressure",pressure,numberofnodes);
243 newresults->AddObject(newresult);
244
245 /*do some cleanup: */
246 xfree((void**)&p_g_serial);
247 xfree((void**)&partition);
[1964]248 xfree((void**)&pressure);
[1943]249 VecFree(&p_g);
[643]250 }
[694]251 else if(strcmp(result->GetFieldName(),"t_g")==0){
252 /*easy, t_g is of size numberofnodes, on 1 dof, just repartition: */
253 result->GetField(&t_g);
254 VecToMPISerial(&t_g_serial,t_g);
255 fem_t->parameters->FindParam((void*)&numberofnodes,"numberofnodes");
[2316]256 VecToMPISerial(&partition,fem_t->partition->vector);
[694]257
258 temperature=(double*)xmalloc(numberofnodes*sizeof(double));
259
260 for(i=0;i<numberofnodes;i++){
261 temperature[i]=t_g_serial[(int)partition[i]];
262 }
263
264 /*Ok, add pressure to newresults: */
265 newresult=new Result(newresults->Size()+1,result->GetTime(),result->GetStep(),"temperature",temperature,numberofnodes);
266 newresults->AddObject(newresult);
267
268 /*do some cleanup: */
269 xfree((void**)&t_g_serial);
270 xfree((void**)&partition);
[1964]271 xfree((void**)&temperature);
[1943]272 VecFree(&t_g);
[694]273 }
274 else if(strcmp(result->GetFieldName(),"m_g")==0){
275 /*easy, m_g is of size numberofnodes, on 1 dof, just repartition: */
276 result->GetField(&m_g);
277 VecToMPISerial(&m_g_serial,m_g);
278 fem_t->parameters->FindParam((void*)&numberofnodes,"numberofnodes");
[695]279 fem_t->parameters->FindParam((void*)&yts,"yts");
[2316]280 VecToMPISerial(&partition,fem_t->partition->vector);
[694]281
282 melting=(double*)xmalloc(numberofnodes*sizeof(double));
283
284 for(i=0;i<numberofnodes;i++){
[695]285 melting[i]=m_g_serial[(int)partition[i]]*yts;
[694]286 }
287
288 /*Ok, add pressure to newresults: */
289 newresult=new Result(newresults->Size()+1,result->GetTime(),result->GetStep(),"melting",melting,numberofnodes);
290 newresults->AddObject(newresult);
291
292 /*do some cleanup: */
293 xfree((void**)&m_g_serial);
294 xfree((void**)&partition);
[1964]295 xfree((void**)&melting);
[1943]296 VecFree(&m_g);
[694]297 }
[667]298 else if(strcmp(result->GetFieldName(),"h_g")==0){
299 /*easy, h_g is of size numberofnodes, on 1 dof, just repartition: */
300 result->GetField(&h_g);
301 VecToMPISerial(&h_g_serial,h_g);
302 fem_p->parameters->FindParam((void*)&numberofnodes,"numberofnodes");
[2316]303 VecToMPISerial(&partition,fem_p->partition->vector);
[667]304
305 thickness=(double*)xmalloc(numberofnodes*sizeof(double));
306
307 for(i=0;i<numberofnodes;i++){
308 thickness[i]=h_g_serial[(int)partition[i]];
309 }
310
311 /*Ok, add pressure to newresults: */
312 newresult=new Result(newresults->Size()+1,result->GetTime(),result->GetStep(),"thickness",thickness,numberofnodes);
313 newresults->AddObject(newresult);
314
315 /*do some cleanup: */
316 xfree((void**)&h_g_serial);
[1964]317 xfree((void**)&thickness);
[667]318 xfree((void**)&partition);
[1943]319 VecFree(&h_g);
[667]320 }
[853]321 else if(strcmp(result->GetFieldName(),"s_g")==0){
322 /*easy, s_g is of size numberofnodes, on 1 dof, just repartition: */
323 result->GetField(&s_g);
324 VecToMPISerial(&s_g_serial,s_g);
325 fem_p->parameters->FindParam((void*)&numberofnodes,"numberofnodes");
[2316]326 VecToMPISerial(&partition,fem_p->partition->vector);
[853]327
328 surface=(double*)xmalloc(numberofnodes*sizeof(double));
329
330 for(i=0;i<numberofnodes;i++){
331 surface[i]=s_g_serial[(int)partition[i]];
332 }
333
334 /*Ok, add pressure to newresults: */
335 newresult=new Result(newresults->Size()+1,result->GetTime(),result->GetStep(),"surface",surface,numberofnodes);
336 newresults->AddObject(newresult);
337
338 /*do some cleanup: */
339 xfree((void**)&s_g_serial);
340 xfree((void**)&partition);
[1964]341 xfree((void**)&surface);
[1943]342 VecFree(&s_g);
[853]343 }
344 else if(strcmp(result->GetFieldName(),"b_g")==0){
345 /*easy, b_g is of size numberofnodes, on 1 dof, just repartition: */
346 result->GetField(&b_g);
347 VecToMPISerial(&b_g_serial,b_g);
348 fem_p->parameters->FindParam((void*)&numberofnodes,"numberofnodes");
[2316]349 VecToMPISerial(&partition,fem_p->partition->vector);
[853]350
351 bed=(double*)xmalloc(numberofnodes*sizeof(double));
352
353 for(i=0;i<numberofnodes;i++){
354 bed[i]=b_g_serial[(int)partition[i]];
355 }
356
357 /*Ok, add pressure to newresults: */
358 newresult=new Result(newresults->Size()+1,result->GetTime(),result->GetStep(),"bed",bed,numberofnodes);
359 newresults->AddObject(newresult);
360
361 /*do some cleanup: */
362 xfree((void**)&b_g_serial);
363 xfree((void**)&partition);
[1964]364 xfree((void**)&bed);
[1943]365 VecFree(&b_g);
[853]366 }
[669]367 else if(strcmp(result->GetFieldName(),"param_g")==0){
368 /*easy, param_g is of size numberofnodes, on 1 dof, just repartition: */
[758]369 result->GetField(&param_g);
370 fem_dh->parameters->FindParam((void*)&numberofnodes,"numberofnodes");
[2316]371 VecToMPISerial(&partition,fem_dh->partition->vector);
[669]372
373 parameter=(double*)xmalloc(numberofnodes*sizeof(double));
374
375 for(i=0;i<numberofnodes;i++){
[1046]376 parameter[i]=param_g[(int)partition[i]];
[669]377 }
378
379 /*Ok, add parameter to newresults: */
380 newresult=new Result(newresults->Size()+1,result->GetTime(),result->GetStep(),"parameter",parameter,numberofnodes);
381 newresults->AddObject(newresult);
382
383 /*do some cleanup: */
384 xfree((void**)&partition);
[1943]385 xfree((void**)&param_g);
[1964]386 xfree((void**)&parameter);
[669]387 }
[1805]388 else if(strcmp(result->GetFieldName(),"riftproperties")==0){
389 result->GetField(&riftproperties);
390 fem_dh->parameters->FindParam((void*)&numrifts,"numrifts");
391 VecToMPISerial(&riftproperties_serial,riftproperties);
392
393 /*Ok, add parameter to newresults: */
394 newresult=new Result(newresults->Size()+1,result->GetTime(),result->GetStep(),"riftproperties",riftproperties_serial,numrifts);
395 newresults->AddObject(newresult);
[1943]396 xfree((void**)&riftproperties);
[1805]397
398 }
[669]399 else{
400 /*Just copy the result into the new results dataset: */
[758]401 newresults->AddObject(result->copy());
[669]402 }
[643]403 }
404
405 /*Assign output pointers:*/
[2112]406 *pnewresults=newresults;
[643]407}
Note: See TracBrowser for help on using the repository browser.