Changeset 20719

Timestamp:

06/09/16 11:23:46 (9 years ago)

Author:

Mathieu Morlighem

Message:

NEW: finished implementation of calving criterion based on minimum thickness

Location:

issm/trunk-jpl/src/c

Files:

• analyses/LevelsetAnalysis.cpp (modified) (3 diffs)
• classes/Elements/Element.cpp (modified) (1 diff)
• classes/Elements/Element.h (modified) (1 diff)
• classes/Elements/Tria.cpp (modified) (1 diff)
• classes/Elements/Tria.h (modified) (1 diff)
• classes/FemModel.cpp (modified) (1 diff)
• classes/FemModel.h (modified) (1 diff)
• modules/Calvingx/Calvingx.cpp (modified) (1 diff)

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

@@ -198 +198 @@
                 case DefaultCalvingEnum:
                 case CalvingDevEnum:
+                case CalvingMinthicknessEnum:
                         lsf_slopex_input  = basalelement->GetInput(LevelsetfunctionSlopeXEnum); _assert_(lsf_slopex_input);
                         if(dim==2) lsf_slopey_input  = basalelement->GetInput(LevelsetfunctionSlopeYEnum); _assert_(lsf_slopey_input);
@@ -251 +252 @@
                 switch(calvinglaw){
                         case DefaultCalvingEnum:
+                        case CalvingDevEnum:
+                        case CalvingMinthicknessEnum:
                                 lsf_slopex_input->GetInputValue(&dlsf[0],gauss);
                                 if(dim==2) lsf_slopey_input->GetInputValue(&dlsf[1],gauss);
@@ -271 +274 @@
 
                         case CalvingLevermannEnum:
-                        case CalvingDevEnum:
                                 calvingratex_input->GetInputValue(&c[0],gauss);
                                 if(dim==2) calvingratey_input->GetInputValue(&c[1],gauss);

issm/trunk-jpl/src/c/classes/Elements/Element.cpp

@@ -2077 +2077 @@
                                         this->CalvingRateDev();
                                         break;
+                                case CalvingMinthicknessEnum:
+                                        this->CalvingRateMinthickness();
+                                        break;
                                 default:
                                         _error_("Calving law "<<EnumToStringx(calvinglaw)<<" not supported yet");

issm/trunk-jpl/src/c/classes/Elements/Element.h

@@ -176 +176 @@
                 virtual void        CalvingRateLevermann(void)=0;
                 virtual void       CalvingRateDev(void){_error_("not implemented yet");};
+                virtual void       CalvingRateMinthickness(void){_error_("not implemented yet");};
                 virtual void       WriteLevelsetSegment(DataSet* segments){_error_("not implemented yet");};
                 virtual void       ResetLevelsetFromSegmentlist(IssmDouble* segments,int numsegments){_error_("not implemented yet");};

issm/trunk-jpl/src/c/classes/Elements/Tria.cpp

@@ -327 +327 @@
                 calvingratey[iv]=vy*sigma_vm/sigma_max;
                 calvingrate[iv]=sqrt(calvingratex[iv]*calvingratex[iv] + calvingratey[iv]*calvingratey[iv]);
+        }
+
+        /*Add input*/
+        this->inputs->AddInput(new TriaInput(CalvingratexEnum,&calvingratex[0],P1Enum));
+        this->inputs->AddInput(new TriaInput(CalvingrateyEnum,&calvingratey[0],P1Enum));
+        this->inputs->AddInput(new TriaInput(CalvingCalvingrateEnum,&calvingrate[0],P1Enum));
+
+        /*Clean up and return*/
+        delete gauss;
+}
+/*}}}*/
+void       Tria::CalvingRateMinthickness(){/*{{{*/
+
+        IssmDouble  calvingratex[NUMVERTICES];
+        IssmDouble  calvingratey[NUMVERTICES];
+        IssmDouble  calvingrate[NUMVERTICES];
+        IssmDouble  vx,vy,H,minH;
+
+        /*Retrieve all inputs and parameters we will need*/
+        Input* vx_input = inputs->GetInput(VxEnum); _assert_(vx_input);
+        Input* vy_input = inputs->GetInput(VyEnum); _assert_(vy_input);
+        Input* H_input  = inputs->GetInput(ThicknessEnum); _assert_(H_input);
+        parameters->FindParam(&minH,CalvingMinthicknessEnum);
+
+        /* Start looping on the number of vertices: */
+        GaussTria* gauss=new GaussTria();
+        for(int iv=0;iv<NUMVERTICES;iv++){
+                gauss->GaussVertex(iv);
+
+                /*Get velocity components and thickness*/
+                vx_input->GetInputValue(&vx,gauss);
+                vy_input->GetInputValue(&vy,gauss);
+                H_input->GetInputValue(&H,gauss);
+
+                /*Assign values*/
+                if(H>minH){
+                        calvingratex[iv]=0.;
+                        calvingratey[iv]=0.;
+                        calvingrate[iv]=0.;
+                }
+                else{
+                        calvingratex[iv]=vx+1e-2; //counter balance advance + add some retreat
+                        calvingratey[iv]=vy+1e-2;
+                        calvingrate[iv]=sqrt(calvingratex[iv]*calvingratex[iv] + calvingratey[iv]*calvingratey[iv]);
+                }
         }
 

issm/trunk-jpl/src/c/classes/Elements/Tria.h

@@ -54 +54 @@
                 void                    CalvingRateLevermann();
                 void                    CalvingRateDev();
+                void                    CalvingRateMinthickness();
                 void                    WriteLevelsetSegment(DataSet* segments);
                 void        ResetLevelsetFromSegmentlist(IssmDouble* segments,int numsegments);

issm/trunk-jpl/src/c/classes/FemModel.cpp

@@ -2166 +2166 @@
 }
 /*}}}*/
+void FemModel::CalvingRateMinthicknessx(){/*{{{*/
+
+        for(int i=0;i<elements->Size();i++){
+                Element* element=dynamic_cast<Element*>(this->elements->GetObjectByOffset(i));
+                element->CalvingRateMinthickness();
+        }
+}
+/*}}}*/
 void FemModel::StrainRateparallelx(){/*{{{*/
 

issm/trunk-jpl/src/c/classes/FemModel.h

@@ -95 +95 @@
                 void CalvingRateLevermannx();
                 void CalvingRateDevx();
+                void CalvingRateMinthicknessx();
                 void ResetLevelset();
                 #ifdef  _HAVE_DAKOTA_

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

@@ -27 +27 @@
                         femmodel->ElementOperationx(&Element::CalvingRateDev);
                         break;
+                case CalvingMinthicknessEnum:
+                        femmodel->CalvingRateMinthicknessx();
+                        break;
                 default:
                         _error_("Caving law "<<EnumToStringx(calvinglaw)<<" not supported yet");