Changeset 26973

Timestamp:

04/27/22 11:22:23 (3 years ago)

Author:

Cheng Gong

Message:

CHG: add a special test case for killberg with only one single node of ice and has ocean_levelset crossing the neighbour elements, select the initial mask more strict to fix this bug

Location:

issm/trunk-jpl

Files:

• src/c/cores/movingfront_core.cpp (modified) (2 diffs)
• src/c/modules/KillIcebergsx/KillIcebergsx.cpp (modified) (9 diffs)
• src/c/modules/KillIcebergsx/KillIcebergsx.h (modified) (1 diff)
• src/c/shared/Enum/EnumDefinitions.h (modified) (1 diff)
• src/m/parameterization/killberg.m (modified) (1 diff)
• test/Archives/Archive813.arch (added)
• test/NightlyRun/test813.m (added)

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

@@ -16 +16 @@
 
         /* intermediaries */
-        bool save_results,isstressbalance,ismasstransport,isthermal,isenthalpy,islevelset,ismovingfront,killicebergs,haskillberg;
+        bool save_results,isstressbalance,ismasstransport,isthermal,isenthalpy,islevelset,ismovingfront,killicebergs;
         int  domaintype, num_extrapol_vars, index,reinit_frequency,step;
         int* extrapol_vars=NULL;
@@ -120 +120 @@
         /*Kill ice berg to avoid free body motion*/
         if(killicebergs){
+                int killberg = 0;
                 if(VerboseSolution()) _printf0_("   looking for icebergs to kill\n");
-                KillIcebergsx(femmodel);
-                femmodel->parameters->FindParam(&haskillberg,LevelsetHasKillIcebergsEnum);
-                if (haskillberg) {
-                        if(VerboseSolution()) _printf0_("   reinitializing level set after killing icebergs\n");
+                killberg = KillIcebergsx(femmodel);
+                /*wait for all cores*/
+                int totalkill;
+                ISSM_MPI_Reduce(&killberg,&totalkill,1,ISSM_MPI_INT,ISSM_MPI_SUM,0,IssmComm::GetComm());
+                ISSM_MPI_Bcast(&totalkill,1,ISSM_MPI_INT,0,IssmComm::GetComm());
+
+                if (totalkill > 0) {
+                        if(VerboseSolution()) _printf0_("   reinitializing level set after killing " << totalkill << " icebergs\n");
                         femmodel->ResetLevelset();
                         ResetBoundaryConditions(femmodel,LevelsetAnalysisEnum);

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

@@ -9 +9 @@
 #include "../InputUpdateFromVectorx/InputUpdateFromVectorx.h"
 
-void KillIcebergsx(FemModel* femmodel){
+int KillIcebergsx(FemModel* femmodel){
 
         /*Intermediaries*/
@@ -18 +18 @@
         IssmDouble *local_mask   = xNewZeroInit<IssmDouble>(nbv_local);
         bool       *element_flag = xNewZeroInit<bool>(femmodel->elements->Size());
+        IssmDouble ice;
 
         /*Step 1, go through all elements and put 1 in local_mask if the element is grounded*/
@@ -27 +28 @@
                         element_flag[i] = true;
                 }
-                else{
+                else {
                         if(element->IsAllGrounded()){
+                                /* only look at element with ice but not fully grounded */
                                 int numvertices = element->GetNumberOfVertices();
-                                for(int v=0;v<numvertices;v++) local_mask[element->vertices[v]->Lid()] = 1.;
+                                Gauss* gauss=element->NewGauss();
+                                Input* icelevel_input = element->GetInput(MaskIceLevelsetEnum);                         _assert_(icelevel_input);
+
+                                for(int v=0;v<numvertices;v++) {
+                                        gauss->GaussVertex(v);
+                                        icelevel_input->GetInputValue(&ice,gauss);
+                                        /* The initial mask is very strict, we look at all grounded elements and set the mask for ice nodes only. */
+                                        if (ice < 0) local_mask[element->vertices[v]->Lid()] = 1.;
+                                }
+                                delete gauss;
                         }
                 }
@@ -48 +59 @@
                 /*Local iterations on partition*/
                 bool keepgoing    = true;
-                int  didsomething = 0;
                 int  iter         = 1;
                 while(keepgoing){
@@ -54 +64 @@
 
                         keepgoing    = false;
-                        didsomething = 0;
                         int i=0;
                         for(Object* & object : femmodel->elements->objects){
@@ -62 +71 @@
                                         int numvertices = element->GetNumberOfVertices();
                                         bool found1 = false;
+                                        IssmDouble sumlocalmask = 0.;
+
                                         for(int j=0;j<numvertices;j++){
                                                 lid = element->vertices[j]->Lid();
-                                                if(local_mask[lid]>0.){
+                                                /*we need to have at least a sharing edge, to extend the mask*/
+                                                sumlocalmask += local_mask[lid];
+                                                if(sumlocalmask > 1.5){
                                                         found1 = true;
                                                         break;
@@ -73 +86 @@
                                                 for(int j=0;j<numvertices;j++){
                                                         lid = element->vertices[j]->Lid();
-                                                        if(local_mask[lid]==0.){
-                                                                local_mask[lid]=1.;
-                                                                keepgoing = true;
-                                                                didsomething = 1;
-                                                        }
+                                                        local_mask[lid]=1.;
                                                 }
+                                                keepgoing = true;
                                         }
                                 }
@@ -99 +109 @@
         xDelete<bool>(element_flag);
 
-        bool killbergReinit = false;
+        int killbergReinit = 0;
         /*OK, now deactivate iceberg and count the number of deactivated vertices*/
         for(Object* & object : femmodel->elements->objects){
@@ -121 +131 @@
                                 element->AddInput(MaskIceLevelsetEnum,values,P1Enum);
                                 xDelete<IssmDouble>(values);
-                                killbergReinit = true;
+                                killbergReinit += 1;
                         }
                 }
         }
-        /*Recompute the sign distance for the levelset function*/
-   femmodel->parameters->SetParam(killbergReinit,LevelsetHasKillIcebergsEnum);
-
         /*cleanup*/
         xDelete<IssmDouble>(local_mask);
+
+        /*Recompute the sign distance for the levelset function*/
+        return killbergReinit;
 }

issm/trunk-jpl/src/c/modules/KillIcebergsx/KillIcebergsx.h ¶

@@ -6 +6 @@
 
 /* local prototypes: */
-void KillIcebergsx(FemModel* femmodel);
+int KillIcebergsx(FemModel* femmodel);
 
 #endif

issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h ¶

@@ -270 +270 @@
         LevelsetReinitFrequencyEnum,
         LevelsetStabilizationEnum,
-        LevelsetHasKillIcebergsEnum,
         LockFileNameEnum,
         LoveAllowLayerDeletionEnum,

issm/trunk-jpl/src/m/parameterization/killberg.m ¶

@@ -28 +28 @@
 %isgrounded = max(md.mask.ocean_levelset(md.mesh.elements),[],2)>0;
 pos = find(isgrounded);
-element_flag(pos) = 1;
+%element_flag(pos) = 1;
 mask(md.mesh.elements(pos,:)) = 1;
+mask(md.mask.ice_levelset>=0) = 0;
 
 iter = 1;