Changeset 27789

Timestamp:

06/15/23 12:27:15 (21 months ago)

Author:

Cheng Gong

Message:

ADD: a new calving law for calvingMIP experiments

Location:

issm/trunk-jpl/src

Files:

• c/analyses/LevelsetAnalysis.cpp (modified) (4 diffs)
• c/classes/Elements/Element.cpp (modified) (1 diff)
• c/classes/Elements/Element.h (modified) (1 diff)
• c/classes/Elements/Tria.cpp (modified) (3 diffs)
• c/classes/Elements/Tria.h (modified) (1 diff)
• c/modules/Calvingx/Calvingx.cpp (modified) (1 diff)
• c/shared/Enum/EnumDefinitions.h (modified) (1 diff)
• c/shared/io/Marshalling/IoCodeConversions.cpp (modified) (1 diff)
• m/classes/calvingcalvingmip.m (added)
• m/contrib/chenggong/dataprocessing/expxy2shpll.m (moved) (moved from issm/trunk-jpl/src/m/contrib/chenggong/expxy2shpll.m )
• m/contrib/chenggong/modeldata/interpFromAtlasDEM.m (added)
• m/contrib/chenggong/modeldata/interpFromMEaSUREsGeotiff.m (moved) (moved from issm/trunk-jpl/src/m/contrib/chenggong/interpFromMEaSUREsGeotiff.m )
• m/contrib/chenggong/visualization/plotCompareTransientFlowline.m (modified) (2 diffs)
• m/exp/isoline.m (modified) (1 diff)
• m/modeldata/interpBedmachineGreenland.m (modified) (1 diff)

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

@@ -120 +120 @@
                         iomodel->FetchDataToInput(inputs,elements,"md.geometry.bed",BedEnum);
                         break;
+                case CalvingCalvingMIPEnum:
+                        break;
 
                 /*"Discrete" calving laws (need to specify rate as 0 so that we can still solve the level set equation)*/
@@ -150 +152 @@
                 case FrontalForcingsDefaultEnum:
                         iomodel->FetchDataToInput(inputs,elements,"md.frontalforcings.meltingrate",CalvingMeltingrateEnum);
-                        if (calvinglaw == CalvingParameterizationEnum) {
+                        if ((calvinglaw == CalvingParameterizationEnum) || (calvinglaw == CalvingCalvingMIPEnum)) {
                                 iomodel->FetchDataToInput(inputs,elements,"md.frontalforcings.ablationrate",CalvingAblationrateEnum);
                         }
@@ -256 +258 @@
                 case CalvingPollardEnum:
                         parameters->AddObject(iomodel->CopyConstantObject("md.calving.rc",CalvingRcEnum));
+                        break;
+                case CalvingCalvingMIPEnum:
+                        parameters->AddObject(iomodel->CopyConstantObject("md.calving.experiment",CalvingUseParamEnum));
+                        parameters->AddObject(iomodel->CopyConstantObject("md.calving.min_thickness",CalvingMinthicknessEnum));
                         break;
                 default:
@@ -748 +754 @@
 
         /*Apply minimum thickness criterion*/
-        if(calvinglaw==CalvingMinthicknessEnum || calvinglaw==CalvingVonmisesEnum || calvinglaw==CalvingParameterizationEnum || calvinglaw==CalvingVonmisesADEnum){
+        if(calvinglaw==CalvingMinthicknessEnum || calvinglaw==CalvingVonmisesEnum || calvinglaw==CalvingParameterizationEnum || calvinglaw==CalvingVonmisesADEnum || calvinglaw==CalvingCalvingMIPEnum){
 
                 IssmDouble mig_max = femmodel->parameters->FindParam(MigrationMaxEnum);

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

@@ -3793 +3793 @@
                                                                                                           case CalvingParameterizationEnum:
                                                                                                                   this->CalvingRateParameterization();
+                                                                                                                  break;
+                                                                                                          case CalvingCalvingMIPEnum:
+                                                                                                                  this->CalvingRateCalvingMIP();
                                                                                                                   break;
                                                                                                           case CalvingTestEnum:

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

@@ -238 +238 @@
                 virtual void             BasalNodeIndices(int* pnumindices,int** pindices,int finiteelement){_error_("not implemented yet");};
                 virtual void       CalvingRateParameterization(void){_error_("not implemented yet");};
+                virtual void       CalvingRateCalvingMIP(void){_error_("not implemented yet");};
                 virtual void       CalvingRateVonmises(void){_error_("not implemented yet");};
                 virtual void       CalvingRateVonmisesAD(void){_error_("not implemented yet");};

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

@@ -1115 +1115 @@
         this->AddInput(SigmaVMEnum,&sigma_vm[0],P1DGEnum);
         this->CalvingRateToVector();
+
+        /*Clean up and return*/
+        delete gauss;
+}
+/*}}}*/
+void       Tria::CalvingRateCalvingMIP(){/*{{{*/
+
+        IssmDouble  calvingrate[NUMVERTICES];
+        IssmDouble  calvingratex[NUMVERTICES];
+        IssmDouble  calvingratey[NUMVERTICES];
+        int                     experiment = 1;  /* exp:1 by default */
+        int         dim, domaintype;
+        IssmDouble      vx, vy, vel, c, wrate;
+
+        /*Get problem dimension and whether there is moving front or not*/
+        this->FindParam(&domaintype,DomainTypeEnum);
+
+        switch(domaintype){
+                case Domain2DverticalEnum:   dim = 1; break;
+                case Domain2DhorizontalEnum: dim = 2; break;
+                case Domain3DEnum:           dim = 2; break;
+                default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
+        }
+        if(dim==1) _error_("not implemented in 1D...");
+
+        /*Retrieve all inputs and parameters we will need*/
+        Input *vx_input      = this->GetInput(VxEnum);                                _assert_(vx_input);
+        Input *vy_input      = this->GetInput(VyEnum);                                _assert_(vy_input);
+        Input *wrate_input   = this->GetInput(CalvingAblationrateEnum);               _assert_(wrate_input); 
+
+        /* Use which experiment: use existing Enum */
+        this->FindParam(&experiment, CalvingUseParamEnum);
+
+        /* Start looping on the number of vertices: */
+        GaussTria* gauss=new GaussTria();
+        for(int iv=0;iv<NUMVERTICES;iv++){
+                gauss->GaussVertex(iv);
+
+                /*Get velocity components */
+                vx_input->GetInputValue(&vx,gauss);
+                vy_input->GetInputValue(&vy,gauss);
+                vel=sqrt(vx*vx+vy*vy)+1.e-14;
+
+
+                switch (experiment) { 
+                        case 1:
+                        case 3:
+                                /* Exp 1 and 3: set c=v-wrate, wrate=0, so that w=0 */
+                                wrate = 0.0;
+                                break;
+                        case 2:
+                        case 4:
+                                /* Exp 2 and 4: set c=v-wrate(given)*/
+                                /*Get wrate*/
+                                wrate_input->GetInputValue(&wrate,gauss);
+                                break;
+                        default:
+                                _error_("The experiment is not supported yet!");
+                }
+
+                calvingrate[iv] = vel - wrate;
+                calvingratex[iv] = vx - wrate*vx/vel;
+                calvingratey[iv] = vy - wrate*vy/vel;
+        }
+        /*Add input*/
+        this->AddInput(CalvingCalvingrateEnum,&calvingrate[0],P1DGEnum);
+        this->AddInput(CalvingratexEnum,&calvingratex[0],P1DGEnum);
+        this->AddInput(CalvingrateyEnum,&calvingratey[0],P1DGEnum);
 
         /*Clean up and return*/
@@ -4483 +4551 @@
                 case CalvingTestEnum:
                 case CalvingParameterizationEnum:
+                case CalvingCalvingMIPEnum:
                         calvingratex_input=this->GetInput(CalvingratexEnum); _assert_(calvingratex_input);
                         calvingratey_input=this->GetInput(CalvingrateyEnum); _assert_(calvingratey_input);
@@ -4528 +4597 @@
                         case CalvingLevermannEnum:
                         case CalvingPollardEnum:
+                        case CalvingCalvingMIPEnum:
                                 calvingratex_input->GetInputValue(&c[0],&gauss);
                                 calvingratey_input->GetInputValue(&c[1],&gauss);

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

@@ -63 +63 @@
                 void                    CalvingMeltingFluxLevelset();
                 void                    CalvingRateParameterization();
+                void                    CalvingRateCalvingMIP();
                 IssmDouble  CharacteristicLength(void);
                 void        ComputeBasalStress(void);

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

@@ -51 +51 @@
                         femmodel->ElementOperationx(&Element::CalvingPollard);
                         break;
+                case CalvingCalvingMIPEnum:
+                        femmodel->ElementOperationx(&Element::CalvingRateCalvingMIP);
+                        break;
                 default:
                         _error_("Caving law "<<EnumToStringx(calvinglaw)<<" not supported yet");

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

@@ -1391 +1391 @@
         CalvingTestEnum,
         CalvingParameterizationEnum,
+        CalvingCalvingMIPEnum,
         CalvingVonmisesEnum,
         CalvingVonmisesADEnum,

issm/trunk-jpl/src/c/shared/io/Marshalling/IoCodeConversions.cpp

@@ -293 +293 @@
                 case 10: return CalvingPollardEnum;
                 case 11: return CalvingVonmisesADEnum;
+                case 12:  return CalvingCalvingMIPEnum;
                 default: _error_("Marshalled Calving law code \""<<enum_in<<"\" not supported yet");
         }

issm/trunk-jpl/src/m/contrib/chenggong/visualization/plotCompareTransientFlowline.m

@@ -16 +16 @@
 %
 %   Author: Cheng Gong
-%   Date: 2023-04-10
+%   Date: 2021-12-06
 
 N = length(velList);
@@ -46 +46 @@
     h=colorbar;
     title(h,'m/a')
-         % add a vertical line to indicate initial ice front position
-         id = max(find(~isnan(cumsum(vel_flowline{p}(:,1)))));
-         icefront = flowline.Xmain(id);
-         plot([icefront, icefront], [2007,2020], '-r')
 end

issm/trunk-jpl/src/m/exp/isoline.m

@@ -258 +258 @@
 elseif strcmp(outputformat,'struct')
         %nothing to do, this is the default
+elseif strcmp(outputformat,'longest')
+        [~, mId] = max([contours.nods]);
+        contours = contours(mId);
 else
         disp('output format not supported, returning struct');

issm/trunk-jpl/src/m/modeldata/interpBedmachineGreenland.m

@@ -32 +32 @@
         ncdate='2021-08-27';
         ncdate='2022-03-17';
+        ncdate='2022-05-18';
         ncdate='2022-07-28';
+        ncdate='v6.0';
+        ncdate='v6.1';
 end
 basename = 'BedMachineGreenland';