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

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

bug typo stokes_reconditionning -> reconditionning_number
File size: 13.1 KB
Line 
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
20#include "./Model.h"
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: */
34
35
36 model->repository=NULL;
37 model->meshtype=NULL;
38 model->analysis_type=NULL;
39 model->numberofelements=0;
40 model->numberofnodes=0;
41 model->x=NULL;
42 model->y=NULL;
43 model->z=NULL;
44 model->elements=NULL;
45 model->elements_type=NULL;
46 model->numberofnodes2d=0;
47 model->elements2d=NULL;
48 model->deadgrids=NULL;
49 model->numlayers=0;
50 model->uppernodes=NULL;
51 model->gridonhutter=NULL;
52
53 model->vx_obs=NULL;
54 model->vy_obs=NULL;
55 model->vx=NULL;
56 model->vy=NULL;
57 model->vz=NULL;
58 model->pressure=NULL;
59 model->elementonbed=NULL;
60 model->elementonsurface=NULL;
61 model->gridonbed=NULL;
62 model->gridonsurface=NULL;
63 model->gridonstokes=NULL;
64 model->borderstokes=NULL;
65 model->thickness=NULL;
66 model->surface=NULL;
67 model->bed=NULL;
68 model->elementoniceshelf=NULL;
69
70 model->drag_type=0;
71 model->drag=NULL;
72 model->p=NULL;
73
74
75 model->numberofsegs_diag=0;
76 model->segmentonneumann_diag=NULL;
77 model-> neumannvalues_diag=NULL;
78 model-> gridondirichlet_diag=NULL;
79 model-> dirichletvalues_diag=NULL;
80 //prognostic
81 model-> segmentonneumann_prog=NULL;
82 model-> neumannvalues_prog=NULL;
83 model-> gridondirichlet_prog=NULL;
84 model-> dirichletvalues_prog=NULL;
85 //prognostic2
86 model-> segmentonneumann_prog2=NULL;
87 model-> neumannvalues_prog2=NULL;
88 //thermal
89 model-> gridondirichlet_thermal=NULL;
90 model-> dirichletvalues_thermal=NULL;
91 model-> geothermalflux=NULL;
92
93 /*!materials: */
94 model->rho_water=0;
95 model->rho_ice=0;
96 model->g=0;
97 model->n=0;
98 model->B=NULL;
99
100 /*!control methods: */
101 model->control_type=NULL;
102
103 /*!solution parameters: */
104 model->fit=NULL;
105 model->meanvel=0;
106 model->epsvel=0;
107 model->artificial_diffusivity=0;
108 model->nsteps=0;
109 model->tolx=0;
110 model->maxiter=NULL;
111 model->mincontrolconstraint=0;
112 model->maxcontrolconstraint=0;
113 model->debug=0;
114 model->plot=0;
115 model->eps_rel=0;
116 model->eps_abs=0;
117 model->dt=0;
118 model->ndt=0;
119 model->penalty_offset=0;
120 model->penalty_melting=0;
121 model->penalty_lock=0;
122 model->sparsity=0;
123 model->connectivity=0;
124 model->lowmem=0;
125 model->optscal=NULL;
126 model->yts=0;
127 model->viscosity_overshoot=0;
128 model->reconditioning_number=0;
129 model->waitonlock=0;
130
131 /*!thermal parameters: */
132 model->beta=0;
133 model->meltingpoint=0;
134 model->latentheat=0;
135 model->heatcapacity=0;
136 model->thermalconductivity=0;
137 model->min_thermal_constraints=0;
138 model->mixed_layer_capacity=0;
139 model->thermal_exchange_velocity=0;
140
141
142 model->numrifts=0;
143 model->riftsnumpenaltypairs=NULL;
144 model->riftspenaltypairs=NULL;
145 model->riftsfill=NULL;
146 model->riftsfriction=NULL;
147
148 /*!penalties: */
149 model->numpenalties=0;
150 model->penalties=NULL;
151 model->penaltypartitioning=NULL;
152
153 /*!basal: */
154 model->melting=NULL;
155 model->accumulation=NULL;
156
157 /*elements type: */
158 model->ishutter=0;
159 model->ismacayealpattyn=0;
160 model->isstokes=0;
161
162 return model;
163}
164
165
166/*!--------------------------------------------------
167 DeleteModel
168 --------------------------------------------------*/
169
170void DeleteModel(Model** pmodel){
171
172 /*!Recover structure: */
173 Model* model = *pmodel;
174
175 int i;
176
177 /*!Two cases here:
178 * - serial mode: matlab's memory manager will take care of delete model when returning from Imp. Do nothing here, so as not to confuse
179 * the memory manager.
180 * - in parallel, anything the io layer does (FetchData) did needs to be erased explicitely in the model.
181 */
182
183 #ifdef _PARALLEL_
184 xfree((void**)&model->x);
185 xfree((void**)&model->y);
186 xfree((void**)&model->z);
187 xfree((void**)&model->elements);
188 xfree((void**)&model->elements_type);
189 xfree((void**)&model->gridonhutter);
190 if (strcmp(model->meshtype,"3d")==0){
191 xfree((void**)&model->elements2d);
192 xfree((void**)&model->deadgrids);
193 xfree((void**)&model->uppernodes);
194 }
195 xfree((void**)&model->solverstring);
196 xfree((void**)&model->elementonbed);
197 xfree((void**)&model->elementonsurface);
198 xfree((void**)&model->gridonbed);
199 xfree((void**)&model->gridonsurface);
200 xfree((void**)&model->gridonstokes);
201 xfree((void**)&model->borderstokes);
202 xfree((void**)&model->thickness);
203 xfree((void**)&model->surface);
204 xfree((void**)&model->bed);
205 xfree((void**)&model->vx_obs);
206 xfree((void**)&model->vy_obs);
207 xfree((void**)&model->vx);
208 xfree((void**)&model->vy);
209 xfree((void**)&model->vz);
210 xfree((void**)&model->pressure);
211 xfree((void**)&model->drag);
212 xfree((void**)&model->p);
213 xfree((void**)&model->q);
214 xfree((void**)&model->elementoniceshelf);
215 xfree((void**)&model->segmentonneumann_diag);
216 xfree((void**)&model->neumannvalues_diag);
217 xfree((void**)&model->gridondirichlet_diag);
218 xfree((void**)&model->dirichletvalues_diag);
219 xfree((void**)&model->segmentonneumann_prog);
220 xfree((void**)&model->neumannvalues_prog);
221 xfree((void**)&model->gridondirichlet_prog);
222 xfree((void**)&model->dirichletvalues_prog);
223 xfree((void**)&model->segmentonneumann_prog2);
224 xfree((void**)&model->neumannvalues_prog2);
225 xfree((void**)&model->gridondirichlet_thermal);
226 xfree((void**)&model->dirichletvalues_thermal);
227 xfree((void**)&model->geothermalflux);
228 xfree((void**)&model->melting);
229 xfree((void**)&model->accumulation);
230 xfree((void**)&model->B);
231 xfree((void**)&model->n);
232 xfree((void**)&model->fit);
233 xfree((void**)&model->optscal);
234 xfree((void**)&model->maxiter);
235
236
237 /*!Delete structure fields: */
238 xfree((void**)&model->repository);
239 xfree((void**)&model->meshtype);
240 xfree((void**)&model->analysis_type);
241
242 if(model->numrifts){
243 for(i=0;i<model->numrifts;i++){
244 double* riftpenaltypairs=model->riftspenaltypairs[i];
245 xfree((void**)&riftpenaltypairs);
246 }
247 xfree((void**)&model->riftspenaltypairs);
248 xfree((void**)&model->riftsnumpenaltypairs);
249 xfree((void**)&model->riftsfill);
250 xfree((void**)&model->riftsfriction);
251 }
252
253 xfree((void**)&model->penalties);
254 xfree((void**)&model->penaltypartitioning);
255
256 xfree((void**)&model->control_type);
257
258 /*!Delete entire structure: */
259 xfree((void**)pmodel);
260 #endif
261}
262
263/*!--------------------------------------------------
264 ModelInit
265 --------------------------------------------------*/
266
267#undef __FUNCT__
268#define __FUNCT__ "ModelInit"
269
270int ModelInit(Model** pmodel,ConstDataHandle model_handle){
271
272 int i,j;
273
274 /*output: */
275 Model* model=NULL;
276
277 /*Allocate model: */
278 model=NewModel();
279
280 /*In ModelInit, we get all the data that is not difficult to get, and that is small: */
281 ModelFetchData((void**)&model->analysis_type,NULL,NULL,model_handle,"analysis_type","String",NULL);
282
283 ModelFetchData((void**)&model->meshtype,NULL,NULL,model_handle,"type","String",NULL);
284 /*!Get numberofelements and numberofnodes: */
285 ModelFetchData((void**)&model->numberofnodes,NULL,NULL,model_handle,"numberofgrids","Integer",NULL);
286 ModelFetchData((void**)&model->numberofelements,NULL,NULL,model_handle,"numberofelements","Integer",NULL);
287
288 /*!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: */
289 if (strcmp(model->meshtype,"3d")==0){
290
291 /*!Deal with 2d mesh: */
292 ModelFetchData((void**)&model->numberofelements2d,NULL,NULL,model_handle,"numberofelements2d","Integer",NULL);
293 ModelFetchData((void**)&model->numberofnodes2d,NULL,NULL,model_handle,"numberofgrids2d","Integer",NULL);
294 ModelFetchData((void**)&model->numlayers,NULL,NULL,model_handle,"numlayers","Integer",NULL);
295 }
296
297 /*elements type: */
298 ModelFetchData((void**)&model->ishutter,NULL,NULL,model_handle,"ishutter","Integer",NULL);
299 ModelFetchData((void**)&model->ismacayealpattyn,NULL,NULL,model_handle,"ismacayealpattyn","Integer",NULL);
300
301 /*!Get drag_type, drag and p,q: */
302 ModelFetchData((void**)&model->drag_type,NULL,NULL,model_handle,"drag_type","Integer",NULL);
303
304 /*!Get materials: */
305 ModelFetchData((void**)&model->rho_water,NULL,NULL,model_handle,"rho_water","Scalar",NULL);
306 ModelFetchData((void**)&model->rho_ice,NULL,NULL,model_handle,"rho_ice","Scalar",NULL);
307 ModelFetchData((void**)&model->g,NULL,NULL,model_handle,"g","Scalar",NULL);
308
309 /*Get control parameters: */
310 ModelFetchData((void**)&model->control_type,NULL,NULL,model_handle,"control_type","String",NULL);
311
312 /*!Get solution parameters: */
313 ModelFetchData((void**)&model->yts,NULL,NULL,model_handle,"yts","Scalar",NULL);
314 ModelFetchData((void**)&model->meanvel,NULL,NULL,model_handle,"meanvel","Scalar",NULL);
315 ModelFetchData((void**)&model->epsvel,NULL,NULL,model_handle,"epsvel","Scalar",NULL);
316 ModelFetchData((void**)&model->debug,NULL,NULL,model_handle,"debug","Integer",NULL);
317 ModelFetchData((void**)&model->plot,NULL,NULL,model_handle,"plot","Integer",NULL);
318 ModelFetchData((void**)&model->artificial_diffusivity,NULL,NULL,model_handle,"artificial_diffusivity","Integer",NULL);
319 ModelFetchData((void**)&model->nsteps,NULL,NULL,model_handle,"nsteps","Integer",NULL);
320 ModelFetchData((void**)&model->tolx,NULL,NULL,model_handle,"tolx","Scalar",NULL);
321 ModelFetchData((void**)&model->mincontrolconstraint,NULL,NULL,model_handle,"mincontrolconstraint","Scalar",NULL);
322 ModelFetchData((void**)&model->maxcontrolconstraint,NULL,NULL,model_handle,"maxcontrolconstraint","Scalar",NULL);
323 ModelFetchData((void**)&model->eps_rel,NULL,NULL,model_handle,"eps_rel","Scalar",NULL);
324 ModelFetchData((void**)&model->eps_abs,NULL,NULL,model_handle,"eps_abs","Scalar",NULL);
325 ModelFetchData((void**)&model->dt,NULL,NULL,model_handle,"dt","Scalar",NULL);
326 ModelFetchData((void**)&model->ndt,NULL,NULL,model_handle,"ndt","Scalar",NULL);
327 ModelFetchData((void**)&model->penalty_offset,NULL,NULL,model_handle,"penalty_offset","Scalar",NULL);
328 ModelFetchData((void**)&model->penalty_melting,NULL,NULL,model_handle,"penalty_melting","Scalar",NULL);
329 ModelFetchData((void**)&model->penalty_lock,NULL,NULL,model_handle,"penalty_lock","Integer",NULL);
330 ModelFetchData((void**)&model->sparsity,NULL,NULL,model_handle,"sparsity","Scalar",NULL);
331 ModelFetchData((void**)&model->connectivity,NULL,NULL,model_handle,"connectivity","Integer",NULL);
332 ModelFetchData((void**)&model->lowmem,NULL,NULL,model_handle,"lowmem","Integer",NULL);
333 ModelFetchData((void**)&model->solverstring,NULL,NULL,model_handle,"solverstring","String",NULL);
334 ModelFetchData((void**)&model->viscosity_overshoot,NULL,NULL,model_handle,"viscosity_overshoot","Scalar",NULL);
335 ModelFetchData((void**)&model->reconditioning_number,NULL,NULL,model_handle,"reconditioning_number","Scalar",NULL);
336 ModelFetchData((void**)&model->waitonlock,NULL,NULL,model_handle,"waitonlock","Integer",NULL);
337
338 /*!Get thermal parameters: */
339 ModelFetchData((void**)&model->beta,NULL,NULL,model_handle,"beta","Scalar",NULL);
340 ModelFetchData((void**)&model->meltingpoint,NULL,NULL,model_handle,"meltingpoint","Scalar",NULL);
341 ModelFetchData((void**)&model->latentheat,NULL,NULL,model_handle,"latentheat","Scalar",NULL);
342 ModelFetchData((void**)&model->heatcapacity,NULL,NULL,model_handle,"heatcapacity","Scalar",NULL);
343 ModelFetchData((void**)&model->thermalconductivity,NULL,NULL,model_handle,"thermalconductivity","Scalar",NULL);
344 ModelFetchData((void**)&model->min_thermal_constraints,NULL,NULL,model_handle,"min_thermal_constraints","Integer",NULL);
345 ModelFetchData((void**)&model->mixed_layer_capacity,NULL,NULL,model_handle,"mixed_layer_capacity","Scalar",NULL);
346 ModelFetchData((void**)&model->thermal_exchange_velocity,NULL,NULL,model_handle,"thermal_exchange_velocity","Scalar",NULL);
347
348 /*rifts: */
349 ModelFetchData((void**)&model->numrifts,NULL,NULL,model_handle,"numrifts","Integer",NULL);
350
351 /*Assign output pointers: */
352 *pmodel=model;
353
354 return 1;
355}
356
357/*!--------------------------------------------------
358 ModelEcho
359 --------------------------------------------------*/
360void ModelEcho(Model* model,int which_part,int rank) {
361
362 //which_part determines what gets echoed, otherwise, we'll get too much output.
363 //1-> penalties
364
365 int i,j;
366
367 if(which_part==1 && my_rank==rank && (strcmp(model->meshtype,"3d")==0)){
368 printf("Model penalties: \n");
369 printf(" number of penalties: %i\n",model->numpenalties);
370 printf(" grids: \n");
371
372 for(i=0;i<model->numpenalties;i++){
373 for(j=0;j<model->numlayers;j++){
374 printf("%i ",(int)*(model->penalties+model->numlayers*i+j));
375 }
376 printf("\n");
377 }
378 }
379
380 if(which_part==2 && my_rank==rank){
381 printf("Model rifts: \n");
382 printf(" number of rifts: %i\n",model->numrifts);
383 for(i=0;i<model->numrifts;i++){
384 double* penaltypairs=model->riftspenaltypairs[i];
385 printf(" rift #%i\n",i);
386 for (j=0;j<model->riftsnumpenaltypairs[i];j++){
387 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),
388 *(penaltypairs+7*j+4),*(penaltypairs+7*j+5),*(penaltypairs+7*j+6));
389 }
390 printf(" friction %g fill %i\n",model->riftsfriction[i],model->riftsfill[i]);
391 }
392 }
393 cleanup_and_return:
394 return;
395}
Note: See TracBrowser for help on using the repository browser.