Changeset 23867

Timestamp:

04/18/19 15:48:42 (6 years ago)

Author:

Mathieu Morlighem

Message:

NEW: added 508 cost function

Location:

issm/trunk-jpl/src/c

Files:

• analyses/AdjointHorizAnalysis.cpp (modified) (5 diffs)
• analyses/AdjointHorizAnalysis.h (modified) (1 diff)
• classes/FemModel.cpp (modified) (1 diff)
• modules/ModelProcessorx/Control/UpdateElementsAndMaterialsControl.cpp (modified) (1 diff)
• modules/RheologyBAbsGradientx/RheologyBAbsGradientx.cpp (modified) (1 diff)
• modules/RheologyBAbsGradientx/RheologyBAbsGradientx.h (modified) (1 diff)
• shared/Enum/EnumDefinitions.h (modified) (1 diff)
• shared/Enum/EnumToStringx.cpp (modified) (1 diff)
• shared/Enum/StringToEnumx.cpp (modified) (7 diffs)

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

@@ -553 +553 @@
                                         /*Nothing in P vector*/
                                         break;
+                                case RheologyBInitialguessMisfitEnum:
+                                        /*Nothing in P vector*/
+                                        break;
                                 default:
                                         _error_("response " << EnumToStringx(responses[resp]) << " not supported yet");
@@ -788 +791 @@
                                         /*Nothing in P vector*/
                                         break;
+                                case RheologyBInitialguessMisfitEnum:
+                                        /*Nothing in P vector*/
+                                        break;
                                 default:
                                         _error_("response " << EnumToStringx(responses[resp]) << " not supported yet");
@@ -1036 +1042 @@
                                         /*Nothing in P vector*/
                                         break;
+                                case RheologyBInitialguessMisfitEnum:
+                                        /*Nothing in P vector*/
+                                        break;
                                 default:
                                         _error_("response " << EnumToStringx(responses[resp]) << " not supported yet");
@@ -1097 +1106 @@
                 case RheologyBbarAbsGradientEnum:    GradientJBbarGradient(element,gradient,control_index); break;
                 case RheologyBAbsGradientEnum:       GradientJBGradient(element,gradient,control_index);    break;
+                case RheologyBInitialguessMisfitEnum:  GradientJBinitial(element,gradient,control_index);    break;
                 default: _error_("response " << EnumToStringx(responses[resp]) << " not supported yet");
         }
@@ -1398 +1408 @@
         delete gauss;
 }/*}}}*/
+void           AdjointHorizAnalysis::GradientJBinitial(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
+
+        /*Intermediaries*/
+        int      domaintype;
+
+        /*Get basal element*/
+        element->FindParam(&domaintype,DomainTypeEnum);
+        switch(domaintype){
+                case Domain2DhorizontalEnum:
+                        break;
+                case Domain2DverticalEnum:
+                        break;
+                case Domain3DEnum:
+                        break;
+                default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
+        }
+
+        /*Intermediaries*/
+        IssmDouble Jdet,weight;
+        IssmDouble B,B0; 
+        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);
+        int*        vertexpidlist = xNew<int>(numvertices);
+
+        /*Retrieve all inputs we will be needing: */
+        element->GetVerticesCoordinates(&xyz_list);
+        element->GradientIndexing(&vertexpidlist[0],control_index);
+        Input* rheology_input  = element->GetInput(MaterialsRheologyBbarEnum);              _assert_(rheology_input);
+        Input* rheology0_input = element->GetInput(RheologyBInitialguessEnum);              _assert_(rheology0_input);
+        Input* weights_input   = element->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
+
+        /* Start  looping on the number of gaussian points: */
+        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);
+                weights_input->GetInputValue(&weight,gauss,RheologyBInitialguessMisfitEnum);
+
+                /*Build alpha_complement_list: */
+                rheology_input->GetInputValue(&B,gauss);
+                rheology0_input->GetInputValue(&B0,gauss);
+
+                /*Build gradje_g_gaussian vector (actually -dJ/ddrag): */
+                for(int i=0;i<numvertices;i++){
+                        ge[i]+=-weight*Jdet*gauss->weight*basis[i]*(B-B0);
+                        _assert_(!xIsNan<IssmDouble>(ge[i]));
+                }
+        }
+        gradient->SetValues(numvertices,vertexpidlist,ge,ADD_VAL);
+
+        /*Clean up and return*/
+        xDelete<IssmDouble>(xyz_list);
+        xDelete<IssmDouble>(basis);
+        xDelete<IssmDouble>(ge);
+        xDelete<int>(vertexpidlist);
+        delete gauss;
+}/*}}}*/
 void           AdjointHorizAnalysis::GradientJBHO(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
         /*Intermediaries*/

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

@@ -38 +38 @@
                 void           GradientJBbarFS(Element* element,Vector<IssmDouble>* gradient,int control_index);
                 void           GradientJBbarGradient(Element* element,Vector<IssmDouble>* gradient,int control_index);
+                void           GradientJBinitial(Element* element,Vector<IssmDouble>* gradient,int control_index);
                 void           GradientJBbarL1L2(Element* element,Vector<IssmDouble>* gradient,int control_index);
                 void           GradientJBbarHO(Element* element,Vector<IssmDouble>* gradient,int control_index);

issm/trunk-jpl/src/c/classes/FemModel.cpp ¶

@@ -2223 +2223 @@
                                 case RheologyBbarAbsGradientEnum:   RheologyBbarAbsGradientx(&double_result,elements,nodes,vertices,loads,materials,parameters);    break;
                                 case RheologyBAbsGradientEnum:      RheologyBAbsGradientx(&double_result,elements,nodes,vertices,loads,materials,parameters);       break;
+                                case RheologyBInitialguessMisfitEnum:  RheologyBInitialguessMisfitx(&double_result,elements,nodes,vertices,loads,materials,parameters);       break;
                                 case DragCoefficientAbsGradientEnum:DragCoefficientAbsGradientx(&double_result,elements,nodes,vertices,loads,materials,parameters); break;
                                 case BalancethicknessMisfitEnum:    BalancethicknessMisfitx(&double_result);                                                        break;

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

@@ -60 +60 @@
                 if(     cost_function==ThicknessAbsMisfitEnum) iomodel->FetchDataToInput(elements,"md.inversion.thickness_obs",InversionThicknessObsEnum);
                 else if(cost_function==SurfaceAbsMisfitEnum)   iomodel->FetchDataToInput(elements,"md.inversion.surface_obs",InversionSurfaceObsEnum);
+                else if(cost_function==RheologyBInitialguessMisfitEnum) iomodel->FetchDataToInput(elements,"md.materials.rheology_B",RheologyBInitialguessEnum);
                 else if(cost_function==SurfaceAbsVelMisfitEnum
                           || cost_function==SurfaceRelVelMisfitEnum

issm/trunk-jpl/src/c/modules/RheologyBAbsGradientx/RheologyBAbsGradientx.cpp ¶

@@ -78 +78 @@
         return Jelem;
 }
+
+void RheologyBInitialguessMisfitx( IssmDouble* pJ, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials,Parameters* parameters){
+
+        /*output: */
+        IssmDouble J=0.;
+        IssmDouble J_sum;
+
+        /*Compute Misfit: */
+        for(int i=0;i<elements->Size();i++){
+                Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
+                J+=RheologyBInitialguessMisfit(element);
+        }
+
+        /*Sum all J from all cpus of the cluster:*/
+        ISSM_MPI_Reduce (&J,&J_sum,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
+        ISSM_MPI_Bcast(&J_sum,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
+        J=J_sum;
+
+        /*Assign output pointers: */
+        *pJ=J;
+}
+
+IssmDouble RheologyBInitialguessMisfit(Element* element){
+
+        int         domaintype,numcomponents;
+        IssmDouble  Jelem=0.;
+        IssmDouble  misfit,Jdet;
+        IssmDouble  B,B0,weight;
+        IssmDouble* xyz_list      = NULL;
+
+        /*If on water, return 0: */
+        if(!element->IsIceInElement()) return 0.;
+
+        /*Get problem dimension*/
+        element->FindParam(&domaintype,DomainTypeEnum);
+        switch(domaintype){
+                case Domain2DverticalEnum:   numcomponents   = 1; break;
+                case Domain3DEnum:           numcomponents   = 2; break;
+                case Domain2DhorizontalEnum: numcomponents   = 2; break;
+                default: _error_("not supported yet");
+        }
+
+        /* Get node coordinates*/
+        element->GetVerticesCoordinates(&xyz_list);
+
+        /*Retrieve all inputs we will be needing: */
+        Input* weights_input=element->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
+        Input* rheologyb_input=element->GetInput(MaterialsRheologyBbarEnum);            _assert_(rheologyb_input);
+        Input* rheologyb0_input=element->GetInput(RheologyBInitialguessEnum);           _assert_(rheologyb0_input);
+
+        /* Start  looping on the number of gaussian points: */
+        Gauss* gauss=element->NewGauss(2);
+        for(int ig=gauss->begin();ig<gauss->end();ig++){
+
+                gauss->GaussPoint(ig);
+
+                /* Get Jacobian determinant: */
+                element->JacobianDeterminant(&Jdet,xyz_list,gauss);
+
+                /*Get all parameters at gaussian point*/
+                weights_input->GetInputValue(&weight,gauss,RheologyBInitialguessMisfitEnum);
+                rheologyb_input->GetInputValue(&B,gauss);
+                rheologyb0_input->GetInputValue(&B0,gauss);
+
+                /*Tikhonov regularization: J = 1/2 (B-B0)^2 */ 
+                Jelem+=weight*1./2.*(B-B0)*(B-B0)*Jdet*gauss->weight;
+        }
+
+        /*clean up and Return: */
+        xDelete<IssmDouble>(xyz_list);
+        delete gauss;
+        return Jelem;
+}

issm/trunk-jpl/src/c/modules/RheologyBAbsGradientx/RheologyBAbsGradientx.h ¶

@@ -12 +12 @@
 IssmDouble RheologyBAbsGradient(Element* element);
 
+void RheologyBInitialguessMisfitx( IssmDouble* pJ, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters);
+IssmDouble RheologyBInitialguessMisfit(Element* element);
+
 #endif

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

@@ -601 +601 @@
         RheologyBAbsGradientEnum,
         RheologyBbarAbsGradientEnum,
+        RheologyBInitialguessMisfitEnum,
+        RheologyBInitialguessEnum,
         SealevelEnum,
         SealevelriseCumDeltathicknessEnum,

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

@@ -607 +607 @@
                 case RheologyBAbsGradientEnum : return "RheologyBAbsGradient";
                 case RheologyBbarAbsGradientEnum : return "RheologyBbarAbsGradient";
+                case RheologyBInitialguessMisfitEnum : return "RheologyBInitialguessMisfit";
+                case RheologyBInitialguessEnum : return "RheologyBInitialguess";
                 case SealevelEnum : return "Sealevel";
                 case SealevelriseCumDeltathicknessEnum : return "SealevelriseCumDeltathickness";

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

@@ -619 +619 @@
               else if (strcmp(name,"RheologyBAbsGradient")==0) return RheologyBAbsGradientEnum;
               else if (strcmp(name,"RheologyBbarAbsGradient")==0) return RheologyBbarAbsGradientEnum;
+              else if (strcmp(name,"RheologyBInitialguessMisfit")==0) return RheologyBInitialguessMisfitEnum;
+              else if (strcmp(name,"RheologyBInitialguess")==0) return RheologyBInitialguessEnum;
               else if (strcmp(name,"Sealevel")==0) return SealevelEnum;
               else if (strcmp(name,"SealevelriseCumDeltathickness")==0) return SealevelriseCumDeltathicknessEnum;
@@ -627 +629 @@
               else if (strcmp(name,"SealevelUEsa")==0) return SealevelUEsaEnum;
               else if (strcmp(name,"SealevelRSLEustaticRate")==0) return SealevelRSLEustaticRateEnum;
-              else if (strcmp(name,"SealevelriseSpcthickness")==0) return SealevelriseSpcthicknessEnum;
-              else if (strcmp(name,"SealevelriseStericRate")==0) return SealevelriseStericRateEnum;
          else stage=6;
    }
    if(stage==6){
-              if (strcmp(name,"SealevelNEsa")==0) return SealevelNEsaEnum;
+              if (strcmp(name,"SealevelriseSpcthickness")==0) return SealevelriseSpcthicknessEnum;
+              else if (strcmp(name,"SealevelriseStericRate")==0) return SealevelriseStericRateEnum;
+              else if (strcmp(name,"SealevelNEsa")==0) return SealevelNEsaEnum;
               else if (strcmp(name,"SealevelUGia")==0) return SealevelUGiaEnum;
               else if (strcmp(name,"SealevelNGia")==0) return SealevelNGiaEnum;
@@ -750 +752 @@
               else if (strcmp(name,"StressTensoryy")==0) return StressTensoryyEnum;
               else if (strcmp(name,"StressTensoryz")==0) return StressTensoryzEnum;
-              else if (strcmp(name,"StressTensorzz")==0) return StressTensorzzEnum;
-              else if (strcmp(name,"SurfaceAbsMisfit")==0) return SurfaceAbsMisfitEnum;
          else stage=7;
    }
    if(stage==7){
-              if (strcmp(name,"SurfaceAbsVelMisfit")==0) return SurfaceAbsVelMisfitEnum;
+              if (strcmp(name,"StressTensorzz")==0) return StressTensorzzEnum;
+              else if (strcmp(name,"SurfaceAbsMisfit")==0) return SurfaceAbsMisfitEnum;
+              else if (strcmp(name,"SurfaceAbsVelMisfit")==0) return SurfaceAbsVelMisfitEnum;
               else if (strcmp(name,"SurfaceArea")==0) return SurfaceAreaEnum;
               else if (strcmp(name,"SurfaceAverageVelMisfit")==0) return SurfaceAverageVelMisfitEnum;
@@ -873 +875 @@
               else if (strcmp(name,"Outputdefinition77")==0) return Outputdefinition77Enum;
               else if (strcmp(name,"Outputdefinition78")==0) return Outputdefinition78Enum;
-              else if (strcmp(name,"Outputdefinition79")==0) return Outputdefinition79Enum;
-              else if (strcmp(name,"Outputdefinition7")==0) return Outputdefinition7Enum;
          else stage=8;
    }
    if(stage==8){
-              if (strcmp(name,"Outputdefinition80")==0) return Outputdefinition80Enum;
+              if (strcmp(name,"Outputdefinition79")==0) return Outputdefinition79Enum;
+              else if (strcmp(name,"Outputdefinition7")==0) return Outputdefinition7Enum;
+              else if (strcmp(name,"Outputdefinition80")==0) return Outputdefinition80Enum;
               else if (strcmp(name,"Outputdefinition81")==0) return Outputdefinition81Enum;
               else if (strcmp(name,"Outputdefinition82")==0) return Outputdefinition82Enum;
@@ -996 +998 @@
               else if (strcmp(name,"FrontalForcingsRignot")==0) return FrontalForcingsRignotEnum;
               else if (strcmp(name,"FSApproximation")==0) return FSApproximationEnum;
-              else if (strcmp(name,"Fset")==0) return FsetEnum;
-              else if (strcmp(name,"FSpressure")==0) return FSpressureEnum;
          else stage=9;
    }
    if(stage==9){
-              if (strcmp(name,"FSSolver")==0) return FSSolverEnum;
+              if (strcmp(name,"Fset")==0) return FsetEnum;
+              else if (strcmp(name,"FSpressure")==0) return FSpressureEnum;
+              else if (strcmp(name,"FSSolver")==0) return FSSolverEnum;
               else if (strcmp(name,"FSvelocity")==0) return FSvelocityEnum;
               else if (strcmp(name,"FullMeltOnPartiallyFloating")==0) return FullMeltOnPartiallyFloatingEnum;
@@ -1119 +1121 @@
               else if (strcmp(name,"NoMeltOnPartiallyFloating")==0) return NoMeltOnPartiallyFloatingEnum;
               else if (strcmp(name,"NoneApproximation")==0) return NoneApproximationEnum;
-              else if (strcmp(name,"None")==0) return NoneEnum;
-              else if (strcmp(name,"Numberedcostfunction")==0) return NumberedcostfunctionEnum;
          else stage=10;
    }
    if(stage==10){
-              if (strcmp(name,"Numericalflux")==0) return NumericalfluxEnum;
+              if (strcmp(name,"None")==0) return NoneEnum;
+              else if (strcmp(name,"Numberedcostfunction")==0) return NumberedcostfunctionEnum;
+              else if (strcmp(name,"Numericalflux")==0) return NumericalfluxEnum;
               else if (strcmp(name,"OldGradient")==0) return OldGradientEnum;
               else if (strcmp(name,"OneLayerP4z")==0) return OneLayerP4zEnum;
@@ -1242 +1244 @@
               else if (strcmp(name,"XY")==0) return XYEnum;
               else if (strcmp(name,"XYZ")==0) return XYZEnum;
-              else if (strcmp(name,"BalancethicknessD0")==0) return BalancethicknessD0Enum;
-              else if (strcmp(name,"BalancethicknessDiffusionCoefficient")==0) return BalancethicknessDiffusionCoefficientEnum;
          else stage=11;
    }
    if(stage==11){
-              if (strcmp(name,"BilinearInterp")==0) return BilinearInterpEnum;
+              if (strcmp(name,"BalancethicknessD0")==0) return BalancethicknessD0Enum;
+              else if (strcmp(name,"BalancethicknessDiffusionCoefficient")==0) return BalancethicknessDiffusionCoefficientEnum;
+              else if (strcmp(name,"BilinearInterp")==0) return BilinearInterpEnum;
               else if (strcmp(name,"CalvingdevCoeff")==0) return CalvingdevCoeffEnum;
               else if (strcmp(name,"DeviatoricStress")==0) return DeviatoricStressEnum;