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

Timestamp:

05/04/22 05:52:26 (3 years ago)

Author:

Mathieu Morlighem

Message:

CHG: moving more calving laws away from MovingFront

File:

• issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (modified) (6 diffs)

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

@@ -273 +273 @@
         this->ComputeSigmaVM();
 
-        IssmDouble  calvingratex[NUMVERTICES];
-        IssmDouble  calvingratey[NUMVERTICES];
         IssmDouble  calvingrate[NUMVERTICES];
         IssmDouble  vx,vy;
@@ -317 +315 @@
                 /*Assign values*/
                 if(bed>sealevel){
-                        calvingratex[iv]=0.;
-                        calvingratey[iv]=0.;
+                        calvingrate[iv] = 0.;
                 }
                 else{
-                        calvingratex[iv]=vx*sigma_vm/sigma_max;
-                        calvingratey[iv]=vy*sigma_vm/sigma_max;
-                }
-                calvingrate[iv] =sqrt(calvingratex[iv]*calvingratex[iv] + calvingratey[iv]*calvingratey[iv]);
+                        calvingrate[iv] = sqrt(vx*vx+vy*vy)*sigma_vm/sigma_max;
+                }
         }
 
         /*Add input*/
-        this->AddBasalInput(CalvingratexEnum,&calvingratex[0],P1DGEnum);
-        this->AddBasalInput(CalvingrateyEnum,&calvingratey[0],P1DGEnum);
         this->AddBasalInput(CalvingCalvingrateEnum,&calvingrate[0],P1DGEnum);
-
+        this->CalvingFromRate();
+
+        /*Extrude*/
+        this->InputExtrude(CalvingCalvingrateEnum,-1);
         this->InputExtrude(CalvingratexEnum,-1);
         this->InputExtrude(CalvingrateyEnum,-1);
-        this->InputExtrude(CalvingCalvingrateEnum,-1);
 }
 /*}}}*/
@@ -2871 +2866 @@
 }
 /*}}}*/
-void       Penta::MovingFrontalVelocity(void){/*{{{*/
+void          Penta::MovingFrontalVelocity(void){/*{{{*/
+
         if(!this->IsOnBase()) return;
         int        domaintype, calvinglaw, i;
@@ -2902 +2898 @@
                 case DefaultCalvingEnum:
                 case CalvingVonmisesEnum:
-                        calvingrate_input = this->GetInput(CalvingCalvingrateEnum);     _assert_(calvingrate_input);
-                        meltingrate_input = this->GetInput(CalvingMeltingrateEnum);     _assert_(meltingrate_input);
-                        break;
                 case CalvingLevermannEnum:
                         calvingratex_input=this->GetInput(CalvingratexEnum); _assert_(calvingratex_input);
@@ -2941 +2934 @@
                         case DefaultCalvingEnum:
                         case CalvingVonmisesEnum:
-                                calvingrate_input->GetInputValue(&calvingrate,&gauss);
-                                meltingrate_input->GetInputValue(&meltingrate,&gauss);
-
-                                if(groundedice<0) meltingrate = 0.;
-
-                                if(norm_dlsf>1.e-10)
-                                 for(i=0;i<dim;i++){
-                                         c[i]=calvingrate*dlsf[i]/norm_dlsf; m[i]=meltingrate*dlsf[i]/norm_dlsf;
-                                 }
-                                else
-                                 for(i=0;i<dim;i++){
-                                         c[i]=0.; m[i]=0.;
-                                 }
-                                break;
-
                         case CalvingLevermannEnum:
                                 calvingratex_input->GetInputValue(&c[0],&gauss);
                                 calvingratey_input->GetInputValue(&c[1],&gauss);
                                 meltingrate_input->GetInputValue(&meltingrate,&gauss);
-                                norm_calving=0.;
-                                for(i=0;i<dim;i++) norm_calving+=pow(c[i],2);
-                                norm_calving=sqrt(norm_calving)+1.e-14;
-                                for(i=0;i<dim;i++) m[i]=meltingrate*c[i]/norm_calving;
+                                if(groundedice<0) meltingrate = 0.;
+                                m[0]=meltingrate*dlsf[0]/norm_dlsf;
+                                m[1]=meltingrate*dlsf[1]/norm_dlsf;
                                 break;
 
@@ -3036 +3013 @@
         this->AddInput(MovingFrontalVxEnum,&movingfrontvx[0],P1DGEnum);
         this->AddInput(MovingFrontalVyEnum,&movingfrontvy[0],P1DGEnum);
-
         this->InputExtrude(MovingFrontalVxEnum,-1);
         this->InputExtrude(MovingFrontalVyEnum,-1);