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

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

Heavy debugging of Mass flux capability
File size: 12.9 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** pnewresults, DataSet* results,Model* model,int analysis_type){
24
25 int i,n;
26 Result* result=NULL;
27 Result* newresult=NULL;
28
29 /*output: */
30 DataSet* newresults=NULL;
31
32 /*fem diagnostic models: */
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
39 /*fem thermal models: */
40 FemModel* fem_t=NULL;
41
42 /*fem prognostic models: */
43 FemModel* fem_p=NULL;
44
45 /*fem control models: */
46 FemModel* fem_c=NULL;
47
48 /*some parameters*/
49 int ishutter;
50 int ismacayealpattyn;
51 int isstokes;
52 int dim;
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;
60 double* vel=NULL;
61 Vec p_g=NULL;
62 double* p_g_serial=NULL;
63 double* pressure=NULL;
64 double* partition=NULL;
65 double yts;
66
67 double* param_g=NULL;
68 double* parameter=NULL;
69 Vec riftproperties=NULL;
70 double* riftproperties_serial=NULL;
71 int numrifts=0;
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 /*Initialize new results: */
96 newresults=new DataSet(ResultsEnum());
97
98 /*some flags needed: */
99 model->FindParam(&dim,"dim");
100 model->FindParam(&ishutter,"ishutter");
101 model->FindParam(&isstokes,"isstokes");
102 model->FindParam(&ismacayealpattyn,"ismacayealpattyn");
103
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());
113
114 for(n=0;n<results->Size();n++){
115 result=(Result*)results->GetObjectByOffset(n);
116
117 if(strcmp(result->GetFieldName(),"u_g")==0){
118
119 /*Ok, are we dealing with velocities coming from MacAyeal, Pattyin, Hutter, on 2,3 dofs or
120 *Stokes on 4 dofs: */
121 result->GetField(&u_g);
122 VecToMPISerial(&u_g_serial,u_g);
123
124 //2d results -> 2 dofs per node
125 if (dim==2){
126 /*ok, 2 dofs, on number of nodes: */
127 if(ismacayealpattyn){
128 fem_dh->parameters->FindParam((void*)&numberofnodes,"numberofnodes");
129 VecToMPISerial(&partition,fem_dh->partition);
130 fem_dh->parameters->FindParam((void*)&yts,"yts");
131 }
132 else{
133 fem_dhu->parameters->FindParam((void*)&numberofnodes,"numberofnodes");
134 VecToMPISerial(&partition,fem_dhu->partition);
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));
140 vel=(double*)xmalloc(numberofnodes*sizeof(double));
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;
146 vel[i]=sqrt(pow(vx[i],2)+pow(vy[i],2)+pow(vz[i],2));
147 }
148 }
149 //3d results -> 3 or 4 (stokes) dofs per node
150 else{
151 if(!isstokes){
152 /*ok, 3 dofs, on number of nodes: */
153 if(ismacayealpattyn){
154 fem_dh->parameters->FindParam((void*)&numberofnodes,"numberofnodes");
155 VecToMPISerial(&partition,fem_dh->partition);
156 fem_dh->parameters->FindParam((void*)&yts,"yts");
157 }
158 else{
159 fem_dhu->parameters->FindParam((void*)&numberofnodes,"numberofnodes");
160 VecToMPISerial(&partition,fem_dhu->partition);
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 }
174 }
175 else{
176 /* 4 dofs on number of nodes. discard pressure: */
177 fem_ds->parameters->FindParam((void*)&numberofnodes,"numberofnodes");
178 VecToMPISerial(&partition,fem_ds->partition);
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 }
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
203 newresult=new Result(newresults->Size()+1,result->GetTime(),result->GetStep(),"vel",vel,numberofnodes);
204 newresults->AddObject(newresult);
205
206 /*do some cleanup: */
207 xfree((void**)&u_g_serial);
208 xfree((void**)&partition);
209 xfree((void**)&vx);
210 xfree((void**)&vy);
211 xfree((void**)&vz);
212 xfree((void**)&vel);
213 VecFree(&u_g);
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");
223 VecToMPISerial(&partition,fem_dh->partition);
224 }
225 else{
226 fem_dhu->parameters->FindParam((void*)&numberofnodes,"numberofnodes");
227 VecToMPISerial(&partition,fem_dhu->partition);
228 }
229 }
230 else{
231 fem_ds->parameters->FindParam((void*)&numberofnodes,"numberofnodes");
232 VecToMPISerial(&partition,fem_ds->partition);
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);
248 xfree((void**)&pressure);
249 VecFree(&p_g);
250 }
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");
256 VecToMPISerial(&partition,fem_t->partition);
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);
271 xfree((void**)&temperature);
272 VecFree(&t_g);
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");
279 fem_t->parameters->FindParam((void*)&yts,"yts");
280 VecToMPISerial(&partition,fem_t->partition);
281
282 melting=(double*)xmalloc(numberofnodes*sizeof(double));
283
284 for(i=0;i<numberofnodes;i++){
285 melting[i]=m_g_serial[(int)partition[i]]*yts;
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);
295 xfree((void**)&melting);
296 VecFree(&m_g);
297 }
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");
303 VecToMPISerial(&partition,fem_p->partition);
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);
317 xfree((void**)&thickness);
318 xfree((void**)&partition);
319 VecFree(&h_g);
320 }
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");
326 VecToMPISerial(&partition,fem_p->partition);
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);
341 xfree((void**)&surface);
342 VecFree(&s_g);
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");
349 VecToMPISerial(&partition,fem_p->partition);
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);
364 xfree((void**)&bed);
365 VecFree(&b_g);
366 }
367 else if(strcmp(result->GetFieldName(),"param_g")==0){
368 /*easy, param_g is of size numberofnodes, on 1 dof, just repartition: */
369 result->GetField(&param_g);
370 fem_dh->parameters->FindParam((void*)&numberofnodes,"numberofnodes");
371 VecToMPISerial(&partition,fem_dh->partition);
372
373 parameter=(double*)xmalloc(numberofnodes*sizeof(double));
374
375 for(i=0;i<numberofnodes;i++){
376 parameter[i]=param_g[(int)partition[i]];
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);
385 xfree((void**)&param_g);
386 xfree((void**)&parameter);
387 }
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);
396 xfree((void**)&riftproperties);
397
398 }
399 else{
400 /*Just copy the result into the new results dataset: */
401 newresults->AddObject(result->copy());
402 }
403 }
404
405 /*Assign output pointers:*/
406 *pnewresults=newresults;
407}
Note: See TracBrowser for help on using the repository browser.