source: issm/trunk/src/c/modules/GroundinglineMigrationx/GroundinglineMigrationxUtils.cpp@ 10300

Last change on this file since 10300 was 10300, checked in by seroussi, 13 years ago

minor grounding line
File size: 11.2 KB
Line 
1/*!\file GroundinglineMigrationx
2 * \brief: migration grounding line position.
3 */
4
5#include "./GroundinglineMigrationx.h"
6#include "./GroundinglineMigrationxLocal.h"
7#include "../../shared/shared.h"
8#include "../../io/io.h"
9#include "../../io/io.h"
10#include "../../include/include.h"
11#include "../../toolkits/toolkits.h"
12#include "../../EnumDefinitions/EnumDefinitions.h"
13
14/*FUNCTION AgressiveMigration{{{1*/
15void AgressiveMigration(Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads,Materials* materials, Parameters* parameters){
16 /*Here, whatever nodes inside the ice sheet want to unground, we allow -> instantaneous transmission of water through the bedrock: */
17
18 int i;
19 Element* element=NULL;
20
21 _printf_(VerboseModule()," Migrating grounding line\n");
22
23 /*Carry out grounding line migration for those elements: */
24 for(i=0;i<elements->Size();i++){
25 element=(Element*)elements->GetObjectByOffset(i);
26 element->AgressiveMigration();
27 }
28
29 /*Synchronize shelf status: */
30 for(i=0;i<elements->Size();i++){
31 element=(Element*)elements->GetObjectByOffset(i);
32 element->ShelfSync();
33 }
34
35}
36/*}}}*/
37/*FUNCTION SoftMigration {{{1*/
38void SoftMigration(Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads,Materials* materials, Parameters* parameters){
39
40 /*intermediary: */
41 int i;
42 double* potential_sheet_ungrounding=NULL;
43 double* sheet_ungrounding=NULL;
44 Element* element=NULL;
45
46 _printf_(VerboseModule()," Migrating grounding line\n");
47
48 /*First, figure out which nodes that are on the ice sheet want to unground. Build a flags vector for this (size number of nodes: */
49 potential_sheet_ungrounding=PotentialSheetUngrounding(elements,nodes,parameters);
50
51 /*Now, propagate ice shelf into connex areas of the ice sheet that potentially want to unground: */
52 sheet_ungrounding=PropagateShelfIntoConnexIceSheet(elements,nodes,parameters,potential_sheet_ungrounding);
53
54 /*Now, use the sheet_ungrounding flags to unground the ice sheet (at the same time, take care of grounding elements of the ice shelf
55 * that want to: */
56 for(i=0;i<elements->Size();i++){
57 element=(Element*)elements->GetObjectByOffset(i);
58 element->SoftMigration(sheet_ungrounding);
59 }
60
61 /*Synchronize shelf status: */
62 for(i=0;i<elements->Size();i++){
63 element=(Element*)elements->GetObjectByOffset(i);
64 element->ShelfSync();
65 }
66
67 /*free ressouces: */
68 xfree((void**)&potential_sheet_ungrounding);
69 xfree((void**)&sheet_ungrounding);
70
71}
72/*}}}*/
73/*FUNCTION PotentialSheetUngrounding {{{1*/
74double* PotentialSheetUngrounding(Elements* elements,Nodes* nodes,Parameters* parameters){
75
76 int i;
77 Element* element=NULL;
78 Vec vec_potential_sheet_ungrounding=NULL;
79 double* potential_sheet_ungrounding=NULL;
80 int numnods;
81 int configuration_type;
82
83 /*recover parameters: */
84 parameters->FindParam(&configuration_type,ConfigurationTypeEnum);
85
86 /*First, initialize vec_new_shelf_nodes, which will track which nodes have changed status: */
87 numnods=nodes->NumberOfNodes(configuration_type);
88 vec_potential_sheet_ungrounding=NewVec(numnods);
89
90 /*Loop through elements, and fill vec_potential_sheet_ungrounding: */
91 for(i=0;i<elements->Size();i++){
92 element=(Element*)elements->GetObjectByOffset(i);
93 element->PotentialSheetUngrounding(vec_potential_sheet_ungrounding);
94 }
95
96 /*Assemble vector: */
97 VecAssemblyBegin(vec_potential_sheet_ungrounding);
98 VecAssemblyEnd(vec_potential_sheet_ungrounding);
99
100 /*Serialize vector: */
101 VecToMPISerial(&potential_sheet_ungrounding,vec_potential_sheet_ungrounding);
102
103 /*free ressouces: */
104 VecFree(&vec_potential_sheet_ungrounding);
105
106 return potential_sheet_ungrounding;
107}
108/*}}}*/
109/*FUNCTION PropagateShelfIntoConnexIceSheet {{{1*/
110double* PropagateShelfIntoConnexIceSheet(Elements* elements,Nodes* nodes,Parameters* parameters,double* potential_sheet_ungrounding){
111
112 int i;
113 Element* element=NULL;
114 int numnods;
115 int analysis_type;
116 Vec vec_nodes_on_iceshelf=NULL;
117 double* nodes_on_iceshelf=NULL;
118 int nflipped,local_nflipped;
119 double* elements_touching_iceshelf=NULL;
120
121 /*recover parameters: */
122 parameters->FindParam(&analysis_type,AnalysisTypeEnum);
123
124 /*recover vec_nodes_on_iceshelf: */
125 vec_nodes_on_iceshelf=CreateNodesOnIceShelf(nodes,analysis_type);
126
127 nflipped=1; //bootstrap
128 while(nflipped){
129
130 /*get a list of potential elements that have nodes on ice shelf: */
131 elements_touching_iceshelf=CreateElementTouchingIceShelf(elements,vec_nodes_on_iceshelf);
132
133 /*now, go through elements_touching_iceshelf, and if they have nodes inside potential_sheet_ungrounding,
134 * flag it: */
135 local_nflipped=0;
136
137 /*serialize vec_nodes_on_iceshelf, needed by element's UpdatePotentialSheetUngrounding routine, to figure out if
138 * nodes have flipped from grounded to ungrounded: */
139 VecToMPISerial(&nodes_on_iceshelf,vec_nodes_on_iceshelf);
140
141 for(i=0;i<elements->Size();i++){
142 element=(Element*)elements->GetObjectByOffset(i);
143 if(elements_touching_iceshelf[element->Sid()]){
144 local_nflipped+=element->UpdatePotentialSheetUngrounding(potential_sheet_ungrounding,vec_nodes_on_iceshelf,nodes_on_iceshelf);
145 }
146 }
147
148 MPI_Allreduce(&local_nflipped,&nflipped,1,MPI_INT,MPI_SUM,MPI_COMM_WORLD);
149 _printf_(VerboseConvergence()," number of grounded vertices connected to grounding line: %i\n",nflipped);
150
151 /*Avoid leaks: */
152 xfree((void**)&nodes_on_iceshelf);
153 xfree((void**)&elements_touching_iceshelf);
154
155 /*Assemble:*/
156 VecAssemblyBegin(vec_nodes_on_iceshelf);
157 VecAssemblyEnd(vec_nodes_on_iceshelf);
158 }
159
160 /*Serialize: */
161 VecToMPISerial(&nodes_on_iceshelf,vec_nodes_on_iceshelf);
162
163 /*Free ressources:*/
164 VecFree(&vec_nodes_on_iceshelf);
165 xfree((void**)&elements_touching_iceshelf);
166
167 return nodes_on_iceshelf;
168}
169/*}}}*/
170/*FUNCTION CreateElementOnGroundingnine {{{1*/
171bool* CreateElementOnGroundingline(Elements* elements,double* element_on_iceshelf){
172
173 int i;
174 int j;
175 Element *element = NULL;
176 bool *element_on_gl = NULL;
177 bool ongl=false;
178 int *neighboorsids = NULL;
179 int sid;
180 int shelf;
181
182 /*Go through elements, and look for elements that can possibly have a grounding line migration. These
183 * are elements that touch the grounding line: */
184 element_on_gl=(bool*)xmalloc(elements->Size()*sizeof(bool));
185
186 for(i=0;i<elements->Size();i++){
187 element=(Element*)elements->GetObjectByOffset(i);
188
189 sid=element->Sid();
190 neighboorsids=element->GetHorizontalNeighboorSids();
191
192 shelf=(int)element_on_iceshelf[sid];
193 ongl=false;
194 for(j=0;j<3;j++){
195 if (neighboorsids[j]<0)continue; //this neighboor does not exist
196 if ((shelf==1) & (element_on_iceshelf[neighboorsids[j]]==1))continue; //both this element and this neighboor are on th ice shelf
197 if ((shelf==0) & (element_on_iceshelf[neighboorsids[j]]==0))continue; //both this element and this neighboor are on the ice sheet
198 ongl=true; //neighboor j is on a different location than us, ie we are touching the grounding line.
199 }
200 element_on_gl[i]=ongl;
201 }
202
203 return element_on_gl;
204}
205/*}}}*/
206/*FUNCTION CreateElementOnIceShelf {{{1*/
207double* CreateElementOnIceShelf(Elements* elements){
208
209 int i;
210 Element* element=NULL;
211 Vec vec_element_on_iceshelf=NULL;
212 double* element_on_iceshelf=NULL;
213
214 /*Create vector holding all the elements IsFloating flags: */
215 vec_element_on_iceshelf=NewVec(elements->NumberOfElements(),true);
216
217 /*Loop through elements, and fill vec_element_on_iceshelf: */
218 for(i=0;i<elements->Size();i++){
219 element=(Element*)elements->GetObjectByOffset(i);
220 VecSetValue(vec_element_on_iceshelf,element->Sid(),(int)element->IsFloating(),INSERT_VALUES);
221 }
222
223 /*Assemble vector: */
224 VecAssemblyBegin(vec_element_on_iceshelf);
225 VecAssemblyEnd(vec_element_on_iceshelf);
226
227 /*Serialize vector: */
228 VecToMPISerial(&element_on_iceshelf,vec_element_on_iceshelf);
229
230 /*free ressouces: */
231 VecFree(&vec_element_on_iceshelf);
232
233 return element_on_iceshelf;
234}
235/*}}}*/
236/*FUNCTION CreateElementTouchingIceShelf {{{1*/
237double* CreateElementTouchingIceShelf(Elements* elements,Vec vec_nodes_on_iceshelf){
238
239 int i;
240 Element* element=NULL;
241 Vec vec_element_touching_iceshelf=NULL;
242 double* nodes_on_iceshelf=NULL;
243
244 /*output: */
245 double* element_touching_iceshelf=NULL;
246
247 /*serialize: */
248 VecToMPISerial(&nodes_on_iceshelf,vec_nodes_on_iceshelf);
249
250 /*Create vector holding all the elements IsFloating flags: */
251 vec_element_touching_iceshelf=NewVec(elements->NumberOfElements(),true);
252
253 /*Loop through elements, and fill vec_element_touching_iceshelf: */
254 for(i=0;i<elements->Size();i++){
255 element=(Element*)elements->GetObjectByOffset(i);
256 VecSetValue(vec_element_touching_iceshelf,element->Sid(),element->IsNodeOnShelfFromFlags(nodes_on_iceshelf)?1.0:0.0,INSERT_VALUES);
257 }
258
259 /*Assemble vector: */
260 VecAssemblyBegin(vec_element_touching_iceshelf);
261 VecAssemblyEnd(vec_element_touching_iceshelf);
262
263 /*Serialize vector: */
264 VecToMPISerial(&element_touching_iceshelf,vec_element_touching_iceshelf);
265
266 /*free ressouces: */
267 VecFree(&vec_element_touching_iceshelf);
268 xfree((void**)&nodes_on_iceshelf);
269
270 return element_touching_iceshelf;
271}
272/*}}}*/
273/*FUNCTION UpdateShelfStatus {{{1*/
274int UpdateShelfStatus(Elements* elements,Nodes* nodes,Parameters* parameters,double* element_touching_iceshelf){
275
276 int i;
277 Element* element=NULL;
278 int configuration_type;
279 int numnods;
280 Vec vec_new_shelf_nodes=NULL;
281 double* new_shelf_nodes=NULL;
282
283 /*output: */
284 int local_nflipped=0;
285 int nflipped=0;
286
287 /*recover parameters: */
288 parameters->FindParam(&configuration_type,ConfigurationTypeEnum);
289
290 /*First, initialize vec_new_shelf_nodes, which will track which nodes have changed status: */
291 numnods=nodes->NumberOfNodes(configuration_type);
292 vec_new_shelf_nodes=NewVec(numnods);
293
294 /*Ok, now go through the elements that have gone through grounding line migration,
295 * and update their flags: */
296 for(i=0;i<elements->Size();i++){
297 element=(Element*)elements->GetObjectByOffset(i);
298 if(element_touching_iceshelf[element->Sid()]){
299 local_nflipped+=element->UpdateShelfStatus(vec_new_shelf_nodes);
300 }
301 }
302 MPI_Allreduce(&local_nflipped,&nflipped,1,MPI_INT,MPI_SUM,MPI_COMM_WORLD);
303
304 /*Serialize vec_new_shelf_nodes: */
305 VecToMPISerial(&new_shelf_nodes,vec_new_shelf_nodes);
306
307 /*Now, go through ALL elements, and update the status of the nodes, to propagate what happened at the grounding line: */
308 /*Carry out grounding line migration for those elements: */
309 for(i=0;i<elements->Size();i++){
310 element=(Element*)elements->GetObjectByOffset(i);
311 element->UpdateShelfFlags(new_shelf_nodes);
312 }
313
314 /*Free ressources: */
315 VecFree(&vec_new_shelf_nodes);
316 xfree((void**)&new_shelf_nodes);
317
318 return nflipped;
319}
320/*}}}*/
321/*FUNCTION CreateNodesOnIceShelf {{{1*/
322Vec CreateNodesOnIceShelf(Nodes* nodes,int configuration_type){
323
324 /*output: */
325 Vec nodes_on_iceshelf=NULL;
326
327 /*intermediary: */
328 int numnods;
329 int i;
330 Node* node=NULL;
331
332
333 /*First, initialize nodes_on_iceshelf, which will track which nodes have changed status: */
334 numnods=nodes->NumberOfNodes(configuration_type);
335 nodes_on_iceshelf=NewVec(numnods);
336
337 /*Loop through nodes, and fill nodes_on_iceshelf: */
338 for(i=0;i<nodes->Size();i++){
339 node=(Node*)nodes->GetObjectByOffset(i);
340 if(node->InAnalysis(configuration_type)){
341 if(node->IsFloating()){
342 VecSetValue(nodes_on_iceshelf,node->Sid(),1.0,INSERT_VALUES);
343 }
344 }
345 }
346
347 /*Assemble vector: */
348 VecAssemblyBegin(nodes_on_iceshelf);
349 VecAssemblyEnd(nodes_on_iceshelf);
350
351 return nodes_on_iceshelf;
352}
353/*}}}*/
Note: See TracBrowser for help on using the repository browser.