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

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

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