Changeset 26148

Timestamp:

03/24/21 11:44:25 (4 years ago)

Author:

tsantos

Message:

CHG: working on mono layer HO

File:

• issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (modified) (12 diffs)

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

@@ -2707 +2707 @@
 /*MLHO*/
 ElementMatrix* StressbalanceAnalysis::CreateKMatrixMLHO(Element* element){/*{{{*/
-        _error_("not implemented yet");
-
+_error_("not implemented yet");
         /* Check if ice in element */
         if(!element->IsIceInElement()) return NULL;
@@ -2723 +2722 @@
 }/*}}}*/
 ElementMatrix* StressbalanceAnalysis::CreateKMatrixMLHOFriction(Element* element){/*{{{*/
-        _error_("not implemented yet");
 
         if(!element->IsOnBase() || element->IsFloating()) return NULL;
@@ -2746 +2744 @@
 }/*}}}*/
 ElementMatrix* StressbalanceAnalysis::CreateKMatrixMLHOViscous(Element* element){/*{{{*/
-        _error_("not implemented yet");
 
         /*Intermediaries*/
@@ -2861 +2858 @@
 }/*}}}*/
 ElementVector* StressbalanceAnalysis::CreatePVectorMLHO(Element* element){/*{{{*/
-        _error_("not implemented yet");
 
         /* Check if ice in element */
@@ -2895 +2891 @@
 }/*}}}*/
 ElementVector* StressbalanceAnalysis::CreatePVectorMLHODrivingStress(Element* element){/*{{{*/
-        _error_("not implemented yet");
 
         /*Intermediaries */
@@ -2943 +2938 @@
 }/*}}}*/
 ElementVector* StressbalanceAnalysis::CreatePVectorMLHOFront(Element* element){/*{{{*/
-        _error_("not implemented yet");
 
         /*If no front, return NULL*/
@@ -2992 +2986 @@
                 pressure = ice_pressure + water_pressure;
                 /*Vertically integrated pressure - HO type*/
-                s=max(0.,thickness+base-sealevel);
-                b=min(0.,base-sealevel); 
+                b=min(0.,base-sealevel); // 0 if the bottom of the glacier is above water level
+                s=thickness+b; // ice surface regardless of whether the top of the glacier is above water level or not
                 water_pressure_sh = gravity*rho_water*(-b*b/2 + pow(-b,n+2)*( -s/(n+2) -b/(n+3) )/pow(thickness,n+1));
                 ice_pressure_sh   = gravity*rho_ice*thickness*thickness*(n+1)/(2*(n+3));
@@ -3017 +3011 @@
 }/*}}}*/
 void           StressbalanceAnalysis::InputUpdateFromSolutionMLHO(IssmDouble* solution,Element* element){/*{{{*/
-
-        _error_("not implemented yet");
 
         int         i,dim,domaintype;
@@ -3056 +3048 @@
                         break;
                 case Domain3DEnum:
-                        if(!element->IsOnBase()){xDelete<IssmDouble>(xyz_list); return;}
-                        basalelement=element->SpawnBasalElement();
-                        break;
+                        _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
+                        //if(!element->IsOnBase()){xDelete<IssmDouble>(xyz_list); return;}
+                        //basalelement=element->SpawnBasalElement();
+                        //break;
                 default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
         }
@@ -3064 +3057 @@
         /*Fetch number of nodes and dof for this finite element*/
         int numnodes = basalelement->GetNumberOfNodes();
-        int numdof   = numnodes*2;
+        int numdof   = numnodes*2*2; //4 DOFs per node
 
         /*Fetch dof list and allocate solution vectors*/
-        basalelement->GetDofListLocal(&doflist,NoneApproximationEnum,GsetEnum);
+        basalelement->GetDofListLocal(&doflist,MLHOApproximationEnum,GsetEnum); //itapopo check this
+        //basalelement->GetDofListLocal(&doflist,NoneApproximationEnum,GsetEnum);
         IssmDouble* values    = xNew<IssmDouble>(numdof);
         IssmDouble* vx        = xNew<IssmDouble>(numnodes);
         IssmDouble* vy        = xNew<IssmDouble>(numnodes);
+   IssmDouble* vshx      = xNew<IssmDouble>(numnodes); 
+   IssmDouble* vshy      = xNew<IssmDouble>(numnodes);
         IssmDouble* vz        = xNew<IssmDouble>(numnodes);
         IssmDouble* vel       = xNew<IssmDouble>(numnodes);
-
+        
         /*Use the dof list to index into the solution vector: */
         for(i=0;i<numdof;i++) values[i]=solution[doflist[i]];
@@ -3079 +3075 @@
         /*Transform solution in Cartesian Space*/
         basalelement->TransformSolutionCoord(&values[0],XYEnum);
-        basalelement->FindParam(&domaintype,DomainTypeEnum);
-
-        /*Ok, we have vx and vy in values, fill in vx and vy arrays: */
-        for(i=0;i<numnodes;i++){
-                vx[i]=values[i*2+0];
-                vy[i]=values[i*2+1];
-
-                /*Check solution*/
-                if(xIsNan<IssmDouble>(vx[i])) _error_("NaN found in solution vector");
-                if(xIsInf<IssmDouble>(vx[i])) _error_("Inf found in solution vector");
-                if(xIsNan<IssmDouble>(vy[i])) _error_("NaN found in solution vector");
-                if(xIsInf<IssmDouble>(vy[i])) _error_("Inf found in solution vector");
+        //basalelement->FindParam(&domaintype,DomainTypeEnum);
+
+   /*Ok, we have vx and vy in values, fill in vx and vy arrays: */
+   for(i=0;i<numnodes;i++){ //numnodes of the 2D mesh in which the MLHO is written
+      vx[i]  =values[i+0*numnodes]; //basal vx
+      vshx[i]=values[i+1*numnodes];
+                if(xIsNan<IssmDouble>(vx[i]))           _error_("NaN found in solution vector");
+      if(xIsInf<IssmDouble>(vx[i]))             _error_("Inf found in solution vector");
+                if(xIsNan<IssmDouble>(vshx[i])) _error_("NaN found in solution vector");
+      if(xIsInf<IssmDouble>(vshx[i]))   _error_("Inf found in solution vector");
+      //if(dim==3){
+         vy[i] =values[i+2*numnodes];
+         vshy[i]=values[i+3*numnodes];
+         if(xIsNan<IssmDouble>(vy[i]))          _error_("NaN found in solution vector");
+         if(xIsInf<IssmDouble>(vy[i]))          _error_("Inf found in solution vector");
+         if(xIsNan<IssmDouble>(vshy[i]))        _error_("NaN found in solution vector");
+         if(xIsInf<IssmDouble>(vshy[i]))        _error_("Inf found in solution vector");
+      //}
+   }
+
+        /*Compute suface velocities vx and vy*/
+        if(domaintype==Domain2DhorizontalEnum) {
+                for(i=0;i<numnodes;i++) vx[i]=vx[i]+vshx[i];
+                for(i=0;i<numnodes;i++) vy[i]=vy[i]+vshy[i];
         }
 
         /*Get Vz and compute vel*/
         basalelement->GetInputListOnNodes(&vz[0],VzEnum,0.);
-        for(i=0;i<numnodes;i++) vel[i]=sqrt(vx[i]*vx[i] + vy[i]*vy[i] + vz[i]*vz[i]);
+        for(i=0;i<numnodes;i++) vel[i]=sqrt(vx[i]*vx[i] + vy[i]*vy[i] + vz[i]*vz[i]); 
 
         /*Add vx and vy as inputs to the tria element: */
@@ -3107 +3115 @@
         xDelete<IssmDouble>(vy);
         xDelete<IssmDouble>(vx);
+        xDelete<IssmDouble>(vshy);
+        xDelete<IssmDouble>(vshx);
         xDelete<IssmDouble>(values);
         xDelete<IssmDouble>(xyz_list);