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

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

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