source: issm/trunk/src/c/objects/Pengrid.cpp@ 2223

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

Added new steady state control test.
Added reset of stabilization in thermal_core_nonlinear, so that frozen penalties are unfrozen at each
step of a thermal analysis.
File size: 17.3 KB
Line 
1/*!\file Pengrid.c
2 * \brief: implementation of the Pengrid object
3 */
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 "stdio.h"
13#include "./Pengrid.h"
14#include <string.h>
15#include "../EnumDefinitions/EnumDefinitions.h"
16#include "../shared/shared.h"
17#include "../include/typedefs.h"
18
19
20Pengrid::Pengrid(){
21 return;
22}
23
24Pengrid::Pengrid(int pengrid_id, int pengrid_node_id,int pengrid_mparid, int pengrid_dof, int pengrid_active, double pengrid_penalty_offset,int pengrid_thermal_steadystate,int pengrid_stabilize_constraints){
25
26 id=pengrid_id;
27 node_id=pengrid_node_id;
28 mparid=pengrid_mparid;
29 dof=pengrid_dof;
30 active=pengrid_active;
31 penalty_offset =pengrid_penalty_offset;
32 thermal_steadystate=pengrid_thermal_steadystate;
33 stabilize_constraints=pengrid_stabilize_constraints;
34
35 node_offset=UNDEF;
36 node=NULL;
37 matpar=NULL;
38 matpar_offset=UNDEF;
39
40 zigzag_counter=0;
41
42 return;
43}
44
45Pengrid::~Pengrid(){
46 return;
47}
48
49void Pengrid::Echo(void){
50
51 printf("Pengrid:\n");
52 printf(" id: %i\n",id);
53 printf(" mparid: %i\n",mparid);
54 printf(" dof: %i\n",dof);
55 printf(" active: %i\n",active);
56 printf(" penalty_offset: %g\n",penalty_offset);
57 printf(" thermal_steadystate: %i\n",thermal_steadystate);
58 printf(" node_id: [%i]\n",node_id);
59 printf(" node_offset: [%i]\n",node_offset);
60 printf(" matpar_offset=%i\n",matpar_offset);
61
62 return;
63}
64void Pengrid::DeepEcho(void){
65
66 printf("Pengrid:\n");
67 printf(" id: %i\n",id);
68 printf(" mparid: %i\n",mparid);
69 printf(" dof: %i\n",dof);
70 printf(" active: %i\n",active);
71 printf(" penalty_offset: %g\n",penalty_offset);
72 printf(" thermal_steadystate: %i\n",thermal_steadystate);
73 printf(" node_id: [%i]\n",node_id);
74 printf(" node_offset: [%i]\n",node_offset);
75 printf(" matpar_offset=%i\n",matpar_offset);
76
77 if(node)node->Echo();
78 if(matpar)matpar->Echo();
79 return;
80}
81void Pengrid::Marshall(char** pmarshalled_dataset){
82
83 char* marshalled_dataset=NULL;
84 int enum_type=0;
85
86 /*recover marshalled_dataset: */
87 marshalled_dataset=*pmarshalled_dataset;
88
89 /*get enum type of Pengrid: */
90 enum_type=PengridEnum();
91
92 /*marshall enum: */
93 memcpy(marshalled_dataset,&enum_type,sizeof(enum_type));marshalled_dataset+=sizeof(enum_type);
94
95 /*marshall Pengrid data: */
96 memcpy(marshalled_dataset,&id,sizeof(id));marshalled_dataset+=sizeof(id);
97 memcpy(marshalled_dataset,&mparid,sizeof(mparid));marshalled_dataset+=sizeof(mparid);
98 memcpy(marshalled_dataset,&dof,sizeof(dof));marshalled_dataset+=sizeof(dof);
99 memcpy(marshalled_dataset,&active,sizeof(active));marshalled_dataset+=sizeof(active);
100 memcpy(marshalled_dataset,&penalty_offset,sizeof(penalty_offset));marshalled_dataset+=sizeof(penalty_offset);
101 memcpy(marshalled_dataset,&thermal_steadystate,sizeof(thermal_steadystate));marshalled_dataset+=sizeof(thermal_steadystate);
102 memcpy(marshalled_dataset,&node_id,sizeof(node_id));marshalled_dataset+=sizeof(node_id);
103 memcpy(marshalled_dataset,&node_offset,sizeof(node_offset));marshalled_dataset+=sizeof(node_offset);
104 memcpy(marshalled_dataset,&matpar,sizeof(matpar));marshalled_dataset+=sizeof(matpar);
105 memcpy(marshalled_dataset,&matpar_offset,sizeof(matpar_offset));marshalled_dataset+=sizeof(matpar_offset);
106
107 *pmarshalled_dataset=marshalled_dataset;
108 return;
109}
110
111int Pengrid::MarshallSize(){
112
113 return sizeof(id)+
114 sizeof(mparid)+
115 sizeof(dof)+
116 sizeof(active)+
117 sizeof(penalty_offset)+
118 sizeof(thermal_steadystate)+
119 sizeof(node_id)+
120 sizeof(node_offset)+
121 sizeof(matpar)+
122 sizeof(matpar_offset)+
123 sizeof(int); //sizeof(int) for enum type
124}
125
126char* Pengrid::GetName(void){
127 return "pengrid";
128}
129
130
131void Pengrid::Demarshall(char** pmarshalled_dataset){
132
133 char* marshalled_dataset=NULL;
134
135 /*recover marshalled_dataset: */
136 marshalled_dataset=*pmarshalled_dataset;
137
138 /*this time, no need to get enum type, the pointer directly points to the beginning of the
139 *object data (thanks to DataSet::Demarshall):*/
140
141 memcpy(&id,marshalled_dataset,sizeof(id));marshalled_dataset+=sizeof(id);
142 memcpy(&mparid,marshalled_dataset,sizeof(mparid));marshalled_dataset+=sizeof(mparid);
143 memcpy(&dof,marshalled_dataset,sizeof(dof));marshalled_dataset+=sizeof(dof);
144 memcpy(&active,marshalled_dataset,sizeof(active));marshalled_dataset+=sizeof(active);
145 memcpy(&penalty_offset,marshalled_dataset,sizeof(penalty_offset));marshalled_dataset+=sizeof(penalty_offset);
146 memcpy(&thermal_steadystate,marshalled_dataset,sizeof(thermal_steadystate));marshalled_dataset+=sizeof(thermal_steadystate);
147 memcpy(&node_id,marshalled_dataset,sizeof(node_id));marshalled_dataset+=sizeof(node_id);
148 memcpy(&node_offset,marshalled_dataset,sizeof(node_offset));marshalled_dataset+=sizeof(node_offset);
149 memcpy(&matpar,marshalled_dataset,sizeof(matpar));marshalled_dataset+=sizeof(matpar);
150 memcpy(&matpar_offset,marshalled_dataset,sizeof(matpar_offset));marshalled_dataset+=sizeof(matpar_offset);
151
152
153 node=NULL;
154 matpar=NULL;
155
156 /*return: */
157 *pmarshalled_dataset=marshalled_dataset;
158 return;
159}
160int Pengrid::Enum(void){
161
162 return PengridEnum();
163}
164
165int Pengrid::GetId(void){ return id; }
166
167int Pengrid::MyRank(void){
168 extern int my_rank;
169 return my_rank;
170}
171void Pengrid::DistributeNumDofs(int* numdofpernode,int analysis_type,int sub_analysis_type){return;}
172
173#undef __FUNCT__
174#define __FUNCT__ "Pengrid::Configure"
175
176void Pengrid::Configure(void* pelementsin,void* pnodesin,void* pmaterialsin){
177
178 DataSet* nodesin=NULL;
179 DataSet* materialsin=NULL;
180
181 /*Recover pointers :*/
182 nodesin=(DataSet*)pnodesin;
183 materialsin=(DataSet*)pmaterialsin;
184
185 /*Link this load with its nodes: */
186 ResolvePointers((Object**)&node,&node_id,&node_offset,1,nodesin);
187 ResolvePointers((Object**)&matpar,&mparid,&matpar_offset,1,materialsin);
188}
189
190
191#undef __FUNCT__
192#define __FUNCT__ "Pengrid::CreateKMatrix"
193
194void Pengrid::CreateKMatrix(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type){
195
196 /*No loads applied, do nothing: */
197 return;
198
199}
200
201#undef __FUNCT__
202#define __FUNCT__ "Pengrid::CreatePVector"
203void Pengrid::CreatePVector(Vec pg, void* inputs, int analysis_type,int sub_analysis_type){
204
205 /*No loads applied, do nothing: */
206 return;
207
208}
209#undef __FUNCT__
210#define __FUNCT__ "Pengrid::UpdateFromInputs"
211void Pengrid::UpdateFromInputs(void* inputs){
212
213}
214
215#undef __FUNCT__
216#define __FUNCT__ "Pengrid::PenaltyCreateKMatrix"
217void Pengrid::PenaltyCreateKMatrix(Mat Kgg,void* inputs,double kmax,int analysis_type,int sub_analysis_type){
218
219 if ((analysis_type==DiagnosticAnalysisEnum()) && ((sub_analysis_type==StokesAnalysisEnum()))){
220
221 PenaltyCreateKMatrixDiagnosticStokes( Kgg,inputs,kmax,analysis_type,sub_analysis_type);
222 }
223 else if (analysis_type==ThermalAnalysisEnum()){
224
225 PenaltyCreateKMatrixThermal( Kgg,inputs,kmax,analysis_type,sub_analysis_type);
226
227 }
228 else if (analysis_type==MeltingAnalysisEnum()){
229
230 PenaltyCreateKMatrixMelting( Kgg,inputs,kmax,analysis_type,sub_analysis_type);
231
232 }
233 else{
234 throw ErrorException(__FUNCT__,exprintf("%s%i%s%i%s","analysis: ",analysis_type," and sub_analysis_type: ",sub_analysis_type," not supported yet"));
235 }
236
237}
238
239#undef __FUNCT__
240#define __FUNCT__ "Pengrid::PenaltyCreateKMatrixDiagnosticStokes"
241void Pengrid::PenaltyCreateKMatrixDiagnosticStokes(Mat Kgg,void* vinputs,double kmax,int analysis_type,int sub_analysis_type){
242
243 const int numgrids=1;
244 const int NDOF4=4;
245 const int numdof=numgrids*NDOF4;
246 int doflist[numdof];
247 int numberofdofspernode;
248
249 int dofs1[1]={0};
250 int dofs2[1]={1};
251 double slope[2];
252 int found=0;
253 double Ke[4][4]={0.0};
254
255 ParameterInputs* inputs=NULL;
256
257 /*recover pointers: */
258 inputs=(ParameterInputs*)vinputs;
259
260 /*Get dof list: */
261 GetDofList(&doflist[0],&numberofdofspernode);
262
263 /*recover slope: */
264 found=inputs->Recover("bedslopex",&slope[0],1,dofs1,numgrids,(void**)&node);
265 if(!found)throw ErrorException(__FUNCT__," bedslopex needed in inputs!");
266 found=inputs->Recover("bedslopey",&slope[1],1,dofs2,numgrids,(void**)&node);
267 if(!found)throw ErrorException(__FUNCT__," bedslopey needed in inputs!");
268
269 //Create elementary matrix: add penalty to contrain wb (wb=ub*db/dx+vb*db/dy)
270 Ke[2][0]=-slope[0]*kmax*pow((double)10.0,penalty_offset);
271 Ke[2][1]=-slope[1]*kmax*pow((double)10.0,penalty_offset);
272 Ke[2][2]=kmax*pow((double)10,penalty_offset);
273
274 /*Add Ke to global matrix Kgg: */
275 MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke,ADD_VALUES);
276}
277
278#undef __FUNCT__
279#define __FUNCT__ "Pengrid::PenaltyCreateKMatrixThermal"
280void Pengrid::PenaltyCreateKMatrixThermal(Mat Kgg,void* vinputs,double kmax,int analysis_type,int sub_analysis_type){
281
282 int found=0;
283
284 const int numgrids=1;
285 const int NDOF1=1;
286 const int numdof=numgrids*NDOF1;
287 double Ke[numdof][numdof];
288 int doflist[numdof];
289 int numberofdofspernode;
290
291 ParameterInputs* inputs=NULL;
292
293 /*recover pointers: */
294 inputs=(ParameterInputs*)vinputs;
295
296
297 if(!active)return;
298
299 /*Get dof list: */
300 GetDofList(&doflist[0],&numberofdofspernode);
301
302 Ke[0][0]=kmax*pow((double)10,penalty_offset);
303
304 /*Add Ke to global matrix Kgg: */
305 MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke,ADD_VALUES);
306}
307
308#undef __FUNCT__
309#define __FUNCT__ "Pengrid::PenaltyCreateKMatrixMelting"
310void Pengrid::PenaltyCreateKMatrixMelting(Mat Kgg,void* vinputs,double kmax,int analysis_type,int sub_analysis_type){
311
312
313 int found=0;
314 const int numgrids=1;
315 const int NDOF1=1;
316 const int numdof=numgrids*NDOF1;
317 double Ke[numdof][numdof]={0.0};
318 int dofs1[1]={0};
319 int doflist[numdof];
320 int numberofdofspernode;
321 double meltingpoint;
322
323 double pressure;
324 double temperature;
325 double beta,t_pmp;
326
327 ParameterInputs* inputs=NULL;
328
329 /*check that pengrid is not a clone (penalty to be added only once)*/
330 if (node->IsClone()) return;
331
332 /*recover pointers: */
333 inputs=(ParameterInputs*)vinputs;
334
335 found=inputs->Recover("pressure",&pressure,1,dofs1,numgrids,(void**)&node);
336 if(!found)throw ErrorException(__FUNCT__," could not find pressure in inputs!");
337
338 found=inputs->Recover("temperature",&temperature,1,dofs1,numgrids,(void**)&node);
339 if(!found)throw ErrorException(__FUNCT__," could not find temperature in inputs!");
340
341 /*Get dof list: */
342 GetDofList(&doflist[0],&numberofdofspernode);
343
344 //Compute pressure melting point
345 meltingpoint=matpar->GetMeltingPoint();
346 beta=matpar->GetBeta();
347 t_pmp=meltingpoint-beta*pressure;
348
349 //Add penalty load
350 if (temperature<t_pmp){ //If T<Tpmp, there must be no melting. Therefore, melting should be constrained to 0 when T<Tpmp, instead of using spcs, use penalties
351 Ke[0][0]=kmax*pow((double)10,penalty_offset);
352 }
353
354 MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke,ADD_VALUES);
355}
356
357#undef __FUNCT__
358#define __FUNCT__ "Pengrid::PenaltyCreatePVector"
359void Pengrid::PenaltyCreatePVector(Vec pg,void* inputs,double kmax,int analysis_type,int sub_analysis_type){
360
361 if (analysis_type==ThermalAnalysisEnum()){
362
363 PenaltyCreatePVectorThermal( pg,inputs,kmax,analysis_type,sub_analysis_type);
364
365 }
366 else if (analysis_type==MeltingAnalysisEnum()){
367
368 PenaltyCreatePVectorMelting( pg,inputs,kmax,analysis_type,sub_analysis_type);
369
370 }
371 else if (analysis_type==DiagnosticAnalysisEnum()){
372
373 /*No loads applied, do nothing: */
374 return;
375
376 }
377 else{
378 throw ErrorException(__FUNCT__,exprintf("%s%i%s%i%s","analysis: ",analysis_type," and sub_analysis_type: ",sub_analysis_type," not supported yet"));
379 }
380
381}
382
383Object* Pengrid::copy() {
384 return new Pengrid(*this);
385}
386
387
388void Pengrid::GetDofList(int* doflist,int* pnumberofdofspernode){
389
390 int j;
391 int doflist_per_node[MAXDOFSPERNODE];
392 int numberofdofspernode;
393
394 node->GetDofList(&doflist_per_node[0],&numberofdofspernode);
395 for(j=0;j<numberofdofspernode;j++){
396 doflist[j]=doflist_per_node[j];
397 }
398
399 /*Assign output pointers:*/
400 *pnumberofdofspernode=numberofdofspernode;
401}
402
403void Pengrid::PenaltyCreatePVectorThermal(Vec pg, void* vinputs, double kmax,int analysis_type,int sub_analysis_type){
404
405 const int numgrids=1;
406 const int NDOF1=1;
407 const int numdof=numgrids*NDOF1;
408 int doflist[numdof];
409 double P_terms[numdof]={0.0};
410 int numberofdofspernode;
411 int found=0;
412 double pressure;
413 int dofs1[1]={0};
414 double meltingpoint;
415 double beta;
416 double t_pmp;
417
418 ParameterInputs* inputs=NULL;
419
420 /*recover pointers: */
421 inputs=(ParameterInputs*)vinputs;
422
423 if(!active)return;
424
425 /*Get dof list: */
426 GetDofList(&doflist[0],&numberofdofspernode);
427
428 //First recover pressure
429 found=inputs->Recover("pressure",&pressure,1,dofs1,numgrids,(void**)&node);
430 if(!found)throw ErrorException(__FUNCT__," could not find pressure in inputs!");
431
432 //Compute pressure melting point
433 meltingpoint=matpar->GetMeltingPoint();
434 beta=matpar->GetBeta();
435 t_pmp=meltingpoint-beta*pressure;
436
437 //Add penalty load
438 P_terms[0]=kmax*pow((double)10,penalty_offset)*t_pmp;
439
440 /*Add P_terms to global vector pg: */
441 VecSetValues(pg,numdof,doflist,(const double*)P_terms,ADD_VALUES);
442}
443
444void Pengrid::PenaltyCreatePVectorMelting(Vec pg, void* vinputs, double kmax,int analysis_type,int sub_analysis_type){
445
446 const int numgrids=1;
447 const int NDOF1=1;
448 const int numdof=numgrids*NDOF1;
449 int doflist[numdof];
450 double P_terms[numdof]={0.0};
451 int numberofdofspernode;
452 int found=0;
453 int dofs1[1]={0};
454 double pressure;
455 double temperature;
456 double melting_offset;
457 double meltingpoint;
458 double beta, heatcapacity;
459 double latentheat;
460 double t_pmp;
461 double dt;
462
463 ParameterInputs* inputs=NULL;
464
465 /*check that pengrid is not a clone (penalty to be added only once)*/
466 if (node->IsClone()) return;
467
468 /*recover pointers: */
469 inputs=(ParameterInputs*)vinputs;
470
471 /*Get dof list: */
472 GetDofList(&doflist[0],&numberofdofspernode);
473
474 //First recover pressure,melting offset and temperature vectors
475 found=inputs->Recover("pressure",&pressure,1,dofs1,numgrids,(void**)&node);
476 if(!found)throw ErrorException(__FUNCT__," could not find pressure in inputs!");
477
478 found=inputs->Recover("temperature",&temperature,1,dofs1,numgrids,(void**)&node);
479 if(!found)throw ErrorException(__FUNCT__," could not find temperature in inputs!");
480
481 found=inputs->Recover("melting_offset",&melting_offset);
482 if(!found)throw ErrorException(__FUNCT__," could not find melting_offset in inputs!");
483
484 found=inputs->Recover("dt",&dt);
485 if(!found)throw ErrorException(__FUNCT__," could not find dt in inputs!");
486
487 meltingpoint=matpar->GetMeltingPoint();
488 beta=matpar->GetBeta();
489 heatcapacity=matpar->GetHeatCapacity();
490 latentheat=matpar->GetLatentHeat();
491
492 //Compute pressure melting point
493 t_pmp=meltingpoint-beta*pressure;
494
495 //Add penalty load
496 //This time, the penalty must have the same value as the one used for the thermal computation
497 //so that the corresponding melting can be computed correctly
498 //In the thermal computation, we used kmax=melting_offset, and the same penalty_offset
499 if (temperature<t_pmp){ //%no melting
500 P_terms[0]=0;
501 }
502 else{
503 if (dt){
504 P_terms[0]=melting_offset*pow((double)10,penalty_offset)*(temperature-t_pmp)/dt;
505 }
506 else{
507 P_terms[0]=melting_offset*pow((double)10,penalty_offset)*(temperature-t_pmp);
508 }
509 }
510
511 /*Add P_terms to global vector pg: */
512 VecSetValues(pg,numdof,doflist,(const double*)P_terms,ADD_VALUES);
513}
514
515
516#undef __FUNCT__
517#define __FUNCT__ "Pengrid::PenaltyConstrain"
518void Pengrid::PenaltyConstrain(int* punstable,void* vinputs,int analysis_type,int sub_analysis_type){
519
520 if ((analysis_type==DiagnosticAnalysisEnum()) && ((sub_analysis_type==StokesAnalysisEnum()))){
521
522 /*No penalty to check*/
523 return;
524
525 }
526 else if (analysis_type==ThermalAnalysisEnum()){
527
528 PenaltyConstrainThermal(punstable,vinputs,analysis_type,sub_analysis_type);
529
530 }
531 else if (analysis_type==MeltingAnalysisEnum()){
532
533 /*No penalty to check*/
534 return;
535
536 }
537 else{
538 throw ErrorException(__FUNCT__,exprintf("%s%i%s%i%s","analysis: ",analysis_type," and sub_analysis_type: ",sub_analysis_type," not supported yet"));
539 }
540
541}
542
543#undef __FUNCT__
544#define __FUNCT__ "Pengrid::PenaltyConstrainThermal"
545void Pengrid::PenaltyConstrainThermal(int* punstable,void* vinputs,int analysis_type,int sub_analysis_type){
546
547 // The penalty is stable if it doesn't change during to successive iterations.
548
549 int found=0;
550 const int numgrids=1;
551
552
553 double pressure;
554 double temperature;
555 double beta,t_pmp;
556 double meltingpoint;
557 int new_active;
558 int dofs1[1]={0};
559 int unstable=0;
560 int reset_penalties=0;
561
562 ParameterInputs* inputs=NULL;
563
564 /*check that pengrid is not a clone (penalty to be added only once)*/
565 if (node->IsClone()){
566 unstable=0;
567 *punstable=unstable;
568 return;
569 }
570
571 /*recover pointers: */
572 inputs=(ParameterInputs*)vinputs;
573
574 //First recover beta, pressure and temperature vectors;
575 found=inputs->Recover("pressure",&pressure,1,dofs1,numgrids,(void**)&node);
576 if(!found)throw ErrorException(__FUNCT__," could not find pressure in inputs!");
577
578 found=inputs->Recover("temperature",&temperature,1,dofs1,numgrids,(void**)&node);
579 if(!found)throw ErrorException(__FUNCT__," could not find temperature in inputs!");
580
581 found=inputs->Recover("reset_penalties",&reset_penalties);
582
583 if(reset_penalties)zigzag_counter=0;
584
585 //Compute pressure melting point
586 meltingpoint=matpar->GetMeltingPoint();
587 beta=matpar->GetBeta();
588
589 t_pmp=meltingpoint-beta*pressure;
590
591 //Figure out if temperature is over melting_point, in which case, this penalty needs to be activated.
592
593 if (temperature>t_pmp){
594 new_active=1;
595 }
596 else{
597 new_active=0;
598 }
599
600
601 //Figure out stability of this penalty
602 if (active==new_active){
603 unstable=0;
604 }
605 else{
606 unstable=1;
607 if(stabilize_constraints)zigzag_counter++;
608 }
609
610 /*If penalty keeps zigzagging more than 5 times: */
611 if(stabilize_constraints){
612 if(zigzag_counter>stabilize_constraints){
613 unstable=0;
614 active=1;
615 }
616 }
617
618 //Set penalty flag
619 active=new_active;
620
621 //*Assign output pointers:*/
622 *punstable=unstable;
623}
Note: See TracBrowser for help on using the repository browser.