Changeset 7288

Timestamp:

02/03/11 06:26:36 (14 years ago)

Author:

Eric.Larour

Message:

Grounding line migration: trying to debug.

Location:

issm/trunk/src/c

Files:

• Makefile.am (modified) (2 diffs)
• modules/GroundingLineMigrationx/GroundingLineMigrationx.cpp (modified) (2 diffs)
• modules/GroundingLineMigrationx/GroundingLineMigrationxLocal.h (modified) (2 diffs)
• modules/GroundingLineMigrationx/UpdateShelfStatus.cpp (modified) (3 diffs)
• objects/Elements/Element.h (modified) (1 diff)
• objects/Elements/Penta.cpp (modified) (2 diffs)
• objects/Elements/Penta.h (modified) (1 diff)
• objects/Elements/Tria.cpp (modified) (6 diffs)
• objects/Elements/Tria.h (modified) (1 diff)
• solutions/control_core.cpp (modified) (1 diff)
• solutions/transient2d_core.cpp (modified) (1 diff)
• solutions/transient3d_core.cpp (modified) (1 diff)

issm/trunk/src/c/Makefile.am

@@ -359 +359 @@
                                         ./modules/GroundingLineMigrationx/CreateElementTouchingIceShelf.cpp\
                                         ./modules/GroundingLineMigrationx/CreateElementOnGroundingLine.cpp\
+                                        ./modules/GroundingLineMigrationx/UpdateShelfStatus.cpp\
                                         ./modules/ModelProcessorx/ModelProcessorx.h\
                                         ./modules/ModelProcessorx/ModelProcessorx.cpp\
@@ -937 +938 @@
                                         ./modules/GroundingLineMigrationx/CreateElementTouchingIceShelf.cpp\
                                         ./modules/GroundingLineMigrationx/CreateElementOnGroundingLine.cpp\
+                                        ./modules/GroundingLineMigrationx/UpdateShelfStatus.cpp\
                                         ./modules/ModelProcessorx/ModelProcessorx.h\
                                         ./modules/ModelProcessorx/ModelProcessorx.cpp\

issm/trunk/src/c/modules/GroundingLineMigrationx/GroundingLineMigrationx.cpp

@@ -14 +14 @@
         int i,j;
         Element* element=NULL;
-        double* element_touching_iceshelf=NULL;
+        double*  element_touching_iceshelf=NULL;
+        int      nflipped;
 
         _printf_(VerboseModule(),"   Migrating grounding line\n");
@@ -20 +21 @@
         /*Loop until no more nodes and elements  change from grounded to floating: */
         nflipped=1; //arbitrary, get things started
+
         while(nflipped){
 

issm/trunk/src/c/modules/GroundingLineMigrationx/GroundingLineMigrationxLocal.h

@@ -7 +7 @@
 
 class Elements;
+class Nodes;
+class Parameters;
 
 /* local prototypes: */
@@ -12 +14 @@
 bool*   CreateElementOnGroundingLine(Elements* elements,double* element_on_iceshelf);
 double* CreateElementTouchingIceShelf(Elements* elements);
-int UpdateShelfStatus(Elements* elements,Nodes* nodes,Parameters* parameters);
 int UpdateShelfStatus(Elements* elements,Nodes* nodes,Parameters* parameters,double* element_touching_iceshelf);
 

issm/trunk/src/c/modules/GroundingLineMigrationx/UpdateShelfStatus.cpp

@@ -1 +1 @@
 /*!\file UpdateShelfStatus
- * \brief: update elements and nodes shelf status. Plus return how many nodes have changed status. */
+ * \brief: update elements and nodes shelf status. Plus return how many nodes have changed status.
  */
 
@@ -16 +16 @@
         int numnods;
         Vec vec_new_shelf_nodes=NULL;
+        double* new_shelf_nodes=NULL;
 
         /*output: */
+        int local_nflipped=0;
         int nflipped=0;
 
@@ -31 +33 @@
         for(i=0;i<elements->Size();i++){
                 element=(Element*)elements->GetObjectByOffset(i);
-                nflipped+=element->UpdateShelfStatus(vec_new_shelf_nodes);
+                if(element_touching_iceshelf[element->Sid()]){
+                        local_nflipped+=element->UpdateShelfStatus(vec_new_shelf_nodes);
+                }
         }
+        MPI_Allreduce(&local_nflipped,&nflipped,1,MPI_INT,MPI_SUM,MPI_COMM_WORLD);
 
         /*Serialize vec_new_shelf_nodes: */

issm/trunk/src/c/objects/Elements/Element.h

@@ -84 +84 @@
                 virtual double TimeAdapt()=0;
                 virtual void   MigrateGroundingLine()=0;
-                virtual void   UpdateShelfStatus()=0;
+                virtual int    UpdateShelfStatus(Vec new_shelf_nodes)=0;
+                virtual void   UpdateShelfFlags(double* new_shelf_nodes)=0;
 
                 /*Implementation: */

issm/trunk/src/c/objects/Elements/Penta.cpp

@@ -5040 +5040 @@
                 /*Check solution*/
                 if(isnan(values[i])) _error_("NaN found in solution vector");
-                if(values[i]<0)      printf("temperature < 0°K found in solution vector\n");
-                if(values[i]>275)    printf("temperature > 275°K found in solution vector (Paterson's rheology associated is negative)\n");
+                //if(values[i]<0)      printf("temperature < 0°K found in solution vector\n");
+                //if(values[i]>275)    printf("temperature > 275°K found in solution vector (Paterson's rheology associated is negative)\n");
         }
 
@@ -6131 +6131 @@
 /*}}}*/
 /*FUNCTION Penta::UpdateShelfStatus{{{1*/
-void Penta::UpdateShelfStatus(void){
+int Penta::UpdateShelfStatus(Vec new_shelf_nodes){
         _error_("Not implemented yet");
 }
 /*}}}*/
+/*FUNCTION Penta::UpdateShelfFlags{{{1*/
+void Penta::UpdateShelfFlags(double* new_shelf_nodes){
+        _error_("Not implemented yet");
+}
+/*}}}*/

issm/trunk/src/c/objects/Elements/Penta.h

@@ -125 +125 @@
                 double SurfaceArea(void);
                 void   Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type);
-                void   UpdateShelfStatus(void);
+                int    UpdateShelfStatus(Vec new_shelf_nodes);
+                void   UpdateShelfFlags(double* new_shelf_nodes);
                 double TimeAdapt();
                 int*   GetHorizontalNeighboorSids(void);

issm/trunk/src/c/objects/Elements/Tria.cpp

@@ -53 +53 @@
                 /*Build horizontalneighborsids list: */
                 _assert_(iomodel->elementconnectivity);
-                for (i=0;i<3;i++) this->horizontalneighborsids[i]=(int)iomodel->elementconnectivity[3*index+i]-1;
+                //for (i=0;i<3;i++) this->horizontalneighborsids[i]=(int)iomodel->elementconnectivity[3*index+i]-1;
 
                 /*intialize inputs and results: */
@@ -4437 +4437 @@
                                 b[i]=ba[i];
                                 s[i]=b[i]+h[i];
+                                //printf("Node id %i is starting to ground\n",nodes[i]->Id());
                         }
                         else{
@@ -4450 +4451 @@
                                 s[i]=(1-density)*h[i];
                                 b[i]=-density*h[i];
+                                //printf("Node id %i is starting to float\n",nodes[i]->Id());
                         }
                         else{
@@ -5327 +5329 @@
 /*}}}*/
 /*FUNCTION Tria::UpdateShelfStatus{{{1*/
-void  Tria::UpdateShelfStatus(void){
+int  Tria::UpdateShelfStatus(Vec new_shelf_nodes){
 
         int     i;
@@ -5341 +5343 @@
         double  density;
         bool    elementonshelf = false;
+        int     flipped=0;
+        bool    shelfstatus[3];
 
 
@@ -5358 +5362 @@
         density=rho_ice/rho_water;
 
+        /*Initialize current status of nodes: */
+        for(i=0;i<3;i++){
+                shelfstatus[i]=nodes[i]->inputs->GetInput(NodeOnIceShelfEnum);
+        }
+        
         /*go through vertices, and figure out if they are grounded or not, then update their status: */
+        flipped=0;
         for(i=0;i<3;i++){
                 if(b[i]<=ba[i]){
                         nodes[i]->inputs->AddInput(new BoolInput(NodeOnIceShelfEnum,false));
                         nodes[i]->inputs->AddInput(new BoolInput(NodeOnIceSheetEnum,true));
+                        if(shelfstatus[i])flipped++;
+                        //printf("node %i flipping to sheet\n",nodes[i]->Sid());
+                        VecSetValue(new_shelf_nodes,nodes[i]->Sid(),0,INSERT_VALUES);
                 }
                 else{
                         nodes[i]->inputs->AddInput(new BoolInput(NodeOnIceShelfEnum,true));
                         nodes[i]->inputs->AddInput(new BoolInput(NodeOnIceSheetEnum,false));
+                        if(!shelfstatus[i])flipped++;
+                        //printf("node %i flipping to shelf\n",nodes[i]->Sid());
+                        VecSetValue(new_shelf_nodes,nodes[i]->Sid(),1,INSERT_VALUES);
                 }
         }
 
         /*Now, update  shelf status of element. An element can only be on shelf if all its nodes are on shelf: */
-        elementonshelf=true;
+        elementonshelf=false;
         for(i=0;i<3;i++){
-                if(nodes[i]->IsOnSheet()){
-                        elementonshelf=false;
+                if(nodes[i]->IsOnShelf()){
+                        elementonshelf=true;
                         break;
                 }
         }
     this->inputs->AddInput(new BoolInput(ElementOnIceShelfEnum,elementonshelf));
-}
-/*}}}*/
+
+         return flipped;
+}
+/*}}}*/
+/*FUNCTION Tria::UpdateShelfFlags{{{1*/
+void Tria::UpdateShelfFlags(double* new_shelf_nodes){
+
+        /*go through vertices, and update the status of NodeOnIceShelfEnum and NodeOnIceSheetEnum flags: */
+        bool flag;
+        int  i;
+
+        for(i=0;i<3;i++){
+                flag=(bool)new_shelf_nodes[nodes[i]->Sid()];
+                nodes[i]->inputs->AddInput(new BoolInput(NodeOnIceShelfEnum,flag));
+                nodes[i]->inputs->AddInput(new BoolInput(NodeOnIceSheetEnum,!flag));
+        }
+}
+/*}}}*/

issm/trunk/src/c/objects/Elements/Tria.h

@@ -129 +129 @@
                 double SurfaceArea(void);
                 void   Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type);
-                void   UpdateShelfStatus(void);
+                int    UpdateShelfStatus(Vec new_shelf_nodes);
+                void   UpdateShelfFlags(double* new_shelf_nodes);
                 double TimeAdapt();
                 int*   GetHorizontalNeighboorSids(void);

issm/trunk/src/c/solutions/control_core.cpp

@@ -108 +108 @@
 
                 /*Temporary saving every 5 control steps: */
-                if (((n+1)%5)==0){
+                /*if (((n+1)%5)==0){
                         _printf_(VerboseControl(),"%s\n","   saving temporary results");
                         controlrestart(femmodel,J);
-                }
+                }*/
         }
 

issm/trunk/src/c/solutions/transient2d_core.cpp

@@ -53 +53 @@
                 step+=1;
 
-                _printf_(VerboseSolution(),"%s%g%s%i%s%g%s%g\n","time [yr]: ",time/yts,"    iteration number: ",step,"/",floor(finaltime/dt)," dt [yr]: ",dt/yts);
+                _printf_(VerboseSolution(),"%s%g%s%i%s%g%s%g\n","time [yr]: ",time/yts,"    iteration number: ",step,"/",floor((finaltime-time)/dt)," dt [yr]: ",dt/yts);
 
                 _printf_(VerboseSolution(),"%s\n","   computing new velocity");

issm/trunk/src/c/solutions/transient3d_core.cpp

@@ -51 +51 @@
                 time+=dt;
 
-                _printf_(VerboseSolution(),"%s%g%s%i%s%g%s%g\n","time [yr]: ",time/yts,"    iteration number: ",step,"/",floor(finaltime/dt)," dt [yr]: ",dt/yts);
+                _printf_(VerboseSolution(),"%s%g%s%i%s%g%s%g\n","time [yr]: ",time/yts,"    iteration number: ",step,"/",floor((finaltime-time)/dt)," dt [yr]: ",dt/yts);
 
                 _printf_(VerboseSolution(),"   computing temperatures:\n");