Changeset 18604

Timestamp:

10/08/14 14:28:10 (10 years ago)

Author:

Mathieu Morlighem

Message:

CHG: working on balance thickness2

Location:

issm/trunk-jpl/src

Files:

• c/analyses/AdjointBalancethickness2Analysis.cpp (modified) (1 diff)
• c/analyses/Balancethickness2Analysis.cpp (modified) (5 diffs)
• c/analyses/SmoothAnalysis.cpp (modified) (2 diffs)
• c/cores/balancethickness2_core.cpp (modified) (2 diffs)
• c/cores/balancevelocity_core.cpp (modified) (1 diff)
• c/cores/controlvalidation_core.cpp (modified) (4 diffs)
• c/modules/ModelProcessorx/Control/UpdateElementsAndMaterialsControl.cpp (modified) (1 diff)
• 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/enum/EnumDefinitions.py (modified) (1 diff)
• m/enum/SmoothThicknessMultiplierEnum.m (added)
• m/enum/SmoothedSurfaceSlopeXAnalysisEnum.m (deleted)
• m/enum/SmoothedSurfaceSlopeYAnalysisEnum.m (deleted)

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

@@ -195 +195 @@
 void AdjointBalancethickness2Analysis::GradientJdHdt(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
 
+        /*Intermediaries*/
+        IssmDouble lambda,Jdet; 
+        IssmDouble *xyz_list= NULL;
+
         /*Fetch number of vertices for this finite element*/
         int numvertices = element->GetNumberOfVertices();
 
         /*Initialize some vectors*/
+        IssmDouble* basis         = xNew<IssmDouble>(numvertices);
         IssmDouble* ge            = xNewZeroInit<IssmDouble>(numvertices);
-        IssmDouble* lambda        = xNew<IssmDouble>(numvertices);
         int*        vertexpidlist = xNew<int>(numvertices);
 
         /*Retrieve all inputs we will be needing: */
+        element->GetVerticesCoordinates(&xyz_list);
         element->GradientIndexing(&vertexpidlist[0],control_index);
-        element->GetInputListOnVertices(lambda,AdjointEnum);
-        for(int i=0;i<numvertices;i++){
-                ge[i]= - lambda[i];
-                _assert_(!xIsNan<IssmDouble>(ge[i]));
-        }
-        gradient->SetValues(numvertices,vertexpidlist,ge,INS_VAL);
+        Input* adjoint_input = element->GetInput(AdjointEnum);            _assert_(adjoint_input);
+
+        Gauss* gauss=element->NewGauss(2);
+        for(int ig=gauss->begin();ig<gauss->end();ig++){
+                gauss->GaussPoint(ig);
+
+                element->JacobianDeterminant(&Jdet,xyz_list,gauss);
+                element->NodalFunctionsP1(basis,gauss);
+
+                adjoint_input->GetInputValue(&lambda,gauss);
+
+                /*Build gradient vector (actually -dJ/da): */
+                for(int i=0;i<numvertices;i++){
+                        ge[i]+= -Jdet*gauss->weight*basis[i]*lambda;
+                        _assert_(!xIsNan<IssmDouble>(ge[i]));
+                }
+        }
+        gradient->SetValues(numvertices,vertexpidlist,ge,ADD_VAL);
 
         /*Clean up and return*/
         xDelete<IssmDouble>(ge);
-        xDelete<IssmDouble>(lambda);
+        xDelete<IssmDouble>(xyz_list);
+        xDelete<IssmDouble>(basis);
         xDelete<int>(vertexpidlist);
+        delete gauss;
 }/*}}}*/
 void AdjointBalancethickness2Analysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/

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

@@ -82 +82 @@
         Input* omega_input = element->GetInput(BalancethicknessOmegaEnum); _assert_(omega_input);
         Input* D0_input    = element->GetInput(BalancethicknessD0Enum);
-        if(!D0_input){
+        //if(!D0_input){
                 this->CreateD0(element);
                 D0_input = element->GetInput(BalancethicknessD0Enum); _assert_(D0_input);
-        }
+        //}
 
         /* Start  looping on the number of gaussian points: */
@@ -220 +220 @@
 
         /*Intermediaries */
-        IssmDouble       Gamma,h,mu0,ds[2],Cmu,B,k,s,b;
+        IssmDouble       Gamma,h,mu0,ds[2],Cmu,B,k,s,b,normds;
         const int        n = 3;
         const IssmDouble Hstar = 500.;
@@ -232 +232 @@
         /*retrieve what we need: */
         element->GetVerticesCoordinates(&xyz_list);
-        Input* surface_input    = element->GetInput(SurfaceEnum);            _assert_(surface_input);
-        Input* B_input          = element->GetInput(MaterialsRheologyBEnum); _assert_(B_input);
-        IssmDouble rhog = element->GetMaterialParameter(MaterialsRhoIceEnum)*element->GetMaterialParameter(ConstantsGEnum);
+        Input* surfaceslopex_input = element->GetInput(SurfaceSlopeXEnum);
+        Input* surfaceslopey_input = element->GetInput(SurfaceSlopeYEnum);
+        Input* surface_input       = element->GetInput(SurfaceEnum);            _assert_(surface_input);
+        Input* B_input             = element->GetInput(MaterialsRheologyBEnum); _assert_(B_input);
+        IssmDouble rhog            = element->GetMaterialParameter(MaterialsRhoIceEnum)*element->GetMaterialParameter(ConstantsGEnum);
 
         /*Calculate damage evolution source term: */
@@ -242 +244 @@
                 
                 B_input->GetInputValue(&B,gauss);
-                surface_input->GetInputDerivativeValue(&ds[0],xyz_list,gauss);
+                if(surfaceslopex_input && surfaceslopey_input){
+                        surfaceslopex_input->GetInputValue(&ds[0],gauss);
+                        surfaceslopey_input->GetInputValue(&ds[1],gauss);
+                }
+                else{
+                        surface_input->GetInputDerivativeValue(&ds[0],xyz_list,gauss);
+                }
+
+                /*check slopes*/
+                normds = sqrt(ds[0]*ds[0]+ds[1]*ds[1]);
+                if (normds==0.){
+                        _error_("surface slope is zero");
+                }
+                if(normds<1.e-5){
+                        ds[0] = ds[0]/normds*1.e+5;
+                        ds[1] = ds[1]/normds*1.e+5;
+                        normds = 1.e-5;
+                }
 
                 mu0   = pow(2.,(1-3*n)/(2.*n)) * B;
@@ -252 +271 @@
         /*Add input*/
         element->AddInput(BalancethicknessD0Enum,D0,element->GetElementType());
+        //if(surfaceslopex_input && surfaceslopey_input){
+        //      element->DeleteInput(SurfaceSlopeXEnum);
+        //      element->DeleteInput(SurfaceSlopeYEnum);
+        //}
         
         /*Clean up and return*/

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

@@ -48 +48 @@
         /* Intermediaries */
         int         domaintype;
-        IssmDouble  Jdet,thickness,l=8.;
+        IssmDouble  Jdet,thickness,l;
         IssmDouble *xyz_list = NULL;
 
@@ -68 +68 @@
 
         /*Retrieve all inputs and parameters*/
+        element->FindParam(&l,SmoothThicknessMultiplierEnum); _assert_(l>0.);
         element->GetVerticesCoordinates(&xyz_list);
         Input* thickness_input = element->GetInput(ThicknessEnum); _assert_(thickness_input);

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

@@ -15 +15 @@
         bool save_results;
 
-        /*activate formulation: */
-        femmodel->SetCurrentConfiguration(Balancethickness2AnalysisEnum);
-
         /*recover parameters: */
         femmodel->parameters->FindParam(&save_results,SaveResultsEnum);
 
+        if(VerboseSolution()) _printf0_("computing smooth surface slopes:\n");
+        femmodel->parameters->SetParam(8.,SmoothThicknessMultiplierEnum);
+        femmodel->SetCurrentConfiguration(SmoothAnalysisEnum);
+        femmodel->parameters->SetParam(SurfaceSlopeXEnum,InputToSmoothEnum);
+        solutionsequence_linear(femmodel);
+        femmodel->parameters->SetParam(SurfaceSlopeYEnum,InputToSmoothEnum);
+        solutionsequence_linear(femmodel);
+
         if(VerboseSolution()) _printf0_("call computational core:\n");
+        femmodel->SetCurrentConfiguration(Balancethickness2AnalysisEnum);
         solutionsequence_linear(femmodel);
         //solutionsequence_nonlinear(femmodel,false);
@@ -27 +33 @@
         if(save_results){
                 if(VerboseSolution()) _printf0_("   saving results\n");
-                const int numoutputs = 4;
-                int outputs[numoutputs] = {SurfaceEnum,VxEnum,VyEnum,VelEnum};
+                const int numoutputs = 6;
+                int outputs[numoutputs] = {SurfaceEnum,SurfaceSlopeXEnum,SurfaceSlopeYEnum,VxEnum,VyEnum,VelEnum};
+                //const int numoutputs = 4;
+                //int outputs[numoutputs] = {SurfaceEnum,VxEnum,VyEnum,VelEnum};
                 femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],numoutputs);
         }

issm/trunk-jpl/src/c/cores/balancevelocity_core.cpp

@@ -18 +18 @@
         femmodel->parameters->FindParam(&save_results,SaveResultsEnum);
 
-        if(VerboseSolution()) _printf0_("computing smoothed slopes:\n");
+        if(VerboseSolution()) _printf0_("computing smooth driving stress:\n");
+        femmodel->parameters->SetParam(8.,SmoothThicknessMultiplierEnum);
         femmodel->SetCurrentConfiguration(SmoothAnalysisEnum);
         femmodel->parameters->SetParam(DrivingStressXEnum,InputToSmoothEnum);

issm/trunk-jpl/src/c/cores/controlvalidation_core.cpp

@@ -13 +13 @@
         int         solution_type,n;
         int         num_responses;
-        IssmDouble  j0,j;
+        IssmDouble  j0,j,yts;
         IssmDouble  Ialpha,exponent,alpha;
         IssmDouble* jlist = NULL;
@@ -28 +28 @@
         femmodel->parameters->SetParam(false,SaveResultsEnum);
         femmodel->parameters->FindParam(&num_responses,InversionNumCostFunctionsEnum);
+        femmodel->parameters->FindParam(&yts,ConstantsYtsEnum);
 
         /*Get initial guess*/
@@ -64 +65 @@
         IssmDouble exp0 = 0.;
         IssmDouble incr = -0.2;
-        IssmDouble exp1 = -10.;
+        IssmDouble exp1 = -18.;
         int        num  = reCast<int,IssmDouble>((exp1-exp0)/incr);
 
@@ -88 +89 @@
 
                 IssmDouble Den = 0.;
-                for(int i=0;i<n;i++) Den += alpha * G[i] * 1.;
+                for(int i=0;i<n;i++) Den += alpha* G[i] * 1.;
                 Ialpha = fabs((j - j0)/Den - 1.);
 

issm/trunk-jpl/src/c/modules/ModelProcessorx/Control/UpdateElementsAndMaterialsControl.cpp

@@ -29 +29 @@
         iomodel->FetchDataToInput(elements,InversionVyObsEnum,0.); 
         iomodel->FetchDataToInput(elements,InversionThicknessObsEnum,0.);
-        //iomodel->FetchDataToInput(elements,InversionSurfaceObsEnum,0.);
+        iomodel->FetchDataToInput(elements,InversionSurfaceObsEnum,0.);
 
         iomodel->FetchData(5,InversionControlParametersEnum,InversionCostFunctionsEnum,InversionCostFunctionsCoefficientsEnum,InversionMinParametersEnum,InversionMaxParametersEnum);

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

@@ -450 +450 @@
         InputToDepthaverageEnum,
         InputToSmoothEnum,
+        SmoothThicknessMultiplierEnum,
         IntParamEnum,
         IntVecParamEnum,

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

@@ -447 +447 @@
                 case InputToDepthaverageEnum : return "InputToDepthaverage";
                 case InputToSmoothEnum : return "InputToSmooth";
+                case SmoothThicknessMultiplierEnum : return "SmoothThicknessMultiplier";
                 case IntParamEnum : return "IntParam";
                 case IntVecParamEnum : return "IntVecParam";

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

@@ -456 +456 @@
               else if (strcmp(name,"InputToDepthaverage")==0) return InputToDepthaverageEnum;
               else if (strcmp(name,"InputToSmooth")==0) return InputToSmoothEnum;
+              else if (strcmp(name,"SmoothThicknessMultiplier")==0) return SmoothThicknessMultiplierEnum;
               else if (strcmp(name,"IntParam")==0) return IntParamEnum;
               else if (strcmp(name,"IntVecParam")==0) return IntVecParamEnum;
@@ -505 +506 @@
               else if (strcmp(name,"BalancethicknessMisfit")==0) return BalancethicknessMisfitEnum;
               else if (strcmp(name,"BedSlopeX")==0) return BedSlopeXEnum;
-              else if (strcmp(name,"BedSlopeY")==0) return BedSlopeYEnum;
          else stage=5;
    }
    if(stage==5){
-              if (strcmp(name,"Boundary")==0) return BoundaryEnum;
+              if (strcmp(name,"BedSlopeY")==0) return BedSlopeYEnum;
+              else if (strcmp(name,"Boundary")==0) return BoundaryEnum;
               else if (strcmp(name,"Converged")==0) return ConvergedEnum;
               else if (strcmp(name,"Fill")==0) return FillEnum;
@@ -628 +629 @@
               else if (strcmp(name,"Time")==0) return TimeEnum;
               else if (strcmp(name,"WaterColumnOld")==0) return WaterColumnOldEnum;
-              else if (strcmp(name,"Outputdefinition")==0) return OutputdefinitionEnum;
          else stage=6;
    }
    if(stage==6){
-              if (strcmp(name,"OutputdefinitionList")==0) return OutputdefinitionListEnum;
+              if (strcmp(name,"Outputdefinition")==0) return OutputdefinitionEnum;
+              else if (strcmp(name,"OutputdefinitionList")==0) return OutputdefinitionListEnum;
               else if (strcmp(name,"Massfluxatgate")==0) return MassfluxatgateEnum;
               else if (strcmp(name,"MassfluxatgateName")==0) return MassfluxatgateNameEnum;
@@ -751 +752 @@
               else if (strcmp(name,"BasalforcingsOceanVy")==0) return BasalforcingsOceanVyEnum;
               else if (strcmp(name,"SurfaceforcingsRhoAir")==0) return SurfaceforcingsRhoAirEnum;
-              else if (strcmp(name,"SurfaceforcingsAirCoef")==0) return SurfaceforcingsAirCoefEnum;
          else stage=7;
    }
    if(stage==7){
-              if (strcmp(name,"SurfaceforcingsAirLinDragCoef")==0) return SurfaceforcingsAirLinDragCoefEnum;
+              if (strcmp(name,"SurfaceforcingsAirCoef")==0) return SurfaceforcingsAirCoefEnum;
+              else if (strcmp(name,"SurfaceforcingsAirLinDragCoef")==0) return SurfaceforcingsAirLinDragCoefEnum;
               else if (strcmp(name,"SurfaceforcingsAirQuadDragCoef")==0) return SurfaceforcingsAirQuadDragCoefEnum;
               else if (strcmp(name,"SurfaceforcingsWindVx")==0) return SurfaceforcingsWindVxEnum;

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

@@ -439 +439 @@
 def InputToDepthaverageEnum(): return StringToEnum("InputToDepthaverage")[0]
 def InputToSmoothEnum(): return StringToEnum("InputToSmooth")[0]
+def SmoothThicknessMultiplierEnum(): return StringToEnum("SmoothThicknessMultiplier")[0]
 def IntParamEnum(): return StringToEnum("IntParam")[0]
 def IntVecParamEnum(): return StringToEnum("IntVecParam")[0]