Changeset 26450

Timestamp:

09/21/21 12:06:54 (3 years ago)

Author:

Mathieu Morlighem

Message:

CHG: rewriting a bit

File:

• issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp (modified) (11 diffs)

issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp

@@ -125 +125 @@
                 case CalvingCrevasseDepthEnum:
                         parameters->AddObject(iomodel->CopyConstantObject("md.calving.crevasse_opening_stress",CalvingCrevasseDepthEnum));
+                        parameters->AddObject(iomodel->CopyConstantObject("md.calving.crevasse_threshold",CalvingCrevasseThresholdEnum));
                         break;
                 case CalvingDev2Enum:
@@ -472 +473 @@
 
         /*Intermediaries*/
-        int         calvinglaw;
-        IssmDouble  min_thickness,thickness,hab_fraction;
-        IssmDouble      crevassedepth,surface_crevasse,surface,critical_fraction;
-        IssmDouble  rho_ice,rho_water;
-        IssmDouble  bed,water_depth;
-        IssmDouble  levelset,sealevel;
-
+        int  calvinglaw;
         femmodel->parameters->FindParam(&calvinglaw,CalvingLawEnum);
 
         if(calvinglaw==CalvingMinthicknessEnum || calvinglaw==CalvingVonmisesEnum){
 
-                /*Get minimum thickness threshold*/
-                femmodel->parameters->FindParam(&min_thickness,CalvingMinthicknessEnum);
+                /*Intermediaries*/
+                IssmDouble thickness,bed,sealevel;
+                IssmDouble min_thickness = femmodel->parameters->FindParam(CalvingMinthicknessEnum);
 
                 /*Loop over all elements of this partition*/
@@ -518 +514 @@
    else if(calvinglaw==CalvingHabEnum){
 
+                /*Intermediaries*/
+                IssmDouble  thickness,water_depth,levelset,hab_fraction;
+
                 /*Get the fraction of the flotation thickness at the terminus*/
                 InputDuplicatex(femmodel,MaskIceLevelsetEnum,DistanceToCalvingfrontEnum);
@@ -526 +525 @@
                         Element* element  = xDynamicCast<Element*>(object);
 
-                        rho_ice = element->FindParam(MaterialsRhoIceEnum);
-                        rho_water = element->FindParam(MaterialsRhoSeawaterEnum);
+                        IssmDouble rho_ice   = element->FindParam(MaterialsRhoIceEnum);
+                        IssmDouble rho_water = element->FindParam(MaterialsRhoSeawaterEnum);
 
                         int      numnodes           = element->GetNumberOfNodes();
@@ -560 +559 @@
    else if(calvinglaw==CalvingCrevasseDepthEnum){
 
-                int   nflipped,local_nflipped;
+                /*Intermediaries*/
+                IssmDouble  levelset,crevassedepth,bed,surface_crevasse,thickness,surface;
+                int         nflipped,local_nflipped;
                 IssmDouble* constraint_nodes = NULL;
 
@@ -571 +572 @@
       int localmasters  = femmodel->nodes->NumberOfNodesLocal();
       Vector<IssmDouble>* vec_constraint_nodes = vec_constraint_nodes=new Vector<IssmDouble>(localmasters,numnodes);
+
+      IssmDouble crevasse_threshold = femmodel->parameters->FindParam(CalvingCrevasseThresholdEnum);
 
                 for(Object* & object : femmodel->elements->objects){
@@ -596 +599 @@
                                         surface_input->GetInputValue(&surface,gauss);
 
-                                        if((surface_crevasse-surface>0. || crevassedepth-thickness>0.) && bed<0.){
+                                        if((surface_crevasse>surface || crevassedepth>crevasse_threshold*thickness) && bed<0.){
                                                 vec_constraint_nodes->SetValue(node->Pid(),1.0,INS_VAL);
                                         }
@@ -616 +619 @@
                                 Gauss*   gauss    = element->NewGauss();
 
-                                Input*   levelset_input         = element->GetInput(DistanceToCalvingfrontEnum); _assert_(levelset_input);
-                                Input*   crevassedepth_input    = element->GetInput(CrevasseDepthEnum); _assert_(crevassedepth_input);
-                                Input*   bed_input              = element->GetInput(BedEnum); _assert_(bed_input);
-                                Input*   surface_crevasse_input = element->GetInput(SurfaceCrevasseEnum); _assert_(surface_crevasse_input);
-                                Input*   thickness_input        = element->GetInput(ThicknessEnum); _assert_(thickness_input);
-                                Input*   surface_input          = element->GetInput(SurfaceEnum); _assert_(surface_input);
+                                Input *levelset_input         = element->GetInput(DistanceToCalvingfrontEnum); _assert_(levelset_input);
+                                Input *crevassedepth_input    = element->GetInput(CrevasseDepthEnum);          _assert_(crevassedepth_input);
+                                Input *bed_input              = element->GetInput(BedEnum);                    _assert_(bed_input);
+                                Input *surface_crevasse_input = element->GetInput(SurfaceCrevasseEnum);        _assert_(surface_crevasse_input);
+                                Input *thickness_input        = element->GetInput(ThicknessEnum);              _assert_(thickness_input);
+                                Input *surface_input          = element->GetInput(SurfaceEnum);                _assert_(surface_input);
 
                                 /*Is this element connected to a node that should be calved*/
@@ -645 +648 @@
                                                 surface_input->GetInputValue(&surface,gauss);
 
-                                                if((surface_crevasse-surface>0. || crevassedepth-thickness>0.) && bed<0. && levelset>-300. && levelset<0. && constraint_nodes[node->Lid()]==0.){
+                  if((surface_crevasse>surface || crevassedepth>crevasse_threshold*thickness) && bed<0. && levelset>-300. && levelset<0. && constraint_nodes[node->Lid()]==0.){
                                                         local_nflipped++;
                                                         vec_constraint_nodes->SetValue(node->Pid(),1.0,INS_VAL);
@@ -655 +658 @@
                         /*Count how many new nodes were found*/
                         ISSM_MPI_Allreduce(&local_nflipped,&nflipped,1,ISSM_MPI_INT,ISSM_MPI_SUM,IssmComm::GetComm());
-                        //_printf0_("Found "<<nflipped<<" to flip\n");
+                        _printf0_("Found "<<nflipped<<" to flip\n");
 
                         /*Assemble and serialize flag vector*/
@@ -662 +665 @@
          femmodel->GetLocalVectorWithClonesNodes(&constraint_nodes,vec_constraint_nodes);
                 }
+
                 /*Free ressources:*/
                 delete vec_constraint_nodes;