Changeset 27905

Timestamp:

09/19/23 12:09:44 (18 months ago)

Author:

Mathieu Morlighem

Message:

BUG: wrong PETSc version

Location:

issm/trunk-jpl/src/c

Files:

• classes/Elements/Tria.cpp (modified) (6 diffs)
• toolkits/petsc/patches/PetscOptionsDetermineSolverType.cpp (modified) (1 diff)

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

@@ -29 +29 @@
 #define NUMVERTICES   3
 #define NUMVERTICES1D 2
-//#define MICI          1 //1 = DeConto & Pollard, 2 = DOMINOS
+#define MICI          2 //1 = DeConto & Pollard, 2 = Anna Crawford DOMINOS
 
 /*Constructors/destructor/copy*/
@@ -2772 +2772 @@
 IssmDouble Tria::GetIcefrontArea(){/*{{{*/
 
-        IssmDouble  bed[NUMVERTICES]; //basinId[NUMVERTICES];
+        IssmDouble  bed[NUMVERTICES];
         IssmDouble      Haverage,frontarea;
         IssmDouble  x1,y1,x2,y2,distance;
         IssmDouble lsf[NUMVERTICES], Haux[NUMVERTICES], surfaces[NUMVERTICES], bases[NUMVERTICES];
         int* indices=NULL;
-        IssmDouble* H=NULL;;
-        int nrfrontbed,numiceverts;
-
-        if(!IsZeroLevelset(MaskIceLevelsetEnum)) return 0;
+
+        /*Return if no ice front present*/
+        if(!this->IsIcefront()) return 0.;
 
         /*Retrieve all inputs and parameters*/
@@ -2788 +2787 @@
         Element::GetInputListOnVertices(&lsf[0],MaskIceLevelsetEnum);
 
-        nrfrontbed=0;
-        for(int i=0;i<NUMVERTICES;i++){
-                /*Find if bed<0*/
-                if(bed[i]<0.) nrfrontbed++;
-        }
-
-        if(nrfrontbed==3){
-                /*2. Find coordinates of where levelset crosses 0*/
-                int         numiceverts;
-                IssmDouble  s[2],x[2],y[2];
-                this->GetLevelsetIntersection(&indices, &numiceverts,&s[0],MaskIceLevelsetEnum,0.);
-                _assert_(numiceverts);
-
-                /*3 Write coordinates*/
-                IssmDouble  xyz_list[NUMVERTICES][3];
-                ::GetVerticesCoordinates(&xyz_list[0][0],this->vertices,NUMVERTICES);
-                int counter = 0;
-                if((numiceverts>0) && (numiceverts<NUMVERTICES)){
-                        for(int i=0;i<numiceverts;i++){
-                                for(int n=numiceverts;n<NUMVERTICES;n++){ // iterate over no-ice vertices
-                                        x[counter] = xyz_list[indices[i]][0]+s[counter]*(xyz_list[indices[n]][0]-xyz_list[indices[i]][0]);
-                                        y[counter] = xyz_list[indices[i]][1]+s[counter]*(xyz_list[indices[n]][1]-xyz_list[indices[i]][1]);
-                                        counter++;
+        /*Only continue if all 3 vertices are below sea level*/
+        for(int i=0;i<NUMVERTICES;i++) if(bed[i]>=0.) return 0.;
+
+        /*2. Find coordinates of where levelset crosses 0*/
+        int         numiceverts;
+        IssmDouble  s[2],x[2],y[2];
+        this->GetLevelsetIntersection(&indices, &numiceverts, &s[0],MaskIceLevelsetEnum,0.);
+        _assert_(numiceverts);
+        if(numiceverts>2){
+                Input* ls_input = this->GetInput(MaskIceLevelsetEnum);
+                ls_input->Echo();
+        }
+
+        /*3 Write coordinates*/
+        IssmDouble  xyz_list[NUMVERTICES][3];
+        ::GetVerticesCoordinates(&xyz_list[0][0],this->vertices,NUMVERTICES);
+        int counter = 0;
+        if((numiceverts>0) && (numiceverts<NUMVERTICES)){
+                for(int i=0;i<numiceverts;i++){
+                        for(int n=numiceverts;n<NUMVERTICES;n++){ // iterate over no-ice vertices
+                                x[counter] = xyz_list[indices[i]][0]+s[counter]*(xyz_list[indices[n]][0]-xyz_list[indices[i]][0]);
+                                y[counter] = xyz_list[indices[i]][1]+s[counter]*(xyz_list[indices[n]][1]-xyz_list[indices[i]][1]);
+                                counter++;
+                        }
+                }
+        }
+        else if(numiceverts==NUMVERTICES){ //NUMVERTICES ice vertices: calving front lies on element edge
+
+                for(int i=0;i<NUMVERTICES;i++){
+                        if(lsf[indices[i]]==0.){
+                                x[counter]=xyz_list[indices[i]][0];
+                                y[counter]=xyz_list[indices[i]][1];
+                                counter++;
+                        }
+                        if(counter==2) break;
+                }
+                if(counter==1){
+                        /*We actually have only 1 vertex on levelset, write a single point as a segment*/
+                        x[counter]=x[0];
+                        y[counter]=y[0];
+                        counter++;
+                }
+        }
+        else{
+                _error_("not sure what's going on here...");
+        }
+        x1=x[0]; y1=y[0]; x2=x[1]; y2=y[1];
+        distance=sqrt(pow((x1-x2),2)+pow((y1-y2),2));
+        if(distance<1e-3) return 0.;
+
+        IssmDouble H[4];
+        for(int iv=0;iv<NUMVERTICES;iv++) Haux[iv]=-bed[indices[iv]]; //sort bed in ice/noice
+        xDelete<int>(indices);
+
+        switch(numiceverts){
+                case 1: // average over triangle
+                        H[0]=Haux[0];
+                        H[1]=Haux[0]+s[0]*(Haux[1]-Haux[0]);
+                        H[2]=Haux[0]+s[1]*(Haux[2]-Haux[0]);
+                        Haverage=(H[1]+H[2])/2;
+                        break;
+                case 2: // average over quadrangle
+                        H[0]=Haux[0];
+                        H[1]=Haux[1];
+                        H[2]=Haux[0]+s[0]*(Haux[2]-Haux[0]);
+                        H[3]=Haux[1]+s[1]*(Haux[2]-Haux[1]);
+                        Haverage=(H[2]+H[3])/2;
+                        break;
+                case 3:
+                        if(counter==1) distance = 0; //front has 0 width on this element because levelset is 0 at a single vertex
+                        else if(counter==2){ //two vertices with levelset=0: averaging ice front depth over both
+                                Haverage = 0;
+                                for(int i=0;i<NUMVERTICES;i++){
+                                        if(lsf[indices[i]]==0.) Haverage -= Haux[indices[i]]/2;
+                                        if(Haverage<Haux[indices[i]]/2-1e-3) break; //done with the two vertices
                                 }
                         }
-                }
-                else if(numiceverts==NUMVERTICES){ //NUMVERTICES ice vertices: calving front lies on element edge
-
-                        for(int i=0;i<NUMVERTICES;i++){
-                                if(lsf[indices[i]]==0.){
-                                        x[counter]=xyz_list[indices[i]][0];
-                                        y[counter]=xyz_list[indices[i]][1];
-                                        counter++;
-                                }
-                                if(counter==2) break;
-                        }
-                        if(counter==1){
-                                /*We actually have only 1 vertex on levelset, write a single point as a segment*/
-                                x[counter]=x[0];
-                                y[counter]=y[0];
-                                counter++;
-                        }
-                }
-                else{
-                        _error_("not sure what's going on here...");
-                }
-                x1=x[0]; y1=y[0]; x2=x[1]; y2=y[1];
-                distance=sqrt(pow((x1-x2),2)+pow((y1-y2),2));
-
-                int numthk=numiceverts+2;
-                H=xNew<IssmDouble>(numthk);
-                for(int iv=0;iv<NUMVERTICES;iv++) Haux[iv]=-bed[indices[iv]]; //sort bed in ice/noice
-
-                switch(numiceverts){
-                        case 1: // average over triangle
-                                H[0]=Haux[0];
-                                H[1]=Haux[0]+s[0]*(Haux[1]-Haux[0]);
-                                H[2]=Haux[0]+s[1]*(Haux[2]-Haux[0]);
-                                Haverage=(H[1]+H[2])/2;
-                                break;
-                        case 2: // average over quadrangle
-                                H[0]=Haux[0];
-                                H[1]=Haux[1];
-                                H[2]=Haux[0]+s[0]*(Haux[2]-Haux[0]);
-                                H[3]=Haux[1]+s[1]*(Haux[2]-Haux[1]);
-                                Haverage=(H[2]+H[3])/2;
-                                break;
-                        case 3:
-                                if(counter==1) distance = 0; //front has 0 width on this element because levelset is 0 at a single vertex
-                                else if(counter==2){ //two vertices with levelset=0: averaging ice front depth over both
-                                        Haverage = 0;
-                                        for(int i=0;i<NUMVERTICES;i++){
-                                                if(lsf[indices[i]]==0.) Haverage -= Haux[indices[i]]/2;
-                                                if(Haverage<Haux[indices[i]]/2-1e-3) break; //done with the two vertices
-                                        }
-                                }
-                                break;
-                        default:
-                                _error_("Number of ice covered vertices wrong in Tria::GetIceFrontArea(void)");
-                                break;
-                }
-                frontarea=distance*Haverage;
-        }
-        else return 0;
-
-        xDelete<int>(indices);
-        xDelete<IssmDouble>(H);
+                        break;
+                default:
+                        _error_("Number of ice covered vertices wrong in Tria::GetIceFrontArea(void)");
+                        break;
+        }
+        frontarea=distance*Haverage;
 
         _assert_(frontarea>0);
@@ -4733 +4728 @@
 
                 /*Do we assume that the calving front does not move if MICI is not engaged?*/
-                movingfrontvx[iv] = 0.;
-                movingfrontvy[iv] = 0.;
+                bool regrowth = false;
+                bool apply_as_retreat = true;
+                if(!regrowth){
+                        movingfrontvx[iv] = 0.;
+                        movingfrontvy[iv] = 0.;
+                }
+
                 //movingfrontvx[iv] = -2000./(365*24*3600.)*dlsf[0]/norm_dlsf;
                 //movingfrontvy[iv] = -2000./(365*24*3600.)*dlsf[1]/norm_dlsf;
@@ -4746 +4746 @@
                 }
                 else if (MICI==2 && Hc>135. && bed<0. && fabs(ls)<100.e3){ // Crawford et all
+
+                        /*if 1: RETREAT rate
+                         *if 0: calving rate*/
+                        if(0) v[0]=0.; v[1]=0.;
 
                         /*5C Bn (worst case scenario)*/
@@ -4755 +4759 @@
                         movingfrontvx[iv] = v[0] -C*dlsf[0]/norm_dlsf;
                         movingfrontvy[iv] = v[1] -C*dlsf[1]/norm_dlsf;
+
+                        /*disable regrowth if calving rate is too low*/
+                        if(!regrowth && C<vel){
+                                movingfrontvx[iv] = 0.;
+                                movingfrontvy[iv] = 0.;
+                        }
                 }
         }

issm/trunk-jpl/src/c/toolkits/petsc/patches/PetscOptionsDetermineSolverType.cpp

@@ -30 +30 @@
         /*retrieve mat_type option: */
         #if PETSC_VERSION_LT(3,7,0)
+        PetscOptionsGetString(PETSC_NULL,"-mat_type",&option[0],100,&flag);
+        #elif PETSC_VERSION_LT(3,19,0)
         PetscOptionsGetString(NULL,PETSC_NULL,"-mat_type",&option[0],100,&flag);
-        #elif PETSC_VERSION_LT(3,19,0)
-        PetscOptionsGetString(PETSC_NULL,"-mat_type",&option[0],100,&flag);
-        #else
+        #else /*newest version*/
         PetscOptionsGetString(NULL,PETSC_NULLPTR,"-mat_type",&option[0],100,&flag);
         #endif