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

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

Model and FemModel are now classes in their own right.
This changes the cores quite a bit.

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