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

Last change on this file since 758 was 758, checked in by Mathieu Morlighem, 16 years ago

fixed control method output
File size: 10.5 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,FemModel* fems,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
73 Vec h_g=NULL;
74 double* h_g_serial=NULL;
75 double* thickness=NULL;
76
77 Vec t_g=NULL;
78 double* t_g_serial=NULL;
79 double* temperature=NULL;
80
81 Vec m_g=NULL;
82 double* m_g_serial=NULL;
83 double* melting=NULL;
84
85 int numberofnodes;
86
87 /*recover input results: */
88 results=*presults;
89
90 /*Initialize new results: */
91 newresults=new DataSet(ResultsEnum());
92
93 /*Recover femmodels first: */
94 if(analysis_type==DiagnosticAnalysisEnum()){
95
96 fem_dh=fems+0;
97 fem_dv=fems+1;
98 fem_ds=fems+2;
99 fem_dhu=fems+3;
100 fem_sl=fems+4;
101
102 /*some flags needed: */
103 fem_dh->parameters->FindParam((void*)&dim,"dim");
104 fem_dhu->parameters->FindParam((void*)&ishutter,"ishutter");
105 fem_ds->parameters->FindParam((void*)&isstokes,"isstokes");
106 fem_dh->parameters->FindParam((void*)&ismacayealpattyn,"ismacayealpattyn");
107 }
108
109 if(analysis_type==PrognosticAnalysisEnum()){
110 fem_p=fems+0;
111 }
112
113 if(analysis_type==ThermalAnalysisEnum()){
114 fem_t=fems+0;
115 }
116
117 if(analysis_type==ControlAnalysisEnum()){
118 fem_dh=fems+0; //load u_g
119 fem_c=fems+0; //load param_g
120 }
121
122 for(n=0;n<results->Size();n++){
123
124 result=(Result*)results->GetObjectByOffset(n);
125
126 if(strcmp(result->GetFieldName(),"u_g")==0){
127 /*Ok, are we dealing with velocities coming from MacAyeal, Pattyin, Hutter, on 2,3 dofs or
128 *Stokes on 4 dofs: */
129 result->GetField(&u_g);
130 VecToMPISerial(&u_g_serial,u_g);
131
132 //2d results -> 2 dofs per node
133 if (dim==2){
134 /*ok, 2 dofs, on number of nodes: */
135 if(ismacayealpattyn){
136 fem_dh->parameters->FindParam((void*)&numberofnodes,"numberofnodes");
137 VecToMPISerial(&partition,fem_dh->partition);
138 fem_dh->parameters->FindParam((void*)&yts,"yts");
139 }
140 else{
141 fem_dhu->parameters->FindParam((void*)&numberofnodes,"numberofnodes");
142 VecToMPISerial(&partition,fem_dhu->partition);
143 fem_dhu->parameters->FindParam((void*)&yts,"yts");
144 }
145 vx=(double*)xmalloc(numberofnodes*sizeof(double));
146 vy=(double*)xmalloc(numberofnodes*sizeof(double));
147 vz=(double*)xmalloc(numberofnodes*sizeof(double));
148 vel=(double*)xmalloc(numberofnodes*sizeof(double));
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;
154 vel[i]=sqrt(pow(vx[i],2)+pow(vy[i],2)+pow(vz[i],2));
155 }
156 }
157
158 //3d results -> 3 or 4 (stokes) dofs per node
159 else{
160
161 if(!isstokes){
162 /*ok, 3 dofs, on number of nodes: */
163 if(ismacayealpattyn){
164 fem_dh->parameters->FindParam((void*)&numberofnodes,"numberofnodes");
165 VecToMPISerial(&partition,fem_dh->partition);
166 fem_dh->parameters->FindParam((void*)&yts,"yts");
167 }
168 else{
169 fem_dhu->parameters->FindParam((void*)&numberofnodes,"numberofnodes");
170 VecToMPISerial(&partition,fem_dhu->partition);
171 fem_dhu->parameters->FindParam((void*)&yts,"yts");
172 }
173 vx=(double*)xmalloc(numberofnodes*sizeof(double));
174 vy=(double*)xmalloc(numberofnodes*sizeof(double));
175 vz=(double*)xmalloc(numberofnodes*sizeof(double));
176 vel=(double*)xmalloc(numberofnodes*sizeof(double));
177
178 for(i=0;i<numberofnodes;i++){
179 vx[i]=u_g_serial[3*(int)partition[i]+0]*yts;
180 vy[i]=u_g_serial[3*(int)partition[i]+1]*yts;
181 vz[i]=u_g_serial[3*(int)partition[i]+2]*yts;
182 vel[i]=sqrt(pow(vx[i],2)+pow(vy[i],2)+pow(vz[i],2));
183 }
184 }
185
186 else{
187 /* 4 dofs on number of nodes. discard pressure: */
188 fem_ds->parameters->FindParam((void*)&numberofnodes,"numberofnodes");
189 VecToMPISerial(&partition,fem_ds->partition);
190 fem_ds->parameters->FindParam((void*)&yts,"yts");
191 vx=(double*)xmalloc(numberofnodes*sizeof(double));
192 vy=(double*)xmalloc(numberofnodes*sizeof(double));
193 vz=(double*)xmalloc(numberofnodes*sizeof(double));
194 vel=(double*)xmalloc(numberofnodes*sizeof(double));
195 for(i=0;i<numberofnodes;i++){
196 vx[i]=u_g_serial[4*(int)partition[i]+0]*yts;
197 vy[i]=u_g_serial[4*(int)partition[i]+1]*yts;
198 vz[i]=u_g_serial[4*(int)partition[i]+2]*yts;
199 vel[i]=sqrt(pow(vx[i],2)+pow(vy[i],2)+pow(vz[i],2));
200 }
201 }
202 }
203
204 /*Ok, add vx,vy and vz to newresults: */
205 newresult=new Result(newresults->Size()+1,result->GetTime(),result->GetStep(),"vx",vx,numberofnodes);
206 newresults->AddObject(newresult);
207
208 newresult=new Result(newresults->Size()+1,result->GetTime(),result->GetStep(),"vy",vy,numberofnodes);
209 newresults->AddObject(newresult);
210
211 newresult=new Result(newresults->Size()+1,result->GetTime(),result->GetStep(),"vz",vz,numberofnodes);
212 newresults->AddObject(newresult);
213
214 newresult=new Result(newresults->Size()+1,result->GetTime(),result->GetStep(),"vel",vel,numberofnodes);
215 newresults->AddObject(newresult);
216
217 /*do some cleanup: */
218 xfree((void**)&u_g_serial);
219 xfree((void**)&partition);
220 }
221 else if(strcmp(result->GetFieldName(),"p_g")==0){
222 /*easy, p_g is of size numberofnodes, on 1 dof, just repartition: */
223 result->GetField(&p_g);
224 VecToMPISerial(&p_g_serial,p_g);
225
226 if(!isstokes){
227 if(ismacayealpattyn){
228 fem_dh->parameters->FindParam((void*)&numberofnodes,"numberofnodes");
229 VecToMPISerial(&partition,fem_dh->partition);
230 }
231 else{
232 fem_dhu->parameters->FindParam((void*)&numberofnodes,"numberofnodes");
233 VecToMPISerial(&partition,fem_dhu->partition);
234 }
235 }
236 else{
237 fem_ds->parameters->FindParam((void*)&numberofnodes,"numberofnodes");
238 VecToMPISerial(&partition,fem_ds->partition);
239 }
240
241 pressure=(double*)xmalloc(numberofnodes*sizeof(double));
242
243 for(i=0;i<numberofnodes;i++){
244 pressure[i]=p_g_serial[(int)partition[i]];
245 }
246
247 /*Ok, add pressure,vy and vz to newresults: */
248 newresult=new Result(newresults->Size()+1,result->GetTime(),result->GetStep(),"pressure",pressure,numberofnodes);
249 newresults->AddObject(newresult);
250
251 /*do some cleanup: */
252 xfree((void**)&p_g_serial);
253 xfree((void**)&partition);
254 }
255 else if(strcmp(result->GetFieldName(),"t_g")==0){
256 /*easy, t_g is of size numberofnodes, on 1 dof, just repartition: */
257 result->GetField(&t_g);
258 VecToMPISerial(&t_g_serial,t_g);
259 fem_t->parameters->FindParam((void*)&numberofnodes,"numberofnodes");
260 VecToMPISerial(&partition,fem_t->partition);
261
262 temperature=(double*)xmalloc(numberofnodes*sizeof(double));
263
264 for(i=0;i<numberofnodes;i++){
265 temperature[i]=t_g_serial[(int)partition[i]];
266 }
267
268 /*Ok, add pressure to newresults: */
269 newresult=new Result(newresults->Size()+1,result->GetTime(),result->GetStep(),"temperature",temperature,numberofnodes);
270 newresults->AddObject(newresult);
271
272 /*do some cleanup: */
273 xfree((void**)&t_g_serial);
274 xfree((void**)&partition);
275 }
276 else if(strcmp(result->GetFieldName(),"m_g")==0){
277 /*easy, m_g is of size numberofnodes, on 1 dof, just repartition: */
278 result->GetField(&m_g);
279 VecToMPISerial(&m_g_serial,m_g);
280 fem_t->parameters->FindParam((void*)&numberofnodes,"numberofnodes");
281 fem_t->parameters->FindParam((void*)&yts,"yts");
282 VecToMPISerial(&partition,fem_t->partition);
283
284 melting=(double*)xmalloc(numberofnodes*sizeof(double));
285
286 for(i=0;i<numberofnodes;i++){
287 melting[i]=m_g_serial[(int)partition[i]]*yts;
288 }
289
290 /*Ok, add pressure to newresults: */
291 newresult=new Result(newresults->Size()+1,result->GetTime(),result->GetStep(),"melting",melting,numberofnodes);
292 newresults->AddObject(newresult);
293
294 /*do some cleanup: */
295 xfree((void**)&m_g_serial);
296 xfree((void**)&partition);
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**)&partition);
318 }
319 else if(strcmp(result->GetFieldName(),"param_g")==0){
320 /*easy, param_g is of size numberofnodes, on 1 dof, just repartition: */
321 result->GetField(&param_g);
322 fem_dh->parameters->FindParam((void*)&numberofnodes,"numberofnodes");
323 VecToMPISerial(&partition,fem_dh->partition);
324
325 parameter=(double*)xmalloc(numberofnodes*sizeof(double));
326
327 for(i=0;i<numberofnodes;i++){
328 parameter[i]=param_g[2*(int)partition[i]];
329 }
330
331 /*Ok, add parameter to newresults: */
332 newresult=new Result(newresults->Size()+1,result->GetTime(),result->GetStep(),"parameter",parameter,numberofnodes);
333 newresults->AddObject(newresult);
334
335 /*do some cleanup: */
336 xfree((void**)&partition);
337 }
338 else{
339 /*Just copy the result into the new results dataset: */
340 newresults->AddObject(result->copy());
341 }
342 }
343
344 /*Delete results: */
345 delete results;
346
347
348 /*Assign output pointers:*/
349 *presults=newresults;
350}
Note: See TracBrowser for help on using the repository browser.