Changeset 24036

Timestamp:

06/24/19 02:21:51 (6 years ago)

Author:

Mathieu Morlighem

Message:

CHG: first implementation of kill berg, seems to be workingmodules/KillIcebergsx/KillIcebergsx.cpp

Location:

issm/trunk-jpl/src/c

Files:

• cores/movingfront_core.cpp (modified) (1 diff)
• modules/KillIcebergsx/KillIcebergsx.cpp (modified) (6 diffs)

issm/trunk-jpl/src/c/cores/movingfront_core.cpp

@@ -88 +88 @@
         delete analysis;
 
+        /*Kill ice berg to avoid free body motion*/
+        if(VerboseSolution()) _printf0_("   looking for icebergs to kill\n");
+        KillIcebergsx(femmodel);
+
         /*Reset levelset if needed*/
         if(reinit_frequency && (step%reinit_frequency==0)){

issm/trunk-jpl/src/c/modules/KillIcebergsx/KillIcebergsx.cpp

@@ -7 +7 @@
 #include "../../toolkits/toolkits.h"
 
+#include "../InputUpdateFromVectorx/InputUpdateFromVectorx.h"
+
 void KillIcebergsx(FemModel* femmodel){
 
         /*Intermediaries*/
-        IssmDouble* local_mask = NULL;
-        const int MAXVERTICES = 6;
-        bool      found1;
-        int       sidlist[MAXVERTICES];
-        int       lidlist[MAXVERTICES];
+        int lid;
 
-        /*retrieve vertex info*/
-        int nbv_global = femmodel->vertices->NumberOfVertices();
-        int nbv_local  = femmodel->vertices->NumberOfVerticesLocal();
-        if(nbv_global==0)  return;
-        Vector<IssmDouble>* vec_connected_to_land=new Vector<IssmDouble>(nbv_local,nbv_global);
+        /*retrieve vertex info and prepare element flag to speed up process*/
+        int         nbv_local    = femmodel->vertices->Size();
+        IssmDouble *local_mask   = xNewZeroInit<IssmDouble>(nbv_local);
+        bool       *element_flag = xNewZeroInit<bool>(femmodel->elements->Size());
 
-        /*Prepare element flag to speed up process*/
-        bool* element_flag = xNewZeroInit<bool>(femmodel->elements->Size());
-
-        /*Fill vector with 1 once for all*/
-        IssmDouble eflags[MAXVERTICES];
-        for(int i=0;i<MAXVERTICES;i++) eflags[i] = 1.;
-
-        /*Step 1, go through all elements and put 1 in vec_connected_to_land if the element is grounded*/
+        /*Step 1, go through all elements and put 1 in local_mask if the element is grounded*/
         for(int i=0;i<femmodel->elements->Size();i++){
                 Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
-
                 if(!element->IsIceInElement()){
                         /*Nothing to do, just flag element to speed up the computation*/
@@ -39 +28 @@
                 else{
                         if(element->IsGrounded()){
-                                int  numvertices = element->GetNumberOfVertices();
-                                if(numvertices>MAXVERTICES) _error_("need to increase MAXVERTICES");
-                                element->GetVerticesSidList(&sidlist[0]);
-                                vec_connected_to_land->SetValues(numvertices,&sidlist[0],&eflags[0],ADD_VAL);
+                                int numvertices = element->GetNumberOfVertices();
+                                for(int v=0;v<numvertices;v++) local_mask[element->vertices[v]->Lid()] = 1.;
                         }
                 }
         }
-        vec_connected_to_land->Assemble();
 
-        /*Now we have 2 loops, one across cpus, and one for each cpus: we are going to propagate the mask if an element
-         * is connected to a positive mask already.
-         * We then communicate to the other partitions. We stop when the mask stops changing*/
+        /*Now we have 2 loops, one across cpus, and one for each cpus: we are going
+         * to propagate the mask if an element is connected to a positive mask
+         * already.  We then communicate to the other partitions. We stop when the
+         * mask stops changing*/
         bool keepsyncing = true;
         while(keepsyncing){
 
                 /*Get local mask from parallel vector*/
-                if(local_mask) xDelete<IssmDouble>(local_mask);
-                femmodel->GetLocalVectorWithClonesVertices(&local_mask,vec_connected_to_land);
+                femmodel->SyncLocalVectorWithClonesVertices(local_mask);
 
                 /*Local iterations on partition*/
-                bool keepgoing = true;
+                bool keepgoing    = true;
                 int  didsomething = 0;
-                int  iter      = 1;
+                int  iter         = 1;
                 while(keepgoing){
-                        _printf0_("   -- Kill icebergs: iteration "<<iter<<"\n");
+                        _printf0_("   -- Kill icebergs: local iteration "<<iter<<"\n");
 
                         keepgoing    = false;
@@ -72 +58 @@
                                 if(!element_flag[i]){
                                         int numvertices = element->GetNumberOfVertices();
-                                        element->GetVerticesLidList(&lidlist[0]);
                                         bool found1 = false;
                                         for(int j=0;j<numvertices;j++){
-                                                if(local_mask[lidlist[j]]>0.){
+                                                lid = element->vertices[j]->Lid();
+                                                if(local_mask[lid]>0.){
                                                         found1 = true;
                                                         break;
@@ -83 +69 @@
                                                 element_flag[i] = true;
                                                 for(int j=0;j<numvertices;j++){
-                                                        if(local_mask[lidlist[j]]==0.){
-                                                                local_mask[lidlist[j]]=1.;
+                                                        lid = element->vertices[j]->Lid();
+                                                        if(local_mask[lid]==0.){
+                                                                local_mask[lid]=1.;
                                                                 keepgoing = true;
                                                                 didsomething = 1;
@@ -96 +83 @@
 
                 /*Check how many iterations all cpus did*/
-                int didsomething_max;
-                ISSM_MPI_Reduce(&didsomething,&didsomething_max,1,ISSM_MPI_INT,ISSM_MPI_MAX,0,IssmComm::GetComm());
-                ISSM_MPI_Bcast(&didsomething_max,1,ISSM_MPI_INT,0,IssmComm::GetComm());
-                if(didsomething_max==0){
+                int iter_max;
+                ISSM_MPI_Reduce(&iter,&iter_max,1,ISSM_MPI_INT,ISSM_MPI_MAX,0,IssmComm::GetComm());
+                ISSM_MPI_Bcast(&iter_max,1,ISSM_MPI_INT,0,IssmComm::GetComm());
+                if(iter_max==2){
+                        /*If iter is only 2, nothing else was changed in the while loop above (iter is initialized as 1 and then ++)*/
                         keepsyncing = false;
-                }
-                else{
-                        for(int i=0;i<femmodel->elements->Size();i++){
-                                Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
-
-                                if(element->IsIceInElement()){
-                                                int  numvertices = element->GetNumberOfVertices();
-                                                element->GetVerticesSidList(&sidlist[0]);
-                                                element->GetVerticesLidList(&lidlist[0]);
-                                                for(int j=0;j<numvertices;j++) eflags[j] = local_mask[lidlist[j]];
-                                                vec_connected_to_land->SetValues(numvertices,&sidlist[0],&eflags[0],ADD_VAL);
-                                        }
-                        }
-                        vec_connected_to_land->Assemble();
                 }
         }
 
+        InputUpdateFromVectorx(femmodel,local_mask,PressureEnum,VertexLIdEnum);
         /*Cleanup*/
         xDelete<bool>(element_flag);
-        delete vec_connected_to_land;
 
         /*OK, now deactivate iceberg and count the number of deactivated vertices*/
@@ -128 +102 @@
                 if(element->IsIceInElement()){
                         int  numvertices = element->GetNumberOfVertices();
-                        element->GetVerticesLidList(&lidlist[0]);
                         bool deactivate = false;
                         for(int j=0;j<numvertices;j++){
-                                if(local_mask[lidlist[j]]==0.){
+                                lid = element->vertices[j]->Lid();
+                                if(local_mask[lid]==0.){
                                         deactivate = true;
                                         break;