Changeset 18061


Ignore:
Timestamp:
05/26/14 19:42:33 (11 years ago)
Author:
Mathieu Morlighem
Message:

NEW: readded adjoint for Balancethickness2

Location:
issm/trunk-jpl/src/c/analyses
Files:
3 edited

Legend:

Unmodified
Added
Removed

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

    r18057 r18061  
    124124}/*}}}*/
    125125void AdjointBalancethickness2Analysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index){/*{{{*/
    126         _error_("Not implemented yet");
     126        /*The gradient of the cost function is calculated in 2 parts.
     127         *
     128         * dJ    \partial J   \partial lambda^T(KU-F)
     129         * --  = ---------- + ------------------------
     130         * dk    \partial k   \parial k                 
     131         *
     132         * */
     133
     134        /*If on water, grad = 0: */
     135        if(!element->IsIceInElement()) return;
     136
     137        /*Get list of cost functions*/
     138        int *responses = NULL;
     139        int num_responses,resp;
     140        element->FindParam(&num_responses,InversionNumCostFunctionsEnum);
     141        element->FindParam(&responses,NULL,InversionCostFunctionsEnum);
     142
     143        /*Check that control_type is supported*/
     144        if(control_type!=BalancethicknessApparentMassbalanceEnum){
     145                _error_("Control "<<EnumToStringx(control_type)<<" not supported");
     146        }
     147
     148        /*Deal with first part (partial derivative a J with respect to k)*/
     149        for(resp=0;resp<num_responses;resp++) switch(responses[resp]){
     150                case Balancethickness2MisfitEnum: /*Nothing, \partial J/\partial k = 0*/ break;
     151                default: _error_("response " << EnumToStringx(responses[resp]) << " not supported yet");
     152        }
     153
     154        /*Deal with second term*/
     155        switch(control_type){
     156                case BalancethicknessApparentMassbalanceEnum: GradientJAdot(element,gradient,control_index); break;
     157                default: _error_("control type not supported yet: " << EnumToStringx(control_type));
     158        }
     159
     160        /*Clean up and return*/
     161        xDelete<int>(responses);
     162
     163}/*}}}*/
     164void AdjointBalancethickness2Analysis::GradientJAdot(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
     165
     166        /*Intermediaries*/
     167        IssmDouble Jdet,weight;
     168        IssmDouble lambda;
     169        IssmDouble *xyz_list= NULL;
     170
     171        /*Fetch number of vertices for this finite element*/
     172        int numvertices = element->GetNumberOfVertices();
     173
     174        /*Initialize some vectors*/
     175        IssmDouble* basis         = xNew<IssmDouble>(numvertices);
     176        IssmDouble* ge            = xNewZeroInit<IssmDouble>(numvertices);
     177        int*        vertexpidlist = xNew<int>(numvertices);
     178
     179        /*Retrieve all inputs we will be needing: */
     180        element->GetVerticesCoordinates(&xyz_list);
     181        element->GradientIndexing(&vertexpidlist[0],control_index);
     182        Input* adjoint_input = element->GetInput(AdjointEnum);                            _assert_(adjoint_input);
     183        Input* weights_input = element->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
     184
     185        Gauss* gauss=element->NewGauss(2);
     186        for(int ig=gauss->begin();ig<gauss->end();ig++){
     187                gauss->GaussPoint(ig);
     188
     189                element->JacobianDeterminant(&Jdet,xyz_list,gauss);
     190                element->NodalFunctionsP1(basis,gauss);
     191                weights_input->GetInputValue(&weight,gauss,Balancethickness2MisfitEnum);
     192                adjoint_input->GetInputValue(&lambda,gauss);
     193
     194                /*Build gradient vector (actually -dJ/da): */
     195                for(int i=0;i<numvertices;i++){
     196                        ge[i]+= - weight*Jdet*gauss->weight*basis[i]*lambda;
     197                        _assert_(!xIsNan<IssmDouble>(ge[i]));
     198                }
     199        }
     200        gradient->SetValues(numvertices,vertexpidlist,ge,ADD_VAL);
     201
     202        /*Clean up and return*/
     203        xDelete<IssmDouble>(ge);
     204        xDelete<IssmDouble>(xyz_list);
     205        xDelete<IssmDouble>(basis);
     206        xDelete<int>(vertexpidlist);
     207        delete gauss;
    127208}/*}}}*/
    128209void AdjointBalancethickness2Analysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/

  • issm/trunk-jpl/src/c/analyses/AdjointBalancethickness2Analysis.h

    r18057 r18061  
    2828                void GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element);
    2929                void GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index);
     30                void GradientJAdot(Element* element,Vector<IssmDouble>* gradient,int control_index);
    3031                void InputUpdateFromSolution(IssmDouble* solution,Element* element);
    3132                void UpdateConstraints(FemModel* femmodel);

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

    r18060 r18061  
    10351035                dragcoefficient_input->GetInputDerivativeValue(&dk[0],xyz_list,gauss);
    10361036
     1037                /*Build gradient vector (actually -dJ/ddrag): */
     1038                for(int i=0;i<numvertices;i++){
     1039                        if(dim==2){
     1040                                ge[i]+=-weight*Jdet*gauss->weight*(dbasis[0*numvertices+i]*dk[0]+dbasis[1*numvertices+i]*dk[1]);
     1041                        }
     1042                        else{
     1043                                ge[i]+=-weight*Jdet*gauss->weight*dbasis[0*numvertices+i]*dk[0];
     1044                        }
     1045                        _assert_(!xIsNan<IssmDouble>(ge[i]));
     1046                }
     1047        }
     1048        gradient->SetValues(numvertices,vertexpidlist,ge,ADD_VAL);
     1049
     1050        /*Clean up and return*/
     1051        xDelete<IssmDouble>(xyz_list);
     1052        xDelete<IssmDouble>(dbasis);
     1053        xDelete<IssmDouble>(ge);
     1054        xDelete<int>(vertexpidlist);
     1055        delete gauss;
     1056        if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
     1057
     1058}/*}}}*/
     1059void AdjointHorizAnalysis::GradientJBGradient(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
     1060
     1061        /*Intermediaries*/
     1062        int      domaintype,dim;
     1063        Element* basalelement;
     1064
     1065        /*Get basal element*/
     1066        element->FindParam(&domaintype,DomainTypeEnum);
     1067        switch(domaintype){
     1068                case Domain2DhorizontalEnum:
     1069                        basalelement = element;
     1070                        dim          = 2;
     1071                        break;
     1072                case Domain2DverticalEnum:
     1073                        if(!element->IsOnBase()) return;
     1074                        basalelement = element->SpawnBasalElement();
     1075                        dim          = 1;
     1076                        break;
     1077                case Domain3DEnum:
     1078                        if(!element->IsOnBase()) return;
     1079                        basalelement = element->SpawnBasalElement();
     1080                        dim          = 2;
     1081                        break;
     1082                default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
     1083        }
     1084
     1085        /*Intermediaries*/
     1086        IssmDouble Jdet,weight;
     1087        IssmDouble dk[3];
     1088        IssmDouble *xyz_list= NULL;
     1089
     1090        /*Fetch number of vertices for this finite element*/
     1091        int numvertices = basalelement->GetNumberOfVertices();
     1092
     1093        /*Initialize some vectors*/
     1094        IssmDouble* dbasis        = xNew<IssmDouble>(2*numvertices);
     1095        IssmDouble* ge            = xNewZeroInit<IssmDouble>(numvertices);
     1096        int*        vertexpidlist = xNew<int>(numvertices);
     1097
     1098        /*Retrieve all inputs we will be needing: */
     1099        basalelement->GetVerticesCoordinates(&xyz_list);
     1100        basalelement->GradientIndexing(&vertexpidlist[0],control_index);
     1101        Input* rheologyb_input = basalelement->GetInput(MaterialsRheologyBbarEnum);              _assert_(rheologyb_input);
     1102        Input* weights_input   = basalelement->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
     1103
     1104        /* Start  looping on the number of gaussian points: */
     1105        Gauss* gauss=basalelement->NewGauss(2);
     1106        for(int ig=gauss->begin();ig<gauss->end();ig++){
     1107                gauss->GaussPoint(ig);
     1108
     1109                basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
     1110                basalelement->NodalFunctionsP1Derivatives(dbasis,xyz_list,gauss);
     1111                weights_input->GetInputValue(&weight,gauss,RheologyBbarAbsGradientEnum);
     1112
     1113                /*Build alpha_complement_list: */
     1114                rheologyb_input->GetInputDerivativeValue(&dk[0],xyz_list,gauss);
     1115
    10371116                /*Build gradje_g_gaussian vector (actually -dJ/ddrag): */
    10381117                for(int i=0;i<numvertices;i++){
    … …  
    10551134        delete gauss;
    10561135        if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
    1057 
    1058 }/*}}}*/
    1059 void AdjointHorizAnalysis::GradientJBGradient(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
    1060 
    1061         /*Intermediaries*/
    1062         int      domaintype,dim;
    1063         Element* basalelement;
    1064 
    1065         /*Get basal element*/
    1066         element->FindParam(&domaintype,DomainTypeEnum);
    1067         switch(domaintype){
    1068                 case Domain2DhorizontalEnum:
    1069                         basalelement = element;
    1070                         dim          = 2;
    1071                         break;
    1072                 case Domain2DverticalEnum:
    1073                         if(!element->IsOnBase()) return;
    1074                         basalelement = element->SpawnBasalElement();
    1075                         dim          = 1;
    1076                         break;
    1077                 case Domain3DEnum:
    1078                         if(!element->IsOnBase()) return;
    1079                         basalelement = element->SpawnBasalElement();
    1080                         dim          = 2;
    1081                         break;
    1082                 default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
    1083         }
    1084 
    1085         /*Intermediaries*/
    1086         IssmDouble Jdet,weight;
    1087         IssmDouble dk[3];
    1088         IssmDouble *xyz_list= NULL;
    1089 
    1090         /*Fetch number of vertices for this finite element*/
    1091         int numvertices = basalelement->GetNumberOfVertices();
    1092 
    1093         /*Initialize some vectors*/
    1094         IssmDouble* dbasis        = xNew<IssmDouble>(2*numvertices);
    1095         IssmDouble* ge            = xNewZeroInit<IssmDouble>(numvertices);
    1096         int*        vertexpidlist = xNew<int>(numvertices);
    1097 
    1098         /*Retrieve all inputs we will be needing: */
    1099         basalelement->GetVerticesCoordinates(&xyz_list);
    1100         basalelement->GradientIndexing(&vertexpidlist[0],control_index);
    1101         Input* rheologyb_input = basalelement->GetInput(MaterialsRheologyBbarEnum);              _assert_(rheologyb_input);
    1102         Input* weights_input   = basalelement->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
    1103 
    1104         /* Start  looping on the number of gaussian points: */
    1105         Gauss* gauss=basalelement->NewGauss(2);
    1106         for(int ig=gauss->begin();ig<gauss->end();ig++){
    1107                 gauss->GaussPoint(ig);
    1108 
    1109                 basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
    1110                 basalelement->NodalFunctionsP1Derivatives(dbasis,xyz_list,gauss);
    1111                 weights_input->GetInputValue(&weight,gauss,RheologyBbarAbsGradientEnum);
    1112 
    1113                 /*Build alpha_complement_list: */
    1114                 rheologyb_input->GetInputDerivativeValue(&dk[0],xyz_list,gauss);
    1115 
    1116                 /*Build gradje_g_gaussian vector (actually -dJ/ddrag): */
    1117                 for(int i=0;i<numvertices;i++){
    1118                         if(dim==2){
    1119                                 ge[i]+=-weight*Jdet*gauss->weight*(dbasis[0*numvertices+i]*dk[0]+dbasis[1*numvertices+i]*dk[1]);
    1120                         }
    1121                         else{
    1122                                 ge[i]+=-weight*Jdet*gauss->weight*dbasis[0*numvertices+i]*dk[0];
    1123                         }
    1124                         _assert_(!xIsNan<IssmDouble>(ge[i]));
    1125                 }
    1126         }
    1127         gradient->SetValues(numvertices,vertexpidlist,ge,ADD_VAL);
    1128 
    1129         /*Clean up and return*/
    1130         xDelete<IssmDouble>(xyz_list);
    1131         xDelete<IssmDouble>(dbasis);
    1132         xDelete<IssmDouble>(ge);
    1133         xDelete<int>(vertexpidlist);
    1134         delete gauss;
    1135         if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
    11361136}/*}}}*/
    11371137void AdjointHorizAnalysis::GradientJDragSSA(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
Note: See TracChangeset for help on using the changeset viewer.