Changeset 18060 for issm/trunk-jpl/src/c/analyses/AdjointBalancethicknessAnalysis.cpp

Timestamp:

05/26/14 11:22:52 (11 years ago)

Author:

Mathieu Morlighem

Message:

CHG: moving gradients from element to adjoint analyses

File:

• issm/trunk-jpl/src/c/analyses/AdjointBalancethicknessAnalysis.cpp (modified) (1 diff)

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

@@ -151 +151 @@
 }/*}}}*/
 void AdjointBalancethicknessAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index){/*{{{*/
-        _error_("Not implemented yet");
+        /*The gradient of the cost function is calculated in 2 parts.
+         *
+         * dJ    \partial J   \partial lambda^T(KU-F)
+         * --  = ---------- + ------------------------
+         * dk    \partial k   \parial k                  
+         *
+         * */
+
+        /*If on water, grad = 0: */
+        if(!element->IsIceInElement()) return;
+
+        /*Get list of cost functions*/
+        int *responses = NULL;
+        int num_responses,resp;
+        element->FindParam(&num_responses,InversionNumCostFunctionsEnum);
+        element->FindParam(&responses,NULL,InversionCostFunctionsEnum);
+
+        /*Check that control_type is supported*/
+        if(control_type!=VxEnum && 
+                control_type!=VyEnum && 
+                control_type!=BalancethicknessThickeningRateEnum){
+                _error_("Control "<<EnumToStringx(control_type)<<" not supported");
+        }
+
+        /*Deal with first part (partial derivative a J with respect to k)*/
+        for(resp=0;resp<num_responses;resp++) switch(responses[resp]){
+                case ThicknessAbsMisfitEnum:      /*Nothing, \partial J/\partial k = 0*/ break;
+                case ThicknessAbsGradientEnum:    /*Nothing, \partial J/\partial k = 0*/ break;
+                case ThicknessAlongGradientEnum:  /*Nothing, \partial J/\partial k = 0*/ break;
+                case ThicknessAcrossGradientEnum: /*Nothing, \partial J/\partial k = 0*/ break;
+                default: _error_("response " << EnumToStringx(responses[resp]) << " not supported yet");
+        }
+
+        /*Deal with second term*/
+        switch(control_type){
+                case BalancethicknessThickeningRateEnum: GradientJDhDt(element,gradient,control_index); break;
+                case VxEnum:                             GradientJVx(  element,gradient,control_index); break;
+                case VyEnum:                             GradientJVy(  element,gradient,control_index); break;
+                default: _error_("control type not supported yet: " << EnumToStringx(control_type));
+        }
+
+        /*Clean up and return*/
+        xDelete<int>(responses);
+
+}/*}}}*/
+void AdjointBalancethicknessAnalysis::GradientJVx(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
+
+        /*Intermediaries*/
+        IssmDouble Jdet,weight;
+        IssmDouble thickness,Dlambda[3],dp[3];
+        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* thickness_input = element->GetInput(ThicknessEnum); _assert_(thickness_input);
+        Input* adjoint_input   = element->GetInput(AdjointEnum);   _assert_(adjoint_input);
+
+        /* Start  looping on the number of gaussian points: */
+        Gauss* gauss=element->NewGauss(4);
+        for(int ig=gauss->begin();ig<gauss->end();ig++){
+                gauss->GaussPoint(ig);
+
+                adjoint_input->GetInputDerivativeValue(&Dlambda[0],xyz_list,gauss);
+                thickness_input->GetInputValue(&thickness, gauss);
+                thickness_input->GetInputDerivativeValue(&dp[0],xyz_list,gauss);
+
+                element->JacobianDeterminant(&Jdet,xyz_list,gauss);
+                element->NodalFunctionsP1(basis,gauss);
+
+                /*Build gradient vector (actually -dJ/dD): */
+                for(int i=0;i<numvertices;i++){
+                        ge[i]+=thickness*Dlambda[0]*Jdet*gauss->weight*basis[i];
+                        _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 AdjointBalancethicknessAnalysis::GradientJVy(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
+
+        /*Intermediaries*/
+        IssmDouble Jdet,weight;
+        IssmDouble thickness,Dlambda[3],dp[3];
+        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* thickness_input = element->GetInput(ThicknessEnum); _assert_(thickness_input);
+        Input* adjoint_input   = element->GetInput(AdjointEnum);   _assert_(adjoint_input);
+
+        /* Start  looping on the number of gaussian points: */
+        Gauss* gauss=element->NewGauss(4);
+        for(int ig=gauss->begin();ig<gauss->end();ig++){
+                gauss->GaussPoint(ig);
+
+                adjoint_input->GetInputDerivativeValue(&Dlambda[0],xyz_list,gauss);
+                thickness_input->GetInputValue(&thickness, gauss);
+                thickness_input->GetInputDerivativeValue(&dp[0],xyz_list,gauss);
+
+                element->JacobianDeterminant(&Jdet,xyz_list,gauss);
+                element->NodalFunctionsP1(basis,gauss);
+
+                /*Build gradient vector (actually -dJ/dvy): */
+                for(int i=0;i<numvertices;i++){
+                        ge[i]+=thickness*Dlambda[1]*Jdet*gauss->weight*basis[i];
+                        _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 AdjointBalancethicknessAnalysis::GradientJDhDt(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
+
+        /*Fetch number of vertices for this finite element*/
+        int numvertices = element->GetNumberOfVertices();
+
+        /*Initialize some vectors*/
+        IssmDouble* ge            = xNewZeroInit<IssmDouble>(numvertices);
+        IssmDouble* lambda        = xNew<IssmDouble>(numvertices);
+        int*        vertexpidlist = xNew<int>(numvertices);
+
+        /*Retrieve all inputs we will be needing: */
+        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);
+
+        /*Clean up and return*/
+        xDelete<IssmDouble>(ge);
+        xDelete<IssmDouble>(lambda);
+        xDelete<int>(vertexpidlist);
 }/*}}}*/
 void AdjointBalancethicknessAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/