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

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

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