Changeset 18060


Ignore:
Timestamp:
05/26/14 11:22:52 (11 years ago)
Author:
Mathieu Morlighem
Message:

CHG: moving gradients from element to adjoint analyses

Location:
issm/trunk-jpl/src/c
Files:
15 edited

Legend:

Unmodified
Added
Removed

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

    r18057 r18060  
    151151}/*}}}*/
    152152void AdjointBalancethicknessAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index){/*{{{*/
    153         _error_("Not implemented yet");
     153        /*The gradient of the cost function is calculated in 2 parts.
     154         *
     155         * dJ    \partial J   \partial lambda^T(KU-F)
     156         * --  = ---------- + ------------------------
     157         * dk    \partial k   \parial k                 
     158         *
     159         * */
     160
     161        /*If on water, grad = 0: */
     162        if(!element->IsIceInElement()) return;
     163
     164        /*Get list of cost functions*/
     165        int *responses = NULL;
     166        int num_responses,resp;
     167        element->FindParam(&num_responses,InversionNumCostFunctionsEnum);
     168        element->FindParam(&responses,NULL,InversionCostFunctionsEnum);
     169
     170        /*Check that control_type is supported*/
     171        if(control_type!=VxEnum &&
     172                control_type!=VyEnum &&
     173                control_type!=BalancethicknessThickeningRateEnum){
     174                _error_("Control "<<EnumToStringx(control_type)<<" not supported");
     175        }
     176
     177        /*Deal with first part (partial derivative a J with respect to k)*/
     178        for(resp=0;resp<num_responses;resp++) switch(responses[resp]){
     179                case ThicknessAbsMisfitEnum:      /*Nothing, \partial J/\partial k = 0*/ break;
     180                case ThicknessAbsGradientEnum:    /*Nothing, \partial J/\partial k = 0*/ break;
     181                case ThicknessAlongGradientEnum:  /*Nothing, \partial J/\partial k = 0*/ break;
     182                case ThicknessAcrossGradientEnum: /*Nothing, \partial J/\partial k = 0*/ break;
     183                default: _error_("response " << EnumToStringx(responses[resp]) << " not supported yet");
     184        }
     185
     186        /*Deal with second term*/
     187        switch(control_type){
     188                case BalancethicknessThickeningRateEnum: GradientJDhDt(element,gradient,control_index); break;
     189                case VxEnum:                             GradientJVx(  element,gradient,control_index); break;
     190                case VyEnum:                             GradientJVy(  element,gradient,control_index); break;
     191                default: _error_("control type not supported yet: " << EnumToStringx(control_type));
     192        }
     193
     194        /*Clean up and return*/
     195        xDelete<int>(responses);
     196
     197}/*}}}*/
     198void AdjointBalancethicknessAnalysis::GradientJVx(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
     199
     200        /*Intermediaries*/
     201        IssmDouble Jdet,weight;
     202        IssmDouble thickness,Dlambda[3],dp[3];
     203        IssmDouble *xyz_list= NULL;
     204
     205        /*Fetch number of vertices for this finite element*/
     206        int numvertices = element->GetNumberOfVertices();
     207
     208        /*Initialize some vectors*/
     209        IssmDouble* basis         = xNew<IssmDouble>(numvertices);
     210        IssmDouble* ge            = xNewZeroInit<IssmDouble>(numvertices);
     211        int*        vertexpidlist = xNew<int>(numvertices);
     212
     213        /*Retrieve all inputs we will be needing: */
     214        element->GetVerticesCoordinates(&xyz_list);
     215        element->GradientIndexing(&vertexpidlist[0],control_index);
     216        Input* thickness_input = element->GetInput(ThicknessEnum); _assert_(thickness_input);
     217        Input* adjoint_input   = element->GetInput(AdjointEnum);   _assert_(adjoint_input);
     218
     219        /* Start  looping on the number of gaussian points: */
     220        Gauss* gauss=element->NewGauss(4);
     221        for(int ig=gauss->begin();ig<gauss->end();ig++){
     222                gauss->GaussPoint(ig);
     223
     224                adjoint_input->GetInputDerivativeValue(&Dlambda[0],xyz_list,gauss);
     225                thickness_input->GetInputValue(&thickness, gauss);
     226                thickness_input->GetInputDerivativeValue(&dp[0],xyz_list,gauss);
     227
     228                element->JacobianDeterminant(&Jdet,xyz_list,gauss);
     229                element->NodalFunctionsP1(basis,gauss);
     230
     231                /*Build gradient vector (actually -dJ/dD): */
     232                for(int i=0;i<numvertices;i++){
     233                        ge[i]+=thickness*Dlambda[0]*Jdet*gauss->weight*basis[i];
     234                        _assert_(!xIsNan<IssmDouble>(ge[i]));
     235                }
     236        }
     237        gradient->SetValues(numvertices,vertexpidlist,ge,ADD_VAL);
     238
     239        /*Clean up and return*/
     240        xDelete<IssmDouble>(xyz_list);
     241        xDelete<IssmDouble>(basis);
     242        xDelete<IssmDouble>(ge);
     243        xDelete<int>(vertexpidlist);
     244        delete gauss;
     245}/*}}}*/
     246void AdjointBalancethicknessAnalysis::GradientJVy(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
     247
     248        /*Intermediaries*/
     249        IssmDouble Jdet,weight;
     250        IssmDouble thickness,Dlambda[3],dp[3];
     251        IssmDouble *xyz_list= NULL;
     252
     253        /*Fetch number of vertices for this finite element*/
     254        int numvertices = element->GetNumberOfVertices();
     255
     256        /*Initialize some vectors*/
     257        IssmDouble* basis         = xNew<IssmDouble>(numvertices);
     258        IssmDouble* ge            = xNewZeroInit<IssmDouble>(numvertices);
     259        int*        vertexpidlist = xNew<int>(numvertices);
     260
     261        /*Retrieve all inputs we will be needing: */
     262        element->GetVerticesCoordinates(&xyz_list);
     263        element->GradientIndexing(&vertexpidlist[0],control_index);
     264        Input* thickness_input = element->GetInput(ThicknessEnum); _assert_(thickness_input);
     265        Input* adjoint_input   = element->GetInput(AdjointEnum);   _assert_(adjoint_input);
     266
     267        /* Start  looping on the number of gaussian points: */
     268        Gauss* gauss=element->NewGauss(4);
     269        for(int ig=gauss->begin();ig<gauss->end();ig++){
     270                gauss->GaussPoint(ig);
     271
     272                adjoint_input->GetInputDerivativeValue(&Dlambda[0],xyz_list,gauss);
     273                thickness_input->GetInputValue(&thickness, gauss);
     274                thickness_input->GetInputDerivativeValue(&dp[0],xyz_list,gauss);
     275
     276                element->JacobianDeterminant(&Jdet,xyz_list,gauss);
     277                element->NodalFunctionsP1(basis,gauss);
     278
     279                /*Build gradient vector (actually -dJ/dvy): */
     280                for(int i=0;i<numvertices;i++){
     281                        ge[i]+=thickness*Dlambda[1]*Jdet*gauss->weight*basis[i];
     282                        _assert_(!xIsNan<IssmDouble>(ge[i]));
     283                }
     284        }
     285        gradient->SetValues(numvertices,vertexpidlist,ge,ADD_VAL);
     286
     287        /*Clean up and return*/
     288        xDelete<IssmDouble>(xyz_list);
     289        xDelete<IssmDouble>(basis);
     290        xDelete<IssmDouble>(ge);
     291        xDelete<int>(vertexpidlist);
     292        delete gauss;
     293}/*}}}*/
     294void AdjointBalancethicknessAnalysis::GradientJDhDt(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
     295
     296        /*Fetch number of vertices for this finite element*/
     297        int numvertices = element->GetNumberOfVertices();
     298
     299        /*Initialize some vectors*/
     300        IssmDouble* ge            = xNewZeroInit<IssmDouble>(numvertices);
     301        IssmDouble* lambda        = xNew<IssmDouble>(numvertices);
     302        int*        vertexpidlist = xNew<int>(numvertices);
     303
     304        /*Retrieve all inputs we will be needing: */
     305        element->GradientIndexing(&vertexpidlist[0],control_index);
     306        element->GetInputListOnVertices(lambda,AdjointEnum);
     307        for(int i=0;i<numvertices;i++){
     308                ge[i]=-lambda[i];
     309                _assert_(!xIsNan<IssmDouble>(ge[i]));
     310        }
     311        gradient->SetValues(numvertices,vertexpidlist,ge,INS_VAL);
     312
     313        /*Clean up and return*/
     314        xDelete<IssmDouble>(ge);
     315        xDelete<IssmDouble>(lambda);
     316        xDelete<int>(vertexpidlist);
    154317}/*}}}*/
    155318void AdjointBalancethicknessAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/

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

    r18057 r18060  
    2828                void GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element);
    2929                void GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index);
     30                void GradientJVx(Element* element,Vector<IssmDouble>* gradient,int control_index);
     31                void GradientJVy(Element* element,Vector<IssmDouble>* gradient,int control_index);
     32                void GradientJDhDt(Element* element,Vector<IssmDouble>* gradient,int control_index);
    3033                void InputUpdateFromSolution(IssmDouble* solution,Element* element);
    3134                void UpdateConstraints(FemModel* femmodel);

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

    r18059 r18060  
    923923
    924924        /*Check that control_type is supported*/
    925         if(control_type!=MaterialsRheologyBbarEnum || control_type!=FrictionCoefficientEnum || control_type!=DamageDbarEnum){
     925        if(control_type!=MaterialsRheologyBbarEnum &&
     926                control_type!=FrictionCoefficientEnum   &&
     927                control_type!=DamageDbarEnum){
    926928                _error_("Control "<<EnumToStringx(control_type)<<" not supported");
    927929        }
    … …  
    975977void AdjointHorizAnalysis::GradientJDragGradient(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
    976978
    977         /*return if floating*/
     979        /*return if floating (gradient is 0)*/
    978980        if(element->IsFloating()) return;
    979981
    … …  
    10121014        /*Initialize some vectors*/
    10131015        IssmDouble* dbasis        = xNew<IssmDouble>(2*numvertices);
    1014         IssmDouble* ge            = xNew<IssmDouble>(numvertices);
     1016        IssmDouble* ge            = xNewZeroInit<IssmDouble>(numvertices);
    10151017        int*        vertexpidlist = xNew<int>(numvertices);
    10161018
    … …  
    10561058}/*}}}*/
    10571059void AdjointHorizAnalysis::GradientJBGradient(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
    1058         _error_("not implemented yet");
    1059 }/*}}}*/
    1060 void AdjointHorizAnalysis::GradientJDragSSA(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
    1061         _error_("not implemented yet");
    1062 }/*}}}*/
    1063 void AdjointHorizAnalysis::GradientJDragHO(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
    1064         _error_("not implemented yet");
    1065 }/*}}}*/
    1066 void AdjointHorizAnalysis::GradientJDragFS(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
    1067         _error_("not implemented yet");
    1068 }/*}}}*/
    1069 void AdjointHorizAnalysis::GradientJBbarSSA(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
    10701060
    10711061        /*Intermediaries*/
    … …  
    10951085        /*Intermediaries*/
    10961086        IssmDouble Jdet,weight;
    1097         IssmDouble thickness,dmudB;
    1098         IssmDouble dvx[3],dvy[3],dadjx[3],dadjy[3];
     1087        IssmDouble dk[3];
    10991088        IssmDouble *xyz_list= NULL;
    11001089
    … …  
    11031092
    11041093        /*Initialize some vectors*/
    1105         IssmDouble* basis         = xNew<IssmDouble>(numvertices);
    1106         IssmDouble* ge            = xNew<IssmDouble>(numvertices);
     1094        IssmDouble* dbasis        = xNew<IssmDouble>(2*numvertices);
     1095        IssmDouble* ge            = xNewZeroInit<IssmDouble>(numvertices);
    11071096        int*        vertexpidlist = xNew<int>(numvertices);
    11081097
    … …  
    11101099        basalelement->GetVerticesCoordinates(&xyz_list);
    11111100        basalelement->GradientIndexing(&vertexpidlist[0],control_index);
    1112         Input* thickness_input = element->GetInput(ThicknessEnum);             _assert_(thickness_input);
    1113         Input* vx_input        = element->GetInput(VxEnum);                    _assert_(vx_input);
    1114         Input* vy_input        = element->GetInput(VyEnum);                    _assert_(vy_input);
    1115         Input* adjointx_input  = element->GetInput(AdjointxEnum);              _assert_(adjointx_input);
    1116         Input* adjointy_input  = element->GetInput(AdjointyEnum);              _assert_(adjointy_input);
    1117         Input* rheologyb_input = element->GetInput(MaterialsRheologyBbarEnum); _assert_(rheologyb_input);
     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;};
     1136}/*}}}*/
     1137void AdjointHorizAnalysis::GradientJDragSSA(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
     1138
     1139        /*return if floating (gradient is 0)*/
     1140        if(element->IsFloating()) return;
     1141
     1142        /*Intermediaries*/
     1143        int      domaintype,dim;
     1144        Element* basalelement;
     1145
     1146        /*Get basal element*/
     1147        element->FindParam(&domaintype,DomainTypeEnum);
     1148        switch(domaintype){
     1149                case Domain2DhorizontalEnum:
     1150                        basalelement = element;
     1151                        dim          = 2;
     1152                        break;
     1153                case Domain2DverticalEnum:
     1154                        if(!element->IsOnBase()) return;
     1155                        basalelement = element->SpawnBasalElement();
     1156                        dim          = 1;
     1157                        break;
     1158                case Domain3DEnum:
     1159                        if(!element->IsOnBase()) return;
     1160                        basalelement = element->SpawnBasalElement();
     1161                        dim          = 2;
     1162                        break;
     1163                default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
     1164        }
     1165
     1166        /*Intermediaries*/
     1167        IssmDouble Jdet,weight;
     1168        IssmDouble drag,dalpha2dk;
     1169        IssmDouble vx,vy,lambda,mu;
     1170        IssmDouble *xyz_list= NULL;
     1171
     1172        /*Fetch number of vertices for this finite element*/
     1173        int numvertices = basalelement->GetNumberOfVertices();
     1174
     1175        /*Initialize some vectors*/
     1176        IssmDouble* basis         = xNew<IssmDouble>(numvertices);
     1177        IssmDouble* ge            = xNewZeroInit<IssmDouble>(numvertices);
     1178        int*        vertexpidlist = xNew<int>(numvertices);
     1179
     1180        /*Build friction element, needed later: */
     1181        Friction* friction=new Friction(basalelement,dim);
     1182
     1183        /*Retrieve all inputs we will be needing: */
     1184        basalelement->GetVerticesCoordinates(&xyz_list);
     1185        basalelement->GradientIndexing(&vertexpidlist[0],control_index);
     1186        Input* vx_input        = basalelement->GetInput(VxEnum);                   _assert_(vx_input);
     1187        Input* vy_input        = basalelement->GetInput(VyEnum);                   _assert_(vy_input);
     1188        Input* adjointx_input  = basalelement->GetInput(AdjointxEnum);             _assert_(adjointx_input);
     1189        Input* adjointy_input  = basalelement->GetInput(AdjointyEnum);             _assert_(adjointy_input);
     1190        Input* dragcoeff_input = basalelement->GetInput(FrictionCoefficientEnum);  _assert_(dragcoeff_input);
    11181191
    11191192        /* Start  looping on the number of gaussian points: */
    … …  
    11221195                gauss->GaussPoint(ig);
    11231196
     1197                adjointx_input->GetInputValue(&lambda, gauss);
     1198                adjointy_input->GetInputValue(&mu, gauss);
     1199                vx_input->GetInputValue(&vx,gauss);
     1200                vy_input->GetInputValue(&vy,gauss);
     1201                dragcoeff_input->GetInputValue(&drag, gauss);
     1202
     1203                friction->GetAlphaComplement(&dalpha2dk,gauss);
     1204
     1205                basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
     1206                basalelement->NodalFunctionsP1(basis,gauss);
     1207
     1208                /*Build gradient vector (actually -dJ/dD): */
     1209                for(int i=0;i<numvertices;i++){
     1210                        ge[i]+=-2.*drag*dalpha2dk*((lambda*vx+mu*vy))*Jdet*gauss->weight*basis[i];
     1211                        _assert_(!xIsNan<IssmDouble>(ge[i]));
     1212                }
     1213        }
     1214        gradient->SetValues(numvertices,vertexpidlist,ge,ADD_VAL);
     1215
     1216        /*Clean up and return*/
     1217        xDelete<IssmDouble>(xyz_list);
     1218        xDelete<IssmDouble>(basis);
     1219        xDelete<IssmDouble>(ge);
     1220        xDelete<int>(vertexpidlist);
     1221        delete gauss;
     1222        delete friction;
     1223        if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
     1224}/*}}}*/
     1225void AdjointHorizAnalysis::GradientJDragHO(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
     1226
     1227        /*return if floating or not on bed (gradient is 0)*/
     1228        if(element->IsFloating()) return;
     1229        if(!element->IsOnBase()) return;
     1230
     1231        /*Intermediaries*/
     1232        int        dim=3;
     1233        IssmDouble Jdet,weight;
     1234        IssmDouble drag,dalpha2dk;
     1235        IssmDouble vx,vy,lambda,mu;
     1236        IssmDouble *xyz_list_base= NULL;
     1237
     1238        /*Fetch number of vertices for this finite element*/
     1239        int numvertices = element->GetNumberOfVertices();
     1240
     1241        /*Initialize some vectors*/
     1242        IssmDouble* basis         = xNew<IssmDouble>(numvertices);
     1243        IssmDouble* ge            = xNewZeroInit<IssmDouble>(numvertices);
     1244        int*        vertexpidlist = xNew<int>(numvertices);
     1245
     1246        /*Build friction element, needed later: */
     1247        Friction* friction=new Friction(element,dim);
     1248
     1249        /*Retrieve all inputs we will be needing: */
     1250        element->GetVerticesCoordinatesBase(&xyz_list_base);
     1251        element->GradientIndexing(&vertexpidlist[0],control_index);
     1252        Input* vx_input        = element->GetInput(VxEnum);                   _assert_(vx_input);
     1253        Input* vy_input        = element->GetInput(VyEnum);                   _assert_(vy_input);
     1254        Input* adjointx_input  = element->GetInput(AdjointxEnum);             _assert_(adjointx_input);
     1255        Input* adjointy_input  = element->GetInput(AdjointyEnum);             _assert_(adjointy_input);
     1256        Input* dragcoeff_input = element->GetInput(FrictionCoefficientEnum);  _assert_(dragcoeff_input);
     1257
     1258        /* Start  looping on the number of gaussian points: */
     1259        Gauss* gauss=element->NewGaussBase(4);
     1260        for(int ig=gauss->begin();ig<gauss->end();ig++){
     1261                gauss->GaussPoint(ig);
     1262
     1263                adjointx_input->GetInputValue(&lambda, gauss);
     1264                adjointy_input->GetInputValue(&mu, gauss);
     1265                vx_input->GetInputValue(&vx,gauss);
     1266                vy_input->GetInputValue(&vy,gauss);
     1267                dragcoeff_input->GetInputValue(&drag, gauss);
     1268
     1269                friction->GetAlphaComplement(&dalpha2dk,gauss);
     1270
     1271                element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss);
     1272                element->NodalFunctionsP1(basis,gauss);
     1273
     1274                /*Build gradient vector (actually -dJ/dD): */
     1275                for(int i=0;i<numvertices;i++){
     1276                        ge[i]+=-2.*drag*dalpha2dk*((lambda*vx+mu*vy))*Jdet*gauss->weight*basis[i];
     1277                        _assert_(!xIsNan<IssmDouble>(ge[i]));
     1278                }
     1279        }
     1280        gradient->SetValues(numvertices,vertexpidlist,ge,ADD_VAL);
     1281
     1282        /*Clean up and return*/
     1283        xDelete<IssmDouble>(xyz_list_base);
     1284        xDelete<IssmDouble>(basis);
     1285        xDelete<IssmDouble>(ge);
     1286        xDelete<int>(vertexpidlist);
     1287        delete gauss;
     1288        delete friction;
     1289}/*}}}*/
     1290void AdjointHorizAnalysis::GradientJDragFS(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
     1291
     1292        /*return if floating or not on bed (gradient is 0)*/
     1293        if(element->IsFloating()) return;
     1294        if(!element->IsOnBase()) return;
     1295
     1296        /*Intermediaries*/
     1297        int        dim=3;
     1298        IssmDouble Jdet,weight;
     1299        IssmDouble drag,dalpha2dk,normal[3];
     1300        IssmDouble vx,vy,vz,lambda,mu,xi;
     1301        IssmDouble *xyz_list_base= NULL;
     1302
     1303        /*Fetch number of vertices for this finite element*/
     1304        int numvertices = element->GetNumberOfVertices();
     1305
     1306        /*Initialize some vectors*/
     1307        IssmDouble* basis         = xNew<IssmDouble>(numvertices);
     1308        IssmDouble* ge            = xNewZeroInit<IssmDouble>(numvertices);
     1309        int*        vertexpidlist = xNew<int>(numvertices);
     1310
     1311        /*Build friction element, needed later: */
     1312        Friction* friction=new Friction(element,dim);
     1313
     1314        /*Retrieve all inputs we will be needing: */
     1315        element->GetVerticesCoordinatesBase(&xyz_list_base);
     1316        element->GradientIndexing(&vertexpidlist[0],control_index);
     1317        Input* vx_input        = element->GetInput(VxEnum);                   _assert_(vx_input);
     1318        Input* vy_input        = element->GetInput(VyEnum);                   _assert_(vy_input);
     1319        Input* vz_input        = element->GetInput(VzEnum);                   _assert_(vy_input);
     1320        Input* adjointx_input  = element->GetInput(AdjointxEnum);             _assert_(adjointx_input);
     1321        Input* adjointy_input  = element->GetInput(AdjointyEnum);             _assert_(adjointy_input);
     1322        Input* adjointz_input  = element->GetInput(AdjointzEnum);             _assert_(adjointz_input);
     1323        Input* dragcoeff_input = element->GetInput(FrictionCoefficientEnum);  _assert_(dragcoeff_input);
     1324
     1325        /* Start  looping on the number of gaussian points: */
     1326        Gauss* gauss=element->NewGaussBase(4);
     1327        for(int ig=gauss->begin();ig<gauss->end();ig++){
     1328                gauss->GaussPoint(ig);
     1329
     1330                adjointx_input->GetInputValue(&lambda, gauss);
     1331                adjointy_input->GetInputValue(&mu, gauss);
     1332                adjointz_input->GetInputValue(&xi    ,gauss);
     1333                vx_input->GetInputValue(&vx,gauss);
     1334                vy_input->GetInputValue(&vy,gauss);
     1335                vz_input->GetInputValue(&vz,gauss);
     1336                dragcoeff_input->GetInputValue(&drag, gauss);
     1337
     1338                friction->GetAlphaComplement(&dalpha2dk,gauss);
     1339                element->NormalBase(&normal[0],xyz_list_base);
     1340
     1341                element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss);
     1342                element->NodalFunctionsP1(basis,gauss);
     1343
     1344                /*Build gradient vector (actually -dJ/dk): */
     1345                for(int i=0;i<numvertices;i++){
     1346                        ge[i]+=(
     1347                                                -lambda*(2*drag*dalpha2dk*(vx - vz*normal[0]*normal[2]))
     1348                                                -mu    *(2*drag*dalpha2dk*(vy - vz*normal[1]*normal[2]))
     1349                                                -xi    *(2*drag*dalpha2dk*(-vx*normal[0]*normal[2]-vy*normal[1]*normal[2]))
     1350                                                )*Jdet*gauss->weight*basis[i];
     1351                        _assert_(!xIsNan<IssmDouble>(ge[i]));
     1352                }
     1353        }
     1354        gradient->SetValues(numvertices,vertexpidlist,ge,ADD_VAL);
     1355
     1356        /*Clean up and return*/
     1357        xDelete<IssmDouble>(xyz_list_base);
     1358        xDelete<IssmDouble>(basis);
     1359        xDelete<IssmDouble>(ge);
     1360        xDelete<int>(vertexpidlist);
     1361        delete gauss;
     1362        delete friction;
     1363}/*}}}*/
     1364void AdjointHorizAnalysis::GradientJBbarSSA(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
     1365
     1366        /*Intermediaries*/
     1367        int      domaintype,dim;
     1368        Element* basalelement;
     1369
     1370        /*Get basal element*/
     1371        element->FindParam(&domaintype,DomainTypeEnum);
     1372        switch(domaintype){
     1373                case Domain2DhorizontalEnum:
     1374                        basalelement = element;
     1375                        dim          = 2;
     1376                        break;
     1377                case Domain2DverticalEnum:
     1378                        if(!element->IsOnBase()) return;
     1379                        basalelement = element->SpawnBasalElement();
     1380                        dim          = 1;
     1381                        break;
     1382                case Domain3DEnum:
     1383                        if(!element->IsOnBase()) return;
     1384                        basalelement = element->SpawnBasalElement();
     1385                        dim          = 2;
     1386                        break;
     1387                default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
     1388        }
     1389
     1390        /*Intermediaries*/
     1391        IssmDouble Jdet,weight;
     1392        IssmDouble thickness,dmudB;
     1393        IssmDouble dvx[3],dvy[3],dadjx[3],dadjy[3];
     1394        IssmDouble *xyz_list= NULL;
     1395
     1396        /*Fetch number of vertices for this finite element*/
     1397        int numvertices = basalelement->GetNumberOfVertices();
     1398
     1399        /*Initialize some vectors*/
     1400        IssmDouble* basis         = xNew<IssmDouble>(numvertices);
     1401        IssmDouble* ge            = xNewZeroInit<IssmDouble>(numvertices);
     1402        int*        vertexpidlist = xNew<int>(numvertices);
     1403
     1404        /*Retrieve all inputs we will be needing: */
     1405        basalelement->GetVerticesCoordinates(&xyz_list);
     1406        basalelement->GradientIndexing(&vertexpidlist[0],control_index);
     1407        Input* thickness_input = basalelement->GetInput(ThicknessEnum);             _assert_(thickness_input);
     1408        Input* vx_input        = basalelement->GetInput(VxEnum);                    _assert_(vx_input);
     1409        Input* vy_input        = basalelement->GetInput(VyEnum);                    _assert_(vy_input);
     1410        Input* adjointx_input  = basalelement->GetInput(AdjointxEnum);              _assert_(adjointx_input);
     1411        Input* adjointy_input  = basalelement->GetInput(AdjointyEnum);              _assert_(adjointy_input);
     1412        Input* rheologyb_input = basalelement->GetInput(MaterialsRheologyBbarEnum); _assert_(rheologyb_input);
     1413
     1414        /* Start  looping on the number of gaussian points: */
     1415        Gauss* gauss=basalelement->NewGauss(4);
     1416        for(int ig=gauss->begin();ig<gauss->end();ig++){
     1417                gauss->GaussPoint(ig);
    11241418
    11251419                thickness_input->GetInputValue(&thickness,gauss);
    … …  
    11531447}/*}}}*/
    11541448void AdjointHorizAnalysis::GradientJBbarHO(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
    1155         _error_("not implemented yet");
     1449
     1450        /*WARNING: We use SSA as an estimate for now*/
     1451        this->GradientJBbarSSA(element,gradient,control_index);
    11561452}/*}}}*/
    11571453void AdjointHorizAnalysis::GradientJBbarFS(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
    1158         _error_("not implemented yet");
     1454        /*WARNING: We use SSA as an estimate for now*/
     1455        this->GradientJBbarSSA(element,gradient,control_index);
    11591456}/*}}}*/
    11601457void AdjointHorizAnalysis::GradientJDSSA(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
    1161         _error_("not implemented yet");
     1458
     1459        /*Intermediaries*/
     1460        int      domaintype,dim;
     1461        Element* basalelement;
     1462
     1463        /*Get basal element*/
     1464        element->FindParam(&domaintype,DomainTypeEnum);
     1465        switch(domaintype){
     1466                case Domain2DhorizontalEnum:
     1467                        basalelement = element;
     1468                        dim          = 2;
     1469                        break;
     1470                case Domain2DverticalEnum:
     1471                        if(!element->IsOnBase()) return;
     1472                        basalelement = element->SpawnBasalElement();
     1473                        dim          = 1;
     1474                        break;
     1475                case Domain3DEnum:
     1476                        if(!element->IsOnBase()) return;
     1477                        basalelement = element->SpawnBasalElement();
     1478                        dim          = 2;
     1479                        break;
     1480                default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
     1481        }
     1482
     1483        /*Intermediaries*/
     1484        IssmDouble Jdet,weight;
     1485        IssmDouble thickness,dmudD;
     1486        IssmDouble dvx[3],dvy[3],dadjx[3],dadjy[3];
     1487        IssmDouble *xyz_list= NULL;
     1488
     1489        /*Fetch number of vertices for this finite element*/
     1490        int numvertices = basalelement->GetNumberOfVertices();
     1491
     1492        /*Initialize some vectors*/
     1493        IssmDouble* basis         = xNew<IssmDouble>(numvertices);
     1494        IssmDouble* ge            = xNewZeroInit<IssmDouble>(numvertices);
     1495        int*        vertexpidlist = xNew<int>(numvertices);
     1496
     1497        /*Retrieve all inputs we will be needing: */
     1498        basalelement->GetVerticesCoordinates(&xyz_list);
     1499        basalelement->GradientIndexing(&vertexpidlist[0],control_index);
     1500        Input* thickness_input = basalelement->GetInput(ThicknessEnum);             _assert_(thickness_input);
     1501        Input* vx_input        = basalelement->GetInput(VxEnum);                    _assert_(vx_input);
     1502        Input* vy_input        = basalelement->GetInput(VyEnum);                    _assert_(vy_input);
     1503        Input* adjointx_input  = basalelement->GetInput(AdjointxEnum);              _assert_(adjointx_input);
     1504        Input* adjointy_input  = basalelement->GetInput(AdjointyEnum);              _assert_(adjointy_input);
     1505        Input* rheologyb_input = basalelement->GetInput(MaterialsRheologyBbarEnum); _assert_(rheologyb_input);
     1506
     1507        /* Start  looping on the number of gaussian points: */
     1508        Gauss* gauss=basalelement->NewGauss(4);
     1509        for(int ig=gauss->begin();ig<gauss->end();ig++){
     1510                gauss->GaussPoint(ig);
     1511
     1512                thickness_input->GetInputValue(&thickness,gauss);
     1513                vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
     1514                vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss);
     1515                adjointx_input->GetInputDerivativeValue(&dadjx[0],xyz_list,gauss);
     1516                adjointy_input->GetInputDerivativeValue(&dadjy[0],xyz_list,gauss);
     1517
     1518                basalelement->dViscositydDSSA(&dmudD,dim,xyz_list,gauss,vx_input,vy_input);
     1519
     1520                basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
     1521                basalelement->NodalFunctionsP1(basis,gauss);
     1522
     1523                for(int i=0;i<numvertices;i++){
     1524                        ge[i]+=-dmudD*thickness*(
     1525                                                (2*dvx[0]+dvy[1])*2*dadjx[0]+(dvx[1]+dvy[0])*(dadjx[1]+dadjy[0])+(2*dvy[1]+dvx[0])*2*dadjy[1]
     1526                                                )*Jdet*gauss->weight*basis[i];
     1527                        _assert_(!xIsNan<IssmDouble>(ge[i]));
     1528                }
     1529        }
     1530        gradient->SetValues(numvertices,vertexpidlist,ge,ADD_VAL);
     1531
     1532        /*Clean up and return*/
     1533        xDelete<IssmDouble>(xyz_list);
     1534        xDelete<IssmDouble>(basis);
     1535        xDelete<IssmDouble>(ge);
     1536        xDelete<int>(vertexpidlist);
     1537        delete gauss;
     1538        if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
    11621539}/*}}}*/
    11631540void AdjointHorizAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/

  • issm/trunk-jpl/src/c/classes/Elements/Element.cpp

    r18059 r18060  
    234234        /*Get viscosity*/
    235235        material->GetViscosity_B(&dmudB,eps_eff);
     236
     237        /*Assign output pointer*/
     238        *pdmudB=dmudB;
     239
     240}
     241/*}}}*/
     242void       Element::dViscositydDSSA(IssmDouble* pdmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){/*{{{*/
     243
     244        /*Intermediaries*/
     245        IssmDouble dmudB;
     246        IssmDouble epsilon2d[3];/* epsilon=[exx,eyy,exy];    */
     247        IssmDouble epsilon1d;   /* epsilon=[exx];    */
     248        IssmDouble eps_eff;
     249
     250        if(dim==2){
     251                /* eps_eff^2 = exx^2 + eyy^2 + exy^2 + exx*eyy*/
     252                this->StrainRateSSA(&epsilon2d[0],xyz_list,gauss,vx_input,vy_input);
     253                eps_eff = sqrt(epsilon2d[0]*epsilon2d[0] + epsilon2d[1]*epsilon2d[1] + epsilon2d[2]*epsilon2d[2] + epsilon2d[0]*epsilon2d[1]);
     254        }
     255        else{
     256                /* eps_eff^2 = 1/2 exx^2*/
     257                this->StrainRateSSA1d(&epsilon1d,xyz_list,gauss,vx_input);
     258                eps_eff = sqrt(epsilon1d*epsilon1d/2.);
     259        }
     260
     261        /*Get viscosity*/
     262        material->GetViscosity_D(&dmudB,eps_eff);
    236263
    237264        /*Assign output pointer*/

  • issm/trunk-jpl/src/c/classes/Elements/Element.h

    r18059 r18060  
    6161                void       DeleteMaterials(void);
    6262                void       dViscositydBSSA(IssmDouble* pdmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input);
     63                void       dViscositydDSSA(IssmDouble* pdmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input);
    6364                IssmDouble Divergence(void);
    6465                void       ThermalToEnthalpy(IssmDouble* penthalpy,IssmDouble temperature,IssmDouble waterfraction,IssmDouble pressure);

  • issm/trunk-jpl/src/c/classes/Elements/Penta.cpp

    r18058 r18060  
    30003000void       Penta::GradjDragHO(Vector<IssmDouble>* gradient,int control_index){/*{{{*/
    30013001
     3002        printf("-------------- file: Penta.cpp line: %i\n",__LINE__);
    30023003        int        i,j;
    30033004        int        analysis_type;

  • issm/trunk-jpl/src/c/classes/Elements/Tria.cpp

    r18058 r18060  
    14571457}
    14581458/*}}}*/
     1459void       Tria::NodalFunctionsP1(IssmDouble* basis, Gauss* gauss){/*{{{*/
     1460
     1461        _assert_(gauss->Enum()==GaussTriaEnum);
     1462        this->GetNodalFunctions(basis,(GaussTria*)gauss,P1Enum);
     1463
     1464}
     1465/*}}}*/
    14591466void       Tria::NodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
    14601467
    14611468        _assert_(gauss->Enum()==GaussTriaEnum);
    14621469        this->GetNodalFunctionsDerivatives(dbasis,xyz_list,(GaussTria*)gauss);
     1470
     1471}
     1472/*}}}*/
     1473void       Tria::NodalFunctionsP1Derivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
     1474
     1475        _assert_(gauss->Enum()==GaussTriaEnum);
     1476        this->GetNodalFunctionsDerivatives(dbasis,xyz_list,(GaussTria*)gauss,P1Enum);
    14631477
    14641478}

  • issm/trunk-jpl/src/c/classes/Elements/Tria.h

    r18058 r18060  
    191191                Gauss*         NewGaussTop(int order);
    192192                void           NodalFunctions(IssmDouble* basis,Gauss* gauss);
    193                 void           NodalFunctionsP1(IssmDouble* basis,Gauss* gauss){_error_("not implemented yet");};
     193                void           NodalFunctionsP1(IssmDouble* basis,Gauss* gauss);
    194194                void           NodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss);
    195                 void           NodalFunctionsP1Derivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss){_error_("not implemented yet");};
     195                void           NodalFunctionsP1Derivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss);
    196196                void           NodalFunctionsMINIDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss){_error_("not implemented yet");};
    197197                void           NodalFunctionsDerivativesVelocity(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss);

  • issm/trunk-jpl/src/c/classes/Elements/TriaRef.cpp

    r17962 r18060  
    44
    55/*Headers:*/
    6 /*{{{*//*{{{*/
     6/*{{{*/
    77#ifdef HAVE_CONFIG_H
    88        #include <config.h>

  • issm/trunk-jpl/src/c/classes/Materials/Material.h

    r18059 r18060  
    2727                virtual void       GetViscosity(IssmDouble* pviscosity,IssmDouble epseff)=0;
    2828                virtual void       GetViscosity_B(IssmDouble* pviscosity,IssmDouble epseff)=0;
     29                virtual void       GetViscosity_D(IssmDouble* pviscosity,IssmDouble epseff)=0;
    2930                virtual void       GetViscosityBar(IssmDouble* pviscosity,IssmDouble epseff)=0;
    3031                virtual void       GetViscosityComplement(IssmDouble* pviscosity_complement, IssmDouble* pepsilon)=0;

  • issm/trunk-jpl/src/c/classes/Materials/Matice.cpp

    r18059 r18060  
    285285/*}}}*/
    286286/*FUNCTION Matice::GetViscosity_B {{{*/
    287 void  Matice::GetViscosity_B(IssmDouble* pviscosity,IssmDouble eps_eff){
    288         /*From a string tensor and a material object, return viscosity, using Glen's flow law.
    289           (1-D)
    290           viscosity= -----------------------
    291           2 eps_eff ^[(n-1)/n]
    292 
    293           where viscosity is the viscotiy, B the flow law parameter , eps_eff is the effective strain rate
    294           and n the flow law exponent.
    295 
    296           If eps_eff = 0 , it means this is the first time SystemMatrices is being run, and we
    297           return 10^14, initial viscosity.
    298           */
    299 
    300         /*output: */
    301         IssmDouble viscosity;
     287void  Matice::GetViscosity_B(IssmDouble* pdmudB,IssmDouble eps_eff){
     288
     289        /*output: */
     290        IssmDouble dmudB;
    302291
    303292        /*Intermediary: */
    … …  
    311300        }
    312301
    313         if (n==1.){
    314                 /*Linear Viscous behavior (Newtonian fluid) viscosity=B/2: */
    315                 viscosity=(1.-D)/2.;
    316         }
    317         else{
    318 
    319                 /*if no strain rate, return maximum viscosity*/
    320                 if(eps_eff==0.){
    321                         viscosity = 1.e+14/2.;
    322                 }
    323 
    324                 else{
    325                         viscosity=(1.-D)/(2.*pow(eps_eff,(n-1.)/n));
    326                 }
    327         }
    328 
    329         /*Checks in debugging mode*/
    330         if(viscosity<=0) _error_("Negative viscosity");
     302        if(n==1.){
     303                /*Linear Viscous behavior (Newtonian fluid) dmudB=B/2: */
     304                dmudB=(1.-D)/2.;
     305        }
     306        else{
     307                if(eps_eff==0.) dmudB = 0.;
     308                else            dmudB = (1.-D)/(2.*pow(eps_eff,(n-1.)/n));
     309        }
    331310
    332311        /*Return: */
    333         *pviscosity=viscosity;
     312        *pdmudB=dmudB;
     313}
     314/*}}}*/
     315/*FUNCTION Matice::GetViscosity_D {{{*/
     316void  Matice::GetViscosity_D(IssmDouble* pdmudD,IssmDouble eps_eff){
     317
     318        /*output: */
     319        IssmDouble dmudD;
     320
     321        /*Intermediary: */
     322        IssmDouble n,B;
     323
     324        /*Get B and n*/
     325        n=GetN(); _assert_(n>0.);
     326        B=GetBbar();
     327        _assert_(this->isdamaged);
     328
     329        if(n==1.){
     330                /*Linear Viscous behavior (Newtonian fluid) dmudB=B/2: */
     331                dmudD=-B/2.;
     332        }
     333        else{
     334                if(eps_eff==0.) dmudD = 0.;
     335                else            dmudD = -B/(2.*pow(eps_eff,(n-1.)/n));
     336        }
     337
     338        /*Return: */
     339        *pdmudD=dmudD;
    334340}
    335341/*}}}*/

  • issm/trunk-jpl/src/c/classes/Materials/Matice.h

    r18059 r18060  
    5656                void       GetViscosity(IssmDouble* pviscosity, IssmDouble eps_eff);
    5757                void       GetViscosity_B(IssmDouble* pviscosity, IssmDouble eps_eff);
     58                void       GetViscosity_D(IssmDouble* pviscosity, IssmDouble eps_eff);
    5859                void       GetViscosityBar(IssmDouble* pviscosity, IssmDouble eps_eff);
    5960                void       GetViscosityComplement(IssmDouble* pviscosity_complement, IssmDouble* pepsilon);

  • issm/trunk-jpl/src/c/classes/Materials/Matpar.h

    r18059 r18060  
    7777                void       GetViscosity(IssmDouble* pviscosity,IssmDouble eps_eff){_error_("not supported");};
    7878                void       GetViscosity_B(IssmDouble* pviscosity,IssmDouble eps_eff){_error_("not supported");};
     79                void       GetViscosity_D(IssmDouble* pviscosity,IssmDouble eps_eff){_error_("not supported");};
    7980                void       GetViscosityBar(IssmDouble* pviscosity,IssmDouble eps_eff){_error_("not supported");};
    8081                void       GetViscosityComplement(IssmDouble* pviscosity_complement, IssmDouble* pepsilon){_error_("not supported");};

  • issm/trunk-jpl/src/c/modules/Gradjx/Gradjx.cpp

    r17899 r18060  
    99void Gradjx(Vector<IssmDouble>** pgradient,IssmDouble** pnorm_list, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters){
    1010
    11         int     i,j,numberofvertices;
    12         int     num_controls;
     11        int          numberofvertices;
     12        int          num_controls,analysisenum;
    1313        IssmDouble   norm_inf;
    1414        IssmDouble  *norm_list     = NULL;
    … …  
    2222        numberofvertices=vertices->NumberOfVertices();
    2323
     24        /*Get current analysis*/
     25        parameters->FindParam(&analysisenum,AnalysisTypeEnum);
     26        Analysis* analysis = EnumToAnalysis(analysisenum);
     27
    2428        /*Allocate gradient_list */
    2529        gradient_list = xNew<Vector<IssmDouble>*>(num_controls);
    2630        norm_list = xNew<IssmDouble>(num_controls);
    27         for(i=0;i<num_controls;i++){
     31        for(int i=0;i<num_controls;i++){
    2832                gradient_list[i]=new Vector<IssmDouble>(num_controls*numberofvertices);
    2933        }
    … …  
    3135
    3236        /*Compute all gradient_list*/
    33         for(i=0;i<num_controls;i++){
    34 
    35                 for(j=0;j<elements->Size();j++){
    36 
     37        for(int i=0;i<num_controls;i++){
     38                for(int j=0;j<elements->Size();j++){
    3739                        Element* element=(Element*)elements->GetObjectByOffset(j);
    38                         element->Gradj(gradient_list[i],control_type[i],i);
     40                        //element->Gradj(gradient_list[i],control_type[i],i);
     41                        analysis->GradientJ(gradient_list[i],element,control_type[i],i);
    3942                }
    40 
    4143                gradient_list[i]->Assemble();
    42 
    4344                norm_list[i]=gradient_list[i]->Norm(NORM_INF);
    4445        }
    4546
    4647        /*Add all gradient_list together*/
    47         for(i=0;i<num_controls;i++){
     48        for(int i=0;i<num_controls;i++){
    4849                gradient->AXPY(gradient_list[i],1.0);
    4950                delete gradient_list[i];
    … …  
    5253        /*Check that gradient is clean*/
    5354        norm_inf=gradient->Norm(NORM_INF);
    54         if(norm_inf<=0)    _error_("||∂J/∂α||∞ = 0    gradient norm is zero");
    55         if(xIsNan<IssmDouble>(norm_inf))_error_("||∂J/∂α||∞ = NaN  gradient norm is NaN");
     55        if(norm_inf<=0)                 _error_("||dJ/dk|| = 0    gradient norm is zero");
     56        if(xIsNan<IssmDouble>(norm_inf))_error_("||dJ/dk|| = NaN  gradient norm is NaN");
    5657
    5758        /*Clean-up and assign output pointer*/

  • issm/trunk-jpl/src/c/modules/Gradjx/Gradjx.h

    r17899 r18060  
    77
    88#include "../../classes/classes.h"
     9#include "../../analyses/analyses.h"
    910
    1011/* local prototypes: */
Note: See TracChangeset for help on using the changeset viewer.