Changeset 17850

Timestamp:

04/25/14 13:51:31 (11 years ago)

Author:

Mathieu Morlighem

Message:

CHG: added new enums for balancethickness2solution

Location:

issm/trunk-jpl/src

Files:

• c/analyses/Balancethickness2Analysis.cpp (modified) (14 diffs)
• c/analyses/Balancethickness2Analysis.h (modified) (1 diff)
• c/analyses/BalancethicknessAnalysis.cpp (modified) (1 diff)
• c/cores/balancethickness2_core.cpp (modified) (1 diff)
• c/modules/ModelProcessorx/CreateParameters.cpp (modified) (2 diffs)
• c/shared/Enum/EnumDefinitions.h (modified) (1 diff)
• c/shared/Enum/EnumToStringx.cpp (modified) (1 diff)
• c/shared/Enum/StringToEnumx.cpp (modified) (4 diffs)
• m/classes/balancethickness2.m (added)
• m/enum/BalancethicknessApparentMassbalanceEnum.m (added)
• m/enum/BalancethicknessNuxEnum.m (added)
• m/enum/BalancethicknessNuyEnum.m (added)
• m/enum/BalancethicknessSpcpotentialEnum.m (added)
• m/enum/BalancethicknessThicknessObsEnum.m (added)
• m/enum/BalancethicknessVxObsEnum.m (added)
• m/enum/BalancethicknessVyObsEnum.m (added)
• m/enum/EnumDefinitions.py (modified) (1 diff)
• m/enum/NuXEnum.m (deleted)
• m/enum/NuYEnum.m (deleted)

issm/trunk-jpl/src/c/analyses/Balancethickness2Analysis.cpp ¶

@@ -26 +26 @@
         }
 
+        iomodel->FetchDataToInput(elements,BalancethicknessApparentMassbalanceEnum);
+        iomodel->FetchDataToInput(elements,BalancethicknessNuxEnum);
+        iomodel->FetchDataToInput(elements,BalancethicknessNuyEnum);
+        iomodel->FetchDataToInput(elements,BalancethicknessVxObsEnum);
+        iomodel->FetchDataToInput(elements,BalancethicknessVyObsEnum);
+        iomodel->FetchDataToInput(elements,BalancethicknessThicknessObsEnum);
+        iomodel->FetchDataToInput(elements,MeshVertexonboundaryEnum);
         iomodel->FetchDataToInput(elements,MaskIceLevelsetEnum);
-        iomodel->FetchDataToInput(elements,VxEnum);
-        iomodel->FetchDataToInput(elements,VyEnum);
-        iomodel->FetchDataToInput(elements,BasalforcingsMeltingRateEnum);
-        iomodel->FetchDataToInput(elements,SurfaceforcingsMassBalanceEnum);
-        iomodel->FetchDataToInput(elements,BalancethicknessThickeningRateEnum);
-        iomodel->FetchDataToInput(elements,MeshVertexonboundaryEnum);
-        iomodel->FetchDataToInput(elements,ThicknessEnum);
 }/*}}}*/
 void Balancethickness2Analysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/
@@ -42 +42 @@
 void Balancethickness2Analysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/
 
-        /*Do not add constraints in DG*/
-        IoModelToConstraintsx(constraints,iomodel,BalancethicknessSpcthicknessEnum,Balancethickness2AnalysisEnum,P1Enum);
+
+        int finiteelement = P1Enum;
+        IoModelToConstraintsx(constraints,iomodel,BalancethicknessSpcpotentialEnum,Balancethickness2AnalysisEnum,finiteelement);
 
 }/*}}}*/
@@ -72 +73 @@
         /*Initialize Element vector and other vectors*/
         ElementMatrix* Ke  = element->NewElementMatrix();
-        IssmDouble*    B   = xNew<IssmDouble>(2*numnodes);
-        IssmDouble     D[2][2];
 
         /*Retrieve all inputs and parameters*/
@@ -83 +82 @@
                 gauss->GaussPoint(ig);
 
-                element->JacobianDeterminant(&Jdet,xyz_list,gauss);
-
                 IssmDouble* dbasis=xNew<IssmDouble>(2*numnodes);
                 element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
+                element->JacobianDeterminant(&Jdet,xyz_list,gauss);
 
                 for(int i=0;i<numnodes;i++){
@@ -93 +91 @@
                         }
                 }
-                //GetB(B,element,xyz_list,gauss);
-
-                //D_scalar=gauss->weight*Jdet;
-
-                //D[0][0]=D_scalar;
-                //D[0][1]=0.;
-                //D[1][0]=0.;
-                //D[1][1]=D_scalar;
-                //TripleMultiply(B,2,numnodes,1,
-                //                      &D[0][0],2,2,0,
-                //                      B,2,numnodes,0,
-                //                      &Ke->values[0],1);
-        }
-
-        /*Clean up and return*/
-        xDelete<IssmDouble>(xyz_list);
-        xDelete<IssmDouble>(B);
+        }
+
+        /*Clean up and return*/
+        xDelete<IssmDouble>(xyz_list);
         delete gauss;
         return Ke;
@@ -128 +113 @@
 
         /*Intermediaries */
-        IssmDouble  dhdt,mb,ms,Jdet;
+        IssmDouble  adot,Jdet;
         IssmDouble* xyz_list = NULL;
 
@@ -140 +125 @@
         /*Retrieve all inputs and parameters*/
         element->GetVerticesCoordinates(&xyz_list);
-        Input* mb_input   = element->GetInput(BasalforcingsMeltingRateEnum);       _assert_(mb_input);
-        Input* ms_input   = element->GetInput(SurfaceforcingsMassBalanceEnum);     _assert_(ms_input);
-        Input* dhdt_input = element->GetInput(BalancethicknessThickeningRateEnum); _assert_(dhdt_input);
-
-        /*Initialize mb_correction to 0, do not forget!:*/
+        Input* adot_input   = element->GetInput(BalancethicknessApparentMassbalanceEnum); _assert_(adot_input);
+
         /* Start  looping on the number of gaussian points: */
         Gauss* gauss=element->NewGauss(2);
@@ -152 +134 @@
                 element->JacobianDeterminant(&Jdet,xyz_list,gauss);
                 element->NodalFunctions(basis,gauss);
-
-                ms_input->GetInputValue(&ms,gauss);
-                mb_input->GetInputValue(&mb,gauss);
-                dhdt_input->GetInputValue(&dhdt,gauss);
-
-                for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*(ms-mb-dhdt)*basis[i];
+                adot_input->GetInputValue(&adot,gauss);
+
+                for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*adot*basis[i];
         }
 
@@ -185 +164 @@
 
         /*Retrieve all inputs and parameters*/
-        Input* thickness_input = element->GetInput(ThicknessEnum); _assert_(thickness_input);
-        Input* vx_input        = element->GetInput(VxEnum);        _assert_(vx_input);
-        Input* vy_input        = element->GetInput(VyEnum);        _assert_(vy_input);
+        Input* thickness_input = element->GetInput(BalancethicknessThicknessObsEnum); _assert_(thickness_input);
+        Input* vx_input        = element->GetInput(BalancethicknessVxObsEnum);        _assert_(vx_input);
+        Input* vy_input        = element->GetInput(BalancethicknessVyObsEnum);        _assert_(vy_input);
 
         element->GetVerticesCoordinates(&xyz_list);
@@ -204 +183 @@
                 element->NodalFunctions(basis,gauss);
 
-                for(int i=0;i<numnodes;i++) pe->values[i] += Jdet*gauss->weight*thickness*(vx*normal[0] + vy*normal[1])*basis[i];
+                for(int i=0;i<numnodes;i++) pe->values[i] += - Jdet*gauss->weight*thickness*(vx*normal[0] + vy*normal[1])*basis[i];
         }
 
@@ -214 +193 @@
         return pe;
 }/*}}}*/
-void Balancethickness2Analysis::GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
-        /*Compute B'  matrix. B'=[B1' B2' B3'] where Bi' is of size 3*NDOF2. 
-         * For node i, Bi can be expressed in the actual coordinate system
-         * by: 
-         *       Bi=[ dN/dx ]
-         *          [ dN/dy ]
-         * where N is the finiteelement function for node i.
-         *
-         * We assume B' has been allocated already, of size: 3x(NDOF2*numnodes)
-         */
-
-        /*Fetch number of nodes for this finite element*/
-        int numnodes = element->GetNumberOfNodes();
-
-        /*Get nodal functions derivatives*/
-        IssmDouble* dbasis=xNew<IssmDouble>(2*numnodes);
-        element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
-
-        /*Build B': */
-        for(int i=0;i<numnodes;i++){
-                B[numnodes*0+i] = dbasis[0*numnodes+i];
-                B[numnodes*1+i] = dbasis[1*numnodes+i];
-        }
-
-        /*Clean-up*/
-        xDelete<IssmDouble>(dbasis);
-
-}/*}}}*/
 void Balancethickness2Analysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
            _error_("not implemented yet");
@@ -248 +199 @@
 
         /*Intermediaries */
-        IssmDouble  vx,vy,vel,h,normdphi,dphi[2];
+        IssmDouble  vx,vy,vbar,nux,nuy,normdphi,dphi[2];
         IssmDouble* xyz_list = NULL;
         int       * doflist  = NULL;
@@ -271 +222 @@
         }
 
-        element->AddInput(EnthalpyEnum,values,element->GetElementType());//FIXME
-
-        /*Retrieve all inputs and parameters*/
-        element->GetVerticesCoordinates(&xyz_list);
-        Input* potential_input = element->GetInput(EnthalpyEnum); _assert_(potential_input);
-        Input* vx_input        = element->GetInput(VxEnum);       _assert_(vx_input);
-        Input* vy_input        = element->GetInput(VyEnum);       _assert_(vy_input);
+        element->AddInput(PotentialEnum,values,element->GetElementType());
+
+        /*Retrieve all inputs and parameters*/
+        element->GetVerticesCoordinates(&xyz_list);
+        Input* potential_input = element->GetInput(PotentialEnum);             _assert_(potential_input);
+        Input* vx_input        = element->GetInput(BalancethicknessVxObsEnum); _assert_(vx_input);
+        Input* vy_input        = element->GetInput(BalancethicknessVyObsEnum); _assert_(vy_input);
+        Input* nux_input       = element->GetInput(BalancethicknessNuxEnum);   _assert_(nux_input);
+        Input* nuy_input       = element->GetInput(BalancethicknessNuyEnum);   _assert_(nuy_input);
 
         Gauss* gauss=element->NewGauss();
@@ -285 +238 @@
                 vx_input->GetInputValue(&vx,gauss);
                 vy_input->GetInputValue(&vy,gauss);
+                nux_input->GetInputValue(&nux,gauss);
+                nuy_input->GetInputValue(&nuy,gauss);
                 potential_input->GetInputDerivativeValue(&dphi[0],xyz_list,gauss);
 
-                vel = sqrt(vx*vx + vy*vy) + 1.e-10;
+                vx = vx*nux; vy = vy*nuy;
+                vbar = sqrt(vx*vx + vy*vy) + 1.e-10;
                 normdphi = sqrt(dphi[0]*dphi[0] + dphi[1]*dphi[1]);
 
-                thickness_list[iv] = normdphi/vel;
-                vx_list[iv]        = 1./thickness_list[iv] * dphi[0];
-                vy_list[iv]        = 1./thickness_list[iv] * dphi[1];
+                thickness_list[iv] = normdphi/vbar;
+                vx_list[iv]        = -1./thickness_list[iv] * dphi[0];
+                vy_list[iv]        = -1./thickness_list[iv] * dphi[1];
         }
         element->AddInput(ThicknessEnum,thickness_list,P1Enum);

issm/trunk-jpl/src/c/analyses/Balancethickness2Analysis.h ¶

@@ -28 +28 @@
                 ElementVector* CreatePVectorVolume(Element* element);
                 ElementVector* CreatePVectorBoundary(Element* element);
-                void GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss);
                 void GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element);
                 void InputUpdateFromSolution(IssmDouble* solution,Element* element);

issm/trunk-jpl/src/c/analyses/BalancethicknessAnalysis.cpp ¶

@@ -10 +10 @@
 }/*}}}*/
 void BalancethicknessAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
+        parameters->AddObject(iomodel->CopyConstantObject(BalancethicknessStabilizationEnum));
 }/*}}}*/
 void BalancethicknessAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/

issm/trunk-jpl/src/c/cores/balancethickness2_core.cpp ¶

@@ -26 +26 @@
         if(save_results){
                 if(VerboseSolution()) _printf0_("   saving results\n");
-                int outputs[4] = {ThicknessEnum,EnthalpyEnum,VxEnum,VyEnum};
+                int outputs[4] = {ThicknessEnum,PotentialEnum,VxEnum,VyEnum};
                 femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],4);
         }

issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp ¶

@@ -57 +57 @@
         parameters->AddObject(iomodel->CopyConstantObject(MeshAverageVertexConnectivityEnum));
         parameters->AddObject(iomodel->CopyConstantObject(ConstantsReferencetemperatureEnum));
-        parameters->AddObject(iomodel->CopyConstantObject(BalancethicknessStabilizationEnum));
         parameters->AddObject(iomodel->CopyConstantObject(GroundinglineMeltingRateEnum));
         parameters->AddObject(iomodel->CopyConstantObject(SettingsWaitonlockEnum));
@@ -219 +218 @@
         /*}}}*/
 
-        
         /*Before returning, create parameters in case we are running Qmu or control types runs: */
         CreateParametersControl(parameters,iomodel,solution_type);

issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h ¶

@@ -288 +288 @@
         TransientRequestedOutputsEnum,
         PotentialEnum,
+        BalancethicknessSpcpotentialEnum,
+        BalancethicknessApparentMassbalanceEnum,
         Balancethickness2MisfitEnum,
-        NuXEnum,
-        NuYEnum,
+        BalancethicknessNuxEnum,
+        BalancethicknessNuyEnum,
+        BalancethicknessVxObsEnum,
+        BalancethicknessVyObsEnum,
+        BalancethicknessThicknessObsEnum,
         /*}}}*/
         /*Surfaceforcings{{{*/

issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp ¶

@@ -296 +296 @@
                 case TransientRequestedOutputsEnum : return "TransientRequestedOutputs";
                 case PotentialEnum : return "Potential";
+                case BalancethicknessSpcpotentialEnum : return "BalancethicknessSpcpotential";
+                case BalancethicknessApparentMassbalanceEnum : return "BalancethicknessApparentMassbalance";
                 case Balancethickness2MisfitEnum : return "Balancethickness2Misfit";
-                case NuXEnum : return "NuX";
-                case NuYEnum : return "NuY";
+                case BalancethicknessNuxEnum : return "BalancethicknessNux";
+                case BalancethicknessNuyEnum : return "BalancethicknessNuy";
+                case BalancethicknessVxObsEnum : return "BalancethicknessVxObs";
+                case BalancethicknessVyObsEnum : return "BalancethicknessVyObs";
+                case BalancethicknessThicknessObsEnum : return "BalancethicknessThicknessObs";
                 case SurfaceforcingsEnum : return "Surfaceforcings";
                 case SMBEnum : return "SMB";

issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp ¶

@@ -302 +302 @@
               else if (strcmp(name,"TransientRequestedOutputs")==0) return TransientRequestedOutputsEnum;
               else if (strcmp(name,"Potential")==0) return PotentialEnum;
+              else if (strcmp(name,"BalancethicknessSpcpotential")==0) return BalancethicknessSpcpotentialEnum;
+              else if (strcmp(name,"BalancethicknessApparentMassbalance")==0) return BalancethicknessApparentMassbalanceEnum;
               else if (strcmp(name,"Balancethickness2Misfit")==0) return Balancethickness2MisfitEnum;
-              else if (strcmp(name,"NuX")==0) return NuXEnum;
-              else if (strcmp(name,"NuY")==0) return NuYEnum;
+              else if (strcmp(name,"BalancethicknessNux")==0) return BalancethicknessNuxEnum;
+              else if (strcmp(name,"BalancethicknessNuy")==0) return BalancethicknessNuyEnum;
+              else if (strcmp(name,"BalancethicknessVxObs")==0) return BalancethicknessVxObsEnum;
+              else if (strcmp(name,"BalancethicknessVyObs")==0) return BalancethicknessVyObsEnum;
+              else if (strcmp(name,"BalancethicknessThicknessObs")==0) return BalancethicknessThicknessObsEnum;
               else if (strcmp(name,"Surfaceforcings")==0) return SurfaceforcingsEnum;
               else if (strcmp(name,"SMB")==0) return SMBEnum;
@@ -378 +383 @@
               else if (strcmp(name,"MeshdeformationSolution")==0) return MeshdeformationSolutionEnum;
               else if (strcmp(name,"MeshdeformationAnalysis")==0) return MeshdeformationAnalysisEnum;
-              else if (strcmp(name,"LevelsetAnalysis")==0) return LevelsetAnalysisEnum;
+         else stage=4;
+   }
+   if(stage==4){
+              if (strcmp(name,"LevelsetAnalysis")==0) return LevelsetAnalysisEnum;
               else if (strcmp(name,"ExtrapolationAnalysis")==0) return ExtrapolationAnalysisEnum;
               else if (strcmp(name,"LsfReinitializationAnalysis")==0) return LsfReinitializationAnalysisEnum;
               else if (strcmp(name,"Approximation")==0) return ApproximationEnum;
               else if (strcmp(name,"NoneApproximation")==0) return NoneApproximationEnum;
-         else stage=4;
-   }
-   if(stage==4){
-              if (strcmp(name,"SIAApproximation")==0) return SIAApproximationEnum;
+              else if (strcmp(name,"SIAApproximation")==0) return SIAApproximationEnum;
               else if (strcmp(name,"SSAApproximation")==0) return SSAApproximationEnum;
               else if (strcmp(name,"SSAHOApproximation")==0) return SSAHOApproximationEnum;
@@ -501 +506 @@
               else if (strcmp(name,"ResetPenalties")==0) return ResetPenaltiesEnum;
               else if (strcmp(name,"SegmentOnIceShelf")==0) return SegmentOnIceShelfEnum;
-              else if (strcmp(name,"SurfaceAbsVelMisfit")==0) return SurfaceAbsVelMisfitEnum;
+         else stage=5;
+   }
+   if(stage==5){
+              if (strcmp(name,"SurfaceAbsVelMisfit")==0) return SurfaceAbsVelMisfitEnum;
               else if (strcmp(name,"SurfaceArea")==0) return SurfaceAreaEnum;
               else if (strcmp(name,"SurfaceAverageVelMisfit")==0) return SurfaceAverageVelMisfitEnum;
               else if (strcmp(name,"SurfaceLogVelMisfit")==0) return SurfaceLogVelMisfitEnum;
               else if (strcmp(name,"SurfaceLogVxVyMisfit")==0) return SurfaceLogVxVyMisfitEnum;
-         else stage=5;
-   }
-   if(stage==5){
-              if (strcmp(name,"SurfaceRelVelMisfit")==0) return SurfaceRelVelMisfitEnum;
+              else if (strcmp(name,"SurfaceRelVelMisfit")==0) return SurfaceRelVelMisfitEnum;
               else if (strcmp(name,"SurfaceSlopeX")==0) return SurfaceSlopeXEnum;
               else if (strcmp(name,"SurfaceSlopeY")==0) return SurfaceSlopeYEnum;
@@ -624 +629 @@
               else if (strcmp(name,"MinVz")==0) return MinVzEnum;
               else if (strcmp(name,"MaxVz")==0) return MaxVzEnum;
-              else if (strcmp(name,"MaxAbsVz")==0) return MaxAbsVzEnum;
+         else stage=6;
+   }
+   if(stage==6){
+              if (strcmp(name,"MaxAbsVz")==0) return MaxAbsVzEnum;
               else if (strcmp(name,"IceVolume")==0) return IceVolumeEnum;
               else if (strcmp(name,"IceVolumeAboveFloatation")==0) return IceVolumeAboveFloatationEnum;
               else if (strcmp(name,"TotalSmb")==0) return TotalSmbEnum;
               else if (strcmp(name,"Absolute")==0) return AbsoluteEnum;
-         else stage=6;
-   }
-   if(stage==6){
-              if (strcmp(name,"Incremental")==0) return IncrementalEnum;
+              else if (strcmp(name,"Incremental")==0) return IncrementalEnum;
               else if (strcmp(name,"AugmentedLagrangianR")==0) return AugmentedLagrangianREnum;
               else if (strcmp(name,"AugmentedLagrangianTheta")==0) return AugmentedLagrangianThetaEnum;

issm/trunk-jpl/src/m/enum/EnumDefinitions.py ¶

@@ -288 +288 @@
 def TransientRequestedOutputsEnum(): return StringToEnum("TransientRequestedOutputs")[0]
 def PotentialEnum(): return StringToEnum("Potential")[0]
+def BalancethicknessSpcpotentialEnum(): return StringToEnum("BalancethicknessSpcpotential")[0]
+def BalancethicknessApparentMassbalanceEnum(): return StringToEnum("BalancethicknessApparentMassbalance")[0]
 def Balancethickness2MisfitEnum(): return StringToEnum("Balancethickness2Misfit")[0]
-def NuXEnum(): return StringToEnum("NuX")[0]
-def NuYEnum(): return StringToEnum("NuY")[0]
+def BalancethicknessNuxEnum(): return StringToEnum("BalancethicknessNux")[0]
+def BalancethicknessNuyEnum(): return StringToEnum("BalancethicknessNuy")[0]
+def BalancethicknessVxObsEnum(): return StringToEnum("BalancethicknessVxObs")[0]
+def BalancethicknessVyObsEnum(): return StringToEnum("BalancethicknessVyObs")[0]
+def BalancethicknessThicknessObsEnum(): return StringToEnum("BalancethicknessThicknessObs")[0]
 def SurfaceforcingsEnum(): return StringToEnum("Surfaceforcings")[0]
 def SMBEnum(): return StringToEnum("SMB")[0]