Changeset 24373

Timestamp:

11/20/19 17:37:35 (5 years ago)

Author:

Mathieu Morlighem

Message:

BUG: fixed problem with calving rate

Location:

issm/trunk-jpl/src/c

Files:

• analyses/LevelsetAnalysis.cpp (modified) (4 diffs)
• classes/Elements/Penta.cpp (modified) (4 diffs)
• classes/Inputs2/PentaInput2.cpp (modified) (2 diffs)

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

@@ -1 +1 @@
 #ifdef HAVE_CONFIG_H
-   #include <config.h>
+#include <config.h>
 #else
 #error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
@@ -435 +435 @@
 
                         case CalvingDev2Enum:
-                                {
-                                lsf_slopex_input->GetInputValue(&dlsf[0],gauss);
-                                if(dim==2) lsf_slopey_input->GetInputValue(&dlsf[1],gauss);
-                                calvingrate_input->GetInputValue(&calvingrate,gauss);
-                                meltingrate_input->GetInputValue(&meltingrate,gauss);
-                                gr_input->GetInputValue(&groundedice,gauss);
-
-                                //idea: no retreat on ice above critical calving height "calvinghaf" . Limit using regularized Heaviside function.
-                                vel=sqrt(v[0]*v[0] + v[1]*v[1]);
-                                haf_eps=10.;
-                                if(groundedice-calvinghaf<=-haf_eps){
-                                        // ice floats freely below calvinghaf: calve freely
-                                        // undercutting has no effect:
-                                        meltingrate=0.;
-                                }
-                                else if(groundedice-calvinghaf>=haf_eps){
-                                        // ice is well above calvinghaf -> no calving back, i.e. limit calving rate to ice velocity
-                                        calvingrate=min(calvingrate,vel);
-                                        // ice is almost grounded: frontal undercutting has maximum effect (do nothing).
-                                }
-                                else{ // ice is close to calvinghaf: smooth transition between limitation and free calving.
-                                        //heaviside: 0 for floating, 1 for grounded
-                                        heaviside=(groundedice-calvinghaf+haf_eps)/(2.*haf_eps) + sin(PI*(groundedice-calvinghaf)/haf_eps)/(2.*PI);
-                                        calvingrate=heaviside*(min(calvingrate,vel)-calvingrate)+calvingrate;
-                                        meltingrate=heaviside*meltingrate+0.;
-                                }
-
-                                norm_dlsf=0.;
-                                for(i=0;i<dim;i++) norm_dlsf+=pow(dlsf[i],2);
-                                norm_dlsf=sqrt(norm_dlsf);
-
-                                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;
-                                }
+                                  {
+                                        lsf_slopex_input->GetInputValue(&dlsf[0],gauss);
+                                        if(dim==2) lsf_slopey_input->GetInputValue(&dlsf[1],gauss);
+                                        calvingrate_input->GetInputValue(&calvingrate,gauss);
+                                        meltingrate_input->GetInputValue(&meltingrate,gauss);
+                                        gr_input->GetInputValue(&groundedice,gauss);
+
+                                        //idea: no retreat on ice above critical calving height "calvinghaf" . Limit using regularized Heaviside function.
+                                        vel=sqrt(v[0]*v[0] + v[1]*v[1]);
+                                        haf_eps=10.;
+                                        if(groundedice-calvinghaf<=-haf_eps){
+                                                // ice floats freely below calvinghaf: calve freely
+                                                // undercutting has no effect:
+                                                meltingrate=0.;
+                                        }
+                                        else if(groundedice-calvinghaf>=haf_eps){
+                                                // ice is well above calvinghaf -> no calving back, i.e. limit calving rate to ice velocity
+                                                calvingrate=min(calvingrate,vel);
+                                                // ice is almost grounded: frontal undercutting has maximum effect (do nothing).
+                                        }
+                                        else{ // ice is close to calvinghaf: smooth transition between limitation and free calving.
+                                                //heaviside: 0 for floating, 1 for grounded
+                                                heaviside=(groundedice-calvinghaf+haf_eps)/(2.*haf_eps) + sin(PI*(groundedice-calvinghaf)/haf_eps)/(2.*PI);
+                                                calvingrate=heaviside*(min(calvingrate,vel)-calvingrate)+calvingrate;
+                                                meltingrate=heaviside*meltingrate+0.;
+                                        }
+
+                                        norm_dlsf=0.;
+                                        for(i=0;i<dim;i++) norm_dlsf+=pow(dlsf[i],2);
+                                        norm_dlsf=sqrt(norm_dlsf);
+
+                                        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;
+                                  }
 
                         default:
@@ -518 +518 @@
                         case 2:
                                   {
-                                /* Streamline Upwinding */
-                                basalelement->ElementSizes(&hx,&hy,&hz);
-                                h=sqrt( pow(hx*w[0]/vel,2) + pow(hy*w[1]/vel,2) );
-                                IssmDouble D[9];
-                                for(row=0;row<dim;row++)
-                                        for(col=0;col<dim;col++)
-                                                D[row*dim+col] = D_scalar*h/(2.*vel)*w[row]*w[col];
-                                GetBprime(Bprime,basalelement,xyz_list,gauss);
-
-                                TripleMultiply(Bprime,dim,numnodes,1,
-                                                        &D[0],dim,dim,0,
-                                                        Bprime,dim,numnodes,0,
-                                                        &Ke->values[0],1);
+                                        /* Streamline Upwinding */
+                                        basalelement->ElementSizes(&hx,&hy,&hz);
+                                        h=sqrt( pow(hx*w[0]/vel,2) + pow(hy*w[1]/vel,2) );
+                                        IssmDouble D[9];
+                                        for(row=0;row<dim;row++)
+                                         for(col=0;col<dim;col++)
+                                          D[row*dim+col] = D_scalar*h/(2.*vel)*w[row]*w[col];
+                                        GetBprime(Bprime,basalelement,xyz_list,gauss);
+
+                                        TripleMultiply(Bprime,dim,numnodes,1,
+                                                                &D[0],dim,dim,0,
+                                                                Bprime,dim,numnodes,0,
+                                                                &Ke->values[0],1);
                                   }
                                 break;
@@ -669 +669 @@
         norm_b=0.;
         for(i=0;i<dim;i++)
-                norm_b+=b[i]*b[i];
+         norm_b+=b[i]*b[i];
         norm_b=sqrt(norm_b);
         _assert_(norm_b>0.);

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

@@ -170 +170 @@
         if(!IsOnBase()) return;
         else{
-                if(interpolation_enum==P1Enum){
+                if(interpolation_enum==P1Enum || interpolation_enum==P1DGEnum){
                         IssmDouble extrudedvalues[NUMVERTICES];
                         for(int i=0;i<NUMVERTICES2D;i++){
@@ -178 +178 @@
                         Penta* penta=this;
                         for(;;){
-                                penta->AddInput2(input_enum,&extrudedvalues[0],P1Enum);
+                                penta->AddInput2(input_enum,&extrudedvalues[0],interpolation_enum);
                                 if (penta->IsOnSurface()) break;
                                 penta=penta->GetUpperPenta(); _assert_(penta->Id()!=this->id);
@@ -382 +382 @@
 
         /*Add input*/
-        this->AddInput2(CalvingratexEnum,&calvingratex[0],P1DGEnum);
-        this->AddInput2(CalvingrateyEnum,&calvingratey[0],P1DGEnum);
-        this->AddInput2(CalvingCalvingrateEnum,&calvingrate[0],P1DGEnum);
-        this->AddInput2(SigmaVMEnum,&sigma_vm[0],P1DGEnum);
+        this->AddBasalInput2(CalvingratexEnum,&calvingratex[0],P1DGEnum);
+        this->AddBasalInput2(CalvingrateyEnum,&calvingratey[0],P1DGEnum);
+        this->AddBasalInput2(CalvingCalvingrateEnum,&calvingrate[0],P1DGEnum);
+        this->AddBasalInput2(SigmaVMEnum,&sigma_vm[0],P1DGEnum);
 
         this->InputExtrude(CalvingratexEnum,-1);
@@ -441 +441 @@
 
         /*Add input*/
-        this->AddInput2(CalvingratexEnum,&calvingratex[0],P1DGEnum);
-        this->AddInput2(CalvingrateyEnum,&calvingratey[0],P1DGEnum);
-        this->AddInput2(CalvingCalvingrateEnum,&calvingrate[0],P1DGEnum);
+        this->AddBasalInput2(CalvingratexEnum,&calvingratex[0],P1DGEnum);
+        this->AddBasalInput2(CalvingrateyEnum,&calvingratey[0],P1DGEnum);
+        this->AddBasalInput2(CalvingCalvingrateEnum,&calvingrate[0],P1DGEnum);
 
         /*Clean up and return*/

issm/trunk-jpl/src/c/classes/Inputs2/PentaInput2.cpp

@@ -110 +110 @@
         if(isserved){
                 _printf_("   current values:      ");
-                _printf_("[ ");
-                for(int i=0;i<6;i++) _printf_(" "<<this->element_values[i]);
-                _printf_("] ("<<EnumToStringx(this->interpolation)<<")\n");
+                if(isserved_collapsed){
+                        _printf_("[ ");
+                        for(int i=0;i<6;i++) _printf_(" "<<this->element_values[i]);
+                        _printf_("] ("<<EnumToStringx(this->interpolation)<<")\n");
+                }
+                else{
+                        _printf_("[ ");
+                        for(int i=0;i<3;i++) _printf_(" "<<this->element_values[i]);
+                        _printf_("] ("<<EnumToStringx(this->interpolation)<<")\n");
+                }
         }
 }
@@ -240 +247 @@
         if(state==1){
                 for(int i=0;i<3;i++) this->element_values[i] = this->values[row*this->N+i];
+                for(int i=3;i<6;i++) this->element_values[i] = 0.;
         }
         else if(state==2){
                 for(int i=0;i<3;i++) this->element_values[i] = this->values[row*this->N+3+i];
+                for(int i=3;i<6;i++) this->element_values[i] = 0.;
         }
         else{