Changeset 28039

Timestamp:

01/04/24 05:49:30 (15 months ago)

Author:

seroussi

Message:

NEW: extended CalvingMIP to 3d Penta

Location:

issm/trunk-jpl/src/c/classes/Elements

Files:

• Penta.cpp (modified) (4 diffs)
• Penta.h (modified) (1 diff)

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

@@ -286 +286 @@
         IssmDouble  groundedice,bed,sealevel;
 
-        /*Depth average B for stress calculation*/
+        /*Depth average velocity for stress calculation*/
         this->InputDepthAverageAtBase(VxEnum,VxAverageEnum);
         this->InputDepthAverageAtBase(VyEnum,VyAverageEnum);
@@ -625 +625 @@
                 delete gauss;
         }
+}
+/*}}}*/
+void       Penta::CalvingRateCalvingMIP(){/*{{{*/
+
+        if(!this->IsOnBase()) return;
+
+        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;
+        IssmDouble  time, groundedice, yts;
+
+        /*Get problem dimension and whether there is moving front or not*/
+        this->FindParam(&domaintype,DomainTypeEnum);
+        this->FindParam(&time,TimeEnum);
+        this->FindParam(&yts,ConstantsYtsEnum);
+
+        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...");
+
+        /*Depth average velocity for stress calculation*/
+        this->InputDepthAverageAtBase(VxEnum,VxAverageEnum);
+        this->InputDepthAverageAtBase(VyEnum,VyAverageEnum);
+
+        /*Retrieve all inputs and parameters we will need*/
+        Input *vx_input      = this->GetInput(VxAverageEnum);                                _assert_(vx_input);
+        Input *vy_input      = this->GetInput(VyAverageEnum);                                _assert_(vy_input);
+        Input *wrate_input   = this->GetInput(CalvingAblationrateEnum);               _assert_(wrate_input); 
+        Input* gr_input      = this->GetInput(MaskOceanLevelsetEnum);                                           _assert_(gr_input);
+
+        /* Use which experiment: use existing Enum */
+        this->FindParam(&experiment, CalvingUseParamEnum);
+
+        /* Start looping on the number of vertices: */
+        GaussPenta gauss;
+        for(int iv=0;iv<3;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;
+
+                /* no calving for grounded ice in EXP4 */
+                gr_input->GetInputValue(&groundedice,&gauss);
+
+                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:
+                                /* Exp 2: set c=v-wrate(given)*/
+                                wrate_input->GetInputValue(&wrate,&gauss);
+                                break;
+                        case 4:
+                                /* Exp 4: set c=v-wrate(given), for the first 500 years, then c=0 for the second 500 years*/
+                                if((groundedice<0) && (time<=500.0*yts)) {
+                                //      wrate_input->GetInputValue(&wrate,&gauss);
+                                        wrate = -750*sin(2.0*M_PI*time/yts/1000)/yts;  // m/a -> m/s
+                                }
+                                else {
+                                        /* no calving on the grounded ice*/
+                                        wrate = vel;
+                                }
+                                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->AddBasalInput(CalvingCalvingrateEnum,&calvingrate[0],P1DGEnum);
+        this->AddBasalInput(CalvingratexEnum,&calvingratex[0],P1DGEnum);
+        this->AddBasalInput(CalvingrateyEnum,&calvingratey[0],P1DGEnum);
+
+        /*Extrude*/
+        this->InputExtrude(CalvingCalvingrateEnum,-1);
+        this->InputExtrude(CalvingratexEnum,-1);
+        this->InputExtrude(CalvingrateyEnum,-1);
 }
 /*}}}*/
@@ -2810 +2901 @@
                 case CalvingVonmisesEnum:
                 case CalvingLevermannEnum:
+                case CalvingCalvingMIPEnum:
                         calvingratex_input=this->GetInput(CalvingratexEnum); _assert_(calvingratex_input);
                         calvingratey_input=this->GetInput(CalvingrateyEnum); _assert_(calvingratey_input);
@@ -2847 +2939 @@
                         case CalvingVonmisesEnum:
                         case CalvingLevermannEnum:
+                        case CalvingCalvingMIPEnum:
                                 calvingratex_input->GetInputValue(&c[0],&gauss);
                                 calvingratey_input->GetInputValue(&c[1],&gauss);

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

@@ -58 +58 @@
                 void           CalvingFluxLevelset();
                 void           CalvingMeltingFluxLevelset();
+                void           CalvingRateCalvingMIP();
                 IssmDouble     CharacteristicLength(void){_error_("not implemented yet");};
                 void           ComputeBasalStress(void);