source: issm/trunk/src/c/ModelProcessorx/Model.cpp@ 1104

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

Transformed some Tria and Penta fields in booleans
File size: 14.9 KB
RevLine 
[1]1/*! \file Model.cpp
2 * \brief Model structure that mirrors the matlab workspace structure. Servers for the serial
3 * and parallel runs.
4 */
5
6#ifdef HAVE_CONFIG_H
7 #include "config.h"
8#else
9#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
10#endif
11
12#include "../shared/shared.h"
13#include "../io/io.h"
14#include "../include/globals.h"
15#include "../include/macros.h"
16
17#include <string.h>
18#include "stdio.h"
19
[117]20#include "./Model.h"
[1]21
22
23/*!--------------------------------------------------
24 NewModel
25 --------------------------------------------------*/
26
27Model* NewModel(void) {
28 /*! create a new Model object */
29 Model* model=NULL;
30
31 model=(Model*)xmalloc(sizeof(Model));
32
33 /*!initialize all pointers to 0: */
[962]34 model->name=NULL;
[1]35 model->repository=NULL;
36 model->meshtype=NULL;
37 model->analysis_type=NULL;
[465]38 model->sub_analysis_type=NULL;
[586]39 model->qmu_analysis=0;
40 model->solverstring=NULL;
41 model->numberofresponses=0;
[765]42 model->numberofvariables=0;
[962]43 model->qmu_npart=0;
[1]44 model->numberofelements=0;
45 model->numberofnodes=0;
46 model->x=NULL;
47 model->y=NULL;
48 model->z=NULL;
49 model->elements=NULL;
50 model->elements_type=NULL;
51 model->numberofnodes2d=0;
52 model->elements2d=NULL;
[88]53 model->deadgrids=NULL;
[1]54 model->numlayers=0;
[88]55 model->uppernodes=NULL;
[308]56 model->gridonhutter=NULL;
[653]57 model->gridonmacayeal=NULL;
58 model->gridonpattyn=NULL;
[1]59
60 model->vx_obs=NULL;
61 model->vy_obs=NULL;
62 model->vx=NULL;
63 model->vy=NULL;
64 model->vz=NULL;
65 model->pressure=NULL;
[575]66 model->temperature=NULL;
67 model->melting=NULL;
[1]68 model->elementonbed=NULL;
69 model->elementonsurface=NULL;
70 model->gridonbed=NULL;
71 model->gridonsurface=NULL;
[377]72 model->gridonstokes=NULL;
73 model->borderstokes=NULL;
[1]74 model->thickness=NULL;
75 model->surface=NULL;
76 model->bed=NULL;
77 model->elementoniceshelf=NULL;
[1104]78 model->elementonwater=NULL;
[387]79 model->gridonicesheet=NULL;
[816]80 model->gridoniceshelf=NULL;
[1]81
82 model->drag_type=0;
83 model->drag=NULL;
84 model->p=NULL;
[586]85 model->q=NULL;
[1]86
87
88 model->numberofsegs_diag=0;
[516]89 model->numberofsegs_diag_stokes=0;
[1]90 model->segmentonneumann_diag=NULL;
[387]91 model->segmentonneumann_diag_stokes=NULL;
[1]92 model-> neumannvalues_diag=NULL;
93 model-> gridondirichlet_diag=NULL;
94 model-> dirichletvalues_diag=NULL;
95 //prognostic
96 model-> segmentonneumann_prog=NULL;
97 model-> neumannvalues_prog=NULL;
98 model-> gridondirichlet_prog=NULL;
99 model-> dirichletvalues_prog=NULL;
100 //prognostic2
101 model-> segmentonneumann_prog2=NULL;
102 model-> neumannvalues_prog2=NULL;
103 //thermal
104 model-> gridondirichlet_thermal=NULL;
105 model-> dirichletvalues_thermal=NULL;
106 model-> geothermalflux=NULL;
107
108 /*!materials: */
109 model->rho_water=0;
110 model->rho_ice=0;
[586]111 model->g=0;
112 model->n=NULL;
[1]113 model->B=NULL;
114
115 /*!control methods: */
116 model->control_type=NULL;
117
118 /*!solution parameters: */
119 model->fit=NULL;
120 model->meanvel=0;
121 model->epsvel=0;
122 model->artificial_diffusivity=0;
123 model->nsteps=0;
124 model->tolx=0;
125 model->maxiter=NULL;
126 model->mincontrolconstraint=0;
127 model->maxcontrolconstraint=0;
128 model->debug=0;
129 model->plot=0;
130 model->eps_rel=0;
131 model->eps_abs=0;
132 model->dt=0;
133 model->ndt=0;
134 model->penalty_offset=0;
135 model->penalty_melting=0;
136 model->penalty_lock=0;
137 model->sparsity=0;
138 model->connectivity=0;
139 model->lowmem=0;
140 model->optscal=NULL;
141 model->yts=0;
142 model->viscosity_overshoot=0;
[387]143 model->stokesreconditioning=0;
[1]144 model->waitonlock=0;
145
146 /*!thermal parameters: */
147 model->beta=0;
148 model->meltingpoint=0;
149 model->latentheat=0;
150 model->heatcapacity=0;
151 model->thermalconductivity=0;
152 model->min_thermal_constraints=0;
153 model->mixed_layer_capacity=0;
154 model->thermal_exchange_velocity=0;
155
156
157 model->numrifts=0;
158 model->riftsnumpenaltypairs=NULL;
159 model->riftspenaltypairs=NULL;
160 model->riftsfill=NULL;
161 model->riftsfriction=NULL;
162
163 /*!penalties: */
164 model->numpenalties=0;
165 model->penalties=NULL;
166 model->penaltypartitioning=NULL;
167
168 /*!basal: */
169 model->melting=NULL;
170 model->accumulation=NULL;
171
[300]172 /*elements type: */
173 model->ishutter=0;
174 model->ismacayealpattyn=0;
175 model->isstokes=0;
176
[586]177
178 model->epart=NULL;
179 model->npart=NULL;
180 model->my_grids=NULL;
181 model->my_bordergrids=NULL;
182
[1]183 return model;
184}
185
186
187/*!--------------------------------------------------
188 DeleteModel
189 --------------------------------------------------*/
190
191void DeleteModel(Model** pmodel){
192
193 /*!Recover structure: */
194 Model* model = *pmodel;
195
196 int i;
197
198 /*!Two cases here:
199 * - serial mode: matlab's memory manager will take care of delete model when returning from Imp. Do nothing here, so as not to confuse
200 * the memory manager.
201 * - in parallel, anything the io layer does (FetchData) did needs to be erased explicitely in the model.
202 */
203
204 #ifdef _PARALLEL_
205 xfree((void**)&model->x);
206 xfree((void**)&model->y);
207 xfree((void**)&model->z);
208 xfree((void**)&model->elements);
209 xfree((void**)&model->elements_type);
[308]210 xfree((void**)&model->gridonhutter);
[653]211 xfree((void**)&model->gridonmacayeal);
[1]212 if (strcmp(model->meshtype,"3d")==0){
213 xfree((void**)&model->elements2d);
214 xfree((void**)&model->deadgrids);
[88]215 xfree((void**)&model->uppernodes);
[653]216 xfree((void**)&model->gridonpattyn);
[1]217 }
218 xfree((void**)&model->solverstring);
219 xfree((void**)&model->elementonbed);
220 xfree((void**)&model->elementonsurface);
221 xfree((void**)&model->gridonbed);
222 xfree((void**)&model->gridonsurface);
[377]223 xfree((void**)&model->gridonstokes);
224 xfree((void**)&model->borderstokes);
[1]225 xfree((void**)&model->thickness);
226 xfree((void**)&model->surface);
227 xfree((void**)&model->bed);
228 xfree((void**)&model->vx_obs);
229 xfree((void**)&model->vy_obs);
230 xfree((void**)&model->vx);
231 xfree((void**)&model->vy);
232 xfree((void**)&model->vz);
233 xfree((void**)&model->pressure);
[575]234 xfree((void**)&model->temperature);
235 xfree((void**)&model->melting);
[1]236 xfree((void**)&model->drag);
237 xfree((void**)&model->p);
238 xfree((void**)&model->q);
239 xfree((void**)&model->elementoniceshelf);
[1104]240 xfree((void**)&model->elementonwater);
[387]241 xfree((void**)&model->gridonicesheet);
[816]242 xfree((void**)&model->gridoniceshelf);
[1]243 xfree((void**)&model->segmentonneumann_diag);
[387]244 xfree((void**)&model->segmentonneumann_diag_stokes);
[1]245 xfree((void**)&model->neumannvalues_diag);
246 xfree((void**)&model->gridondirichlet_diag);
247 xfree((void**)&model->dirichletvalues_diag);
248 xfree((void**)&model->segmentonneumann_prog);
249 xfree((void**)&model->neumannvalues_prog);
250 xfree((void**)&model->gridondirichlet_prog);
251 xfree((void**)&model->dirichletvalues_prog);
252 xfree((void**)&model->segmentonneumann_prog2);
253 xfree((void**)&model->neumannvalues_prog2);
254 xfree((void**)&model->gridondirichlet_thermal);
255 xfree((void**)&model->dirichletvalues_thermal);
256 xfree((void**)&model->geothermalflux);
257 xfree((void**)&model->melting);
258 xfree((void**)&model->accumulation);
259 xfree((void**)&model->B);
260 xfree((void**)&model->n);
261 xfree((void**)&model->fit);
262 xfree((void**)&model->optscal);
263 xfree((void**)&model->maxiter);
264
265
266 /*!Delete structure fields: */
267 xfree((void**)&model->repository);
268 xfree((void**)&model->meshtype);
[962]269 xfree((void**)&model->name);
[1]270 xfree((void**)&model->analysis_type);
[465]271 xfree((void**)&model->sub_analysis_type);
[1]272
273 if(model->numrifts){
274 for(i=0;i<model->numrifts;i++){
275 double* riftpenaltypairs=model->riftspenaltypairs[i];
276 xfree((void**)&riftpenaltypairs);
277 }
278 xfree((void**)&model->riftspenaltypairs);
279 xfree((void**)&model->riftsnumpenaltypairs);
280 xfree((void**)&model->riftsfill);
281 xfree((void**)&model->riftsfriction);
282 }
283
284 xfree((void**)&model->penalties);
285 xfree((void**)&model->penaltypartitioning);
286
287 xfree((void**)&model->control_type);
288
[387]289 xfree((void**)&model->epart);
290 xfree((void**)&model->npart);
291 xfree((void**)&model->my_grids);
292 xfree((void**)&model->my_bordergrids);
293
[1]294 /*!Delete entire structure: */
295 xfree((void**)pmodel);
296 #endif
297}
298
299/*!--------------------------------------------------
300 ModelInit
301 --------------------------------------------------*/
302
303#undef __FUNCT__
304#define __FUNCT__ "ModelInit"
305
306int ModelInit(Model** pmodel,ConstDataHandle model_handle){
307
308 int i,j;
309
310 /*output: */
311 Model* model=NULL;
312
313 /*Allocate model: */
314 model=NewModel();
315
316 /*In ModelInit, we get all the data that is not difficult to get, and that is small: */
[586]317
[962]318 ModelFetchData((void**)&model->name,NULL,NULL,model_handle,"name","String",NULL);
[1]319 ModelFetchData((void**)&model->analysis_type,NULL,NULL,model_handle,"analysis_type","String",NULL);
[465]320 ModelFetchData((void**)&model->sub_analysis_type,NULL,NULL,model_handle,"sub_analysis_type","String",NULL);
[586]321 ModelFetchData((void**)&model->qmu_analysis,NULL,NULL,model_handle,"qmu_analysis","Integer",NULL);
[1]322
323 ModelFetchData((void**)&model->meshtype,NULL,NULL,model_handle,"type","String",NULL);
324 /*!Get numberofelements and numberofnodes: */
325 ModelFetchData((void**)&model->numberofnodes,NULL,NULL,model_handle,"numberofgrids","Integer",NULL);
326 ModelFetchData((void**)&model->numberofelements,NULL,NULL,model_handle,"numberofelements","Integer",NULL);
327 /*!In case we are running 3d, we are going to need the collapsed and non-collapsed 2d meshes, from which the 3d mesh was extruded: */
328 if (strcmp(model->meshtype,"3d")==0){
329
330 /*!Deal with 2d mesh: */
331 ModelFetchData((void**)&model->numberofelements2d,NULL,NULL,model_handle,"numberofelements2d","Integer",NULL);
332 ModelFetchData((void**)&model->numberofnodes2d,NULL,NULL,model_handle,"numberofgrids2d","Integer",NULL);
333 ModelFetchData((void**)&model->numlayers,NULL,NULL,model_handle,"numlayers","Integer",NULL);
334 }
[300]335
[586]336
[300]337 /*elements type: */
338 ModelFetchData((void**)&model->ishutter,NULL,NULL,model_handle,"ishutter","Integer",NULL);
339 ModelFetchData((void**)&model->ismacayealpattyn,NULL,NULL,model_handle,"ismacayealpattyn","Integer",NULL);
[394]340 ModelFetchData((void**)&model->isstokes,NULL,NULL,model_handle,"isstokes","Integer",NULL);
[300]341
[1]342 /*!Get drag_type, drag and p,q: */
343 ModelFetchData((void**)&model->drag_type,NULL,NULL,model_handle,"drag_type","Integer",NULL);
344
345 /*!Get materials: */
346 ModelFetchData((void**)&model->rho_water,NULL,NULL,model_handle,"rho_water","Scalar",NULL);
347 ModelFetchData((void**)&model->rho_ice,NULL,NULL,model_handle,"rho_ice","Scalar",NULL);
348 ModelFetchData((void**)&model->g,NULL,NULL,model_handle,"g","Scalar",NULL);
349
350 /*Get control parameters: */
351 ModelFetchData((void**)&model->control_type,NULL,NULL,model_handle,"control_type","String",NULL);
352
353 /*!Get solution parameters: */
354 ModelFetchData((void**)&model->yts,NULL,NULL,model_handle,"yts","Scalar",NULL);
355 ModelFetchData((void**)&model->meanvel,NULL,NULL,model_handle,"meanvel","Scalar",NULL);
356 ModelFetchData((void**)&model->epsvel,NULL,NULL,model_handle,"epsvel","Scalar",NULL);
357 ModelFetchData((void**)&model->debug,NULL,NULL,model_handle,"debug","Integer",NULL);
358 ModelFetchData((void**)&model->plot,NULL,NULL,model_handle,"plot","Integer",NULL);
359 ModelFetchData((void**)&model->artificial_diffusivity,NULL,NULL,model_handle,"artificial_diffusivity","Integer",NULL);
360 ModelFetchData((void**)&model->nsteps,NULL,NULL,model_handle,"nsteps","Integer",NULL);
361 ModelFetchData((void**)&model->tolx,NULL,NULL,model_handle,"tolx","Scalar",NULL);
362 ModelFetchData((void**)&model->mincontrolconstraint,NULL,NULL,model_handle,"mincontrolconstraint","Scalar",NULL);
363 ModelFetchData((void**)&model->maxcontrolconstraint,NULL,NULL,model_handle,"maxcontrolconstraint","Scalar",NULL);
364 ModelFetchData((void**)&model->eps_rel,NULL,NULL,model_handle,"eps_rel","Scalar",NULL);
365 ModelFetchData((void**)&model->eps_abs,NULL,NULL,model_handle,"eps_abs","Scalar",NULL);
366 ModelFetchData((void**)&model->dt,NULL,NULL,model_handle,"dt","Scalar",NULL);
367 ModelFetchData((void**)&model->ndt,NULL,NULL,model_handle,"ndt","Scalar",NULL);
368 ModelFetchData((void**)&model->penalty_offset,NULL,NULL,model_handle,"penalty_offset","Scalar",NULL);
369 ModelFetchData((void**)&model->penalty_melting,NULL,NULL,model_handle,"penalty_melting","Scalar",NULL);
370 ModelFetchData((void**)&model->penalty_lock,NULL,NULL,model_handle,"penalty_lock","Integer",NULL);
371 ModelFetchData((void**)&model->sparsity,NULL,NULL,model_handle,"sparsity","Scalar",NULL);
372 ModelFetchData((void**)&model->connectivity,NULL,NULL,model_handle,"connectivity","Integer",NULL);
373 ModelFetchData((void**)&model->lowmem,NULL,NULL,model_handle,"lowmem","Integer",NULL);
374 ModelFetchData((void**)&model->solverstring,NULL,NULL,model_handle,"solverstring","String",NULL);
375 ModelFetchData((void**)&model->viscosity_overshoot,NULL,NULL,model_handle,"viscosity_overshoot","Scalar",NULL);
[387]376 ModelFetchData((void**)&model->stokesreconditioning,NULL,NULL,model_handle,"stokesreconditioning","Scalar",NULL);
[1]377 ModelFetchData((void**)&model->waitonlock,NULL,NULL,model_handle,"waitonlock","Integer",NULL);
378
379 /*!Get thermal parameters: */
380 ModelFetchData((void**)&model->beta,NULL,NULL,model_handle,"beta","Scalar",NULL);
381 ModelFetchData((void**)&model->meltingpoint,NULL,NULL,model_handle,"meltingpoint","Scalar",NULL);
382 ModelFetchData((void**)&model->latentheat,NULL,NULL,model_handle,"latentheat","Scalar",NULL);
383 ModelFetchData((void**)&model->heatcapacity,NULL,NULL,model_handle,"heatcapacity","Scalar",NULL);
384 ModelFetchData((void**)&model->thermalconductivity,NULL,NULL,model_handle,"thermalconductivity","Scalar",NULL);
385 ModelFetchData((void**)&model->min_thermal_constraints,NULL,NULL,model_handle,"min_thermal_constraints","Integer",NULL);
386 ModelFetchData((void**)&model->mixed_layer_capacity,NULL,NULL,model_handle,"mixed_layer_capacity","Scalar",NULL);
387 ModelFetchData((void**)&model->thermal_exchange_velocity,NULL,NULL,model_handle,"thermal_exchange_velocity","Scalar",NULL);
[358]388
[1]389 /*rifts: */
390 ModelFetchData((void**)&model->numrifts,NULL,NULL,model_handle,"numrifts","Integer",NULL);
391
[586]392 /*qmu: */
393 if(model->qmu_analysis){
[765]394 ModelFetchData((void**)&model->numberofvariables,NULL,NULL,model_handle,"numberofvariables","Integer",NULL);
[586]395 ModelFetchData((void**)&model->numberofresponses,NULL,NULL,model_handle,"numberofresponses","Integer",NULL);
396 ModelFetchData((void**)&model->qmu_npart,NULL,NULL,model_handle,"npart","Integer",NULL);
397 }
398
[1]399 /*Assign output pointers: */
400 *pmodel=model;
401
402 return 1;
403}
404
405/*!--------------------------------------------------
406 ModelEcho
407 --------------------------------------------------*/
408void ModelEcho(Model* model,int which_part,int rank) {
409
410 //which_part determines what gets echoed, otherwise, we'll get too much output.
411 //1-> penalties
412
413 int i,j;
414
415 if(which_part==1 && my_rank==rank && (strcmp(model->meshtype,"3d")==0)){
416 printf("Model penalties: \n");
417 printf(" number of penalties: %i\n",model->numpenalties);
418 printf(" grids: \n");
419
420 for(i=0;i<model->numpenalties;i++){
421 for(j=0;j<model->numlayers;j++){
422 printf("%i ",(int)*(model->penalties+model->numlayers*i+j));
423 }
424 printf("\n");
425 }
426 }
427
428 if(which_part==2 && my_rank==rank){
429 printf("Model rifts: \n");
430 printf(" number of rifts: %i\n",model->numrifts);
431 for(i=0;i<model->numrifts;i++){
432 double* penaltypairs=model->riftspenaltypairs[i];
433 printf(" rift #%i\n",i);
434 for (j=0;j<model->riftsnumpenaltypairs[i];j++){
435 printf(" grids %g %g elements %g %g normal [%g,%g] length %g\n",*(penaltypairs+7*j+0),*(penaltypairs+7*j+1),*(penaltypairs+7*j+2),*(penaltypairs+7*j+3),
436 *(penaltypairs+7*j+4),*(penaltypairs+7*j+5),*(penaltypairs+7*j+6));
437 }
438 printf(" friction %g fill %i\n",model->riftsfriction[i],model->riftsfill[i]);
439 }
440 }
441 cleanup_and_return:
442 return;
443}
Note: See TracBrowser for help on using the repository browser.